{"id":926,"date":"2024-05-19T20:00:30","date_gmt":"2024-05-20T00:00:30","guid":{"rendered":"https:\/\/www.clayford.net\/statistics\/?p=926"},"modified":"2024-05-21T13:30:45","modified_gmt":"2024-05-21T17:30:45","slug":"age-adjustment-section-2-4-of-regression-and-other-stories","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/age-adjustment-section-2-4-of-regression-and-other-stories\/","title":{"rendered":"Age Adjustment: section 2.4 of Regression and Other Stories"},"content":{"rendered":"

In section 2.4 of Regression and Other Stories<\/em>, Gelman, et al.\u00a0explain the necessity of age adjustment when investigating mortality rates. The book is freely available as a PDF online<\/a> and the section of interest is on pages 31-33. Upon first reading, I had trouble understanding what they were doing. In particular I didn\u2019t follow Figure 2.12. I\u2019m aware that says more about me and than the authors. Fortunately, there\u2019s a footnote in the book that says all data and code are available in the AgePeriodCohort<\/a> folder on GitHub. Good, I thought, I\u2019ll look at the code and figure out what\u2019s going on. Famous last words.<\/p>\n

The R script that performs the age adjustment is births.R<\/a>. It clocks in at over 400 lines has practically no comments outside of the occasional \u201cSum it up.\u201d As you run the code, you\u2019ll find the script generates several plots not in the book. In addition, the plots that are<\/em> in the book are generated in a different order. Trying to parse the R code to help me understand the exposition was frustrating. But I persisted.<\/p>\n

Reading the bibliographic note at the end of the chapter indicated the age adjustment example was first discussed on Gelman\u2019s blog<\/a>. In the blog post he walks through the process of age adjustment, creating the same plots in the book, and provides the R code<\/a>. This is basically the births.R script. He says at the end, \u201cthe code is ugly. Don\u2019t model your code after my practices! If any of you want to make a statistics lesson out of this episode, I recommend you clean the code.\u201d<\/p>\n

This blog post is my statistics lesson trying to understand and clean this code.<\/p>\n

\n

Fig 2.11 (a)<\/h2>\n

The data apparently come from the CDC<\/a>, but I\u2019m using the data file Gelman provides with his R code. The data shows number of deaths per age per gender per year for white non-hispanics in the US. For example, the first row shows 1291 female deaths (Male = 0) in 1999 for those who were 35 years old. The total population of 35 year old women in 1999 was 1,578,829. The rate is 1291\/1,578,829 x 100,000 = 81.8, or 81 deaths per 100,000.<\/p>\n

data <- read.table("white_nonhisp_death_rates_from_1999_to_2013_by_sex.txt", \r\n                   header=TRUE)\r\nhead(data)<\/code><\/pre>\n
##   Age Male Year Deaths Population Rate\r\n## 1  35    0 1999   1291    1578829 81.8\r\n## 2  35    0 2000   1264    1528463 82.7\r\n## 3  35    0 2001   1186    1377466 86.1\r\n## 4  35    0 2002   1194    1333639 89.5\r\n## 5  35    0 2003   1166    1302188 89.5\r\n## 6  35    0 2004   1166    1325435 88.0<\/code><\/pre>\n
The first plot is mortality rate of the 45-54 age group from 1999 – 2013. We first sum both Deaths and Population by year and then calculate the Mortality Rate by dividing Deaths by Population. This is a nice opportunity to use the base R pipe operator, |><\/code>.<\/p>\n
aggregate(cbind(Deaths, Population) ~ Year, data = data, FUN = sum, \r\n          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  plot(Rate ~ Year, data = _, type = "l", ylab = "Death Rate")<\/code><\/pre>\n
<\/p>\n
This is the third plot in Gelman\u2019s blog post titled \u201cSo take the ratio!\u201d<\/p>\n<\/div>\n
\n
Fig 2.11 (b)<\/h2>\n
The second plot shows the average of the 45-54 age group increasing as the baby boomers move through. First we sum the Population by Age and Year for the 45-54 group. Then we take that data and take the mean age per year weighted by the population.<\/p>\n
aggregate(Population ~ Age + Year, data = data, sum, \r\n          subset = Age %in% 45:54) |> \r\n  aggregate(Population ~ Year, data = _, \r\n            function(x)weighted.mean(45:54, x)) |> \r\n  plot(Population ~ Year, data = _, type = "l")<\/code><\/pre>\n
<\/p>\n
To help make this clear, let\u2019s find the mean age of the 45-54 group in 1999. First find the population for each age in 1999:<\/p>\n
tmp <- aggregate(Population ~ Age + Year, data = data, sum, \r\n          subset = Age %in% 45:54 & Year == 1999)\r\ntmp$Population<\/code><\/pre>\n
##  [1] 3166393 3007083 2986252 2805975 2859406 2868751 2804957 3093631 2148382\r\n## [10] 2254975<\/code><\/pre>\n
To find the mean age of the 45-54 group in 1999, we need to weight each age with the population. We can do that with the weighted.mean()<\/code> function.<\/p>\n
weighted.mean(45:54, tmp$Population)<\/code><\/pre>\n
## [1] 49.25585<\/code><\/pre>\n
The code above does this for 1999-2013. I think it\u2019s worth noting that while the plot looks dramatic, the average age only increases from about 49.2 to 49.7. But I suppose when you\u2019re dealing with millions of people that increase makes a difference.<\/p>\n
This is the fourth plot in Gelman\u2019s blog post titled \u201cBut the average age in this group is going up!\u201d<\/p>\n<\/div>\n
\n
Fig 2.11 (c)<\/h2>\n
This is where I began to struggle when reading the book.<\/p>\n
This figure is titled \u201cThe trend in raw death rates since 2005 can be explained by age-aggregation bias\u201d. This is the eighth plot in the blog post where it has a bit more motivation. Let\u2019s recreate the plots in the blog post leading up to this plot.<\/p>\n
The first plot is the sixth plot. It\u2019s basically the previous plot rescaled as a rate. It\u2019s created by first calculating the death rate in 1999, and then taking the weighted mean of that rate by using the total population for each age group.<\/p>\n
dr1999 <- aggregate(cbind(Deaths, Population) ~ Age, data = data, FUN = sum, \r\n          subset = Age %in% 45:54 & Year == 1999) |> \r\n  transform(Rate = Deaths\/Population) \r\n\r\n# Now create plot\r\naggregate(Population ~ Age + Year, data = data, sum, \r\n          subset = Age %in% 45:54) |> \r\n  aggregate(Population ~ Year, data = _, \r\n            function(x)weighted.mean(dr1999$Rate, x)) |> \r\n  plot(Population ~ Year, data = _, type = "l", ylab = "Reconstructed death rate")<\/code><\/pre>\n
<\/p>\n
Next he combines this plot with the plot of the raw death rate (Fig 2.11 (a)). This is the seventh plot in the blog post.<\/p>\n
years <- 1999:2013\r\n\r\nRaw <- aggregate(cbind(Deaths, Population) ~ Year, data = data, FUN = sum, \r\n          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population)\r\n\r\nExpected <- aggregate(Population ~ Age + Year, data = data, sum, \r\n          subset = Age %in% 45:54) |> \r\n  aggregate(Population ~ Year, data = _, \r\n            function(x)weighted.mean(dr1999$Rate, x))\r\n\r\nplot(years, Raw$Rate, type="l", ylab="Death rate for 45-54 non-Hisp whites")\r\nlines(years, Expected$Population, col="green4")\r\ntext(2002.5, .00404, "Raw death rate", cex=.8)\r\ntext(2009, .00394, "Expected just from\\nage shift", col="green4", cex=.8)<\/code><\/pre>\n
<\/p>\n
Then finally he says, \u201cWe can sharpen this comparison by anchoring the expected-trend-in-death-rate-just-from-changing-age-composition graph at 2013, the end of the time series, instead of 1999.\u201d This means we need to calculate the death rate in 2013, and then take the weighted mean of that rate by using the total population for each age group. This is the dr2013 data frame. Then we create the same plot as above except now using the death rate in 2013.<\/p>\n
dr2013 <- aggregate(cbind(Deaths, Population) ~ Age, data = data, FUN = sum, \r\n                    subset = Age %in% 45:54 & Year == 2013) |> \r\n  transform(Rate = Deaths\/Population) \r\n\r\nRaw <- aggregate(cbind(Deaths, Population) ~ Year, data = data, FUN = sum, \r\n                 subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population)\r\n\r\nExpected <- aggregate(Population ~ Age + Year, data = data, sum, \r\n                      subset = Age %in% 45:54) |> \r\n  aggregate(Population ~ Year, data = _, \r\n            function(x)weighted.mean(dr2013$Rate, x))\r\n\r\nplot(years, Raw$Rate,  type="l", \r\n     ylab="Death rate for 45-54 non-Hisp whites")\r\nlines(years, Expected$Population, col="green4")\r\ntext(2002.5, 0.00395, "Raw death rate", cex=.8)\r\ntext(2002, .00409, "Expected just from\\nage shift", col="green4", cex=.8)<\/code><\/pre>\n
<\/p>\n
Gelman notes, \u201csince 2003, all the changes in raw death rate in this group can be explained by changes in age composition.\u201d<\/p>\n<\/div>\n
\n
Fig 2.12 (a)<\/h2>\n
This is the first plot showing age-adjusted death rates. Gelman explains this as follows in his blog post: \u201cfor each year in time, we take the death rates by year of age and average them, thus computing the death rate that would\u2019ve been observed had the population distribution of 45-54-year-olds been completely flat each year.\u201d The book calls it \u201cthe simplest such adjustment, normalizing each year to a hypothetical uniformly distributed population in which the number of people is equal at each age from 45 through 54.\u201d I found this latter explanation a little confusing.<\/p>\n
To create this plot we first sum Deaths and Populations by age and year for the 45-54 age group, then calculate the death rate, and then simply take the mean rate by year. That\u2019s it. Gelman takes the additional step of rescaling the rate so that the rate is 1 in 1999.<\/p>\n
aggregate(cbind(Deaths, Population) ~ Age + Year, data = data, sum, \r\n          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  aggregate(Rate ~ Year, data = _, mean) |> \r\n  transform(AA_Rate = Rate\/Rate[1]) |>   # relative to 1999\r\n  plot(AA_Rate ~ Year, data = _, type = "l", \r\n       ylab = "age-adjusted death rate, relative to 1999")<\/code><\/pre>\n
<\/p>\n<\/div>\n
\n
Fig 2.12 (b)<\/h2>\n
In the book, this plot shows two different age adjustments, even thought the exposition says there are three. They\u2019re probably referring to the original blog post plot which does show three. I recreate the plot in the blog post, which is the second to last plot.<\/p>\n
This plot shows (1) age-adjustment using the simple mean of rates, (i.e., the plot above), (2) age-adjustment using the distribution of ages in 1999, and (3) age-adjustment using the distribution of ages in 2013.<\/p>\n
This plot requires the most work of all. First we need to get the total population for all ages in 1999 and 2013. These are used to make the age adjustments.<\/p>\n
pop1999 <- aggregate(Population ~ Age, data = data, \r\n                     subset = Year == 1999 & Age %in% 45:54, sum)[["Population"]]\r\npop2013 <- aggregate(Population ~ Age, data = data, \r\n                     subset = Year == 2013 & Age %in% 45:54, sum)[["Population"]]<\/code><\/pre>\n
Next we calculate age-adjusted rates using the population distributions from 1999 and 2013. Again we sum Deaths and Populations by age and year for the 45-54 age group and calculate the death rate. Then we calculate the average rate by year using the population distributions to calculate a weighted mean.<\/p>\n
# age-adjustment from Fig 2.12 (a)\r\naa_rate_uniform <- aggregate(cbind(Deaths, Population) ~ Age + Year, \r\n                             data = data, sum, \r\n                             subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  aggregate(Rate ~ Year, data = _, mean) |> \r\n  transform(AA_Rate = Rate\/Rate[1])\r\n\r\naa_rate_1999 <- aggregate(cbind(Deaths, Population) ~ Age + Year, \r\n                          data = data, sum,\r\n                          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  aggregate(Rate ~ Year, data = _, function(x)weighted.mean(x, pop1999))\r\n\r\naa_rate_2013 <- aggregate(cbind(Deaths, Population) ~ Age + Year, \r\n                          data = data, sum, \r\n                          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  aggregate(Rate ~ Year, data = _, function(x)weighted.mean(x, pop2013))<\/code><\/pre>\n
Now we can make the plot. Notice we find the range of all the data to help set the limits of the y-axis. Also notice we rescale the plot so all lines begin at 1.<\/p>\n
rng <- range(aa_rate_uniform$Rate\/aa_rate_uniform$Rate[1], \r\n             aa_rate_1999$Rate\/aa_rate_1999$Rate[1], \r\n             aa_rate_2013$Rate\/aa_rate_2013$Rate[1])\r\nplot(years, aa_rate_uniform$Rate\/aa_rate_uniform$Rate[1], type = "l", ylim=rng,\r\n     ylab = "age-adjusted death rate, relative to 1999")\r\nlines(years, aa_rate_1999$Rate\/aa_rate_1999$Rate[1], lty=2)\r\nlines(years, aa_rate_2013$Rate\/aa_rate_2013$Rate[1], lty=3)\r\ntext(2003, 1.053, "Using 1999\\nage dist", cex=.8)\r\ntext(2004, 1.032, "Using 2013\\nage dist", cex=.8)<\/code><\/pre>\n
<\/p>\n
The point of this plot is to demonstrate it doesn\u2019t matter how the age-adjustment is done.<\/p>\n<\/div>\n
\n
Fig 2.12 (c)<\/h2>\n
The final plot shows age adjusted death rates broken down by sex. This is basically the same code as Fig 2.12 (a) but with male included in the calls to aggregate()<\/code>. To rescale the y-axis so it starts at 1 we need divide each vector of rates by the respective 1999 value.<\/p>\n
aa_rate_sex <- aggregate(cbind(Deaths, Population) ~ Age + Year + Male, \r\n                         data = data, sum, \r\n          subset = Age %in% 45:54) |> \r\n  transform(Rate = Deaths\/Population) |> \r\n  aggregate(Rate ~ Year + Male, data = _, mean) \r\n\r\nplot(years, aa_rate_sex$Rate[aa_rate_sex$Male == 0]\/\r\n       aa_rate_sex$Rate[aa_rate_sex$Year == 1999 & aa_rate_sex$Male == 0], \r\n     col="red", type = "l",\r\n     ylab = "Death rate relative to 1999")\r\nlines(years, aa_rate_sex$Rate[aa_rate_sex$Male == 1]\/\r\n        aa_rate_sex$Rate[aa_rate_sex$Year == 1999 & aa_rate_sex$Male == 1], \r\n      col="blue", type = "l")\r\ntext(2011.5, 1.075, "Women", col="red")\r\ntext(2010.5, 1.02, "Men", col="blue")<\/code><\/pre>\n
<\/p>\n
Gelman called his code \u201cugly\u201d, but it\u2019s his code and he understands it. I don\u2019t claim my code is any better, but it\u2019s my code and I understand it.<\/p>\n","protected":false},"excerpt":{"rendered":"
In section 2.4 of Regression and Other Stories, Gelman, et al.\u00a0explain the necessity of age adjustment when investigating mortality rates…. Read more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[13],"tags":[81,82],"class_list":["post-926","post","type-post","status-publish","format-standard","hentry","category-using-r","tag-age-adjustment","tag-mortality-rates"],"_links":{"self":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/926","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/comments?post=926"}],"version-history":[{"count":4,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/926\/revisions"}],"predecessor-version":[{"id":930,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/posts\/926\/revisions\/930"}],"wp:attachment":[{"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/media?parent=926"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/categories?post=926"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.clayford.net\/statistics\/wp-json\/wp\/v2\/tags?post=926"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}