<jupyter_start><jupyter_text>SOEP Mortality Results for the Paper "Social differences in mortality and life expectancy. Current situation and trends"  Submitted for Journal of Health Monitoring with an aimed publication date of March 2019**Authors:**<NAME> ([@l_kroll](https://twitter.com/l_kroll)), <NAME>, <NAME>**Date:**17.10.2018**Version:** 1.0  Setup<jupyter_code># Path
setwd("~/wissdaten/Abt2_Computation/P3_Mortality")
#setwd("~/Desktop/P3_Mortality")

# SCRATCH_DIR
scratch <- "~/SCRATCH_NOBAK/"
# scratch <- "data/" # home

# Log
sink("Logged_Output.txt", append =FALSE)

# Stop log
sink()
# Attach libraries
options(tidyverse.quiet = TRUE)
library("tidyverse", warn.conflicts=FALSE, quiet=TRUE)
library("broom", warn.conflicts=FALSE, quiet=TRUE)
library("stringr", warn.conflicts=FALSE, quiet=TRUE)
library("survival", warn.conflicts=FALSE, quiet=TRUE)
library("rkicolors", warn.conflicts=FALSE, quiet=TRUE) 
library("wec", warn.conflicts=FALSE, quiet=TRUE)
library("ggfortify", warn.conflicts=FALSE, quiet=TRUE)
library("haven", warn.conflicts=FALSE, quiet=TRUE)
suppressPackageStartupMessages(library("survey", warn.conflicts=FALSE, quiet=TRUE))
library("srvyr", warn.conflicts=FALSE, quiet=TRUE)
library("foreach", warn.conflicts=FALSE, quiet=TRUE)
library("doMC", warn.conflicts=FALSE, quiet=TRUE)
library("openxlsx", warn.conflicts=FALSE, quiet=TRUE)

# Jupyter Set-up
options(repr.plot.width=6, repr.plot.height=(6/16)*9 , repr.plot.dpi=300)
options(jupyter.plot_mimetypes = c("text/plain", "image/png" ))<jupyter_output><empty_output><jupyter_text>Data ManagementWe use data from GSOEP which can be obtained [here](https://www.diw.de/en/diw_02.c.222518.en/research_data_center_of_the_soep.html) as well as annual death rates based on the official German Life Tables which can be obtained from the [Federal Statistical Office for Germany](www.destatis.de).  SOEP  Load SOEP DataWe used the long format of SOEP and performed custom computations and matching of multiple long format files. The Methods of the data matching and creation of our working dataset called *"soep_arbeitsdaten.dta"* cannot be described in detail here.<jupyter_code># Log to File
sink("Logged_Output.txt", append =TRUE)

# Load recoded SOEP data
dir.create(paste0(scratch,"v33"), showWarnings =FALSE)
invisible(file.copy("data/v33/soep_trend_full.dta",paste0(scratch,"v33/soep_trend_full.dta"),overwrite=F))
invisible(file.copy("data/v33/soep_arbeitsdaten.dta",paste0(scratch,"v33/soep_arbeitsdaten.dta"),overwrite=F))

soep <- read_dta(paste0(scratch,"v33/soep_arbeitsdaten.dta"))

soep.rec <- soep %>% ungroup() %>%
    select(pid,hid,cid,sex,age,year=syear,
           todjahr,todinfo,contains("hrf"),
           SDeqinc,SDincpos,SDcasmin) %>%
mutate(
    gebjahr=year-age,
    SDincpos_grp = ifelse(SDincpos<0.60, "<60%",NA),
    SDincpos_grp = ifelse(SDincpos>=0.60 & SDincpos<0.80," 60%-79%",SDincpos_grp),
    SDincpos_grp = ifelse(SDincpos>=0.80 & SDincpos<1.00, " 80%-99%",SDincpos_grp),
    SDincpos_grp = ifelse(SDincpos>=1.00 & SDincpos<1.50, "100%-149%",SDincpos_grp),
    SDincpos_grp = ifelse(SDincpos>=1.50, "150% u.m.",SDincpos_grp),
    SDincpos_grp = factor(SDincpos_grp, levels=c("<60%"," 60%-79%"," 80%-99%","100%-149%","150% u.m.")),
    sex=ifelse(sex==1,"M","W"),
    SDcasmin=as_factor(SDcasmin)) %>%
rename("incpos"=SDincpos_grp)

andata <- soep.rec %>% filter(year>=1992) %>%
    group_by(pid) %>% arrange(pid,year) %>%
    mutate(
        SurvJahr=ifelse(is.na(todjahr),max(year),todjahr),
        SurvTime=ifelse(is.na(todjahr),max(age),todjahr-gebjahr),
        death=ifelse(is.na(todjahr),0,1)) %>% ungroup()
sink()<jupyter_output><empty_output><jupyter_text>Show Properties of the SOEP Dataset for mortality analysisThe following stept analyzes the SOEP Dataset regarding the available Information on Mortality. We show numbers of deaths, persons and episodes and investigate the association of "Fate Studies" for Nonresponders on the number of recorded deaths per year in reference to the sample size.<jupyter_code>andata %>% summarise(Episodes=n())
as_tibble(andata) %>% group_by(pid) %>% filter(row_number()==1) %>% ungroup()  %>% summarise("Respondents"=n())
as_tibble(andata) %>% group_by(pid) %>% filter(row_number()==1) %>% ungroup() %>% filter(!is.na(todjahr)) %>% as_tibble() %>% summarise("Deaths"=n()) 

plotdata <- left_join(andata %>% filter(tsphrf>0) %>% mutate( year = as.numeric(year)) %>% group_by(year) %>% summarise(Respondents=n(), Mean_Age =mean(age, na.rm=T)) ,
andata %>% group_by(pid) %>% filter(row_number()==1)  %>% mutate(todinfo=as_factor(todinfo))  %>% 
                      filter(!is.na(todjahr)) %>% group_by(todjahr) %>% count() %>% ungroup() %>% 
                      rename(year=todjahr,Deaths=n) %>%
                      mutate(Deaths=ifelse(is.na(Deaths),0,Deaths),
                             year = as.numeric(year)), by="year")

# Plot Respondents
plot.soep.respondents <- ggplot(plotdata,aes(x=year,y=Respondents, fill="x")) + geom_bar(stat="identity") + theme_rki() +
scale_fill_rki("qualitative") + labs(fill="",y="",caption="SOEP Participants",x="") + guides(fill=FALSE) +
theme(legend.position="bottom") 

# Plot Deaths
plot.soep.deaths <- andata %>% group_by(pid) %>% filter(row_number()==1) %>% mutate(info= case_when(todinfo %in% c(5,6,7) ~ "Special fate study",
                                            todinfo %in% c(1,2,3,4,8) ~ "Ongoing Fieldwork"))  %>% filter(!is.na(todjahr)) %>% group_by(todjahr,info) %>% count() %>% ungroup() %>%
ggplot(.,aes(x=todjahr,y=n, fill=info)) + 
geom_bar(stat="identity") + theme_rki() + scale_fill_rki("qualitative") + labs(fill="",y="",caption="Sources of Death Cases",x="") +
theme(legend.position="bottom") 

# Plot Comparison with and without fate studies
plot.soep.deathrates <- ggplot(plotdata %>% filter(!is.na(Deaths)), aes(y=Deaths,x=Respondents/1000)) + geom_smooth(method="gam") + geom_point() + 
facet_grid(. ~ ifelse(year<2009,"Before 2009","2009 and later")) + theme_rki() + scale_fill_rki() + scale_color_rki() +
labs(x="Respondents (1k)", y="Death Cases") + scale_x_continuous(limits=c(10,40))

# Write Info to Excel
write.xlsx(plotdata, file = "Results_Mortality.xlsx", sheetName="SOEP")

# Write Plots to Excel
wb <- loadWorkbook("Results_Mortality.xlsx")
plot.soep.respondents
insertPlot(wb, "SOEP", width = 7.11, height = 4, startRow = 1,  startCol = 6  , fileType = "png", units = "in", dpi = 400)
plot.soep.deaths
insertPlot(wb, "SOEP", width = 7.11, height = 4, startRow = 1,  startCol = 18  , fileType = "png", units = "in", dpi = 400)
plot.soep.deathrates
insertPlot(wb, "SOEP", width = 7.11, height = 4, startRow = 22,  startCol = 6  , fileType = "png", units = "in", dpi = 400)
saveWorkbook(wb, "Results_Mortality.xlsx" , overwrite=TRUE)<jupyter_output><empty_output><jupyter_text>ConclusionEstimates for deaths after 2009 are more unreliable, due to a large number of censored cases with unknown whereabouts as a consequence of missing fate/whereabout data on nonrespondents.  Life Tables Load Data<jupyter_code># load data
load("data/destatis/life_tables_brd.RData")
#lifetables.metadata
#head(lifetables.data)

# recode data
lifetables.rec <- lifetables.data %>% filter(DINSG=="DG") %>% 
    select(sex=GES,year=JAHR,age=ALT577,qx=LEB001_val,LE_orig=LEB007_val) %>%
    mutate(sex=as.character(str_sub(sex, 4, -1)),
           age=as.numeric(str_sub(age, 8, -1)))  %>%
           group_by(sex,year) %>% mutate(qx=ifelse(age==last(age) & qx==last(qx) & last(qx)==0,NA,qx))
# Give all Years the same age range
lifetables.rec <- lifetables.rec %>% group_by(sex,year) %>% arrange(sex,year,age) %>%
mutate(qx=ifelse(age==last(age) & qx==0, NA,qx),lastage=max(age)) %>%
left_join(expand.grid("sex"=unique(lifetables.rec$sex),"year"=unique(lifetables.rec$year),"age"=seq(0:112)) %>% arrange(year,sex,age) %>% mutate(sex=as.character(sex)),.
          , by=c("sex","year","age")) %>%  group_by(sex,year) %>% 
mutate(lastage=max(lastage,na.rm=T),qx=ifelse(is.na(qx),max(qx, na.rm=T)*1.06^(age-lastage-1),qx)) %>% select(-lastage) %>% mutate(age=age-1)<jupyter_output><empty_output><jupyter_text>Function to calculate Life Expectancy<jupyter_code># Function to calculate life expactancies from death rates.
life_expectancy <- function(data,deathrate="qx",age="age",indicator="ex_t"){
   tempdata <- data %>%  rename("deathrate"=deathrate,"age"=age)  %>% group_by(age, add=T) %>% summarise(deathrate=mean(deathrate))
   group_vars <- setdiff(names(tempdata) ,c("age","deathrate")) 
   tempdata %>% group_by_at(group_vars) %>%
   mutate(S_t = ifelse(age==0,1,NA),
           S_t = ifelse(age>0,1-lag(deathrate),S_t),
           S_t = cumprod(S_t),
           H_t = 1-S_t,                 
           ex_t= cumsum(S_t),
           ex_t= last(ex_t)-ex_t+.45) %>% select(-deathrate)   %>%     
    select_(group_vars,age,indicator) %>% left_join(data,.)
}

le.comparison.methods <- left_join(lifetables.rec,lifetables.rec %>% group_by(sex,year) %>% life_expectancy(.)) %>% 
    rename(LE_function = ex_t) %>% filter(age==0) %>% 
    mutate(Difference=LE_orig-LE_function)
le.comparison.methods %>% group_by(sex) %>% summarise(Mean_Diff = mean(Difference))
ggplot(le.comparison.methods ,aes(x=year,color=sex,y=Difference)) + geom_smooth() + theme_rki() + scale_color_rki() + scale_y_continuous(limits=c(-1,1)) + labs(y="Diff. of approx. fcn and off. values")<jupyter_output>Adding missing grouping variables: `year`
Joining, by = c("sex", "year", "age")
Joining, by = c("sex", "year", "age", "qx", "LE_orig")<jupyter_text>The Function to compute life expectancies based on age specific death rates (`life_expectancy()`) produces a high accuracy when compared to the calculations that were provided by DESTATIS (via Indcicator `LEB006`). We conclude that it can be used to approximate group specific life expectancies.  AnalysisThe analytical strategy was first outlined by <NAME> (2009) [PDF](https://www.gesis.org/fileadmin/upload/forschung/publikationen/zeitschriften/mda/Vol.3_Heft_1/02_Kroll.pdf) and is slightly modified for the current analysis. We will analyse the data in the following steps: 1. Analyse the Inequalities in Mortality using an Indicator for the current Income Position in SOEP for two age groups (0-50,51-101) using Cox Proportional Hazard Regression models (function `svycoxph` from package `survey`).2. Apply the resulting Hazard Ratios as well as their upper and lower confidence bands in effect coding ([cf. source](https://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html)) to death rates from official life tables.3. Calculate life expectancies based on calculated rates including CI  FunctionsIn this section, we will define several functions, that will be used to  1. Create Groups for Analysis,  2. Merge Effects and Lifetable Data  3. Filter data based on group definitions for Survival Analysis 4. Predict Life Expectancies for the results of the regression models 5. Get Number of Cases and Observations depending on survival Status 6. Perform regressions using Effect Coding and provide Estimates in Reference to the Grand Mean for all Levels 7. Combine all Subfunctions 8. Compute some metrics for the results (gradients mortality risk as well as in life expectancy)<jupyter_code># Create Groups for Analysis
create_angroups <- function (data,seqyear=3,seqage=51)
{
# Sex
sex <-  c("M","W")
# Time
startyear <-min(data$year)
stopyear <- max(data$year)-1 # tbd dynamically
# Age
startage <- 0
stopage <- 100 
# Define Groups
startyear <-seq(startyear,stopyear,seqyear)
startage<-seq(startage,stopage,seqage)
angroups <- expand.grid(sex=sex,begin=startyear,youngest=startage) %>%   
arrange(sex,begin,youngest)  %>% 
mutate(
  group_number=1:n(),  
  sex=as.character(sex)  ,
  end=begin+seqyear-1,
  oldest=youngest+seqage-1) %>%
  select(group_number,sex,begin,end,youngest,oldest) 
print(paste("Number of Groups to analyze:",nrow(angroups)))
angroups
}

# Merge Effects and Lifetable Data
merge_lifetable <- function (myResult,groupdefs) {
    lifetables.rec %>% filter(sex==groupdefs$sex   &
                      year>=groupdefs$begin & 
                      year<=groupdefs$end &
                      age>=groupdefs$youngest &
                      age<=groupdefs$oldest) %>%
    mutate(begin=groupdefs$begin,end=groupdefs$end) %>% group_by(sex,begin,end,age) %>% summarise(qx=mean(qx,na.rm=T)) %>% ungroup() %>% 
    left_join(groupdefs,.) %>%
    left_join(myResult,.,by=names(groupdefs)) %>%
    mutate(qx.base=qx,
           qx.low =qx*conf.low,
           qx.high=qx*conf.high,
           qx     =qx*estimate) 
}

# Predict LE 
predict_lifexp <- function (myResult) {
myResult %>% 
life_expectancy(.,deathrate="qx.base" ) %>% rename("LE_pop"=ex_t) %>% 
life_expectancy(.,deathrate="qx" ) %>% rename("LE_mean"=ex_t) %>% 
life_expectancy(.,deathrate="qx.low" ) %>% rename("LE_max"=ex_t)  %>% 
life_expectancy(.,deathrate="qx.high" ) %>% rename("LE_min"=ex_t) %>%
life_expectancy(.,deathrate="qx.base" , indicator="S_t") %>% rename("Surv_pop"=S_t) %>% 
life_expectancy(.,deathrate="qx" , indicator="S_t") %>% rename("Surv_mean"=S_t) %>% 
life_expectancy(.,deathrate="qx.low" , indicator="S_t" ) %>% rename("Surv_max"=S_t)  %>% 
life_expectancy(.,deathrate="qx.high" , indicator="S_t" ) %>% rename("Surv_min"=S_t) 
    
}

# Filter data for Survival Analysis
filter_soep <- function (data,groupdefs) {
    data %>% 
    filter(sex==groupdefs$sex &
           year>=groupdefs$begin & 
           year<=groupdefs$end & 
           age>=groupdefs$youngest &
           age<=groupdefs$oldest) %>% 
    group_by(pid) %>% mutate(tsphrf=mean(tsphrf,na.rm=T)) %>% filter(row_number()==1) %>% 
    ungroup()    %>% as_survey_design(ids = pid, weight=tsphrf)
}
# Get Number of Cases and Observations depending on survival Status
cases_soep <-  function (data,groupdefs,determinant) {
    left_join(data$variables %>%
              group_by(pid) %>% filter(row_number()==1) %>% ungroup() %>% group_by_(determinant,"death") %>%
              summarise(Number = n()) %>% ungroup() %>% 
              mutate(death=ifelse(death==1,"Cases Death","Cases Censored")) %>% 
              spread(death,Number),
       data$variables %>% group_by_(determinant)  %>% summarise(Episodes = n()), 
        by=c(determinant)) %>% 
    rename("label"=determinant) %>% filter(!is.na(label)) %>% cbind(groupdefs,., row.names = NULL) %>% mutate(label=as.character(label))
}

# Get Hazard Ratios from Survival Model
survival_reg_soep <- function(data,groupdefs,determinant,controls) {
    contrasts(data$variables[[determinant]]) <- "contr.sum"
    Survdef <- "Surv(SurvTime, death, type=\"right\")"
    formula <- as.formula(paste(Survdef, paste(determinant,"+",paste(controls, collapse=" + ")), sep=" ~ ")) 
    data  %>% svycoxph(formula, data=. , design=.) %>%    broom::tidy(exponentiate=TRUE) 
}

# Use Effect Coding and provide Estimates in Reference to the Grand Mean for all Levels
survival_reg_soep_alllevels <- function(data,groupdefs,determinant,controls) {
    # Store original Labels
    myLabels <- cbind(term=determinant,label=levels(data$variables[[determinant]])) %>% 
              as_tibble() %>% 
              mutate(term=paste0(term,1:n()))    
    # Effect coding with standard Refernce Category
    data$variables[[determinant]] <- relevel(data$variables[[determinant]], ref=levels(data$variables[[determinant]])[1])
    myReg <- invisible(survival_reg_soep(data,groupdefs,determinant,controls))
    # Effect coding with alternative Reference Category
    data$variables[[determinant]] <- relevel(data$variables[[determinant]], ref=levels(data$variables[[determinant]])[length(levels(data$variables[[determinant]]))])
    myReg <- bind_rows(myReg,invisible(survival_reg_soep(data,groupdefs,determinant,controls)) %>%  # filter only last level and rename term accordingly
                      filter(term==paste0(determinant,1)) %>% 
                      mutate(term=paste0(determinant,length(levels(data$variables[[determinant]]))))
                      ) %>% arrange(term) 
    # Merge with label
    myReg <-  left_join(myReg,
              myLabels, by="term") %>% arrange(.,term,label) 
myReg
}

# Combine all Subfunctions
survival_soep <- function(data,groupdefs,determinant,controls,byvars) {
    # Group Summary
    myResult <- data.frame(groupdefs)
    # Data
    data <- filter_soep(data,groupdefs)
    # Case Summary
    myResult <- left_join(myResult,cases_soep(data,groupdefs,determinant),by=names(groupdefs))
    # Regression
    myReg <- survival_reg_soep_alllevels(data,groupdefs,determinant,controls)
    # Merge Effects
    myResult <- tryCatch(left_join(myResult,data.frame(myReg)), error=function (e) myResult)
    myResult                    
    }

# Gradient low vs. high
show_gradient <- function(result) {
result %>% select(group_number,sex,begin,end,youngest,oldest,label,estimate,p.value,contains("conf")) %>%
  group_by(group_number) %>%  filter(age==min(age)) %>% mutate_at(c("estimate","conf.low", "conf.high"), funs(./last(.))) %>% 
    mutate_at(c("estimate","conf.low", "conf.high","p.value"), funs(round(., digits = 4))) %>% filter(label=="<60%")
}
                         
# Show Difference low vs. high
show_lediff <- function(result,agefilter=0) {
    result %>% filter(age==agefilter ) %>%  group_by(sex,begin,end) %>% select(sex,begin,end,label,contains("LE_"))  %>%
arrange(sex,begin,label) %>% summarise(Diff_label=paste0(first(label)," vs. ",last(label)), 
                                      Diff = round(last(LE_mean)-first(LE_mean),digits=1),
                                      Diff_min= round(last(LE_min)-first(LE_max),digits=1),
                                      Diff_max= round(last(LE_max)-first(LE_min),digits=1))
}
                         
# Test Functions
#================
#andata$variables  %>% count(year) %>% arrange(-n) %>% head(.)
#andata$variables  %>% count(pid) %>% count(n) %>% head(.)
# filter_soep(andata,angroups[1,]) %>% head(.)
# survival_reg_soep_alllevels(filter_soep(andata,angroups[1,]),angroups[1,],"incpos","age")
# survival_soep(andata,angroups[1,],"incpos","age")
# survival_soep(andata,angroups[1,],"incpos","age") %>% merge_lifetable(.,angroups[1,]) %>% head(10)<jupyter_output><empty_output><jupyter_text>Results We will now perform our analysis using the predefined functions. We decided to provide results for the whole observational period (1992-2016) currently available in [SOEP](https://www.diw.de/SOEPfdz). Additionally we decided to compare risks for two age groups (0-50, 51+). This was due to limited death cases before the age of 51 (see next plot).<jupyter_code>andata %>% filter(!is.na(todjahr)) %>% group_by(pid) %>% filter(row_number()==1) %>%
group_by(sex) %>% count(SurvTime) %>% ungroup() %>% mutate(sex=ifelse(sex=="M","Men","Women")) %>%
ggplot(aes(x=sex, y=SurvTime, fill=sex)) + geom_boxplot() + 
 theme_rki() + scale_fill_rki() + labs(y="Age at Death", x="") + guides(fill=FALSE)

andata %>% filter(!is.na(todjahr)) %>% group_by(pid, year) %>% filter(row_number()==1) %>%
group_by(sex,year) %>% count(SurvTime) %>% ungroup() %>% mutate(sex=ifelse(sex=="M","Men","Women"),
                                                               agegroup=cut(SurvTime, breaks=c(0, 25, 50, 65, Inf), labels=c("0-25","26-50","51-65","65+")),
                                                               year=floor(year/4)*4+2) %>%
ggplot(aes(x=year, fill=sex)) + geom_bar() + 
 theme_rki() + scale_fill_rki() + labs(y="Death Cases", x="") + guides(fill=FALSE) + facet_grid(agegroup ~ sex)

# Set and show analysis Groups
angroups_all <- create_angroups(andata,seqyear=25, seqage=51)
angroups <- create_angroups(andata,seqyear=4, seqage=51)<jupyter_output>[1] "Number of Groups to analyze: 4"
[1] "Number of Groups to analyze: 24"<jupyter_text>Overall Life Expectacy for Obervational Period<jupyter_code>lifetables.rec %>% filter(year>=1992 & year<=2016) %>% group_by(sex,age)  %>% summarise(qx=mean(qx,na.rm=T)) %>% group_by(sex) %>% life_expectancy(.) %>% filter(age==0 | age==65) %>% select(-qx) %>% spread(age,ex_t)<jupyter_output>Joining, by = c("sex", "age")<jupyter_text>Computation<jupyter_code># Log to File

# Model for all Angroups
#######################
registerDoMC(2)
anresult <- foreach(i=1:nrow(angroups)) %do% { # nrow(angroups)
    myresult <- survival_soep(data=andata,groupdefs=angroups[i,],determinant="incpos",controls=c("age"),byvars="sex") %>%
    merge_lifetable(.,angroups[i,])
    myresult
    }

# Store Result
anresult <- plyr::rbind.fill(anresult) %>% 
         filter(stringr::str_detect(term,"incpos")) %>% mutate(label = ordered(label,levels=levels(andata$incpos))) 

# Predict Lifeexpectancy
anresult.le <-   anresult %>% group_by(sex, begin, label) %>% arrange(sex, begin, label , age) %>%  predict_lifexp(.) 

# Overall Result 
#######################
overallresult <-  foreach(i=1:nrow(angroups_all)) %do% { 
    myresult <- survival_soep(data=andata,groupdefs=angroups_all[i,],determinant="incpos",controls=c("age"),byvars="sex") %>%
    merge_lifetable(.,angroups_all[i,])
    myresult
    }
registerDoMC(1) # Register 1 Core to reduce System load

# Store Result
overallresult <- plyr::rbind.fill(overallresult) %>% 
         filter(stringr::str_detect(term,"incpos")) %>% mutate(label = ordered(label,levels=levels(andata$incpos))) 

# Predict Lifeexpectancy
overallresult.le <- overallresult %>% group_by(sex, begin, label) %>% arrange(sex, begin, label , age) %>%  predict_lifexp(.)<jupyter_output>Independent Sampling design (with replacement)
Called via srvyr
Sampling variables:
 - ids: pid
 - weights: tsphrf
Independent Sampling design (with replacement)
Called via srvyr
Sampling variables:
 - ids: pid
 - weights: tsphrf<jupyter_text>Overall Results for the whole observational period Here we will compute the results for the whole observational period and export them to an excel file. We will show differences in Life Expectancy by income as well  plots by income for:- Hazard Ratios- Survival Rates- Life Expectancy at birth- Mortality risks until age of 65- Remaining Life Expectancy at 65<jupyter_code>myresult <- bind_rows(overallresult.le  %>% show_lediff(.,agefilter=0) %>% mutate(Outcome="LE at birth") ,
          overallresult.le  %>% show_lediff(.,agefilter=65)%>% mutate(Outcome="LE at 65 yrs") )
wb <- loadWorkbook("Results_Mortality.xlsx")
addWorksheet(wb, "Overall")
writeDataTable(wb, "Overall", myresult)
addWorksheet(wb, "Overall (Raw)")
writeDataTable(wb, "Overall (Raw)", overallresult.le)
myresult

# Plot LE at birth
le_plot_all <- overallresult.le %>% filter(age==0) %>% 
select(begin,end,sex,label,contains("LE_")) %>% ungroup() %>%
mutate(sex=ifelse(sex=="M","Men","Women"),
       sex=ordered(sex,levels=c("Men","Women"))) %>%  
       ggplot(.,aes(x=label,y=LE_mean, label=round(LE_mean, digits=1), fill=label)) +  facet_grid(.~sex) + theme_rki(base_size=8) + scale_fill_rki(reverse=F) +
labs(x="",y="Life Expectancy at birth", fill="Income")   + geom_bar(stat="identity") + coord_cartesian(ylim=c(65,85)) +
geom_text(aes(y=LE_mean+1), size=2.2) + guides(fill=FALSE)
le_plot_all
insertPlot(wb, "Overall", width = 7.11, height = 4, startRow = 1,  startCol = length(myresult)+5 , fileType = "png", units = "in", dpi = 400)

# Plot LE at 65
le65_plot_all <- overallresult.le %>% filter(age==65) %>% 
select(begin,end,sex,label,contains("LE_")) %>% ungroup() %>%
mutate(sex=ifelse(sex=="M","Men","Women"),
       sex=ordered(sex,levels=c("Men","Women"))) %>%  
       ggplot(.,aes(x=label,y=LE_mean, label=round(LE_mean, digits=1), fill=label)) +  facet_grid(.~sex) + theme_rki(base_size=8) + scale_fill_rki(reverse=F) +
labs(x="",y="Life Expectancy at 65", fill="Income")   + geom_bar(stat="identity") + coord_cartesian(ylim=c(0,20)) +
geom_text(aes(y=LE_mean+1), size=2.2) + guides(fill=FALSE)
le65_plot_all
insertPlot(wb, "Overall", width = 7.11, height = 4, startRow = 22,  startCol = length(myresult)+5 , fileType = "png", units = "in", dpi = 400)


# Rates at 65
surv65_plot_all <- overallresult.le %>% filter(age==65) %>% 
select(begin,end,sex,label,contains("Surv_")) %>% ungroup() %>%
mutate(Surv_mean=(1-Surv_mean)*100,
       sex=ifelse(sex=="M","Men","Women"),
       sex=ordered(sex,levels=c("Men","Women"))) %>%
       ggplot(.,aes(x=label,y=Surv_mean, label=round(Surv_mean, digits=1), fill=label)) +  facet_grid(.~sex) + theme_rki(base_size=8) + scale_fill_rki(reverse=F) +
labs(x="",y="Mortality Rate before age of 65", fill="Income")   + geom_bar(stat="identity")  + coord_cartesian(ylim=c(0,50)) +
geom_text(aes(y=Surv_mean+1), size=2.2) + guides(fill=FALSE)
surv65_plot_all
insertPlot(wb, "Overall", width = 7.11, height = 4, startRow = 7,  startCol = 1 , fileType = "png", units = "in", dpi = 400)

# Survival Rates
surv_rate_plot_all <- overallresult.le %>%  
select(begin,end,sex,age,label,contains("Surv_")) %>% ungroup() %>%
mutate(sex=ifelse(sex=="M","Men","Women"),
       sex=ordered(sex,levels=c("Men","Women"))) %>%  
       ggplot(.,aes(x=age,y=Surv_mean, label=round(Surv_mean, digits=1), color=label)) +  facet_grid(.~sex) + theme_rki(base_size=8) + scale_color_rki(reverse=F) +
labs(x="",y="Survival Rates", color="Income")  +
geom_line() + theme(legend.position="bottom")
surv_rate_plot_all
insertPlot(wb, "Overall", width = 7.11, height = 4, startRow = 27,  startCol = 1  , fileType = "png", units = "in", dpi = 400)

# Effects
effect_plot <- overallresult.le %>% select(sex,begin,end,label,youngest,oldest,estimate,contains("conf.")) %>% 
group_by(sex,begin,end,label,estimate,conf.low,conf.high) %>% filter(row_number()==1) %>% ungroup() %>%
mutate(Range=ifelse(oldest<max(oldest),paste0("Age: ",youngest,"-",oldest),paste0("Age: ",youngest,"+")),
      sex=ifelse(sex=="M","Men","Women")) %>%
ggplot(., aes(x=label,label=round(estimate,digits=1),y=estimate,ymin=conf.low,ymax=conf.high, color=label)) + geom_point() + geom_linerange()  +  facet_grid(Range~sex) +
 theme_rki(base_size=8) + scale_color_rki(reverse=F) + 
scale_y_log10(breaks=c(.5,1,2),labels=c(".5","Ref.","2"),limits=c(.4,3)) + theme(legend.position="bottom") + geom_text(color="black",size=2,hjust=-.5) +
labs(y="Hazard Ratio (Ref. Grand Mean)", x="Income Position") + guides(color=FALSE)
effect_plot
insertPlot(wb, "Overall", width = 7.11, height = 4, startRow = 42,   startCol = length(myresult)+5  , fileType = "png", units = "in", dpi = 400)
saveWorkbook(wb, "Results_Mortality.xlsx" , overwrite=TRUE)<jupyter_output><empty_output><jupyter_text>Results split by Timeframes four years Here we will replicate the results presented above differenciated by Time Frames of 4 years.<jupyter_code>myresult <- bind_rows(anresult.le  %>% show_lediff(.,agefilter=0) %>% mutate(Outcome="LE at birth") ,
          anresult.le  %>% show_lediff(.,agefilter=65)%>% mutate(Outcome="LE at 65 yrs") )
wb <- loadWorkbook("Results_Mortality.xlsx")
addWorksheet(wb, "By Timeframe")
writeDataTable(wb, "By Timeframe", myresult)
addWorksheet(wb, "By Timeframe (Raw)")
writeDataTable(wb, "By Timeframe (Raw)", anresult.le)
myresult
saveWorkbook(wb, "Results_Mortality.xlsx" , overwrite=TRUE)<jupyter_output><empty_output><jupyter_text>Plots for Life Expectancy at birth<jupyter_code>plotdata <- anresult.le %>% filter(age %in% c(0,65))  %>% select(sex,begin,end,age,label,contains("LE"),contains("Surv")) %>% ungroup() %>% 
            mutate(Geschlecht=ifelse(sex=="M","Men","Women"),
                  Geschlecht=ordered(Geschlecht,levels=c("Men","Women")),
                  Jahr=round((begin+end)/2)) %>% arrange(sex,begin,label) %>% 
            mutate_at(vars(contains("Surv")),funs(.*100))
plotdata_at_birth <- plotdata %>% filter(age==0) %>% select(label,year=Jahr,sex=Geschlecht,LE_mean) %>% spread(label,LE_mean) %>% arrange(sex,year)

# Data
plotdata_at_birth <- plotdata %>% filter(age==0) %>% select(label,year=Jahr,sex=Geschlecht,LE_mean) %>% spread(label,LE_mean) %>% arrange(sex,year) %>%
left_join(.,anresult.le  %>% show_lediff(.,agefilter=0) %>% ungroup() %>% 
          mutate(year=round((begin+end)/2),
                 sex=ifelse(sex=="M","Men","Women"),
                 sex=ordered(sex,levels=c("Men","Women"))) )
wb <- loadWorkbook("Results_Mortality.xlsx")
addWorksheet(wb, "Plot LE at birth")
writeDataTable(wb, "Plot LE at birth", plotdata_at_birth)
plotdata_at_birth

# LE
trendplot <- ggplot(plotdata %>% filter(age==0), aes(x=Jahr, y=LE_mean, ymin=LE_min,ymax=LE_max, color=label)) + geom_point() + geom_line() +
facet_grid(.~ Geschlecht)   + theme_rki(base_size = 10) + scale_color_rki() +
labs(color="Income Position", x="" , y="Life Expectancy at birth") + theme(legend.pos="bottom") 
trendplot
insertPlot(wb, "Plot LE at birth", width = 7.11, height = 4, startRow = 1,  startCol = length(plotdata_at_birth)+3 , fileType = "png", units = "in", dpi = 400)

# Difference
diffplot <- ggplot(plotdata_at_birth, aes(x=year , y=Diff, color=sex)) + geom_point() + theme_rki(base_size = 10) + scale_color_rki() +
facet_grid(.~ sex) + geom_smooth(method="gam")  + labs(y="LE at birth highest vs. lowest") + guides(color=FALSE)
diffplot
insertPlot(wb, "Plot LE at birth", width = 7.11, height = 4, startRow = 22,  startCol = length(plotdata_at_birth)+3 , fileType = "png", units = "in", dpi = 400)

# Effects
effect_plot_trend <- anresult.le %>% select(sex,begin,end,label,youngest,oldest,estimate,contains("conf.")) %>% 
group_by(sex,begin,end,label,estimate,conf.low,conf.high) %>% filter(row_number()==1) %>% ungroup() %>%
mutate(Range=ifelse(oldest<max(oldest),paste0("Age: ",youngest,"-",oldest),paste0("Age: ",youngest,"+")),
      sex=ifelse(sex=="M","Men","Women"),
      Year=round((begin+end)/2))  %>% filter(label %in% c("<60%","150% u.m.")) %>%
ggplot(., aes(x=Year,label=round(estimate,digits=1),y=estimate,group=label,color=label, fill=label)) +  
geom_line() + geom_ribbon(aes(ymin=conf.low,ymax=conf.high), color=NA,alpha=.6) +
 theme_rki(base_size=8) + scale_color_rki(reverse=F) + scale_fill_rki(reverse=F) + 
facet_grid(sex ~ Range) +
scale_y_log10(breaks=c(.5,1,2),labels=c(".5","Ref.","2")) + theme(legend.position="bottom") +  
labs(y="Hazard Ratio (Ref. Grand Mean)",color="Income Position") + guides(fill=FALSE)
effect_plot_trend
insertPlot(wb, "Plot LE at birth", width = 7.11, height = 4, startRow = 22,  startCol = 1 , fileType = "png", units = "in", dpi = 400)

# Excel Export
saveWorkbook(wb, "Results_Mortality.xlsx" , overwrite=TRUE)<jupyter_output>Joining, by = c("year", "sex")<jupyter_text>Life Expectancy at 65<jupyter_code># Data
plotdata_at_65 <- plotdata %>% filter(age==65) %>% select(label,year=Jahr,sex=Geschlecht,LE_mean) %>% spread(label,LE_mean) %>% arrange(sex,year) %>%
left_join(.,anresult.le  %>% show_lediff(.,agefilter=65) %>% ungroup() %>% 
          mutate(year=round((begin+end)/2),
                 sex=ifelse(sex=="M","Men","Women"),
                 sex=ordered(sex,levels=c("Men","Women"))) )
wb <- loadWorkbook("Results_Mortality.xlsx")
addWorksheet(wb, "Plot LE at 65")
writeDataTable(wb, "Plot LE at 65", plotdata_at_65)
plotdata_at_65

# LE
trendplot <- ggplot(plotdata %>% filter(age==65), aes(x=Jahr, y=LE_mean, ymin=LE_min,ymax=LE_max, color=label)) + geom_point() + geom_line() +
facet_grid(.~ Geschlecht)   + theme_rki(base_size = 10) + scale_color_rki() +
labs(color="Income Position", x="" , y="Life Expectancy at 65") + theme(legend.pos="bottom") 
trendplot
insertPlot(wb, "Plot LE at 65", width = 7.11, height = 4, startRow = 1,  startCol = length(plotdata_at_birth)+3 , fileType = "png", units = "in", dpi = 400)

# Surv
trendplot_surv <- ggplot(plotdata %>% filter(age==65), aes(x=Jahr, y=100-Surv_mean, ymin=Surv_min,ymax=Surv_max, color=label)) + geom_point() + geom_line() +
facet_grid(.~ Geschlecht)   + theme_rki(base_size = 10) + scale_color_rki() +
labs(color="Income Position", x="" , y="Mortality before age of 65") + theme(legend.pos="bottom") +  coord_cartesian(ylim=c(0,50))
trendplot_surv
insertPlot(wb, "Plot LE at 65", width = 7.11, height = 4, startRow = 22,  startCol = 1 , fileType = "png", units = "in", dpi = 400)


# Difference
diffplot <- ggplot(plotdata_at_65 , aes(x=year , y=Diff, color=sex)) + geom_point() + theme_rki(base_size = 10) + scale_color_rki() +
facet_grid(.~ sex) + geom_smooth(method="gam")  + labs(y="LE at 65 highest vs. lowest")  + guides(color=FALSE)
diffplot
insertPlot(wb, "Plot LE at 65", width = 7.11, height = 4, startRow = 22,  startCol = length(plotdata_at_birth)+3 , fileType = "png", units = "in", dpi = 400)

# Excel Export
saveWorkbook(wb, "Results_Mortality.xlsx" , overwrite=TRUE)<jupyter_output>Joining, by = c("year", "sex")<jupyter_text>Information on R Session and Packages<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1763969465648855, 'fraction_numerical': 0.019786259541984732, 'mean_word_length': 5.34586320480527, 'pattern_counts': {'":': 0, '<': 150, '<?xml version=': 0, '>': 319, 'https://': 5, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '21437', 'license': 'Apache-2.0', 'n_tokens_mistral': 12719, 'n_tokens_neox': 11626, 'n_words': 2147, 'path': 'files/SOEP_Mortality_JOHM_2019.ipynb', 'repo_name': 'lekroll/R', 'size': '32780'}