Notes

This report was generated on 2018-01-04 14:08:41. Output is available at https://graveja0.github.io/sccs

GitHub

The code and primary data for this project is available at https://github.com/graveja0/sccs.

Preparations

# from https://mran.revolutionanalytics.com/web/packages/checkpoint/vignettes/using-checkpoint-with-knitr.html
# if you don't need a package, remove it from here (commenting is probably not sufficient)
# tidyverse: see https://blog.rstudio.org/2016/09/15/tidyverse-1-0-0/
cat("
library(tidyverse) # ggplot2, dplyr, tidyr, readr, purrr, tibble
library(rlang)
library(MatchIt)
library(Hmisc)
library(rms) # for restricted cubic splines
library(magrittr) # pipes
library(stringr) # string manipulation
library(readxl) # excel
library(scales) # scales for ggplot2
library(jsonlite) # json
library(forcats) # easier factor handling,
library(lintr) # code linting
library(haven) # for loading SAS
library(janitor) # to clean up variable names 
library(stargazer)
library(ggthemes)
library(mice)
library(VIM)
library(DiagrammeR)
library(survminer)
library(directlabels)
library(knitr)
library(cem)
",
file = "manifest.R")

# For .gitignore in ./analysis/
# .httr-oauth
# .google_API_keys
# *.Rproj
# *.Rproj.user
# *.Rhistory
# *.Rprofile
# *.html
# auth
# .~lock*
# manifest.R
# output/ignore/*
# input/ignore/*
# scripts/ignore/*
# # if checkpoint is not yet installed, install it (for people using this
# # system for the first time)
# if (!require(checkpoint)) {
#   if (!require(devtools)) {
#     install.packages("devtools", repos = "http://cran.us.r-project.org")
#     require(devtools)
#   }
#   devtools::install_github("checkpoint",
#                            username = "RevolutionAnalytics",
#                            ref = "v0.3.2", # could be adapted later,
#                            # as of now (beginning of July 2017
#                            # this is the current release on CRAN)
#                            repos = "http://cran.us.r-project.org")
#   require(checkpoint)
# }
# # nolint start
# if (!dir.exists("~/.checkpoint")) {
#   dir.create("~/.checkpoint")
# }
# # nolint end
# # install packages for the specified CRAN snapshot date
# checkpoint(snapshotDate = package_date,
#            project = path_to_wd,
#            verbose = T,
#            scanForPackages = T,
#            use.knitr = F)
# rm(package_date)
source("manifest.R")
unlink("manifest.R")
select <- dplyr::select
sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] tcltk     grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] cem_1.1.17              knitr_1.15.1           
##  [3] directlabels_2017.03.31 survminer_0.4.0        
##  [5] ggpubr_0.1.5            DiagrammeR_0.9.0       
##  [7] VIM_4.7.0               data.table_1.10.4      
##  [9] colorspace_1.2-6        mice_2.30              
## [11] ggthemes_3.0.3          stargazer_5.2          
## [13] janitor_0.2.0           haven_0.2.1            
## [15] lintr_1.0.1             forcats_0.1.0          
## [17] jsonlite_1.5            scales_0.5.0           
## [19] readxl_1.0.0            stringr_1.2.0          
## [21] magrittr_1.5            rms_4.5-0              
## [23] SparseM_1.7             Hmisc_4.0-3            
## [25] Formula_1.2-1           survival_2.41-3        
## [27] lattice_0.20-33         MatchIt_2.4-21         
## [29] MASS_7.3-45             rlang_0.1.2            
## [31] dplyr_0.7.2             purrr_0.2.3            
## [33] readr_1.0.0             tidyr_0.6.3            
## [35] tibble_1.3.4            ggplot2_2.2.1          
## [37] tidyverse_1.0.0        
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.0-7       minqa_1.2.4         class_7.3-14       
##  [4] rprojroot_1.2       htmlTable_1.8       base64enc_0.1-3    
##  [7] rstudioapi_0.6      MatrixModels_0.4-1  mvtnorm_1.0-5      
## [10] codetools_0.2-14    splines_3.3.3       mnormt_1.5-4       
## [13] robustbase_0.92-6   nloptr_1.0.4        pbkrtest_0.4-6     
## [16] broom_0.4.1         km.ci_0.5-2         cluster_2.0.4      
## [19] backports_1.0.5     assertthat_0.2.0    Matrix_1.2-6       
## [22] lazyeval_0.2.0      acepack_1.3-3.3     visNetwork_2.0.0   
## [25] htmltools_0.3.6     quantreg_5.26       tools_3.3.3        
## [28] bindrcpp_0.2        igraph_1.0.1        gtable_0.2.0       
## [31] glue_1.1.1          reshape2_1.4.2      Rcpp_0.12.12       
## [34] cellranger_1.1.0    rgexf_0.15.3        nlme_3.1-128       
## [37] lmtest_0.9-34       psych_1.6.6         laeken_0.4.6       
## [40] lme4_1.1-12         XML_3.98-1.4        polspline_1.1.12   
## [43] DEoptimR_1.0-6      zoo_1.7-13          rex_1.1.1          
## [46] parallel_3.3.3      sandwich_2.3-4      RColorBrewer_1.1-2 
## [49] yaml_2.1.14         gridExtra_2.2.1     KMsurv_0.1-5       
## [52] rpart_4.1-10        latticeExtra_0.6-28 stringi_1.1.5      
## [55] Rook_1.1-1          randomForest_4.6-12 e1071_1.6-8        
## [58] checkmate_1.8.2     boot_1.3-18         pkgconfig_2.0.1    
## [61] evaluate_0.10       bindr_0.1           htmlwidgets_0.9    
## [64] cmprsk_2.2-7        plyr_1.8.4          R6_2.2.2           
## [67] combinat_0.0-8      multcomp_1.4-5      foreign_0.8-66     
## [70] mgcv_1.8-12         sp_1.2-3            nnet_7.3-12        
## [73] car_2.1-2           survMisc_0.5.4      rmarkdown_1.6      
## [76] viridis_0.4.0       influenceR_0.1.0    vcd_1.4-3          
## [79] digest_0.6.12       xtable_1.8-2        brew_1.0-6         
## [82] munsell_0.4.3       viridisLite_0.2.0   quadprog_1.5-5

Data Set Up

Read in and Tidy Data

nr.weight.ver <- "1-0" 
imputation.ver <- "1-0"
impute.data <- FALSE # This only needs to be run once
source("./scripts/read-in-and-tidy-data.R")
source("./scripts/project-functions.R")

Define Baseline, Outcome and other Key Variables

XX <- c(
  "education",
  "sex",
  "maritalstatus",
  "hhsize",
  "employed",
  "raceanalysis",
  "insurancecoverage"
)

