rm(list=ls()) ##################################################### # NOT NECESSARY FOR REPRODUCTION OF OUR PAPER # ##################################################### # Step 13 of the code by Schmertmann & Rau on life expectancy estimates # on the Kreis-level. # Purpose: Create the plot for the website where it is explained how # the TOPALS approach works in principle in our case. # License: GPL-2, https://www.gnu.org/licenses/old-licenses/gpl-2.0.html ##################################################### library(tidyverse) library(cowplot) # Purpose of file: load("../data/dths.expos2015_2017.Rdata") load("../initialdata/geo.RData") # Vorpommern Rügen males VPRM = dths.expos2015_2017 %>% filter(AGS == 13073, Sex==1, Age < 85) %>% mutate(mx = Deaths/Exposures, logmx = log(mx)) # posterior means (from full Stan fit on WiWi server) alpha = c(-0.094, -0.101, 0.063, -0.007, 0.337, 0.107, -0.151) # male std (from full Stan fit on WiWi server) male_std = c(-5.594, -8.119, -8.664, -8.942, -9.052, -9.093, -9.158, -9.27, -9.393, -9.492, -9.53, -9.472, -9.307, -9.066, -8.781, -8.488, -8.219, -8.001, -7.833, -7.709, -7.622, -7.566, -7.534, -7.519, -7.513, -7.51, -7.501, -7.48, -7.445, -7.399, -7.345, -7.289, -7.234, -7.181, -7.129, -7.074, -7.013, -6.944, -6.864, -6.774, -6.676, -6.572, -6.464, -6.354, -6.243, -6.13, -6.016, -5.9, -5.782, -5.664, -5.545, -5.429, -5.315, -5.206, -5.102, -5.003, -4.909, -4.819, -4.732, -4.648, -4.567, -4.487, -4.409, -4.331, -4.253, -4.175, -4.097, -4.018, -3.938, -3.858, -3.776, -3.692, -3.604, -3.51, -3.409, -3.301, -3.187, -3.069, -2.948, -2.828, -2.71, -2.595, -2.481, -2.368, -2.255, -2.142, -2.028, -1.913, -1.799, -1.684) # quantiles of logmx for 0..89, males in Vorpommern-Rügen VPRMfit = structure(c(-5.7753, -8.3115, -8.8305, -9.0841, -9.1722, -9.1915, -9.236, -9.3335, -9.4445, -9.5315, -9.5596, -9.498, -9.334, -9.0928, -8.811, -8.5229, -8.261, -8.0513, -7.8938, -7.7811, -7.7098, -7.632, -7.5786, -7.5424, -7.5164, -7.4917, -7.4617, -7.4211, -7.3652, -7.2994, -7.2271, -7.1516, -7.0777, -7.0052, -6.9351, -6.862, -6.7839, -6.6976, -6.6018, -6.4964, -6.3834, -6.2857, -6.1835, -6.0797, -5.9744, -5.8677, -5.7593, -5.6496, -5.5383, -5.4259, -5.3141, -5.2037, -5.0964, -4.9935, -4.8963, -4.8041, -4.7165, -4.6328, -4.5529, -4.4765, -4.4025, -4.3305, -4.2595, -4.1898, -4.1202, -4.0507, -3.9808, -3.9104, -3.8396, -3.7681, -3.695, -3.6225, -3.5457, -3.4637, -3.3755, -3.2797, -3.1783, -3.0732, -2.9669, -2.861, -2.7573, -2.6565, -2.5578, -2.4602, -2.3632, -2.2658, -2.1679, -2.0693, -1.9699, -1.8714, -5.6871, -8.2201, -8.7476, -9.0068, -9.0988, -9.1213, -9.1687, -9.2622, -9.3666, -9.448, -9.4681, -9.4162, -9.2587, -9.0235, -8.7462, -8.4589, -8.1981, -7.9869, -7.8269, -7.7095, -7.63, -7.5573, -7.5077, -7.4758, -7.4527, -7.4318, -7.4063, -7.3677, -7.315, -7.2509, -7.1808, -7.1074, -7.0347, -6.9645, -6.8941, -6.8217, -6.7447, -6.6579, -6.5605, -6.4537, -6.3383, -6.2423, -6.1421, -6.0398, -5.9361, -5.831, -5.7244, -5.6158, -5.5063, -5.3957, -5.2852, -5.1764, -5.0707, -4.969, -4.8728, -4.7817, -4.6951, -4.6124, -4.5334, -4.457, -4.3831, -4.3108, -4.24, -4.17, -4.1001, -4.0299, -3.9594, -3.8882, -3.8162, -3.7437, -3.6699, -3.5989, -3.5241, -3.4442, -3.357, -3.2624, -3.1616, -3.0569, -2.9504, -2.844, -2.7395, -2.6375, -2.5372, -2.4375, -2.3385, -2.2386, -2.1383, -2.0377, -1.9368, -1.8358, -5.5995, -8.1252, -8.662, -8.9295, -9.0251, -9.0509, -9.0962, -9.1871, -9.2872, -9.3614, -9.3735, -9.3317, -9.1829, -8.9553, -8.6836, -8.399, -8.1371, -7.9239, -7.7586, -7.6344, -7.5466, -7.4772, -7.4322, -7.4035, -7.3846, -7.3684, -7.3464, -7.3117, -7.2622, -7.2023, -7.1346, -7.064, -6.9934, -6.9253, -6.8571, -6.7849, -6.7075, -6.6197, -6.5211, -6.4131, -6.2959, -6.2014, -6.1029, -6.0019, -5.8996, -5.7958, -5.6903, -5.5838, -5.4753, -5.3661, -5.2568, -5.1491, -5.0444, -4.9439, -4.8491, -4.7593, -4.6732, -4.5917, -4.5133, -4.4377, -4.3641, -4.2921, -4.221, -4.1503, -4.0799, -4.0087, -3.9372, -3.8651, -3.7925, -3.7185, -3.6434, -3.5751, -3.5019, -3.4238, -3.3383, -3.2449, -3.1447, -3.0404, -2.9336, -2.8265, -2.7209, -2.6176, -2.5154, -2.4146, -2.3137, -2.212, -2.1101, -2.0069, -1.9035, -1.8), .Dim = c(90L, 3L), .Dimnames = list( 0:89, c("Q10", "Q50", "Q90") )) B = structure(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.889, 0.778, 0.667, 0.556, 0.444, 0.333, 0.222, 0.111, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.111, 0.222, 0.333, 0.444, 0.556, 0.667, 0.778, 0.889, 1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 0.967, 0.933, 0.9, 0.867, 0.833, 0.8, 0.767, 0.733, 0.7, 0.667, 0.633, 0.6, 0.567, 0.533, 0.5, 0.467, 0.433, 0.4, 0.367, 0.333, 0.3, 0.267, 0.233, 0.2, 0.167, 0.133, 0.1, 0.067, 0.033, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.033, 0.067, 0.1, 0.133, 0.167, 0.2, 0.233, 0.267, 0.3, 0.333, 0.367, 0.4, 0.433, 0.467, 0.5, 0.533, 0.567, 0.6, 0.633, 0.667, 0.7, 0.733, 0.767, 0.8, 0.833, 0.867, 0.9, 0.933, 0.967, 1, 0.947, 0.895, 0.842, 0.789, 0.737, 0.684, 0.632, 0.579, 0.526, 0.474, 0.421, 0.368, 0.316, 0.263, 0.211, 0.158, 0.105, 0.053, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.053, 0.105, 0.158, 0.211, 0.263, 0.316, 0.368, 0.421, 0.474, 0.526, 0.579, 0.632, 0.684, 0.737, 0.789, 0.842, 0.895, 0.947, 1), .Dim = c(90L, 7L)) ############################################ B.df = as.data.frame(B) %>% mutate(age = 0:89) %>% gather( key='column', value='val', -age) G1 = ggplot( data=B.df) + aes(x=age, y=val, group=column) + geom_line() + geom_line(data=filter(B.df,column=='V5'), lwd=2) + labs(x='',y='') + scale_y_continuous(breaks=seq(0,1,.50)) + scale_x_continuous(breaks=seq(0,90,10)) + theme_bw() + theme(panel.grid = element_blank(), axis.text = element_text(face='bold', size=11)) a.df = data.frame(age=c(0,1,10,20,40,70,90), alpha = alpha) G2 = ggplot(data=a.df) + aes(x=age, y=alpha) + geom_point(shape=16, size=2) + geom_line() + geom_hline(yintercept = 0, lty='dotted') + scale_y_continuous(limits=c(-0.4,+0.4)) + scale_x_continuous(breaks=seq(0,90,10)) + geom_segment(aes(x=age,y=0,xend=age,yend=alpha), lwd=0.25)+ labs(x='') + theme_bw() + theme(panel.grid = element_blank(), axis.text = element_text(face='bold', size=11)) TOPALS.df = data.frame(age=0:89, VPRMfit, male_std) %>% mutate(m10 = exp(Q10), m50 = exp(Q50), m90 = exp(Q90)) G3 = ggplot(data=TOPALS.df) + aes(x=age, y=Q50) + geom_line(lwd=1.5, color='blue', alpha=0.8) + geom_line(aes(x=age,y=male_std), lwd=1.5, lty='dotted') + geom_ribbon(aes(x=age, ymin=Q10, ymax=Q90), fill='blue', color=NA, alpha=.30) + labs(x='Age', y='Ln(Mortality Rate)') + scale_y_continuous(breaks=seq(-8,0,2), labels=function(x) {sprintf('%.f',x)}) + scale_x_continuous(breaks=seq(0,90,10)) + theme_bw() + theme(axis.text = element_text(face='bold', size=11)) boundaries = c(0,1,seq(5,85,5)) VPRM$L = head(boundaries,-1) # lower end of age groups VPRM$U = tail(boundaries,-1) # upper end of age groups G3 = G3 + geom_segment( data=VPRM, aes(x=L,y=logmx,xend=U,yend=logmx), lwd=1.5, color='blue') + geom_text(data=VPRM, aes(x=(L+U)/2, y=logmx, label=ifelse(Deaths>0,Deaths,'')), nudge_y=0.35,color='blue', size=5) + geom_text(aes(x=7.5, y=-9.6, label='0'), color='blue', size=5, inherit.aes = FALSE) ################ graphics.off() pdf("../plots/TOPALS-Explanation.pdf", height=12, width=12, pointsize=16) plot_grid( G1,G2,G3, ncol=1, rel_heights = c(1,1,3.5), labels=LETTERS[1:3]) dev.off()