rm(list=ls()) library(tidyverse) ##################################################### # Step 12 of the code by Schmertmann & Rau on life expectancy estimates # on the Kreis-level. # Purpose: Plotting of the impact of the change by one standard # deviation on life expectancy for each of the indicators, separately # for women and men as well as for east and west. # License: GPL-2, https://www.gnu.org/licenses/old-licenses/gpl-2.0.html ##################################################### # CHOOSE THE LANGUAGE FOR THE PLOT HERE: # must be either 'German' or 'English' make.figure.4 <- function(language = 'German') { if (!(language %in% c("English", "German"))) { stop("Argument 'language' has to be either 'English' or 'German'") } #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # English terms here English_pop_density <- 'Population\nDensity' English_doctors <- 'Doctors/Person' English_gdp <- 'GDP/Person' English_unemp <- 'Unemployment' English_child_pov <- 'Child\nPoverty' English_housing <- 'Housing\nSubsidy' English_assistance <- 'Public\nAssistance' English_elder_supp <- 'Elder\nSupport' English_x_axis <- 'Expected Change in Life Expectancy (years)\nwhen indicator increases by 1 std deviation' English_sex <- c('F'='FEMALE','M'='MALE') English_EW <- c('E'='E','W'='W') #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # German terms here (only edit the text in quotes) German_pop_density <- 'Bevölkerungsdichte' German_doctors <- 'Ärzte/Person' German_gdp <- 'BIP/Kopf' German_unemp <- 'Arbeitslosigkeit' German_child_pov <- 'Kinderarmut' German_housing <- 'Wohngeld' German_assistance <- 'SGB II / \"Hartz IV\"' German_elder_supp <- 'Grundsicherung im Alter' German_x_axis <- 'Erwartete Veränderung in der Lebenswartung (in Jahren)\n bei Anstieg des Indikators um 1 Standardabweichung' German_sex <- c('F'='Weiblich','M'='Männlich') #<< keep F and M as element names German_EW <- c('E'='O', 'W'='W') #<< keep E and W as element names #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ prev_data <- read.table("../initialdata/inkar-2019-09.csv", header=FALSE, sep=";", skip=1, encoding='latin1') names(prev_data) <- c("AGS", "County", "ALQ", "ALQWOM", "ALQMEN", "PopDensity", "DoctorsPer100000", "NoSchoolWom", "NoSchoolMen", "BIP", "SGBII", "GrundsicherungAlter") ### all data for 2016! ## DoctorsPer100000: Allgemeinärzte per 100,000 Einwohner (General ## physicians per 100,000 inhabitants ## BIP: GDP per inhabitants in 1000 Euros ## SGB II: Anteil der erwerbsfähigen und nicht erwerbsfähigen ## Leistungsberechtigten nach SGB II an den unter 65- ## jährigen Einwohnern in % (-> (long-term) unemployment benfits & ## social welfare) "Hartz IV" => https://en.wikipedia.org/wiki/Hartz_concept#Hartz_IV ## Grundsicherung: Anteil der Bevölkerung mit Grundsicherung im Alter ## an den Einwohnern 65 Jahre und älter in ‰ per 1000!!! ## => SGB XII => https://de.wikipedia.org/wiki/Grundsicherung_im_Alter_und_bei_Erwerbsminderung # more indicators (downloaded 14 Jan 2020) if (.Platform$OS.type=="unix") { ## I (RR) had problems on Linux ## Error in make.names(col.names, unique = TRUE) : ## invalid multibyte string 8 inkar_data = read.csv('../initialdata/inkar-Tuesday-table-14-Jan.csv', header=FALSE, skip=1, stringsAsFactors = FALSE, encoding = 'latin1') the.names <- read.csv('../initialdata/inkar-Tuesday-table-14-Jan.csv', header=FALSE, nrow=1, stringsAsFactors = FALSE, encoding='latin1') names(inkar_data) <- as.character(the.names) names(inkar_data)[names(inkar_data)=="Kennziffer"] <- 'AGS' } else { inkar_data = read.csv('../initialdata/inkar-Tuesday-table-14-Jan.csv', header=TRUE, stringsAsFactors = FALSE, encoding = 'latin1') %>% rename(AGS=Kennziffer) } load("../data/e0-of-our-model.RData") tmp = dplyr::select(e0.our.model, AGS, Sex2, e0=pct50) %>% spread( key='Sex2', value='e0') %>% rename(e0M = M, e0F='F') head(tmp) inkar_data = left_join(inkar_data, tmp, by='AGS') %>% mutate( loc = factor( case_when( AGS < 11000 ~ "West", AGS == 11000 ~ "Berlin", AGS > 11000 ~ "East" ))) %>% left_join(prev_data) %>% pivot_longer(cols=contains('e0'), values_to = 'e0', names_to = 'sex') %>% transform(sex=c('F','M')) %>% dplyr::select(-Raumeinheit, -Aggregat) # names of indicator variables vnames = setdiff( names(inkar_data), c('AGS','sex','e0','loc','County')) # rho 1 contains means of indicators rho1 = inkar_data %>% group_by(sex,loc) %>% summarize_at(vnames, ~mean(.)) %>% gather(key='var', value='mean', -sex, -loc) # rho 2 contains correlations of indicators with e0 rho2 = inkar_data %>% group_by(sex,loc) %>% summarize_at(vnames, ~cor(.,e0, use='pairwise.complete.obs')) %>% gather(key='var', value='correl_e0', -sex, -loc) # rho 3 contains standardized regression coefs for # indicators e0 = a + b * [(indicator - mean) / sd ] b = rep(NA, nrow(rho1)) for (i in seq(b)) { this_sex = rho1$sex[i] this_loc = rho1$loc[i] this_var = rho1$var[i] tmp = filter(inkar_data, sex==this_sex, loc==this_loc) xx = pull(tmp, this_var) if (this_loc =='Berlin') { beta=NA} else { xx = (xx - mean(xx,na.rm=TRUE)) / sd(xx,na.rm=TRUE) yy = pull(tmp, e0) beta = coef(lm(yy~xx))[2] # delta = beta * diff( quantile(xx,c(.05,.95),na.rm=TRUE)) } b[i] = beta } rho3 = rho1 %>% ungroup() %>% dplyr::select(sex:var) %>% mutate(beta=b) rho= full_join(rho1,rho2) %>% full_join(rho3) %>% ungroup() ## bar plot of effects for selected indicators ## selected vars in desired order for plotting # assign English or German labels for (var in c('pop_density','doctors','gdp','unemp', 'child_pov','housing','assistance', 'elder_supp','x_axis','sex','EW')) { assign(x = paste0('label_',var), value = get(paste0(language,'_',var))) } #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## bar plot of effects for selected indicators ## selected vars in desired order for plotting selected = tribble( ~var, ~abb, ~ord, 'PopDensity', label_pop_density, 1, 'DoctorsPer100000', label_doctors, 2, 'BIP', label_gdp, 3, 'ALQ', label_unemp,4, 'Kinderarmut', label_child_pov,5, 'Wohngeldhaushalte', label_housing, 6, 'SGBII', label_assistance,7, 'GrundsicherungAlter', label_elder_supp, 8 ) mydata = rho %>% filter(var %in% selected$var, loc != 'Berlin') %>% left_join(selected) %>% mutate(var = factor(abb, levels=rev(selected$abb)), sexlabel = label_sex[sex], EW = label_EW[substr(loc,1,1)] ) ggplot(data=mydata) + aes(x=var, y=beta, label=EW, fill=loc) + geom_bar(position = 'dodge', stat='identity', size=0.2, color='black') + scale_y_continuous(breaks=seq(-.50,.50,.50), limits=c(-.65,+.65)) + facet_grid(~sexlabel) + guides(fill=FALSE) + geom_text(aes(x=var, y=beta + sign(beta)*.04), position = position_dodge(width = 1)) + theme_bw() + labs(fill='Location', x='',y=label_x_axis) + scale_fill_manual(values=c('white','red')) + theme(axis.text.x = element_text(face='bold', size=14,vjust=0.5), axis.text.y = element_text(face='bold', size=14,vjust=0.5), axis.title.x = element_text(face='bold', size=14,vjust=0.5), strip.text = element_text(face='bold',size=16), panel.spacing.x=unit(0, "lines"), legend.text = element_text(size=12)) + coord_flip() filename = ifelse(language=='German', '../plots/deFigure4.pdf', '../plots/enFigure4.pdf') graphics.off() ggsave(file=filename, width=15, height=8.5, units='in') filename.ps = ifelse(language=='German', '../plots/deFigure4.eps', '../plots/enFigure4.eps') graphics.off() ggsave(file=filename.ps, width=15, height=8.5, units='in') } make.figure.4('German') make.figure.4('English')