XX.spline <- "vitalstatus_age"
XX.admin <- c("idnumber","r_1","r_2","r_3","vitalstatus","enrollment_source","state","agef3")
XX.other <- c("maritalstatusf1","hhincomef1","maritalstatusf2","hhincomef2","maritalstatusf3","hhincomef3","vitalstatus",
              "employmentstatusf1","bmif1")
XX.health <- c("smokenowf1",
               "menopausef1",
               "usualsourceofcaref1",
               "polypsf1",
               "anyaspirinf1",
               "statinusef1",
               "mammof1",
               "cancerf1",
               "bladdercancerf1",
               "bonecancerf1",
               "braincancerf1",
               "breastcancerf1",
               "cervicalcancerf1",
               "coloncancerf1",
               "colorectalcancerf1",
               "connecttisscancerf1",
               "esophagealcancerf1",
               "gallbladdercancerf1",
               "gastriccancerf1",
               "kidneycancerf1",
               "laryngealcancerf1",
               "livercancerf1",
               "lungcancerf1",
               "nonmelskincancerf1",
               "oralcancerf1",
               "ovariancancerf1",
               "pancreaticcancerf1",
               "pharyngealcancerf1",
               "prostatecancerf1",
               "rectalcancerf1",
               "thyroidcancerf1",
               "testicularcancerf1",
               "uterinecancerf1",
               "othercancerf1",     
               "diabetesf1",
               "parkinsonsf1",
               "msf1",
               "fibroidsf1",
               "bphf1",
               "mif1",
               "strokef1",
               "hipfxf1",
               "spinefxf1",
               "hodgkinlymphomaf1",
               "leukemiaf1",
               "lymphomaf1",
               "melanomaf1",
               "multmyelomaf1",
               "breastbxf1",
               "papf1",
               "psaf1",
               "dref1",
               "prostatebxf1",
               "colonoscopyf1",
               "sigmoidoscopyf1",
               "fecaloccultf1",
               "diabetesglucosetestf1",
               "aspirinf1")
YY.health <- grep("sf12v2",names(xx0),value=TRUE)
YY.health <- YY.health[-grep("_nbs$",YY.health)]
YY.insurance <- c("insurancecoveragef1","insurancecoveragef2","insurancecoveragef3")

# Add in all other FU1 variables
x <- grep("f1$",names(xx0),value=TRUE)
y <- c(XX.spline,XX.admin,XX.other,XX.health,YY.health,YY.insurance) %>% unique()

# CHeck for variables that have lots of 8888/7777
XX.FU1.other <- x[!x%in%y]
XX.FU1.other <- XX.FU1.other[-grep("_agef1$|^age|questionnaire|months|_med",XX.FU1.other)]
test <- XX.FU1.other %>% map(~prop.table(table(xx0[,.x])))
names(test) = XX.FU1.other

analytic.vars <- c(XX,XX.spline,XX.other,XX.health,YY.health,YY.insurance) %>% unique()
XX.tot <- c(XX,XX.health,XX.other,XX.spline,XX.admin,YY.health,YY.insurance,"intervention","expansion_state","jan1_2008_to_vs_months",XX.FU1.other)

source("./scripts/consort-diagram-main-samples.R")

Imputation

## Imputation of Missing Values
xx.preimp <- xx0 %>% filter(idnumber %in% c(ids.ss.within.state,ids.ss.across.state.nonelderly))
impute.data = FALSE
source("./scripts/impute-missing-values.R")

# Adding in non-imputed variables 
not.in.xx <- c("idnumber",names(xx.preimp)[-(which(names(xx.preimp) %in% intersect(names(xx.preimp),names(xx.i))))])
xx <- xx.i %>% left_join(xx.preimp[,not.in.xx],"idnumber")

Define Samples

In this section we will define some key samples to use later in the analysis. Note that ALL samples are defined among FU1 responders in the 12 SCCS states, and include both FU3 responders as well as those who were lost to follow-up due to death. The FU3 survey nonresponders will be dealt with via nonresponse weights.

  1. The Nonelderly Across-State sample are individuals aged 64 or younger as of the last vital status update (12/31/15).

  2. The Within-State Medicare sample is broadly comprised of FU1 and FU3 responders (plus those who died). We identify Medicare status as of FU1 to ensure consistency in insurance coverage over time. W

  3. The Within-State Nonelderly Medicare sample is a subset of the within-state Medicare sample: only those who were 64 or less as of the last vital status update (12/31/15).

  4. The Within-State Nonelderly General Population and Clinic Population sample is also a subset of the within-state sample: only individuals 64 and younger who were either recruited at the health clinics, or via the general population mail survey.

  5. We define a few sub-samples for robustness checks by only using those identified with FU1 income <400% FPL.

xx[,"temp"] <- xx[,"medicare"]
    xx[,"tempf1"] <- xx[,paste0("medicare","f1")]
    xx[,"tempf2"] <- xx[,paste0("medicare","f2")]
    xx[,"tempf3"] <- xx[,paste0("medicare","f3")]
    xx <- xx %>% 
      mutate(temp = as.numeric(paste0(temp)),
             tempf1 = as.numeric(paste0(tempf1)),
             tempf2 = as.numeric(paste0(tempf2)),
             tempf3 = as.numeric(paste0(tempf3))) %>%
      mutate(
        temp = ifelse(temp==7777,0,ifelse(temp==8888,NA,temp)),
        tempf1 = pmax(as.integer(temp==1),as.integer(tempf1==1),9999*as.integer(tempf1>7777)),
        tempf2 = pmax(as.integer(tempf1==1),as.integer(tempf2==1),9999*as.integer(tempf2>7777)),
        tempf3 = pmax(as.integer(tempf2==1),as.integer(tempf3==1),9999*as.integer(tempf3>7777))
      )
    xx[,"medicare"] <- xx[,"temp"] 
    xx[,paste0("medicare","f1")] <- xx[,"tempf1"]
    xx[,paste0("medicare","f2")] <- xx[,"tempf2"]
    xx[,paste0("medicare","f3")] <-  xx[,"tempf3"]
    
    
