In these exercises, we’ll look at how per capita death rates are influenced by age structure, taking advantage of the stable population model.
The stable population model allows us to compare the per capita (“crude”) death rates of stable populations with the same age-specific mortality rates and differing growth rates.
We will use a stable population growth rate of +2 per cent as a loose approximation of India and a stable population with a growth rate of 0 per cent as a loose approximation of the United States.
Preliminaries
Upload the data and define some variables. Note we’re going to use “x” for age, although the lectures used “a”.
library(data.table)
dt <- fread("exercise_1_data.csv")
x = dt$x
x.mid = x + .5 ## mid-point estimate for age group.
Lx = dt$Lx
hx = dt$hx
Nx = dt$Nx.usa ## population counts
A. Calculate difference in crude death rates for r = +2% and 0% via simulation.
r.0 = 0 ## resembling USA
r.2 = +0.02 ## resembling India
## stable age structures
cx.0 = exp(-r.0 * x) * Lx / sum( exp(-r.0 * x) * Lx )
cx.2 = exp(-r.2 * x) * Lx / sum( exp(-r.2 * x) * Lx )
plot(cx.2, x)
lines(cx.0, x)

Q1: Which population is the solid line and which is the dots?
Now let’s get the crude death rates (AKA “per-capita” aggregate death rates)
## crude death rates, with different age structures but same age-specific mortality
cdr.0 = sum(cx.0 * hx)
cdr.2 = sum(cx.2 * hx)
print("CDR with zero growth")
[1] "CDR with zero growth"
print(round(cdr.0, 4))
[1] 0.0124
print("CDR with +2% growth")
[1] "CDR with +2% growth"
print(round(cdr.2, 4))
[1] 0.0056
Let’s convert this to a ratio:
print("simulated ratio of per capita mortality in high growth to low growth pops")
[1] "simulated ratio of per capita mortality in high growth to low growth pops"
cdr.ratio.simu = cdr.2/cdr.0
print(cdr.ratio.simu)
[1] 0.451164
## [1] 0.451164, less than half!
Q2. Does it appear that the historical growth rate of a population matters “a little” or “a lot” for per capita Covid mortality?
B. Calculate using comparative statics
Now we’re going to use our comparative statics result to re-estimate the change in per capita death rates that results from changing the growth rate.
Our result was \[
\Delta \log CDR(r) = \Delta r (A_0 - e_0)
\]
So, we’ll start by getting the mean age of the stationary population (\(A_0\)) and life expectancy (\(e_0\)).
## mean age of stationary pop
A0 = sum(x.mid * Lx) / sum(Lx)
print("Mean age of stationary pop")
[1] "Mean age of stationary pop"
print(A0)
[1] 41.34563
## mean age of death in stationary pop
e0 = sum(Lx)
print("Mean age of death in stationary pop")
[1] "Mean age of death in stationary pop"
print(e0)
[1] 79.10057
## Change in growth rate
Delta.r = r.2 - r.0
print("Change in 'r'")
[1] "Change in 'r'"
print(Delta.r)
[1] 0.02
Now we’re ready to use our formula
## Comparative statics estimate of how much log of CDR changes
Delta.log.cdr.hat = (A0 - e0) * Delta.r
print(Delta.log.cdr.hat)
In order to make sense of this number, which tells us the change in the log of the death rate, let’s convert to a ratio. Exponentiating should do the trick, since
\[
\exp[\log(A) - \log(B)] = A/B
\]
## Convert to ratio, since A/B = exp(log(A) - log(B))
cdr.ratio.hat = exp(Delta.log.cdr.hat)
print("comparative statics approxiamtion of same ratio")
print(cdr.ratio.hat)
Let’s compare the two results
result = c("cdr.ratio.simu" = cdr.ratio.simu,
"cdr.ratio.hat" = cdr.ratio.hat)
print(round(result, 3))
Q2. How accurate is the comparative statics approximation? “Super accurate”, “pretty accurate”, or “not very accurate”?
Note: please don’t take the good accuracy of this one example as typical. Often the approximations are much worse, because we’re only using the linear term of Taylor approximation.
C. Adapt for Covid
Until now we’ve been discussing the all-cause mortality rate. But the logic should carry over to the per-capita Covid mortality rate.
We’ll use the standard Covid mortality age schedule estimated by Goldstein and Lee, and then use our comparative statics approach to compare the two stable populations.
Import the covid age-specific death rates:
## Import Goldstein and Lee's schedule of Covid mortality rates
## (Note: these have realistic age profile, but the level is set to produce 1 million deaths in US)
dt <- fread("https://raw.githubusercontent.com/josh-goldstein-git/dempersp_covid_mortality/master/data/cleaned/normalized_nMx_out.csv")
hx.covid = dt[x %in% 0:100]$Mx_covid_scaled
The approximation we want to use for the effect of the growth rate on per capita Covid deaths needs to be adapted so that it uses the mean age of death from Covid in the stationary population. It should be
\[
\Delta \log CDR_{covid} = (A_0 - A_{D(covid)}) \Delta r,
\] where the age of death in the stationary pop from Covid is \[
A_{D(covid)} = {\int x h_x^{covid} \ell(x) \,dx \over \int h_x^{covid}\ell(x), \,dx}
\]
A.D.covid = sum(x *hx.covid * Lx) / sum(hx.covid*Lx)
print("Mean age of death from Covid in stationary pop")
print(A.D.covid)
Q3. Is the mean age of death from Covid higher or lower than e0? Is this what you expected?
Q4. Would you expect per capita death rates for covid to be more or less sensitive to the rates of historical growth that determine age structure? Can your answer to Q3 help here?
Q5. Looking back at the graph from Science, after doing this exercise do you think it’s still a mystery why per capita Covid mortality is lower in India than the U.S.?
Q6. Advanced question that we probably won’t have time to discuss but Josh is happy to do on Piazza: How might you do a comparative statics type calculation for per capita Covid mortality comparing countries with similar growth rates but different longevity? (e.g., France and USA). Or countries with different growth rates and different longevity (e.g., Italy and USA)
Congratulations.
You’ve finished the 1st set of exercises of Day 2 of the workshop. One down and two more to go!
LS0tCnRpdGxlOiAiTW9ydGFsaXR5IEV4ZXJjaXNlcyAxOiBQZXIgQ2FwaXRhIERlYXRoIFJhdGVzIgphdXRob3I6IEpvc2h1YSBSLiBHb2xkc3RlaW4Kb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKSW4gdGhlc2UgZXhlcmNpc2VzLCB3ZSdsbCBsb29rIGF0IGhvdyBwZXIgY2FwaXRhIGRlYXRoIHJhdGVzIGFyZSBpbmZsdWVuY2VkIGJ5IGFnZSBzdHJ1Y3R1cmUsIHRha2luZyBhZHZhbnRhZ2Ugb2YgdGhlIHN0YWJsZSBwb3B1bGF0aW9uIG1vZGVsLiAKClRoZSBzdGFibGUgcG9wdWxhdGlvbiBtb2RlbCBhbGxvd3MgdXMgdG8gY29tcGFyZSB0aGUgcGVyIGNhcGl0YSAoImNydWRlIikgZGVhdGggcmF0ZXMgb2Ygc3RhYmxlIHBvcHVsYXRpb25zIHdpdGggdGhlIHNhbWUgYWdlLXNwZWNpZmljIG1vcnRhbGl0eSByYXRlcyBhbmQgZGlmZmVyaW5nIGdyb3d0aCByYXRlcy4gCgpXZSB3aWxsIHVzZSBhIHN0YWJsZSBwb3B1bGF0aW9uIGdyb3d0aCByYXRlIG9mICsyIHBlciBjZW50IGFzIGEgbG9vc2UgYXBwcm94aW1hdGlvbiBvZiBJbmRpYSBhbmQgYSBzdGFibGUgcG9wdWxhdGlvbiB3aXRoIGEgZ3Jvd3RoIHJhdGUgb2YgMCBwZXIgY2VudCBhcyBhIGxvb3NlIGFwcHJveGltYXRpb24gb2YgdGhlIFVuaXRlZCBTdGF0ZXMuCgojIyBQcmVsaW1pbmFyaWVzCgpVcGxvYWQgdGhlIGRhdGEgYW5kIGRlZmluZSBzb21lIHZhcmlhYmxlcy4gTm90ZSB3ZSdyZSBnb2luZyB0byB1c2UgIngiIGZvciBhZ2UsIGFsdGhvdWdoIHRoZSBsZWN0dXJlcyB1c2VkICJhIi4KCmBgYHtyfQpsaWJyYXJ5KGRhdGEudGFibGUpCmR0IDwtIGZyZWFkKCJleGVyY2lzZV8xX2RhdGEuY3N2IikKeCA9IGR0JHgKeC5taWQgPSB4ICsgLjUgICMjIG1pZC1wb2ludCBlc3RpbWF0ZSBmb3IgYWdlIGdyb3VwLgpMeCA9IGR0JEx4IApoeCA9IGR0JGh4Ck54ID0gZHQkTngudXNhICMjIHBvcHVsYXRpb24gY291bnRzCmBgYAoKIyMgIEEuIENhbGN1bGF0ZSBkaWZmZXJlbmNlIGluIGNydWRlIGRlYXRoIHJhdGVzIGZvciByID0gKzIlIGFuZCAwJSB2aWEgc2ltdWxhdGlvbi4KCmBgYHtyfQpyLjAgPSAwICMjIHJlc2VtYmxpbmcgVVNBCnIuMiA9ICswLjAyICMjIHJlc2VtYmxpbmcgSW5kaWEKCiMjIHN0YWJsZSBhZ2Ugc3RydWN0dXJlcwpjeC4wID0gZXhwKC1yLjAgKiB4KSAqIEx4IC8gc3VtKCBleHAoLXIuMCAqIHgpICogTHggKQpjeC4yID0gZXhwKC1yLjIgKiB4KSAqIEx4IC8gc3VtKCBleHAoLXIuMiAqIHgpICogTHggKQpwbG90KGN4LjIsIHgpCmxpbmVzKGN4LjAsIHgpCmBgYAoKPiBRMTogV2hpY2ggcG9wdWxhdGlvbiBpcyB0aGUgc29saWQgbGluZSBhbmQgd2hpY2ggaXMgdGhlIGRvdHM/CgpOb3cgbGV0J3MgZ2V0IHRoZSBjcnVkZSBkZWF0aCByYXRlcyAoQUtBICJwZXItY2FwaXRhIiBhZ2dyZWdhdGUgZGVhdGggcmF0ZXMpCmBgYHtyfQojIyBjcnVkZSBkZWF0aCByYXRlcywgd2l0aCBkaWZmZXJlbnQgYWdlIHN0cnVjdHVyZXMgYnV0IHNhbWUgYWdlLXNwZWNpZmljIG1vcnRhbGl0eQpjZHIuMCA9IHN1bShjeC4wICogaHgpCmNkci4yID0gc3VtKGN4LjIgKiBoeCkKCnByaW50KCJDRFIgd2l0aCB6ZXJvIGdyb3d0aCIpCnByaW50KHJvdW5kKGNkci4wLCA0KSkKcHJpbnQoIkNEUiB3aXRoICsyJSBncm93dGgiKQpwcmludChyb3VuZChjZHIuMiwgNCkpCmBgYAoKTGV0J3MgY29udmVydCB0aGlzIHRvIGEgcmF0aW86CmBgYHtyfQpwcmludCgic2ltdWxhdGVkIHJhdGlvIG9mIHBlciBjYXBpdGEgbW9ydGFsaXR5IGluIGhpZ2ggZ3Jvd3RoIHRvIGxvdyBncm93dGggcG9wcyIpCmNkci5yYXRpby5zaW11ID0gY2RyLjIvY2RyLjAKcHJpbnQoY2RyLnJhdGlvLnNpbXUpCiMjIFsxXSAwLjQ1MTE2NCwgbGVzcyB0aGFuIGhhbGYhCmBgYAoKPiBRMi4gRG9lcyBpdCBhcHBlYXIgdGhhdCB0aGUgaGlzdG9yaWNhbCBncm93dGggcmF0ZSBvZiBhIHBvcHVsYXRpb24KPiBtYXR0ZXJzICJhIGxpdHRsZSIgb3IgImEgbG90IiAgZm9yIHBlciBjYXBpdGEgQ292aWQgbW9ydGFsaXR5PwoKIyMgQi4gQ2FsY3VsYXRlIHVzaW5nIGNvbXBhcmF0aXZlIHN0YXRpY3MKCk5vdyB3ZSdyZSBnb2luZyB0byB1c2Ugb3VyIGNvbXBhcmF0aXZlIHN0YXRpY3MgcmVzdWx0IHRvIHJlLWVzdGltYXRlIHRoZSBjaGFuZ2UgaW4gcGVyIGNhcGl0YSBkZWF0aCByYXRlcyB0aGF0IHJlc3VsdHMgZnJvbSBjaGFuZ2luZyB0aGUgZ3Jvd3RoIHJhdGUuCgpPdXIgcmVzdWx0IHdhcwokJApcRGVsdGEgXGxvZyBDRFIocikgPSBcRGVsdGEgciAoQV8wIC0gZV8wKQokJAoKU28sIHdlJ2xsIHN0YXJ0IGJ5IGdldHRpbmcgdGhlIG1lYW4gYWdlIG9mIHRoZSBzdGF0aW9uYXJ5IHBvcHVsYXRpb24gKCRBXzAkKSBhbmQgbGlmZSBleHBlY3RhbmN5ICgkZV8wJCkuCgpgYGB7cn0KIyMgbWVhbiBhZ2Ugb2Ygc3RhdGlvbmFyeSBwb3AKQTAgPSBzdW0oeC5taWQgKiBMeCkgLyBzdW0oTHgpCnByaW50KCJNZWFuIGFnZSBvZiBzdGF0aW9uYXJ5IHBvcCIpCnByaW50KEEwKQojIyBtZWFuIGFnZSBvZiBkZWF0aCBpbiBzdGF0aW9uYXJ5IHBvcAplMCA9IHN1bShMeCkKcHJpbnQoIk1lYW4gYWdlIG9mIGRlYXRoIGluIHN0YXRpb25hcnkgcG9wIikKcHJpbnQoZTApCgojIyBDaGFuZ2UgaW4gZ3Jvd3RoIHJhdGUKRGVsdGEuciA9IHIuMiAtIHIuMApwcmludCgiQ2hhbmdlIGluICdyJyIpCnByaW50KERlbHRhLnIpCmBgYAoKTm93IHdlJ3JlIHJlYWR5IHRvIHVzZSBvdXIgZm9ybXVsYQoKYGBge3J9CiMjIENvbXBhcmF0aXZlIHN0YXRpY3MgZXN0aW1hdGUgb2YgaG93IG11Y2ggbG9nIG9mIENEUiBjaGFuZ2VzCkRlbHRhLmxvZy5jZHIuaGF0ID0gKEEwIC0gZTApICogRGVsdGEucgpwcmludChEZWx0YS5sb2cuY2RyLmhhdCkKYGBgCgpJbiBvcmRlciB0byBtYWtlIHNlbnNlIG9mIHRoaXMgbnVtYmVyLCB3aGljaCB0ZWxscyB1cyB0aGUgY2hhbmdlIGluIHRoZSBsb2cgb2YgdGhlIGRlYXRoIHJhdGUsIGxldCdzIGNvbnZlcnQgdG8gYSByYXRpby4gRXhwb25lbnRpYXRpbmcgc2hvdWxkIGRvIHRoZSB0cmljaywgc2luY2UgCgokJApcZXhwW1xsb2coQSkgLSBcbG9nKEIpXSA9IEEvQgokJAoKYGBge3J9CiMjIENvbnZlcnQgdG8gcmF0aW8sIHNpbmNlIEEvQiA9IGV4cChsb2coQSkgLSBsb2coQikpCmNkci5yYXRpby5oYXQgPSBleHAoRGVsdGEubG9nLmNkci5oYXQpCnByaW50KCJjb21wYXJhdGl2ZSBzdGF0aWNzIGFwcHJveGlhbXRpb24gb2Ygc2FtZSByYXRpbyIpCnByaW50KGNkci5yYXRpby5oYXQpCmBgYApMZXQncyBjb21wYXJlIHRoZSB0d28gcmVzdWx0cwpgYGB7cn0KcmVzdWx0ID0gYygiY2RyLnJhdGlvLnNpbXUiID0gY2RyLnJhdGlvLnNpbXUsCiAgICAgICAgICAgImNkci5yYXRpby5oYXQiID0gY2RyLnJhdGlvLmhhdCkKcHJpbnQocm91bmQocmVzdWx0LCAzKSkKYGBgCgo+IFEyLiBIb3cgYWNjdXJhdGUgaXMgdGhlIGNvbXBhcmF0aXZlIHN0YXRpY3MgYXBwcm94aW1hdGlvbj8KPiAiU3VwZXIgYWNjdXJhdGUiLCAicHJldHR5IGFjY3VyYXRlIiwgb3IgIm5vdCB2ZXJ5IGFjY3VyYXRlIj8KCk5vdGU6IHBsZWFzZSBkb24ndCB0YWtlIHRoZSBnb29kIGFjY3VyYWN5IG9mIHRoaXMgb25lIGV4YW1wbGUgYXMKdHlwaWNhbC4gT2Z0ZW4gdGhlIGFwcHJveGltYXRpb25zIGFyZSBtdWNoIHdvcnNlLCBiZWNhdXNlIHdlJ3JlCm9ubHkgdXNpbmcgdGhlIGxpbmVhciB0ZXJtIG9mIFRheWxvciBhcHByb3hpbWF0aW9uLgoKIyMgQy4gQWRhcHQgZm9yIENvdmlkCgpVbnRpbCBub3cgd2UndmUgYmVlbiBkaXNjdXNzaW5nIHRoZSBhbGwtY2F1c2UgbW9ydGFsaXR5IHJhdGUuIEJ1dCB0aGUgbG9naWMgc2hvdWxkIGNhcnJ5IG92ZXIgdG8gdGhlIHBlci1jYXBpdGEgQ292aWQgbW9ydGFsaXR5IHJhdGUuCgpXZSdsbCB1c2UgdGhlIHN0YW5kYXJkIENvdmlkIG1vcnRhbGl0eSBhZ2Ugc2NoZWR1bGUgZXN0aW1hdGVkIGJ5IEdvbGRzdGVpbiBhbmQgTGVlLCBhbmQgdGhlbiB1c2Ugb3VyIGNvbXBhcmF0aXZlIHN0YXRpY3MgYXBwcm9hY2ggdG8gY29tcGFyZSB0aGUgdHdvIHN0YWJsZSBwb3B1bGF0aW9ucy4KCkltcG9ydCB0aGUgY292aWQgYWdlLXNwZWNpZmljIGRlYXRoIHJhdGVzOgpgYGB7cn0KIyMgSW1wb3J0IEdvbGRzdGVpbiBhbmQgTGVlJ3Mgc2NoZWR1bGUgb2YgQ292aWQgbW9ydGFsaXR5IHJhdGVzCiMjIChOb3RlOiB0aGVzZSBoYXZlIHJlYWxpc3RpYyBhZ2UgcHJvZmlsZSwgYnV0IHRoZSBsZXZlbCBpcyBzZXQgdG8gcHJvZHVjZSAxIG1pbGxpb24gZGVhdGhzIGluIFVTKQpkdCA8LSBmcmVhZCgiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL2pvc2gtZ29sZHN0ZWluLWdpdC9kZW1wZXJzcF9jb3ZpZF9tb3J0YWxpdHkvbWFzdGVyL2RhdGEvY2xlYW5lZC9ub3JtYWxpemVkX25NeF9vdXQuY3N2IikKaHguY292aWQgPSBkdFt4ICVpbiUgMDoxMDBdJE14X2NvdmlkX3NjYWxlZApgYGAKClRoZSBhcHByb3hpbWF0aW9uIHdlIHdhbnQgdG8gdXNlIGZvciB0aGUgZWZmZWN0IG9mIHRoZSBncm93dGggcmF0ZSBvbiBwZXIgY2FwaXRhIENvdmlkIGRlYXRocyBuZWVkcyB0byBiZSBhZGFwdGVkIHNvIHRoYXQgaXQgdXNlcyB0aGUgbWVhbiBhZ2Ugb2YgZGVhdGggZnJvbSBDb3ZpZCBpbiB0aGUgc3RhdGlvbmFyeSBwb3B1bGF0aW9uLiBJdCBzaG91bGQgYmUKCiQkClxEZWx0YSBcbG9nIENEUl97Y292aWR9ID0gKEFfMCAtIEFfe0QoY292aWQpfSkgXERlbHRhIHIsCiQkCndoZXJlIHRoZSBhZ2Ugb2YgZGVhdGggaW4gdGhlIHN0YXRpb25hcnkgcG9wIGZyb20gQ292aWQgaXMKJCQKQV97RChjb3ZpZCl9ID0ge1xpbnQgeCBoX3hee2NvdmlkfSBcZWxsKHgpIFwsZHggXG92ZXIgXGludCBoX3hee2NvdmlkfVxlbGwoeCksIFwsZHh9CiQkCgpgYGB7cn0KQS5ELmNvdmlkID0gc3VtKHggKmh4LmNvdmlkICogTHgpIC8gc3VtKGh4LmNvdmlkKkx4KQpwcmludCgiTWVhbiBhZ2Ugb2YgZGVhdGggZnJvbSBDb3ZpZCBpbiBzdGF0aW9uYXJ5IHBvcCIpCnByaW50KEEuRC5jb3ZpZCkKYGBgCgo+IFEzLiBJcyB0aGUgbWVhbiBhZ2Ugb2YgZGVhdGggZnJvbSBDb3ZpZCBoaWdoZXIgb3IgbG93ZXIgdGhhbiBlMD8gCj4gIElzIHRoaXMgd2hhdCB5b3UgZXhwZWN0ZWQ/Cgo+IFE0LiBXb3VsZCB5b3UgZXhwZWN0IHBlciBjYXBpdGEgZGVhdGggcmF0ZXMgZm9yIGNvdmlkIHRvIGJlIG1vcmUgb3IKPiAgbGVzcyBzZW5zaXRpdmUgdG8gdGhlIHJhdGVzIG9mIGhpc3RvcmljYWwgZ3Jvd3RoIHRoYXQgZGV0ZXJtaW5lIGFnZQo+ICBzdHJ1Y3R1cmU/IENhbiB5b3VyIGFuc3dlciB0byBRMyBoZWxwIGhlcmU/Cgo+IFE1LiBMb29raW5nIGJhY2sgYXQgdGhlIGdyYXBoIGZyb20gU2NpZW5jZSwgYWZ0ZXIgZG9pbmcgdGhpcwo+ICBleGVyY2lzZSBkbyB5b3UgdGhpbmsgaXQncyBzdGlsbCBhIG15c3Rlcnkgd2h5IHBlciBjYXBpdGEgQ292aWQKPiAgbW9ydGFsaXR5IGlzIGxvd2VyIGluIEluZGlhIHRoYW4gdGhlIFUuUy4/IAoKCj4gUTYuIEFkdmFuY2VkIHF1ZXN0aW9uIHRoYXQgd2UgcHJvYmFibHkgd29uJ3QgaGF2ZSB0aW1lIHRvIGRpc2N1c3MgCj4gYnV0IEpvc2ggaXMgaGFwcHkgdG8gZG8gb24gUGlhenphOiBIb3cgbWlnaHQgeW91IGRvIGEgY29tcGFyYXRpdmUgCj4gc3RhdGljcyB0eXBlIGNhbGN1bGF0aW9uIGZvciBwZXIgY2FwaXRhIENvdmlkIG1vcnRhbGl0eSBjb21wYXJpbmcgCj4gY291bnRyaWVzIHdpdGggc2ltaWxhciBncm93dGggcmF0ZXMgYnV0IGRpZmZlcmVudCBsb25nZXZpdHk/Cj4gKGUuZy4sIEZyYW5jZSBhbmQgVVNBKS4gT3IgY291bnRyaWVzIHdpdGggZGlmZmVyZW50IGdyb3d0aCByYXRlcyBhbmQKPiBkaWZmZXJlbnQgbG9uZ2V2aXR5IChlLmcuLCBJdGFseSBhbmQgVVNBKQoKIyMgQ29uZ3JhdHVsYXRpb25zLiAKCllvdSd2ZSBmaW5pc2hlZCB0aGUgMXN0IHNldCBvZiBleGVyY2lzZXMgb2YgRGF5IDIgb2YgdGhlIHdvcmtzaG9wLiBPbmUgZG93biBhbmQgdHdvIG1vcmUgdG8gZ28hCgoKCg==