# Define Samples
xx <- xx %>%
  mutate(medicare_baseline_to_fu2 = as.integer(medicare==1 | medicaref1==1 | medicaref2==1)) %>% 
  mutate(medicare_at_baseline = as.integer(medicare==1)) %>% 
  mutate(ss.across.state.nonelderly = as.integer(idnumber %in% ids.ss.across.state.nonelderly)) %>% 
  mutate(ss.within.state = as.integer(idnumber %in% ids.ss.within.state)) %>% 
  mutate(expansion_state = as.integer(state %in% c("05", "21", "54", "22")),
         expansion_state_noLA = as.integer(state %in% c("05", "21", "54")),
         inc_lt_400 = as.integer(
                        (maritalstatusf1=="1" &  hhincomef1 %in% c("1","2","3","4")) | 
                        (maritalstatusf1!="1" & hhincomef1 %in% c("1","2","3"))
                      )) %>% 
  mutate(
         general_population = as.integer(enrollment_source=="M"), 
         general_population400 = as.integer(enrollment_source=="M" & inc_lt_400==0)
         ) %>% 
  mutate(intervention = as.integer(inc_lt_400==1 & medicaref1!=1 & !is.na(medicaref1)  )) %>% 
  # Near-Elderly (58-63) vs. Newly Elderly (66-74)
  mutate(ss.neareld = as.integer(
          ss.within.state ==1  & 
         ( ( inc_lt_400 & age_asof_2014 %in% 60:63) | ( age_asof_2014 %in% 65:69 & inc_lt_400 )) )) %>% 
  # Nonelderly Medicare vs. Nonelderly Intevention
  mutate( 
    ss.noneld.medicare = as.integer(ss.within.state ==1  & age_asof_2014<65 & (intervention==1 | medicaref1==1))
    ) %>% 
  # Nonelderly Intervention vs. Nonelderly High Income General Population Sample
  mutate(
    ss.genpop = as.integer(ss.within.state ==1  & age_asof_2014<65 & medicaref1!=1)
  ) %>% 
  # Pooled Controls
  mutate(ss.pooled = as.integer(ss.neareld == 1 | ss.genpop==1 | ss.noneld.medicare==1))
         
XX.samples <- c("expansion_state","expansion_state_noLA","inc_lt_400","general_population",
                "general_population400","medicare_baseline_to_fu2",
                "AFU1_DFU2_DFU3","AFU1_AFU2_DFU3","AFU1_AFU2_AFU3","DIED_AFTER_FU3",
                "age_asof_2014","medicaref1","medicaref3","ss.neareld","ss.noneld.medicare",
                "ss.genpop","ss.pooled")

XX.tot <- c(XX.tot,XX.samples) %>% unique()
# Further Data Edits

# Merge in the Nonresponse weight
# source("./scripts/nonresponse-weights.R")
load(file.path(fileoutputdir,paste0("nonresponse-weight-",nr.weight.ver,".RData")))
xx <-  xx %>% left_join(nr.weight, "idnumber") %>%
  mutate(w.FU3 = ifelse(is.na(w.FU3), 0, w.FU3))

xx <-
  xx %>% 
  mutate(
    insurancef1 = factor(insurancecoveragef1,labels=c("Uninsured","Insured")),
    insurancef2 = factor(ifelse(AFU1_DFU2_DFU3==1,3,insurancecoveragef2),
                         labels=c("Uninsured","Insured","Death")),
    insurancef3 = factor(ifelse(AFU1_DFU2_DFU3==1|AFU1_AFU2_DFU3==1,3,insurancecoveragef3),
                         labels=c("Uninsured","Insured","Death"))
  ) %>% 
  mutate(education=factor(education),
       sex = factor(sex),
       maritalstatus=factor(maritalstatus),
       employed = factor(employed),
       raceanalysis = factor(raceanalysis),
       insurancecoverage = factor(insurancecoverage),
       sf12v2healthf1 = factor(sf12v2healthf1))  %>% 
  # Survival censored at 2014
  mutate(jan1_2010_to_vs_months_c14 = jan1_2010_to_vs_months) %>% 
  mutate(vs_c14 = vitalstatus) %>% 
  mutate(vs_c14 = ifelse(jan1_2010_to_vs_months>480,"A",vs_c14)) %>% 
  mutate(jan1_2010_to_vs_months_c14 = ifelse(jan1_2010_to_vs_months>480,48,jan1_2010_to_vs_months))

Propensity Scores and Matching

Standardized Difference in Means: Near-Elderly Low-Income vs. Newly-Elderly Low-Income Control Group

Standardized Difference in Means: Near-Elderly Low-Income vs. Newly-Elderly Low-Income Control Group 

par(mfrow=c(2,3))
rms::survplot(f.surv.temp,xlab="Months Since 1/1/2010")
title(main="Unadjusted Sample")
rms::survplot(f.surv.temp.w,xlab="Months Since 1/1/2010")
title(main="Weighted Sample")
rms::survplot(f.surv.temp.m,xlab="Months Since 1/1/2010")
title(main="Matched Sample")
rms::survdiffplot(f.surv.diff.temp,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nUnadjusted Sample")
rms::survdiffplot(f.surv.diff.temp.w,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nWeighted Sample")
rms::survdiffplot(f.surv.diff.temp.m,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nMatched Sample")
Survival Comparison: Near-Elderly Low-Income vs. Newly-Elderly Low-Income Control Group

Survival Comparison: Near-Elderly Low-Income vs. Newly-Elderly Low-Income Control Group 

xx$ss.sample = xx$ss.noneld.medicare
ps.noneld.medicare <- psfit(indata = subset(xx,ss.sample==1 & jan1_2010_to_vs_months>=0  ),
                    Txvar="intervention",
                    idvar = "idnumber", 
                    psvars = psvars,
                    psvars.spline = "bmif1",
                    match_caliper = 0.10)

temp.psfit <- ps.noneld.medicare
xx.survival.temp <-   xx %>% 
filter(ss.sample==1 & jan1_2010_to_vs_months>=0) %>% 
mutate(id = idnumber) %>% 
left_join(temp.psfit$psdata,"id")  %>% 
tbl_df()
dd <- suppressWarnings(xx.survival.temp %>% select(-r_1) %>% datadist())
options(datadist = 'dd')
# Unmatched Sample
f.surv.diff.temp <- rms::npsurv(Surv(jan1_2010_to_vs_months_c14,as.integer(vs_c14=="D"))~intervention,data=xx.survival.temp)
f.surv.temp <- rms::cph(Surv(jan1_2010_to_vs_months_c14,as.integer(vs_c14=="D"))~strat(intervention),data=xx.survival.temp,y=TRUE)
# Weighted Sample
f.surv.diff.temp.w <- rms::npsurv(Surv(jan1_2010_to_vs_months_c14,as.integer(vs_c14=="D"))~intervention,data=subset(xx.survival.temp,ps_weight>0),weights=ps_weight)
f.surv.temp.w <- rms::cph(Surv(jan1_2010_to_vs_months,as.integer(vs_c14=="D"))~strat(intervention),data=subset(xx.survival.temp,ps_weight>0),y=TRUE,weights=ps_weight)
# Matched Sample
dd <- xx.survival.temp %>% select(-r_1) %>% filter(matched==1) %>% datadist()
options(datadist = 'dd')
f.surv.diff.temp.m <- rms::npsurv(Surv(jan1_2010_to_vs_months_c14,as.integer(vs_c14=="D"))~intervention,data=subset(xx.survival.temp,matched==1))
f.surv.temp.m <- rms::cph(Surv(jan1_2010_to_vs_months_c14,as.integer(vs_c14=="D"))~strat(intervention),data=subset(xx.survival.temp,matched==1),y=TRUE)
Standardized Difference in Means: Low-Income Nonelderly vs. Medicare-Insured Nonelderly Comparison Group

Standardized Difference in Means: Low-Income Nonelderly vs. Medicare-Insured Nonelderly Comparison Group

Survival Comparison: Low-Income Nonelderly vs. Medicare-Insured Nonelderly Comparison Group

Survival Comparison: Low-Income Nonelderly vs. Medicare-Insured Nonelderly Comparison Group

Standardized Difference in Means: Low-Income Nonelderly vs. Higher-Income Nonelderly

Standardized Difference in Means: Low-Income Nonelderly vs. Higher-Income Nonelderly 

par(mfrow=c(2,3))
rms::survplot(f.surv.temp,xlab="Months Since 1/1/2010")
title(main="Unadjusted Sample")
rms::survplot(f.surv.temp.w,xlab="Months Since 1/1/2010")
title(main="Weighted Sample")
rms::survplot(f.surv.temp.m,xlab="Months Since 1/1/2010")
title(main="Matched Sample")
rms::survdiffplot(f.surv.diff.temp,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nUnadjusted Sample")
rms::survdiffplot(f.surv.diff.temp.w,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nWeighted Sample")
rms::survdiffplot(f.surv.diff.temp.m,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nMatched Sample")
Survival Comparison: Low-Income Nonelderly vs. Higher-Income Nonelderly

Survival Comparison: Low-Income Nonelderly vs. Higher-Income Nonelderly

Standardized Difference in Means: Low-Income Nonelderly vs. Pooled Comparison Groups

Standardized Difference in Means: Low-Income Nonelderly vs. Pooled Comparison Groups 

par(mfrow=c(2,3))
rms::survplot(f.surv.temp,xlab="Months Since 1/1/2010")
title(main="Unadjusted Sample")
rms::survplot(f.surv.temp.w,xlab="Months Since 1/1/2010")
title(main="Weighted Sample")
rms::survplot(f.surv.temp.m,xlab="Months Since 1/1/2010")
title(main="Matched Sample")
rms::survdiffplot(f.surv.diff.temp,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nUnadjusted Sample")
rms::survdiffplot(f.surv.diff.temp.w,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nWeighted Sample")
rms::survdiffplot(f.surv.diff.temp.m,xlab="Months Since 1/1/2010",ylim=c(-.15,.15))
title(main="Survival Difference:\nMatched Sample")
Survival Comparison: Low-Income Nonelderly vs. Pooled Comparison Groups

Survival Comparison: Low-Income Nonelderly vs. Pooled Comparison Groups 

ps.genpop.data <- ps.genpop$psdata %>% select(idnumber = id, ps_weight_genpop = ps_weight , matched_genpop = matched)
ps.neareld.data <- ps.neareld$psdata %>% select(idnumber = id, ps_weight_neareld = ps_weight , matched_neareld = matched)
ps.noneld.medicare.data <- ps.noneld.medicare$psdata %>% select(idnumber = id, ps_weight_noneld.medicare = ps_weight , matched_noneld.medicare = matched)
ps.pooled.data <- ps.pooled$psdata %>% select(idnumber = id, ps_weight_pooled = ps_weight , matched_pooled = matched)

xx <- 
  xx %>% 
  left_join(ps.genpop.data,"idnumber") %>% 
  left_join(ps.neareld.data,"idnumber") %>% 
  left_join(ps.noneld.medicare.data,"idnumber") %>% 
  left_join(ps.pooled.data,"idnumber")

Composite Health Status

ggplot_missing <- function(x){
  x %>%
    is.na %>%
    reshape2::melt() %>%
    ggplot(data = .,
           aes(x = Var2,
               y = Var1)) +
    geom_raster(aes(fill = value)) +
    scale_fill_grey(name = "",
                    labels = c("Present","Missing")) +
    theme_minimal() +
    theme(axis.text.x  = element_text(angle=45, vjust=0.5)) +
    labs(x = "Variables in Dataset",
         y = "Rows / observations")
}

source("./scripts/composite-sf12-score.R")
xx <-
  xx %>% left_join(sf12f1[,c("idnumber","sf12_physicalf1","sf12_mentalf1")], "idnumber") %>% left_join(sf12f3[,c("idnumber","sf12_physicalf3","sf12_mentalf3")], "idnumber")

sf12f1.unimputed %>% select(starts_with("i")) %>% select(-idnumber) %>% filter(row_number() %in% sample(seq(nrow(sf12f1)),200)) %>% ggplot_missing() -> p.sf12f1.unimp

sf12f1 %>% select(starts_with("i")) %>% select(-idnumber) %>% filter(row_number() %in% sample(seq(nrow(sf12f1)),200)) %>% ggplot_missing() ->
  p.sf12f1.imp

#cowplot::plot_grid(p.sf12f1.unimp,p.sf12f1.imp,labels=c("Raw Data","Imputed Data"))

Comparison of Baseline (FU1) Composite Health Status After Matching: Nonelderly Expansion State (Tx) vs. Nonexpansion State (Cx) 

xx %>% filter(matched_neareld==1) %>% 
  mutate(intervention = factor(intervention,labels=c("Control","Intervention"))) %>% 
  mutate(healthf1 = sf12_physicalf1, 
         healthf3 = sf12_physicalf3,
         mentalf1 = sf12_mentalf1,
         mentalf3 = sf12_mentalf3) %>% 
  filter(!is.na(healthf3) & !(is.na(healthf1))) %>% 
  select(intervention,healthf3,healthf1) %>% 
  mutate(health.cs = healthf3-healthf1) %>% 
  with(histboxp(x=healthf1,group=intervention,sd=TRUE,bins=200,xlab="Composite Physical Health Status (FU1)"))

Comparison of Baseline (FU1) Composite Health Status After Matching: Low-Income Near-Elderly (Tx) vs. Low-Income Newly-Elderly (Cx) 

xx %>% filter(matched_noneld.medicare==1) %>% 
  mutate(intervention = factor(intervention,labels=c("Control","Intervention"))) %>% 
  mutate(healthf1 = sf12_physicalf1, 
         healthf3 = sf12_physicalf3,
         mentalf1 = sf12_mentalf1,
         mentalf3 = sf12_mentalf3) %>% 
  filter(!is.na(healthf3) & !(is.na(healthf1))) %>% 
  select(intervention,healthf3,healthf1) %>% 
  mutate(health.cs = healthf3-healthf1) %>% 
  with(histboxp(x=healthf1,group=intervention,sd=TRUE,bins=200,xlab="Composite Physical Health Status (FU1)"))

Comparison of Baseline (FU1) Composite Health Status After Matching: Within-State Sample: Nonelderly Low-Income (Tx) vs. Nonelderly Medicare (Cx) 

xx %>% filter(matched_genpop==1) %>% 
  mutate(intervention = factor(intervention,labels=c("Control","Intervention"))) %>% 
  mutate(healthf1 = sf12_physicalf1, 
         healthf3 = sf12_physicalf3,
         mentalf1 = sf12_mentalf1,
         mentalf3 = sf12_mentalf3) %>% 
  filter(!is.na(healthf3) & !(is.na(healthf1))) %>% 
  select(intervention,healthf3,healthf1) %>% 
  mutate(health.cs = healthf3-healthf1) %>% 
  with(histboxp(x=healthf1,group=intervention,sd=TRUE,bins=200,xlab="Composite Physical Health Status (FU1)"))

Comparison of Baseline (FU1) Composite Health Status After Matching: Low-Income Nonelderly (Tx) vs. Higher-Income Nonelderly (Cx) 

xx %>% filter(matched_pooled==1) %>% 
  mutate(intervention = factor(intervention,labels=c("Control","Intervention"))) %>% 
  mutate(healthf1 = sf12_physicalf1, 
         healthf3 = sf12_physicalf3,
         mentalf1 = sf12_mentalf1,
         mentalf3 = sf12_mentalf3) %>% 
  filter(!is.na(healthf3) & !(is.na(healthf1))) %>% 
  select(intervention,healthf3,healthf1) %>% 
  mutate(health.cs = healthf3-healthf1) %>% 
  with(histboxp(x=healthf1,group=intervention,sd=TRUE,bins=200,xlab="Composite Physical Health Status (FU1)"))

Comparison of Baseline (FU1) Composite Health Status After Matching: Low-Income Nonelderly vs. Pooled Cx Groups 

source("./scripts/sccs-dd-regression-function.R")
xx <- 
  xx %>% 
  mutate(insurance01f1 = as.integer(insurancecoveragef1==2),
         insurance01f3 = as.integer(insurancecoveragef3==2),
         ex_vgf1  = as.integer(sf12v2healthf1 %in% 1:2),
         ex_vgf3  = as.integer(sf12v2healthf3 %in% 1:2),
         fr_prf1  = as.integer(sf12v2healthf1 %in% 4:5),
         fr_prf3  = as.integer(sf12v2healthf3 %in% 4:5))
source("./scripts/PS_M_weighting_SA.R")
source("./scripts/PS_M_weighting.R")
source("./scripts/PS_SM_weighting.R")

fit.health <- function(ss.ps,yy="ex_vg") {
  ss.naive <- ss.ps %>% filter(vitalstatus=="A")
  fit.naive <- 
    Fit(YY= yy,  # Outcome
      XX=unique(c(XX,XX.health)), # Covariates
      XX.spline = c("vitalstatus_age"), # Covariates that enter as a spline
      PS = FALSE, # Fit the principal stratification analysis?
     Sample = "within", # Across-state or within-state analysis?
     Data=ss.naive,
     quiet=TRUE
  )
  
  # PS Bounds
  Y <- ss.ps[,paste0(yy,"f3")] - ss.ps[,paste0(yy,"f1")]
  S <- ss.ps$vitalstatus
  Z <- ss.ps$intervention
  bounds.est <- try(ps_bounds(Z=Z,S=S,Y=Y,Tx=1,Cx=0,PS="A",length=5))
  
  # PS
  set.seed(45)
  ps.est <- try(
    Fit.PS(YY=yy,
    XX=unique(XX.orig),
    XX.spline=XX.spline,
    Data = ss.ps,n.iter=10
  ))
  
  # PS Monotonicity Sensitivity Analysis
  ps.est.sa <- try(
    Fit.PS(YY=yy,
    XX=XX.orig,
    XX.spline=XX.spline,
    Data = ss.ps,
    n.iter=2,
    n.sa = 5,
    point=FALSE
  ))
  
  out <- list(
      tab  = ss.ps %>% count(intervention,vitalstatus),
      naive.unadj = fit.naive$fit.unadj$coefficients["post_intervention_1"],
      naive.adj  = fit.naive$fit.adj$coefficients["post_intervention_1"],
      ps.unadj = ps.est$point[1],
      ps.adj = ps.est$point[2],
      ps.bounds = bounds.est$formatted,
      ps.sa = ps.est.sa
  )
  return(out)
  }
fit.health <- FALSE 

ss.ps.neareld <- xx %>% filter(ss.neareld==1 & matched_neareld==1) 
ss.ps.genpop <- xx %>% filter(ss.genpop==1 & matched_genpop==1) 
ss.ps.noneld.medicare <- xx %>% filter(ss.noneld.medicare==1 & matched_noneld.medicare==1) 

if (fit.health) {
  yy  = "fr_pr"
  fit.neareld <- fit.health(ss.ps = ss.ps.neareld,yy=yy)
  save(fit.neareld,file=paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
  # fit.genpop <- fit.health(ss.ps = ss.ps.genpop,yy=yy)
  # save(fit.genpop,file=paste0("./output/fit-health/fit-health-genpop-",yy,".RData"))
  fit.noneld.medicare <- fit.health(ss.ps = ss.ps.noneld.medicare,yy=yy)
  save(fit.noneld.medicare,file=paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
  
  yy  = "ex_vg"
  fit.neareld <- fit.health(ss.ps = ss.ps.neareld,yy=yy)
  save(fit.neareld,file=paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
  # fit.genpop <- fit.health(ss.ps = ss.ps.genpop,yy=yy)
  # save(fit.genpop,file=paste0("./output/fit-health/fit-health-genpop-",yy,".RData"))
  fit.noneld.medicare <- fit.health(ss.ps = ss.ps.noneld.medicare,yy=yy)
  save(fit.noneld.medicare,file=paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
  
  yy  = "sf12_physical"
  fit.neareld <- fit.health(ss.ps = ss.ps.neareld,yy=yy)
  save(fit.neareld,file=paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
  # fit.genpop <- fit.health(ss.ps = ss.ps.genpop,yy=yy)
  # save(fit.genpop,file=paste0("./output/fit-health/fit-health-genpop-",yy,".RData"))
  fit.noneld.medicare <- fit.health(ss.ps = ss.ps.noneld.medicare,yy=yy)
  save(fit.noneld.medicare,file=paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
  
  yy  = "sf12_mental"
  fit.neareld <- fit.health(ss.ps = ss.ps.neareld,yy=yy)
  save(fit.neareld,file=paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
  # fit.genpop <- fit.health(ss.ps = ss.ps.genpop,yy=yy)
  # save(fit.genpop,file=paste0("./output/fit-health/fit-health-genpop-",yy,".RData"))
  fit.noneld.medicare <- fit.health(ss.ps = ss.ps.noneld.medicare,yy=yy)
  save(fit.noneld.medicare,file=paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
  
}

Near-Elderly Comparison Group

The following analyses take as their sample the matched group of nearl-elderly intervention (i.e., aged 60-63 as of 2014, <400% FPL at FU1) individuals and their matched control group of newly-elderly individuals (i.e., aged 65-69 as of 2014, <400% FPL at FU1).

The first table simply demonstrates the high mortality rates (all were alive as of 1/1/2010) in both groups, though the intervention group has a slightly lower mortality rate indicating relatively healthier individuals.

yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))

fit.neareld[[1]] %>% 
  spread(vitalstatus,n) %>% 
  mutate(intervention = ifelse(intervention == 1,"Intervention Group","Newly-Elderly Controls")) %>% 
  mutate(tot = A+D) %>% 
  mutate(A = 100*round(A/tot,3),
         D = 100*round(D/tot,3)) %>% 
  select(-tot) %>% 
  rename(Alive = A,
         Dead = D,
         Group = intervention)  %>% kable(caption = "Observed Vital Status: Newly-Elderly Controls")
Observed Vital Status: Newly-Elderly Controls
Group Alive Dead
Newly-Elderly Controls 77.9 22.1
Intervention Group 82.0 18.0 
# Get first stage
ss.naive <- ss.ps.neareld %>% filter(vitalstatus=="A")
fit.naive <- 
  Fit(YY= "insurance01",  # Outcome
      XX=unique(c(XX,XX.health)), # Covariates
      XX.spline = c("vitalstatus_age"), # Covariates that enter as a spline
      PS = FALSE, # Fit the principal stratification analysis?
      Sample = "within", # Across-state or within-state analysis?
      Data=ss.naive,
      quiet=TRUE
  )

first_stage_unadj <- fit.naive$fit.unadj$coefficients[["post_intervention_1"]]
first_stage_adj <- fit.naive$fit.adj$coefficients[["post_intervention_1"]]

yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
phys.neareld <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.neareld[[2]] ,fit.neareld[[3]]),
        rbind(fit.neareld[[2]]/first_stage_unadj ,fit.neareld[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(phys.neareld) <- c("DDCoverage","DDHealth","IV") 
phys.neareld <- 
  phys.neareld %>% 
  mutate(Outcome = c("Composite Physical Health","Composite Physical Health"),
         adj = c("Unadjusted","Adjusted"))

yy = "sf12_mental"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
ment.neareld <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.neareld[[2]] ,fit.neareld[[3]]),
        rbind(fit.neareld[[2]]/first_stage_unadj ,fit.neareld[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(ment.neareld) <- c("DDCoverage","DDHealth","IV") 
ment.neareld <- 
  ment.neareld %>% 
  mutate(Outcome = c("Composite Mental Health","Composite Mental Health"),
         adj = c("Unadjusted","Adjusted"))


yy = "ex_vg"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
ex_vg.neareld <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.neareld[[2]] ,fit.neareld[[3]]),
        rbind(fit.neareld[[2]]/first_stage_unadj ,fit.neareld[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(ex_vg.neareld) <- c("DDCoverage","DDHealth","IV") 
ex_vg.neareld <- 
  ex_vg.neareld %>% 
  mutate(Outcome = c("Excellent or Very Good Health","Excellent or Very Good Health"),
         adj = c("Unadjusted","Adjusted"))

yy = "fr_pr"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
fr_pr.neareld <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.neareld[[2]] ,fit.neareld[[3]]),
        rbind(fit.neareld[[2]]/first_stage_unadj ,fit.neareld[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(fr_pr.neareld) <- c("DDCoverage","DDHealth","IV") 
fr_pr.neareld <- 
  fr_pr.neareld %>% 
  mutate(Outcome = c("Fair or Poor Health","Fair or Poor Health"),
         adj = c("Unadjusted","Adjusted"))


rbind(phys.neareld,ment.neareld,ex_vg.neareld,fr_pr.neareld) %>% 
  gather(key=variable,value=value,-adj,-Outcome) %>% 
  #spread(adj,value) %>% 
  reshape2::dcast(Outcome~adj+variable) %>% 
  select(Outcome,Unadjusted_DDCoverage,Unadjusted_DDHealth,Unadjusted_IV,
         Adjusted_DDCoverage,Adjusted_DDHealth,Adjusted_IV) -> result.neareld

The following 2 tables report regression coefficient \(\beta_3\) from a DD regression of the form: \[ h(E(Y_{itk})) = \beta_0 + \beta_1 Year_t + \beta_2 Intervention + \beta_3 PostExpansion_t \times Intervention_t + \beta_4 State_k + \\ \beta_5 State_k \times Year_t + \beta_6 State_k \times Intervention_i + \beta_7RCS(Age_{it})+\beta_8OtherSociodemo_i+\beta_9Health_i \] For these results we do a “naive” analysis that simply subsets to the living. That is, we ignore the possible sample selection effets of mortality in these analyses. The first table reports unadjusted results (i.e., no age, sociodemographic, or health controls) while the second table reports the adjusted results. We consider four outcomes:

  • Composite Physical and Mental Health: based on the composite aggregation of the SF12 questions.
  • Excellent or Very Good Health: based on the ordinal health status question
  • Fair or Poor Health: based on the ordinal health status question.

The first column of results provides DD estimates from a regression with a binary indicator for insurance coverage as the outcome. The second column then reports the \(\hat \beta_3\) coefficient for the health outcome. If we assume that the entire differential change in health outcomes is mediated by the expansions, then we can obtain an effective IV estimate of the effect of coverage by dividing this number by the differential change in coverage from the first column; that is the result depicted in the last column (marked “IV”).

Health Status Results: Ignoring Death, Unadjusted, Newly-Elderly Comparison Group
Outcome Unadjusted_DDCoverage Unadjusted_DDHealth Unadjusted_IV
Composite Mental Health 0.22187 -0.08169 -0.36817
Composite Physical Health 0.22187 0.51462 2.31945
Excellent or Very Good Health 0.22187 0.01210 0.05454
Fair or Poor Health 0.22187 -0.00237 -0.01069
Health Status Results: Ignoring Death, Adjusted, Newly-Elderly Comparison Group
Outcome Adjusted_DDCoverage Adjusted_DDHealth Adjusted_IV
Composite Mental Health 0.22127 -0.10378 -0.46899
Composite Physical Health 0.22127 0.50165 2.26708
Excellent or Very Good Health 0.22127 0.01148 0.05189
Fair or Poor Health 0.22127 -0.00149 -0.00672 
yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
foo1 <- fit.neareld$ps.bounds
aace1 <- fit.neareld$ps.sa[,c(1,2)]
aace_reg1 <- fit.neareld$ps.sa[,c(1,3)]

yy = "sf12_mental"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
foo2 <- fit.neareld$ps.bounds[,3]
aace2 <- fit.neareld$ps.sa[,c(2)]
aace_reg2 <- fit.neareld$ps.sa[,c(3)]

yy = "ex_vg"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
foo3 <- fit.neareld$ps.bounds[,3]
aace3 <- fit.neareld$ps.sa[,c(3)]
aace_reg3 <- fit.neareld$ps.sa[,c(2)]

yy = "fr_pr"
load(paste0("./output/fit-health/fit-health-neareld-",yy,".RData"))
foo4 <- fit.neareld$ps.bounds[,3]
aace4 <- fit.neareld$ps.sa[,c(3)]
aace_reg4 <- fit.neareld$ps.sa[,c(2)]

The following table provides bounds on the survivor average treatment effects (i.e., ATEs in the principal stratum of always-survivors) under various assumptions. The underlying bounds equations are provided in Zhang and Rubin (2003).

  • Monotonicity in this context refers to the assumption that assignment to insurance can only keep you alive. That is, there are no individuals who would die under treatment, but live if not assigned insurance. This may not hold in this context, since some people who receive insurance could receive care (e.g., surgery) with some death risk. We will examinine the sensitivity of results to the monotonicity assumption directly below.

  • The Ranked Average assumption is slightly weaker. It says that the average health in the always-living principal stratum will be greater than the average health in the other two stratum of individuals who would either die without insurance but live with insurance, or who would die with insurance but live with insurance.

One item to note is that the underlying results reflect a DD estimate based on the more standard \(post + intervention + post \times intervention\) specification. That is because this framework maps nicely to the PS framework in that we can regress an outcome defined by the change score (i.e., FU3 health status minus FU1 health status) on a Treatment indicator; the coefficient on this Tx indicator is the same as would be obtained from the interaction term in the DD set up above.

As can be seen in the results (which are shown for each of the outcomes), imposing one or both of these assumptions can help to sharpen the bounds.

cbind(foo1,foo2,foo3,foo4) %>% 
  rename(physical=bound,
         mental = foo2,
         ex_vg = foo3,
         fr_pr = foo4) %>% kable(cap="Principal Stratification Bounds: AACE, Newly-Elderly Comparison Group ")
Principal Stratification Bounds: AACE, Newly-Elderly Comparison Group
monotonicty ranked_avg physical mental ex_vg fr_pr
N N [-8.5, 9.5] [-9.77, 9.75] [-2, 2] [-2, 2]
Y N [-0.74, 1.76] [-1.39, 1.31] [-0.09, 0.1] [-0.13, 0.15]
N Y [-3.75, 5.28] [-4.49, 5.05] [-0.98, 1] [-1.03, 1.02]
Y Y [0.48, 1.76] [-0.01, 1.31] [0.01, 0.1] [-0.01, 0.15]

The following two tables provide point estimates for the DD estimate for the principal stratum of always-survivors. Each row corresponds to different values of a monotonicity sensitivity parameter (\(\eta\)), which is bounded below by 0 and from above by the value in the last column (the equation to find bounds for \(\eta\) can be found in the Ding & Lu paper).

When \(\eta=0\) there is no monotonicity. By comparison, as values of \(\eta\) increase we assume greater and greater departure from monotonicity. The results make clear that our point estimates are largely invariant to sensitivity here.

The first table below shows the AACE (i.e., Always-Survivor Causal Effect) estimates unadjusted by X’s. The second table shows analogous results while conditioning on all the X’s.

cbind(aace1,aace2,aace3,aace4) %>% 
  rename(physical=aace,
         mental = aace2,
         ex_vg = aace3,
         fr_pr = aace4) %>% 
  kable(cap="Principal Stratification Monotonicity Sensitivity: AACE, Newly-Elderly Comparison Group ")
Principal Stratification Monotonicity Sensitivity: AACE, Newly-Elderly Comparison Group
eta physical mental ex_vg fr_pr
0.00000 0.46147 -0.05682 0.01290 0.00053
0.20313 0.45432 -0.08010 0.01271 0.00107
0.40625 0.43310 -0.10138 0.01171 0.00078
0.60938 0.41616 -0.07439 0.01158 0.00056
0.81250 0.46647 -0.00859 0.01366 0.00045 
cbind(aace_reg1,aace_reg2,aace_reg3,aace_reg4) %>% 
  rename(physical=aace_reg,
         mental = aace_reg2,
         ex_vg = aace_reg3,
         fr_pr = aace_reg4) %>% 
  kable(cap="Principal Stratification Monotonicity Sensitivity: AACE (Adjsuted), Newly-Elderly Comparison Group ")
Principal Stratification Monotonicity Sensitivity: AACE (Adjsuted), Newly-Elderly Comparison Group
eta physical mental ex_vg fr_pr
0.00000 0.25172 -0.03934 0.01496 -0.00462
0.20313 0.24054 -0.04545 0.01463 -0.00391
0.40625 0.20870 -0.04478 0.01325 -0.00370
0.60938 0.19214 -0.06339 0.01306 -0.00472
0.81250 0.22408 -0.08174 0.01573 -0.00704

Nonelderly Medicare Comparison Group

The following tables repeat the analyses above for an alternative comparison group of non-elderly Medicaire beneficiaries.

yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))

fit.noneld.medicare[[1]] %>% 
  spread(vitalstatus,n) %>% 
  mutate(intervention = ifelse(intervention == 1,"Intervention Group","Newly-Elderly Controls")) %>% 
  mutate(tot = A+D) %>% 
  mutate(A = 100*round(A/tot,3),
         D = 100*round(D/tot,3)) %>% 
  select(-tot) %>% 
  rename(Alive = A,
         Dead = D,
         Group = intervention)  %>% kable(caption = "Observed Vital Status: Newly-Elderly Controls")
Observed Vital Status: Newly-Elderly Controls
Group Alive Dead
Newly-Elderly Controls 80.4 19.6
Intervention Group 84.5 15.5 
# Get first stage
ss.naive <- ss.ps.noneld.medicare %>% filter(vitalstatus=="A")
fit.naive <- 
  Fit(YY= "insurance01",  # Outcome
      XX=unique(c(XX,XX.health)), # Covariates
      XX.spline = c("vitalstatus_age"), # Covariates that enter as a spline
      PS = FALSE, # Fit the principal stratification analysis?
      Sample = "within", # Across-state or within-state analysis?
      Data=ss.naive,
      quiet=TRUE
  )

first_stage_unadj <- fit.naive$fit.unadj$coefficients[["post_intervention_1"]]
first_stage_adj <- fit.naive$fit.adj$coefficients[["post_intervention_1"]]

yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
phys.noneld.medicare <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.noneld.medicare[[2]] ,fit.noneld.medicare[[3]]),
        rbind(fit.noneld.medicare[[2]]/first_stage_unadj ,fit.noneld.medicare[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(phys.noneld.medicare) <- c("DDCoverage","DDHealth","IV") 
phys.noneld.medicare <- 
  phys.noneld.medicare %>% 
  mutate(Outcome = c("Composite Physical Health","Composite Physical Health"),
         adj = c("Unadjusted","Adjusted"))

yy = "sf12_mental"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
ment.noneld.medicare <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.noneld.medicare[[2]] ,fit.noneld.medicare[[3]]),
        rbind(fit.noneld.medicare[[2]]/first_stage_unadj ,fit.noneld.medicare[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(ment.noneld.medicare) <- c("DDCoverage","DDHealth","IV") 
ment.noneld.medicare <- 
  ment.noneld.medicare %>% 
  mutate(Outcome = c("Composite Mental Health","Composite Mental Health"),
         adj = c("Unadjusted","Adjusted"))


yy = "ex_vg"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
ex_vg.noneld.medicare <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.noneld.medicare[[2]] ,fit.noneld.medicare[[3]]),
        rbind(fit.noneld.medicare[[2]]/first_stage_unadj ,fit.noneld.medicare[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(ex_vg.noneld.medicare) <- c("DDCoverage","DDHealth","IV") 
ex_vg.noneld.medicare <- 
  ex_vg.noneld.medicare %>% 
  mutate(Outcome = c("Excellent or Very Good Health","Excellent or Very Good Health"),
         adj = c("Unadjusted","Adjusted"))

yy = "fr_pr"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
fr_pr.noneld.medicare <- 
  cbind(rbind(fit.naive$fit.unadj$coefficients[["post_intervention_1"]],
              fit.naive$fit.adj$coefficients[["post_intervention_1"]]),
        rbind(fit.noneld.medicare[[2]] ,fit.noneld.medicare[[3]]),
        rbind(fit.noneld.medicare[[2]]/first_stage_unadj ,fit.noneld.medicare[[3]]/first_stage_adj))  %>% 
  tbl_df() 
names(fr_pr.noneld.medicare) <- c("DDCoverage","DDHealth","IV") 
fr_pr.noneld.medicare <- 
  fr_pr.noneld.medicare %>% 
  mutate(Outcome = c("Fair or Poor Health","Fair or Poor Health"),
         adj = c("Unadjusted","Adjusted"))


rbind(phys.noneld.medicare,ment.noneld.medicare,ex_vg.noneld.medicare,fr_pr.noneld.medicare) %>% 
  gather(key=variable,value=value,-adj,-Outcome) %>% 
  #spread(adj,value) %>% 
  reshape2::dcast(Outcome~adj+variable) %>% 
  select(Outcome,Unadjusted_DDCoverage,Unadjusted_DDHealth,Unadjusted_IV,
         Adjusted_DDCoverage,Adjusted_DDHealth,Adjusted_IV) -> result.noneld.medicare
Health Status Results: Ignoring Death, Unadjusted, Newly-Elderly Comparison Group
Outcome Unadjusted_DDCoverage Unadjusted_DDHealth Unadjusted_IV
Composite Mental Health 0.23272 -0.05247 -0.22549
Composite Physical Health 0.23272 0.43575 1.87245
Excellent or Very Good Health 0.23272 0.00350 0.01504
Fair or Poor Health 0.23272 0.00618 0.02658
Health Status Results: Ignoring Death, Adjusted, Newly-Elderly Comparison Group
Outcome Adjusted_DDCoverage Adjusted_DDHealth Adjusted_IV
Composite Mental Health 0.23332 -0.03744 -0.16048
Composite Physical Health 0.23332 0.43388 1.85964
Excellent or Very Good Health 0.23332 0.00355 0.01522
Fair or Poor Health 0.23332 0.00578 0.02475 
yy = "sf12_physical"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
foo1 <- fit.noneld.medicare$ps.bounds
aace1 <- fit.noneld.medicare$ps.sa[,c(1,2)]
aace_reg1 <- fit.noneld.medicare$ps.sa[,c(1,3)]

yy = "sf12_mental"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
foo2 <- fit.noneld.medicare$ps.bounds[,3]
aace2 <- fit.noneld.medicare$ps.sa[,c(2)]
aace_reg2 <- fit.noneld.medicare$ps.sa[,c(3)]

yy = "ex_vg"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
foo3 <- fit.noneld.medicare$ps.bounds[,3]
aace3 <- fit.noneld.medicare$ps.sa[,c(3)]
aace_reg3 <- fit.noneld.medicare$ps.sa[,c(2)]

yy = "fr_pr"
load(paste0("./output/fit-health/fit-health-noneld-medicare-",yy,".RData"))
foo4 <- fit.noneld.medicare$ps.bounds[,3]
aace4 <- fit.noneld.medicare$ps.sa[,c(3)]
aace_reg4 <- fit.noneld.medicare$ps.sa[,c(2)]
cbind(foo1,foo2,foo3,foo4) %>% 
  rename(physical=bound,
         mental = foo2,
         ex_vg = foo3,
         fr_pr = foo4) %>% kable(cap="Principal Stratification Bounds: AACE, Newly-Elderly Comparison Group ")
Principal Stratification Bounds: AACE, Newly-Elderly Comparison Group
monotonicty ranked_avg physical mental ex_vg fr_pr
N N [-7.35, 8.3] [-9.15, 8.81] [-2, 2] [-2, 2]
Y N [-0.73, 1.65] [-1.5, 1.29] [-0.07, 0.08] [-0.12, 0.17]
N Y [-3.08, 4.75] [-4.16, 4.77] [-1, 1] [-1.04, 1.04]
Y Y [0.45, 1.65] [-0.08, 1.29] [0, 0.08] [0, 0.17] 
cbind(aace1,aace2,aace3,aace4) %>% 
  rename(physical=aace,
         mental = aace2,
         ex_vg = aace3,
         fr_pr = aace4) %>% 
  kable(cap="Principal Stratification Monotonicity Sensitivity: AACE, Newly-Elderly Comparison Group ")
Principal Stratification Monotonicity Sensitivity: AACE, Newly-Elderly Comparison Group
eta physical mental ex_vg fr_pr
0.00000 0.48205 -0.01909 0.00563 0.00415
0.19867 0.47860 -0.02313 0.00530 0.00365
0.39735 0.47376 -0.03766 0.00508 0.00346
0.59602 0.47252 -0.03720 0.00548 0.00349
0.79469 0.49844 0.01376 0.00706 0.00304 
cbind(aace_reg1,aace_reg2,aace_reg3,aace_reg4) %>% 
  rename(physical=aace_reg,
         mental = aace_reg2,
         ex_vg = aace_reg3,
         fr_pr = aace_reg4) %>% 
  kable(cap="Principal Stratification Monotonicity Sensitivity: AACE (Adjsuted), Newly-Elderly Comparison Group ")
Principal Stratification Monotonicity Sensitivity: AACE (Adjsuted), Newly-Elderly Comparison Group
eta physical mental ex_vg fr_pr
0.00000 0.48669 0.04413 0.00530 0.00388
0.19867 0.47648 0.05511 0.00477 0.00379
0.39735 0.46409 0.05705 0.00440 0.00382
0.59602 0.46473 0.05591 0.00471 0.00366
0.79469 0.51757 0.07255 0.00587 0.00210