Living Conditions and Dropout (step 2)

Survival Analyses (Feb 2023)

Show code

require(timereg)

require(MatchIt)
vars_cov<-c("edad_al_ing_sc", "sexo_2", "escolaridad_rec", "sus_ini_mod_mvv", "freq_cons_sus_prin", "dg_trs_cons_sus_or", "sus_principal_mod")

prueba2_imp22$edad_al_ing_sc<-scale(prueba2_imp22$edad_al_ing, scale=F)
prueba2_imp22$edad_al_ing_sc<-as.numeric(prueba2_imp22$edad_al_ing_sc)

m.out1 <- matchit(con_quien_vive_joel_fam_or ~ edad_al_ing + sexo_2+escolaridad_rec+sus_ini_mod_mvv+freq_cons_sus_prin+dg_trs_cons_sus_or+sus_principal_mod+cnt_mod_cie_10_or, 
                  data = prueba2_imp22 %>% dplyr::filter(con_quien_vive_joel_alone!=1), 
                  method = "nearest", 
                  discard = "both",
                  caliper = .01, 
                  standardize = T)

plot(
cobalt::bal.tab(m.out1, un = TRUE, m.threshold = .1, 
        v.threshold = 2))

Show code

m.out2 <- matchit(con_quien_vive_joel_alone ~ edad_al_ing+ sexo_2+escolaridad_rec+sus_ini_mod_mvv+freq_cons_sus_prin+dg_trs_cons_sus_or+sus_principal_mod+cnt_mod_cie_10_or, 
                  data = prueba2_imp22 %>% dplyr::filter(con_quien_vive_joel_fam_or!=1), 
                  method = "nearest", 
                  discard = "both",
                  caliper = .01, 
                  standardize = T)

plot(
cobalt::bal.tab(m.out2, un = TRUE, m.threshold = .1, 
        v.threshold = 2))

Show code

m.out3 <- matchit(con_quien_vive_joel_alone ~ edad_al_ing + sexo_2+escolaridad_rec+sus_ini_mod_mvv+freq_cons_sus_prin+dg_trs_cons_sus_or+sus_principal_mod+cnt_mod_cie_10_or, 
                  data = prueba2_imp22 %>% dplyr::filter(con_quien_vive_joel_coup_chil!=1), 
                  method = "nearest", 
                  discard = "both",
                  caliper = .01, 
                  standardize = T)

plot(
cobalt::bal.tab(m.out3, un = TRUE, m.threshold = .1, 
        v.threshold = 2))

Show code

fam_or_vs_couple_child<-match.data(m.out1)
alone_vs_couple_child<-match.data(m.out2)
alone_vs_fam_orig<-match.data(m.out3)

sub1_statadf_miss<- dplyr::mutate(fam_or_vs_couple_child, yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~0.0001,T~yrs_to_tr_dropout)) %>% 
  #dplyr::mutate(con_quien_vive_joel_rec= ifelse(con_quien_vive_joel=="Alone",1,0)) %>%  
  slice(1:5e5) %>% 
  data.table::data.table()

sub2_statadf_miss<- dplyr::mutate(alone_vs_couple_child, yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~0.0001,T~yrs_to_tr_dropout)) %>% 
  #dplyr::mutate(con_quien_vive_joel_rec= ifelse(con_quien_vive_joel=="Alone",1,0)) %>%  
  slice(1:5e5) %>% 
  data.table::data.table()

sub3_statadf_miss<- dplyr::mutate(alone_vs_fam_orig, yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~0.0001,T~yrs_to_tr_dropout)) %>% 
  #dplyr::mutate(con_quien_vive_joel_rec= ifelse(con_quien_vive_joel=="Alone",1,0)) %>%  
  slice(1:5e5) %>% 
  data.table::data.table()


cat("Living with Family of origin vs. With couple/children at 6 months")

Living with Family of origin vs. With couple/children at 6 months

Show code

survRM2::rmst2(sub1_statadf_miss$yrs_to_tr_dropout, sub1_statadf_miss$event, as.numeric(sub1_statadf_miss$con_quien_vive_joel=="Family of origin"), tau=.5)

The truncation time: tau = 0.5  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 0.394 0.002     0.390     0.397
RMST (arm=0) 0.384 0.002     0.381     0.388


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 0.106 0.002     0.103     0.110
RMTL (arm=0) 0.116 0.002     0.112     0.119


Between-group contrast 
                      Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) 0.010     0.005     0.014 0
RMST (arm=1)/(arm=0) 1.025     1.012     1.037 0
RMTL (arm=1)/(arm=0) 0.918     0.879     0.958 0

Show code

cat("Living with Family of origin vs. With couple/children at 1 year")

Living with Family of origin vs. With couple/children at 1 year

Show code

survRM2::rmst2(sub1_statadf_miss$yrs_to_tr_dropout, sub1_statadf_miss$event, as.numeric(sub1_statadf_miss$con_quien_vive_joel=="Family of origin"), tau=1)

The truncation time: tau = 1  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 0.608 0.004      0.60     0.616
RMST (arm=0) 0.578 0.004      0.57     0.587


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 0.392 0.004     0.384      0.40
RMTL (arm=0) 0.422 0.004     0.413      0.43


Between-group contrast 
                      Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) 0.030     0.018     0.041 0
RMST (arm=1)/(arm=0) 1.051     1.031     1.072 0
RMTL (arm=1)/(arm=0) 0.930     0.903     0.957 0

Show code

cat("Living with Family of origin vs. With couple/children at 3 years")

Living with Family of origin vs. With couple/children at 3 years

Show code

survRM2::rmst2(sub1_statadf_miss$yrs_to_tr_dropout, sub1_statadf_miss$event, as.numeric(sub1_statadf_miss$con_quien_vive_joel=="Family of origin"), tau=3)

The truncation time: tau = 3  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 1.248 0.015     1.218     1.277
RMST (arm=0) 1.127 0.014     1.099     1.155


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 1.752 0.015     1.723     1.782
RMTL (arm=0) 1.873 0.014     1.845     1.901


Between-group contrast 
                      Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) 0.120     0.080     0.161 0
RMST (arm=1)/(arm=0) 1.107     1.069     1.145 0
RMTL (arm=1)/(arm=0) 0.936     0.915     0.957 0

Show code

cat("===============================================================================")

===============================================================================

Show code

cat("Living Alone vs. Family of origin at 6 months")

Living Alone vs. Family of origin at 6 months

Show code

survRM2::rmst2(sub3_statadf_miss$yrs_to_tr_dropout, sub3_statadf_miss$event, sub3_statadf_miss$con_quien_vive_joel_alone, tau=.5)

The truncation time: tau = 0.5  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 0.358 0.004     0.350     0.367
RMST (arm=0) 0.385 0.004     0.377     0.393


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 0.142 0.004     0.133     0.150
RMTL (arm=0) 0.115 0.004     0.107     0.123


Between-group contrast 
                       Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) -0.027    -0.039    -0.015 0
RMST (arm=1)/(arm=0)  0.930     0.901     0.960 0
RMTL (arm=1)/(arm=0)  1.235     1.126     1.355 0

Show code

cat("Living Alone vs. Family of origin at 1 year")

Living Alone vs. Family of origin at 1 year

Show code

survRM2::rmst2(sub3_statadf_miss$yrs_to_tr_dropout, sub3_statadf_miss$event, sub3_statadf_miss$con_quien_vive_joel_alone, tau=1)

The truncation time: tau = 1  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.   se lower .95 upper .95
RMST (arm=1) 0.544 0.01     0.524     0.563
RMST (arm=0) 0.598 0.01     0.578     0.617


Restricted Mean Time Lost (RMTL) by arm 
              Est.   se lower .95 upper .95
RMTL (arm=1) 0.456 0.01     0.437     0.476
RMTL (arm=0) 0.402 0.01     0.383     0.422


Between-group contrast 
                       Est. lower .95 upper .95 p
RMST (arm=1)-(arm=0) -0.054    -0.081    -0.026 0
RMST (arm=1)/(arm=0)  0.910     0.867     0.955 0
RMTL (arm=1)/(arm=0)  1.134     1.063     1.209 0

Show code

cat("Living Alone vs. Family of origin at 3 years")

Living Alone vs. Family of origin at 3 years

Show code

survRM2::rmst2(sub3_statadf_miss$yrs_to_tr_dropout, sub3_statadf_miss$event, sub3_statadf_miss$con_quien_vive_joel_alone, tau=3)

The truncation time: tau = 3  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 1.116 0.034     1.050     1.182
RMST (arm=0) 1.264 0.035     1.196     1.332


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 1.884 0.034     1.818     1.950
RMTL (arm=0) 1.736 0.035     1.668     1.804


Between-group contrast 
                       Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) -0.148    -0.243    -0.053 0.002
RMST (arm=1)/(arm=0)  0.883     0.815     0.957 0.002
RMTL (arm=1)/(arm=0)  1.085     1.030     1.144 0.002

Show code

cat("===============================================================================")

===============================================================================

Show code

cat("Living with With couple/children(1) vs. Alone at 6 months")

Living with With couple/children(1) vs. Alone at 6 months

Show code

survRM2::rmst2(sub2_statadf_miss$yrs_to_tr_dropout, sub2_statadf_miss$event, sub2_statadf_miss$con_quien_vive_joel_coup_chil, tau=.5)

The truncation time: tau = 0.5  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 0.377 0.004     0.369     0.385
RMST (arm=0) 0.359 0.004     0.350     0.367


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 0.123 0.004     0.115     0.131
RMTL (arm=0) 0.141 0.004     0.133     0.150


Between-group contrast 
                      Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) 0.019     0.007     0.030 0.002
RMST (arm=1)/(arm=0) 1.052     1.018     1.086 0.002
RMTL (arm=1)/(arm=0) 0.869     0.794     0.950 0.002

Show code

cat("Living with With couple/children(1) vs. Alone at 1 year")

Living with With couple/children(1) vs. Alone at 1 year

Show code

survRM2::rmst2(sub2_statadf_miss$yrs_to_tr_dropout, sub2_statadf_miss$event, sub2_statadf_miss$con_quien_vive_joel_coup_chil, tau=1)

The truncation time: tau = 1  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.   se lower .95 upper .95
RMST (arm=1) 0.568 0.01     0.549     0.587
RMST (arm=0) 0.544 0.01     0.524     0.564


Restricted Mean Time Lost (RMTL) by arm 
              Est.   se lower .95 upper .95
RMTL (arm=1) 0.432 0.01     0.413     0.451
RMTL (arm=0) 0.456 0.01     0.436     0.476


Between-group contrast 
                      Est. lower .95 upper .95    p
RMST (arm=1)-(arm=0) 0.024    -0.003     0.052 0.08
RMST (arm=1)/(arm=0) 1.045     0.995     1.098 0.08
RMTL (arm=1)/(arm=0) 0.946     0.890     1.007 0.08

Show code

cat("Living with With couple/children(1) vs. Alone at 3 years")

Living with With couple/children(1) vs. Alone at 3 years

Show code

survRM2::rmst2(sub2_statadf_miss$yrs_to_tr_dropout, sub2_statadf_miss$event, sub2_statadf_miss$con_quien_vive_joel_coup_chil, tau=3)

The truncation time: tau = 3  was specified. 

Restricted Mean Survival Time (RMST) by arm 
              Est.    se lower .95 upper .95
RMST (arm=1) 1.111 0.033     1.046     1.175
RMST (arm=0) 1.116 0.034     1.049     1.182


Restricted Mean Time Lost (RMTL) by arm 
              Est.    se lower .95 upper .95
RMTL (arm=1) 1.889 0.033     1.825     1.954
RMTL (arm=0) 1.884 0.034     1.818     1.951


Between-group contrast 
                       Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) -0.005    -0.098     0.087 0.912
RMST (arm=1)/(arm=0)  0.995     0.916     1.082 0.912
RMTL (arm=1)/(arm=0)  1.003     0.955     1.053 0.912

Time for this code chunk to run: 0.2 minutes

Show code

# 
# out<-cox.aalen(Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel, data=prueba2_imp22)+  edad_al_ing_sc+sexo_2+escolaridad_rec+sus_ini_mod_mvv+freq_cons_sus_prin+dg_trs_cons_sus_or+sus_principal_mod+cnt_mod_cie_10_or,data=prueba2_imp22, max.timepoint.sim=NULL, resample.iid=1)
#  
# coxrm <- restricted.residual.mean(out, tau=23, x=rbind(0, 1),iid=1)
# summary(coxrm)
#+ edad_al_ing_sc + sexo_2 + escolaridad_rec + sus_ini_mod_mvv + freq_cons_sus_prin + dg_trs_cons_sus_or + sus_principal_mod + cnt_mod_cie_10_or
out <- cox.aalen(Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel , data=data.table::data.table(prueba2_imp22), max.time=7, n.sim=0, resample.iid=1)# , data=prueba2_imp22)
#Error in model.frame.default(formula = Surv(survs$stop, survs$status) ~  : invalid type (NULL) for variable 'Z'
coxrm <- restricted.residual.mean(out, tau=3, x=rbind(0, 1),iid=1)

Time for this code chunk to run: 0 minutes

Show code

require(rstpm2)

#Define covariates
vars_cov<-c("edad_al_ing", "sexo_2", "escolaridad_rec", "sus_ini_mod_mvv", "freq_cons_sus_prin", "dg_trs_cons_sus_or", "sus_principal_mod", "cnt_mod_cie_10_or")

#Scale age at admission to overcome convergence issues in flexsurv
prueba2_imp22$edad_al_ing_sc<-scale(prueba2_imp22$edad_al_ing, scale=F)

# prueba2_imp22<- 
#   prueba2_imp22 %>% 
#   dplyr::filter()
#select cases with t>=0


warning(paste0("Centered at ",attr(scale(prueba2_imp22$edad_al_ing, scale=F),"scaled:center")," years old"))
##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#Compute frecuency tables to see what is the mode by each covariate to choose an standard covariate profile
tables <- list()
for(i in seq_along(vars_cov)) {
    tables[[i]] <- dim(table(prueba2_imp22[[vars_cov[i]]]))
    attr(tables[[i]], "var_name") <- names(vars_cov)[i]
}


##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
# make database to compute STPM2
rstmp_options<-expand.grid(tvc=1:5,df=1:8)

invisible("Casos completos")
table(complete.cases(prueba2_imp22[,c(vars_cov,"yrs_to_tr_dropout","event","con_quien_vive_joel")]))

 TRUE 
23979 

Show code

rstmp_options$AIC<-0
rstmp_options$BIC<-0
for (i in seq_along(rstmp_options$tvc)){
  dfs<-rstmp_options[i, "df"]
  tvcs <-rstmp_options[i, "tvc"]
  tryCatch({
assign(paste0("rstpm_",dfs,"_tvc_",tvcs),stpm2(Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+  edad_al_ing_sc+ sexo_2+ escolaridad_rec+ sus_ini_mod_mvv+ freq_cons_sus_prin+ dg_trs_cons_sus_or+ sus_principal_mod+ cnt_mod_cie_10_or, data=prueba2_imp22,df=dfs, tvc=list(con_quien_vive_joel_coup_chil=tvcs, con_quien_vive_joel_fam_or=tvcs)), envir=.GlobalEnv)
  }, error= function(e){
    warning("Error: prbabilly starting values")
  })
  rstmp_options$AIC[i]<-if(is.null(AIC(get(paste0("rstpm_",dfs,"_tvc_",tvcs))))){0} else { AIC(get(paste0("rstpm_",dfs,"_tvc_",tvcs)))}
  rstmp_options$BIC[i]<-if(is.null(BIC(get(paste0("rstpm_",dfs,"_tvc_",tvcs))))){
        0
    } else { 
      BIC(get(paste0("rstpm_",dfs,"_tvc_",tvcs)))
      }
}

assign(paste0("rstpm_",5,"_tvc_",0),stpm2(Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+  edad_al_ing_sc+ sexo_2+ escolaridad_rec+ sus_ini_mod_mvv+ freq_cons_sus_prin+ dg_trs_cons_sus_or+ sus_principal_mod+ cnt_mod_cie_10_or, data=prueba2_imp22,df=dfs, envir=.GlobalEnv))


arrange(rstmp_options, AIC, BIC) %>% 
  knitr::kable(size=10, format="markdown",caption= "FIt indices", escape=T)

Table 1: FIt indices

tvc	df	AIC	BIC
2	8	29391.76	29674.68
3	8	29393.82	29692.91
2	7	29395.68	29670.52
5	8	29396.46	29727.88
4	8	29397.47	29712.72
2	6	29399.97	29666.72
5	7	29400.54	29723.88
4	7	29401.53	29708.71
3	6	29402.24	29685.16
5	6	29405.09	29720.35
4	6	29405.92	29705.01
3	5	29419.88	29694.72
4	5	29421.79	29712.80
5	5	29423.44	29730.61
5	4	29436.25	29735.34
3	4	29437.42	29704.18
5	3	29441.66	29732.67
3	7	29453.26	29744.26
2	5	29467.13	29725.80
2	4	29480.97	29731.56
4	4	29485.70	29768.62
1	7	29490.25	29748.93
1	6	29494.99	29745.58
5	2	29504.12	29787.05
4	3	29505.85	29780.69
1	5	29513.93	29756.44
4	2	29522.32	29789.08
1	8	29525.21	29791.97
1	4	29532.67	29767.09
2	3	29672.19	29914.70
3	3	29674.61	29933.29
3	2	29703.39	29953.98
1	3	29769.12	29995.46
2	2	30190.43	30424.85
1	2	30341.16	30559.42
5	1	30531.88	30806.72
4	1	30550.64	30809.31
3	1	30748.58	30991.09
2	1	31218.70	31445.04
1	1	49590.39	49800.56

Show code

anova(rstpm_4_tvc_4, rstpm_5_tvc_5, rstpm_5_tvc_1)

Likelihood Ratio Tests
Model 1: rstpm_4_tvc_4, [negll]: (Intercept)+
          con_quien_vive_joel_coup_chil+con_quien_vive_joel_fam_or+
          edad_al_ing_sc+sexo_2Women+escolaridad_rec.L+
          escolaridad_rec.Q+sus_ini_mod_mvvCocaine hydrochloride+
          sus_ini_mod_mvvMarijuana+sus_ini_mod_mvvOther+
          sus_ini_mod_mvvCocaine paste+freq_cons_sus_prin.L+
          freq_cons_sus_prin.Q+freq_cons_sus_prin.C+
          freq_cons_sus_prin^4+dg_trs_cons_sus_orHazardous
          consumption+sus_principal_modCocaine hydrochloride+
          sus_principal_modMarijuana+sus_principal_modOther+
          sus_principal_modCocaine paste+cnt_mod_cie_10_or.L+
          cnt_mod_cie_10_or.Q+cnt_mod_cie_10_or.C+
          nsx(log(yrs_to_tr_dropout), df = 4)1+
          nsx(log(yrs_to_tr_dropout), df = 4)2+
          nsx(log(yrs_to_tr_dropout), df = 4)3+
          nsx(log(yrs_to_tr_dropout), df = 4)4+
          nsx(log(yrs_to_tr_dropout), df =
          4)1:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          4)2:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          4)3:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          4)4:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          4)1:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          4)2:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          4)3:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          4)4:as.numeric(con_quien_vive_joel_fam_or)
Model 2: rstpm_5_tvc_5, [negll]: (Intercept)+
          con_quien_vive_joel_coup_chil+con_quien_vive_joel_fam_or+
          edad_al_ing_sc+sexo_2Women+escolaridad_rec.L+
          escolaridad_rec.Q+sus_ini_mod_mvvCocaine hydrochloride+
          sus_ini_mod_mvvMarijuana+sus_ini_mod_mvvOther+
          sus_ini_mod_mvvCocaine paste+freq_cons_sus_prin.L+
          freq_cons_sus_prin.Q+freq_cons_sus_prin.C+
          freq_cons_sus_prin^4+dg_trs_cons_sus_orHazardous
          consumption+sus_principal_modCocaine hydrochloride+
          sus_principal_modMarijuana+sus_principal_modOther+
          sus_principal_modCocaine paste+cnt_mod_cie_10_or.L+
          cnt_mod_cie_10_or.Q+cnt_mod_cie_10_or.C+
          nsx(log(yrs_to_tr_dropout), df = 5)1+
          nsx(log(yrs_to_tr_dropout), df = 5)2+
          nsx(log(yrs_to_tr_dropout), df = 5)3+
          nsx(log(yrs_to_tr_dropout), df = 5)4+
          nsx(log(yrs_to_tr_dropout), df = 5)5+
          nsx(log(yrs_to_tr_dropout), df =
          5)1:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          5)2:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          5)3:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          5)4:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          5)5:as.numeric(con_quien_vive_joel_coup_chil)+
          nsx(log(yrs_to_tr_dropout), df =
          5)1:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          5)2:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          5)3:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          5)4:as.numeric(con_quien_vive_joel_fam_or)+
          nsx(log(yrs_to_tr_dropout), df =
          5)5:as.numeric(con_quien_vive_joel_fam_or)
Model 3: rstpm_5_tvc_1, [negll]: (Intercept)+
          con_quien_vive_joel_coup_chil+con_quien_vive_joel_fam_or+
          edad_al_ing_sc+sexo_2Women+escolaridad_rec.L+
          escolaridad_rec.Q+sus_ini_mod_mvvCocaine hydrochloride+
          sus_ini_mod_mvvMarijuana+sus_ini_mod_mvvOther+
          sus_ini_mod_mvvCocaine paste+freq_cons_sus_prin.L+
          freq_cons_sus_prin.Q+freq_cons_sus_prin.C+
          freq_cons_sus_prin^4+dg_trs_cons_sus_orHazardous
          consumption+sus_principal_modCocaine hydrochloride+
          sus_principal_modMarijuana+sus_principal_modOther+
          sus_principal_modCocaine paste+cnt_mod_cie_10_or.L+
          cnt_mod_cie_10_or.Q+cnt_mod_cie_10_or.C+
          nsx(log(yrs_to_tr_dropout), df = 5)1+
          nsx(log(yrs_to_tr_dropout), df = 5)2+
          nsx(log(yrs_to_tr_dropout), df = 5)3+
          nsx(log(yrs_to_tr_dropout), df = 5)4+
          nsx(log(yrs_to_tr_dropout), df = 5)5+
          as.numeric(con_quien_vive_joel_coup_chil):nsx(log(yrs_to_tr_dropout),
          df = 1)+nsx(log(yrs_to_tr_dropout), df =
          1):as.numeric(con_quien_vive_joel_fam_or)
  Tot Df Deviance   Chisq Df Pr(>Chisq)    
1     35    29416                          
2     38    29347  68.261  3  1.006e-14 ***
3     30    29454 106.494  8  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Show code

# predict(rstpm_4_tvc_4,se.fit=TRUE,type="meanhaz")

system.time(
  pred2 <- predict(rstpm_4_tvc_2,
                 newdata = transform(prueba2_imp22, con_quien_vive_joel_fam_or = 0),
                 grid = TRUE, full = TRUE, se.fit = TRUE, type = "rmst",
                 exposed = function(data) transform(prueba2_imp22, con_quien_vive_joel_fam_or = 1))
  )

Error: no se puede ubicar un vector de tamaño 5.3 Gb

Timing stopped at: 53.8 20.22 75.2

Show code

rstpm_4_tvc_4_hrs<-eform(rstpm_4_tvc_4)
rstpm_4_tvc_2_hrs<-eform(rstpm_4_tvc_2)


#https://biostat3.net/download/R/archive/q28.html
plot(rstpm_4_tvc_2,newdata=data.frame(con_quien_vive_joel_fam_or = 0, sexo_2= "Men", escolaridad_rec= "2-Completed high school or less", sus_ini_mod_mvv= "Alcohol", freq_cons_sus_prin= "Daily", dg_trs_cons_sus_or= "Drug dependence", sus_principal_mod= "Cocaine paste", cnt_mod_cie_10_or= 0),#transform(prueba2_imp22, con_quien_vive_joel_fam_or = 0),#
     type="sdiff", rug=FALSE, line.col=1, ci=T,
     exposed=function(data) transform(data, con_quien_vive_joel_fam_or = 1),
     xlab="Time since diagnosis (years)",ylab="Survivals difference")

Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels): variable lengths differ (found for ‘con_quien_vive_joel_fam_or’)

Time for this code chunk to run: 0.5 minutes

Show code

round(rstpm_4_tvc_2_hrs,2) %>% 
  knitr::kable(size=10, format="markdown",caption= "Hazard ratios", escape=T)

Table 2: Hazard ratios

	exp(beta)	2.5 %	97.5 %
(Intercept)	0.00	0.00	0.00
con_quien_vive_joel_coup_chil	0.03	0.03	0.03
con_quien_vive_joel_fam_or	0.15	0.15	0.15
edad_al_ing_sc	0.99	0.99	1.00
sexo_2Women	0.82	0.79	0.85
escolaridad_rec.L	0.78	0.75	0.80
escolaridad_rec.Q	0.97	0.95	0.99
sus_ini_mod_mvvCocaine hydrochloride	1.10	1.02	1.18
sus_ini_mod_mvvMarijuana	1.05	1.02	1.07
sus_ini_mod_mvvOther	1.02	0.90	1.15
sus_ini_mod_mvvCocaine paste	1.09	1.02	1.17
freq_cons_sus_prin.L	1.06	1.02	1.09
freq_cons_sus_prin.Q	0.96	0.93	0.99
freq_cons_sus_prin.C	0.97	0.93	1.01
freq_cons_sus_prin^4	1.02	0.98	1.05
dg_trs_cons_sus_orHazardous consumption	1.02	0.99	1.05
sus_principal_modCocaine hydrochloride	1.18	1.15	1.22
sus_principal_modMarijuana	0.99	0.94	1.04
sus_principal_modOther	0.87	0.74	1.01
sus_principal_modCocaine paste	1.26	1.24	1.29
cnt_mod_cie_10_or.L	0.68	0.66	0.70
cnt_mod_cie_10_or.Q	0.80	0.78	0.83
cnt_mod_cie_10_or.C	1.27	1.24	1.31
nsx(log(yrs_to_tr_dropout), df = 4)1	236.10	231.47	240.48
nsx(log(yrs_to_tr_dropout), df = 4)2	202.63	200.93	204.11
nsx(log(yrs_to_tr_dropout), df = 4)3	5783.82	5623.76	5966.64
nsx(log(yrs_to_tr_dropout), df = 4)4	152.65	151.68	153.79
as.numeric(con_quien_vive_joel_coup_chil):nsx(log(yrs_to_tr_dropout), df = 2)1	373.27	356.10	390.30
as.numeric(con_quien_vive_joel_coup_chil):nsx(log(yrs_to_tr_dropout), df = 2)2	5.03	4.94	5.12
nsx(log(yrs_to_tr_dropout), df = 2)1:as.numeric(con_quien_vive_joel_fam_or)	21.72	21.09	22.30
nsx(log(yrs_to_tr_dropout), df = 2)2:as.numeric(con_quien_vive_joel_fam_or)	2.33	2.31	2.36

Show code

invisible("predict does not work")
# 
# #sexo_2= Men
# #escolaridad_rec= 2-Completed high school or less
# #sus_ini_mod_mvv= Alcohol
# #freq_cons_sus_prin= Daily
# #dg_trs_cons_sus_or= Drug dependence
# #sus_principal_mod= Cocaine paste 
# #cnt_mod_cie_10_or= 0
# 
# predict(rstpm_4_tvc_4, newdata= expand.grid(edad_al_ing_sc= 0, sexo_2= "Men", escolaridad_rec= "2-Completed high school or less", sus_ini_mod_mvv= "Alcohol", freq_cons_sus_prin= "Daily", dg_trs_cons_sus_or= "Drug dependence", sus_principal_mod= "Cocaine paste", con_quien_vive_joel=c("Alone"), yrs_to_dropout=3),#
# type="rmst", se.fit=TRUE, tau=3)
# 
# 
# fiform_intonly<-
# as.formula(paste("Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel")) #rcs(edad_al_ing_1,4) +
# fiform<-
# as.formula(paste("Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel + ", paste(vars_cov, collapse = "+"))) #rcs(edad_al_ing_1,4) +
# 
# #para saber qué era la variable
# round(attr(scale(prueba2_imp22$edad_al_ing_sc,scale=F),"scaled:center"),4)

Time for this code chunk to run: 0 minutes

Show code

require(flexsurv)
n_trans <- 1
fits_wei <- vector(mode = "list", length = n_trans)
fits_weiph <- vector(mode = "list", length = n_trans)
fits_llogis <- vector(mode = "list", length = n_trans)
fits_gomp <- vector(mode = "list", length = n_trans)
fits_logn <- vector(mode = "list", length = n_trans)
fits_exp <- vector(mode = "list", length = n_trans)
fits_gam <- vector(mode = "list", length = n_trans)
fits_ggam <- vector(mode = "list", length = n_trans)
fits_genf <- vector(mode = "list", length = n_trans)
fits_genf_orig <- vector(mode = "list", length = n_trans)
fits_ggam_orig <- vector(mode = "list", length = n_trans)
fits_rp02 <- vector(mode = "list", length = n_trans)
fits_rp03 <- vector(mode = "list", length = n_trans)
fits_rp04 <- vector(mode = "list", length = n_trans)
fits_rp05 <- vector(mode = "list", length = n_trans)
fits_rp06 <- vector(mode = "list", length = n_trans)
fits_rp07 <- vector(mode = "list", length = n_trans)
fits_rp08 <- vector(mode = "list", length = n_trans)
fits_rp09 <- vector(mode = "list", length = n_trans)
fits_rp10 <- vector(mode = "list", length = n_trans)
km.lc <-vector(mode = "list", length = n_trans)


#Specify the formula
fitform <- Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel # ~ 1 #con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or #Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+ 


for (i in 1:n_trans){  
  fits_wei[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "weibull")
}

for (i in 1:n_trans){  
  fits_weiph[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "weibullph")
}

for (i in 1:n_trans){
  fits_llogis[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "llogis")
}

for (i in 1:n_trans){
  fits_gam[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "gamma")
}
#LOS ERRORES OCURREN CUANDO NO HAY COVARIABLES
#In (function (q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,  ... :   NaNs produced
#gamma no funcionó: NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced
for (i in 1:n_trans){
  fits_ggam[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "gengamma")
}

for (i in 1:n_trans){
  fits_gomp[[i]] <- flexsurvreg(formula=fitform,
                                data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                                dist = "gompertz")
}


for (i in 1:n_trans){
  fits_logn[[i]] <- flexsurvreg(formula=fitform,
                                data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                                dist = "lnorm")
}

for (i in 1:n_trans){
  fits_exp[[i]] <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "exp")
}

no_attempts <- 20

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "gengamma.orig")
      )
    }
    fits_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "genf.orig")
      )
    }
    fits_genf_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <-  flexsurvreg(formula=fitform,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), # subset(ms_d_match_surv, trans == i),
                               dist = "genf")
      )
    }
    fits_genf[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=1,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp02[[i]] <- r
}  

Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito
Error in optim(control = list(fnscale = -14000), method = "BFGS", par = c(gamma0 = 0.419784254753456,  : 
  valor inicial en 'vmmin' no es finito

Show code

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=2,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=3,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=4,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=5,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=6,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=7,
                               data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) # subset(ms_d_match_surv, trans == i),
      )
    }
    fits_rp08[[i]] <- r
}

km.lc[[1]] <- survfit(formula= fitform, data = prueba2_imp22%>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame())

Time for this code chunk to run: 0.2 minutes

Show code

get_distinct_hues <- function(ncolor,s=0.5,v=0.95,seed=2125) {
  golden_ratio_conjugate <- 0.618033988749895
  set.seed(seed)
  h <- runif(1)
  H <- vector("numeric",ncolor)
  for(i in seq_len(ncolor)) {
    h <- (h + golden_ratio_conjugate) %% 1
    H[i] <- h
  }
  hsv(H,s=s,v=v)
}

distinct_hues<-get_distinct_hues(18)
#rp02 no converge
plot(km.lc[[1]], col="red", conf.int=F);
   lines(fits_wei[[1]], col=distinct_hues[1], ci=F);
   lines(fits_weiph[[1]], col=distinct_hues[2], ci=F);
   lines(fits_gomp[[1]], col=distinct_hues[3], ci=F);
   lines(fits_llogis[[1]], col=distinct_hues[4], ci=F);#A0A36D
   lines(fits_gam[[1]], col=distinct_hues[5], ci=F);
   lines(fits_ggam[[1]], col=distinct_hues[6], ci=F);
   lines(fits_ggam_orig[[1]], col=distinct_hues[7], ci=F);
   lines(fits_logn[[1]], col=distinct_hues[8], ci=F);
   lines(fits_exp[[1]],col=distinct_hues[9], ci=F);
   lines(fits_genf[[1]],col=distinct_hues[10], ci=F, lty = 2);
   lines(fits_genf_orig[[1]],col=distinct_hues[11], ci=F, lty = 2);  
  # lines(fits_rp02[[1]],col=distinct_hues[12], ci=F, lty = 2);
   lines(fits_rp03[[1]],col=distinct_hues[13], ci=F, lty = 2);
   lines(fits_rp04[[1]],col=distinct_hues[14], ci=F, lty = 2);
   lines(fits_rp05[[1]],col=distinct_hues[15], ci=F, lty = 2);
   lines(fits_rp06[[1]],col=distinct_hues[16], ci=F, lty = 2);
   lines(fits_rp07[[1]],col=distinct_hues[17], ci=F, lty = 3);
   lines(fits_rp08[[1]],col=distinct_hues[18], ci=F, lty = 3)
  
legend("topright", legend = c("Kaplan-Meier", "Weibull (AFT)", "Weibull (PH)","Gompertz", "Log-logistic", "Gamma",                 "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)", "RP3", "RP4", "RP5", "RP6", "RP7", "RP8"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .45, bty = "n", lty=c(rep(1,12),rep(2,4),rep(3,2)),lwd=3)# rep(3,4)

Time for this code chunk to run: 0.1 minutes

Show code

fits_objects = c("fits_wei", "fits_weiph", "fits_gomp", "fits_llogis", "fits_gam", "fits_ggam", "fits_ggam_orig", "fits_logn", "fits_exp",  "fits_genf", "fits_genf_orig", "fits_rp02", "fits_rp03", "fits_rp04", "fits_rp05", "fits_rp06", "fits_rp07", "fits_rp08")

for( i in 1:length(setdiff(fits_objects,"fits_rp02"))){
    print(setdiff(fits_objects,"fits_rp02")[i])
    x<- get(setdiff(fits_objects,"fits_rp02")[i])
    x<-x[[1]]
    print(x$AIC)
}

[1] "fits_wei"
[1] 50183.62
[1] "fits_weiph"
[1] 50183.62
[1] "fits_gomp"
[1] 35043.79
[1] "fits_llogis"
[1] 44145.69
[1] "fits_gam"
[1] 53389.48
[1] "fits_ggam"
[1] 41339.6
[1] "fits_ggam_orig"
[1] 44903.69
[1] "fits_logn"
[1] 44090.38
[1] "fits_exp"
[1] 62918.94
[1] "fits_genf"
[1] 36041.8
[1] "fits_genf_orig"
[1] 36028.8
[1] "fits_rp03"
[1] 251896.7
[1] "fits_rp04"
[1] 197933.6
[1] "fits_rp05"
[1] 207732.2
[1] "fits_rp06"
[1] 201165.6
[1] "fits_rp07"
[1] 201580
[1] "fits_rp08"
[1] 200765.3

Time for this code chunk to run: 0 minutes

Show code

B_inter_flexsurv = 250
limit=500
#detach("package:flexsurv", unload=TRUE)
require(flexsurv)
vector <- c("con_quien_vive_joel_coup_chil", "con_quien_vive_joel_fam_or", "edad_al_ing_sc", "sexo_2", 
            "escolaridad_rec", "sus_ini_mod_mvv", "freq_cons_sus_prin", "dg_trs_cons_sus_or", 
            "sus_principal_mod", "cnt_mod_cie_10_or")
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
set.seed(2125)
fits_c_gomp <- 
flexsurvreg(formula=Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+  edad_al_ing_sc+ sexo_2+ escolaridad_rec+ sus_ini_mod_mvv+ freq_cons_sus_prin+ dg_trs_cons_sus_or+ sus_principal_mod+ cnt_mod_cie_10_or, data=prueba2_imp22, dist = "gompertz")#%>% slice(sample(1:nrow(prueba2_imp22),limit))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
system.time(
ss.gomp.rmst_noboot <- standsurv(fits_c_gomp,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "ratio"
))

   user  system elapsed 
5140.29    2.87 5148.86 

Show code

system.time(
ss.gomp.surv_noboot <- standsurv(fits_c_gomp,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "ratio"
))

   user  system elapsed 
 159.69    5.30  165.14 

Show code

system.time(
ss.gomp.rmst_noboot_diff <- standsurv(fits_c_gomp,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "difference"
))

   user  system elapsed 
5154.82    2.48 5163.52 

Show code

system.time(
ss.gomp.surv_noboot_diff <- standsurv(fits_c_gomp,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "difference"
))

   user  system elapsed 
 173.44    2.13  175.83 

Show code

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
system.time(
ss.gomp.rmst_boot <- standsurv(fits_c_gomp,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "ratio"
))

    user   system  elapsed 
19495.84    10.19 19523.44 

Show code

system.time(
ss.gomp.surv_boot <- standsurv(fits_c_gomp,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "ratio"
))

   user  system elapsed 
 443.10    2.59  445.89 

Show code

system.time(
ss.gomp.rmst_boot_diff <- standsurv(fits_c_gomp,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "difference"
))

    user   system  elapsed 
19310.14     8.58 19333.34 

Show code

system.time(
ss.gomp.surv_boot_diff <- standsurv(fits_c_gomp,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "difference"
))

   user  system elapsed 
 441.14    2.27  443.78 

Show code

#no funciona  Error in eigen(Dmat) : infinite or missing values in 'x'
#flexsurvspline(formula=Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+  edad_al_ing_sc+ sexo_2+ escolaridad_rec+ sus_ini_mod_mvv+ freq_cons_sus_prin+ dg_trs_cons_sus_or+ sus_principal_mod+ cnt_mod_cie_10_or, k=5, data = prueba2_imp22[complete.cases(prueba2_imp22[,c(vector,"yrs_to_tr_dropout","event")]),] %>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000))
#https://github.com/chjackson/flexsurv-dev/issues/69
set.seed(2125)
fits_c_rp04 <-
flexsurvspline(formula=Surv(yrs_to_tr_dropout,  event) ~ con_quien_vive_joel_coup_chil+ con_quien_vive_joel_fam_or+  edad_al_ing_sc+ sexo_2+ escolaridad_rec+ sus_ini_mod_mvv+ freq_cons_sus_prin+ dg_trs_cons_sus_or+ sus_principal_mod+ cnt_mod_cie_10_or, k=3, data = prueba2_imp22[complete.cases(prueba2_imp22[,c(vector,"yrs_to_tr_dropout","event")]),] %>% dplyr::mutate(yrs_to_tr_dropout=ifelse(yrs_to_tr_dropout==0,.001,yrs_to_tr_dropout)) %>% data.frame(), control = list(fnscale = -14000)) #%>% slice(sample(1:nrow(prueba2_imp22),(limit*10)))
#https://github.com/chjackson/flexsurv-dev/issues/69

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
system.time(
ss.rp04.rmst_noboot <- standsurv(fits_c_rp04,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "ratio"
))

    user   system  elapsed 
47882.49    15.24 47933.11 

Show code

system.time(
ss.rp04.surv_noboot <- standsurv(fits_c_rp04,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "ratio"
))

   user  system elapsed 
 210.78    4.01  215.39 

Show code

system.time(
ss.rp04.rmst_noboot_diff <- standsurv(fits_c_rp04,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "difference"
))

    user   system  elapsed 
48085.92    17.08 48172.92 

Show code

system.time(
ss.rp04.surv_noboot_diff <- standsurv(fits_c_rp04,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = FALSE,
                seed = 2125,
                contrast = "difference"
))

   user  system elapsed 
 213.94    1.19  215.38 

Show code

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
system.time(
ss.rp04.rmst_boot <- standsurv(fits_c_rp04,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "ratio"
))

     user    system   elapsed 
169975.91     51.78 170140.53 

Show code

system.time(
ss.rp04.surv_boot <- standsurv(fits_c_rp04,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "ratio"
))

   user  system elapsed 
3016.66    4.43 3023.67 

Show code

system.time(
ss.rp04.rmst_boot_diff <- standsurv(fits_c_rp04,
                type = "rmst",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "difference"
))

     user    system   elapsed 
172055.80     61.48 172289.22 

Show code

system.time(
ss.rp04.surv_boot_diff <- standsurv(fits_c_rp04,
                type = "survival",
                at = list(list(con_quien_vive_joel_coup_chil = 1),
                          list(con_quien_vive_joel_fam_or = 1),
                          list(con_quien_vive_joel_fam_or = 0, con_quien_vive_joel_coup_chil= 0)),
                t = seq(0,3,length=25),
                ci = TRUE,
                boot = T,
                B = B_inter_flexsurv,
                seed = 2125,
                contrast = "difference"
))

   user  system elapsed 
3063.64    4.97 3072.59 

Show code

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

Time for this code chunk to run: 0.1 minutes

RMSTs, flexsurv, Gompertz

Show code

invisible("RMSTs")
ss.gomp.rmst_noboot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Differences, Delta)")

Show code

ss.gomp.rmst_boot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Differences, Bootstrap)")

Show code

ss.gomp.rmst_noboot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Delta)")

Show code

ss.gomp.rmst_boot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
ggtitle("RMSTs of dropout (Bootstrap)")

Time for this code chunk to run: 0.1 minutes

Surv, flexsurv, Gompertz

Show code

ss.gomp.surv_noboot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Differences, Delta)")

Show code

ss.gomp.surv_boot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Differences, Bootstrap)")

Show code

ss.gomp.surv_noboot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Delta)")

Show code

ss.gomp.surv_boot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Bootstrap)")

Time for this code chunk to run: 0 minutes

RMSTs, flexsurv, RP04

Show code

invisible("RMSTs")
ss.gomp.rmst_noboot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Differences, Delta)")

Show code

ss.gomp.rmst_boot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Differences, Bootstrap)")

Show code

ss.gomp.rmst_noboot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("RMSTs of dropout (Delta)")

Show code

ss.gomp.rmst_boot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
ggtitle("RMSTs of dropout (Bootstrap)")

Time for this code chunk to run: 0.1 minutes

Surv, flexsurv, RP04

Show code

ss.gomp.surv_noboot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Differences, Delta)")

Show code

ss.gomp.surv_boot_diff %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=0), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Differences, Bootstrap)")

Show code

ss.gomp.surv_noboot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Delta)")

Show code

ss.gomp.surv_boot %>% 
  dplyr::rename_with(~ paste0(.x, "_est"), .cols =  names(.)[!str_ends(names(.), "time|_lci|_uci")]) %>% 
  dplyr::rename_with(~ sub("2_1","21",.x), .cols =  names(.)[grepl("2_1",names(.))]) %>% 
  dplyr::rename_with(~ sub("3_1","31",.x), .cols =  names(.)[grepl("3_1",names(.))]) %>% 
  pivot_longer(!time, names_to=c("at","est"), names_sep="_", values_to="rmst") %>% 
  pivot_wider(id_cols=c("time", "at"), names_from="est", values_from="rmst") %>% 
  dplyr::mutate(at=factor(at, levels=c(paste0("at",1:3),paste0("contrast",3:2,"1")), labels=c("With couple/children", "Alone", "Family of origin", "Alone vs. Couple/children","Fam of origin vs. Alone"))) %>% 
  ggplot()+
  geom_line(aes(x=time, y=est, fill=at))+
  geom_ribbon(aes(x=time, ymin=lci, ymax=uci, fill=at), alpha=.3)+
  geom_hline(aes(yintercept=1), linewidth=1)+
  theme_classic()+
  ggtitle("Survival probs. of dropout (Bootstrap)")

Time for this code chunk to run: 0.1 minutes

Inverse probability weighted Kaplan-meier

Show code

require(landest)
km_alone_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 

km_alone_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=1, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 

km_alone_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=3, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 




km_famorigin_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 

km_famorigin_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=1, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 

km_famorigin_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=3, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 





km_coupchil_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 

km_coupchil_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=1, ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 

km_coupchil_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=3,  ps.weights = prueba2_imp22$weight_mult_a00[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 







km2_alone_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 

km2_alone_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=1, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 

km2_alone_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Alone")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Alone")], tt=3, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Alone")], conf.int=T) 




km2_famorigin_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 

km2_famorigin_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=1, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 

km2_famorigin_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], tt=3, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="Family of origin")], conf.int=T) 





km2_coupchil_t05<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=.5, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 

km2_coupchil_t1<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=1, ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 

km2_coupchil_t3<-
surv.iptw.km(tl=prueba2_imp22$yrs_to_tr_dropout[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], dl = prueba2_imp22$event[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], tt=3,  ps.weights = prueba2_imp22$weight_mult_a0[which(prueba2_imp22$con_quien_vive_joel=="With couple/children")], conf.int=T) 


get_est <- function(mat_res){
    LCI <- mat_res$conf.int.normal.S[[1]]
  mean <- mat_res$S.estimate
  UCI <- mat_res$conf.int.normal.S[[2]]
  mean_CI <- cbind.data.frame(mean, LCI, UCI)
  return(mean_CI)
}
 
#km2_alone_t05 km2_alone_t1  km2_alone_t3 km2_famorigin_t05 km2_famorigin_t1 km2_famorigin_t3 km2_coupchil_t05 km2_coupchil_t1 km2_coupchil_t3
get_est_km_noboot_05_3<- #cbind.data.frame(group="",time="",mean="",LCI="",UCI=""), 
rbind.data.frame(cbind.data.frame(group="Alone", time="0.5", get_est(km_alone_t05)), 
                 cbind.data.frame(group="Alone", time="1", get_est(km_alone_t1)),
                 cbind.data.frame(group="Alone", time="3", get_est(km_alone_t3)), 
                 cbind.data.frame(group="Family of origin", time="0.5", get_est(km_famorigin_t05)),
                 cbind.data.frame(group="Family of origin", time="1", get_est(km_famorigin_t1)),
                 cbind.data.frame(group="Family of origin", time="3", get_est(km_famorigin_t3)),
                 cbind.data.frame(group="With couple/children", time="0.5", get_est(km_coupchil_t05)),
                 cbind.data.frame(group="With couple/children", time="1", get_est(km_coupchil_t1)),
                 cbind.data.frame(group="With couple/children", time="3", get_est(km_coupchil_t3))
                  )

get_est2_km_noboot_05_3<- #cbind.data.frame(group="",time="",mean="",LCI="",UCI=""), 
rbind.data.frame(cbind.data.frame(group="Alone", time="0.5", get_est(km2_alone_t05)), 
                 cbind.data.frame(group="Alone", time="1", get_est(km2_alone_t1)),
                 cbind.data.frame(group="Alone", time="3", get_est(km2_alone_t3)), 
                 cbind.data.frame(group="Family of origin", time="0.5", get_est(km2_famorigin_t05)),
                 cbind.data.frame(group="Family of origin", time="1", get_est(km2_famorigin_t1)),
                 cbind.data.frame(group="Family of origin", time="3", get_est(km2_famorigin_t3)),
                 cbind.data.frame(group="With couple/children", time="0.5", get_est(km2_coupchil_t05)),
                 cbind.data.frame(group="With couple/children", time="1", get_est(km2_coupchil_t1)),
                 cbind.data.frame(group="With couple/children", time="3", get_est(km2_coupchil_t3))
                  )

get_est_km_noboot_05_3 %>% 
  #dplyr::mutate_if(c("mean", "LCI", "UCI"),~round(readr::parse_number(.,2)) %>% 
  knitr::kable("markdown", digits=2, caption="Kaplan meier (IPW Multinomial logistic), Survival Estimates (95%CI no bootstrap)")

Table 3: Kaplan meier (IPW Multinomial logistic), Survival Estimates (95%CI no bootstrap)

group	time	mean	LCI	UCI
Alone	0.5	0.46	0.43	0.49
Alone	1	0.32	0.29	0.35
Alone	3	0.28	0.26	0.31
Family of origin	0.5	0.53	0.52	0.53
Family of origin	1	0.36	0.35	0.37
Family of origin	3	0.31	0.30	0.32
With couple/children	0.5	0.49	0.48	0.50
With couple/children	1	0.32	0.30	0.33
With couple/children	3	0.26	0.25	0.27

Time for this code chunk to run: 0.1 minutes

Show code

get_est2_km_noboot_05_3 %>% 
  #dplyr::mutate_if(c("mean", "LCI", "UCI"),~round(readr::parse_number(.,2)) %>% 
  knitr::kable("markdown", digits=2, caption="Kaplan meier (IPW BART), Survival Estimates (95%CI no bootstrap)")

Table 4: Kaplan meier (IPW BART), Survival Estimates (95%CI no bootstrap)

group	time	mean	LCI	UCI
Alone	0.5	0.46	0.43	0.49
Alone	1	0.32	0.29	0.35
Alone	3	0.28	0.25	0.31
Family of origin	0.5	0.53	0.52	0.54
Family of origin	1	0.36	0.35	0.37
Family of origin	3	0.31	0.30	0.32
With couple/children	0.5	0.49	0.47	0.50
With couple/children	1	0.31	0.30	0.33
With couple/children	3	0.26	0.25	0.27

Show code

#knitr::kable(digits=2) %>%   kableExtra::kable_classic()

Time for this code chunk to run: 0 minutes

Multinomial logistic

Show code

require(boot) 

B_iter <- 2000

invisible("Alone")
kmfit05alone <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05alone <- boot(data = prueba2_imp22, statistic = kmfit05alone, R = B_iter)

kmfit1alone <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1alone <- boot(data = prueba2_imp22, statistic = kmfit1alone, R = B_iter)

kmfit3alone <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3alone <- boot(data = prueba2_imp22, statistic = kmfit3alone, R = B_iter)


invisible("Family of origin")

kmfit05fam <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05fam <- boot(data = prueba2_imp22, statistic = kmfit05fam, R = B_iter)

kmfit1fam <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1fam <- boot(data = prueba2_imp22, statistic = kmfit1fam, R = B_iter)

kmfit3fam <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3fam <- boot(data = prueba2_imp22, statistic = kmfit3fam, R = B_iter)


invisible("Couple & children")

kmfit05coupchil <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05coupchil <- boot(data = prueba2_imp22, statistic = kmfit05coupchil, R = B_iter)

kmfit1coupchil <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1coupchil <- boot(data = prueba2_imp22, statistic = kmfit1coupchil, R = B_iter)

kmfit3coupchil <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a00)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3coupchil <- boot(data = prueba2_imp22, statistic = kmfit3coupchil, R = B_iter)


#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

#Función para obtener intervalos de confianza al 95% del remuestreo
get_estimation<-function(mat_res){
  LCI <- (boot.ci(mat_res, type = "basic"))[[4]][[4]]
  mean <- as.numeric((boot.ci(mat_res, type = "basic"))[[2]])
  UCI <- (boot.ci(mat_res, type = "basic"))[[4]][[5]]
  mean_CI <- cbind.data.frame(mean, LCI, UCI)
  return(mean_CI)
}

# invisible("replicate times")
# design_weights <- survey::svydesign(ids = ~1, data = dplyr::mutate(prueba2_imp22, yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout)), weights = ~weight_mult_a00)
# z<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, design=design_weights , rho=0, gamma=0)

get_est_km_05_3<-
rbind.data.frame(cbind.data.frame(group="Alone", time="0.5", get_estimation(bootapp05alone)), 
                 cbind.data.frame(group="Alone", time="1", get_estimation(bootapp1alone)),
                 cbind.data.frame(group="Alone", time="3", get_estimation(bootapp3alone)), 
                 cbind.data.frame(group="Family of origin", time="0.5", get_estimation(bootapp05fam)),
                 cbind.data.frame(group="Family of origin", time="1", get_estimation(bootapp1fam)),
                 cbind.data.frame(group="Family of origin", time="3", get_estimation(bootapp3fam)),
                 cbind.data.frame(group="With couple/children", time="0.5", get_estimation(bootapp05coupchil)),
                 cbind.data.frame(group="With couple/children", time="1", get_estimation(bootapp1coupchil)),
                 cbind.data.frame(group="With couple/children", time="3", get_estimation(bootapp3coupchil))
                  )

get_est_km_05_3 %>% 
  knitr::kable("markdown", digits=2, caption="Kaplan meier (IPW Multinomial logistic), Survival Estimates (Bootstrap 95%CI)")

Table 5: Kaplan meier (IPW Multinomial logistic), Survival Estimates (Bootstrap 95%CI)

group	time	mean	LCI	UCI
Alone	0.5	0.46	0.43	0.49
Alone	1	0.32	0.29	0.35
Alone	3	0.28	0.25	0.31
Family of origin	0.5	0.53	0.52	0.54
Family of origin	1	0.36	0.36	0.37
Family of origin	3	0.31	0.30	0.32
With couple/children	0.5	0.49	0.47	0.50
With couple/children	1	0.32	0.30	0.33
With couple/children	3	0.26	0.25	0.27

Time for this code chunk to run: 0.9 minutes

Bart

Show code

invisible("Alone")
kmfit05alone2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05alone2 <- boot(data = prueba2_imp22, statistic = kmfit05alone2, R = B_iter)

kmfit1alone2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1alone2 <- boot(data = prueba2_imp22, statistic = kmfit1alone2, R = B_iter)

kmfit3alone2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$Alone$time,s$Alone$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3alone2 <- boot(data = prueba2_imp22, statistic = kmfit3alone2, R = B_iter)


invisible("Family of origin")

kmfit05fam2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05fam2 <- boot(data = prueba2_imp22, statistic = kmfit05fam2, R = B_iter)

kmfit1fam2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1fam2 <- boot(data = prueba2_imp22, statistic = kmfit1fam2, R = B_iter)

kmfit3fam2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`Family of origin`$time,s$`Family of origin`$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3fam2 <- boot(data = prueba2_imp22, statistic = kmfit3fam2, R = B_iter)


invisible("Couple & children")

kmfit05coupchil2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.49 & V1<0.51) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp05coupchil2 <- boot(data = prueba2_imp22, statistic = kmfit05coupchil2, R = B_iter)

kmfit1coupchil2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>0.9 & V1<1.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp1coupchil2 <- boot(data = prueba2_imp22, statistic = kmfit1coupchil2, R = B_iter)

kmfit3coupchil2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3coupchil2 <- boot(data = prueba2_imp22, statistic = kmfit3coupchil2, R = B_iter)

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

invisible("replicate times")
design_weights <- survey::svydesign(ids = ~1, data = dplyr::mutate(prueba2_imp22, yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout)), weights = ~weight_mult_a0)
z<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, design=design_weights , rho=0, gamma=0)

get_est_km_05_3_2<-
rbind.data.frame(cbind.data.frame(group="Alone", time="0.5", get_estimation(bootapp05alone2)), 
                 cbind.data.frame(group="Alone", time="1", get_estimation(bootapp1alone2)),
                 cbind.data.frame(group="Alone", time="3", get_estimation(bootapp3alone2)), 
                 cbind.data.frame(group="Family of origin", time="0.5", get_estimation(bootapp05fam2)),
                 cbind.data.frame(group="Family of origin", time="1", get_estimation(bootapp1fam2)),
                 cbind.data.frame(group="Family of origin", time="3", get_estimation(bootapp3fam2)),
                 cbind.data.frame(group="With couple/children", time="0.5", get_estimation(bootapp05coupchil2)),
                 cbind.data.frame(group="With couple/children", time="1", get_estimation(bootapp1coupchil2)),
                 cbind.data.frame(group="With couple/children", time="3", get_estimation(bootapp3coupchil2))
                  )

get_est_km_05_3_2 %>% 
  knitr::kable("markdown", digits=2, caption="Kaplan meier (IPW Multinomial logistic), Survival Estimates (95%CI)")

Table 6: Kaplan meier (IPW Multinomial logistic), Survival Estimates (95%CI)

group	time	mean	LCI	UCI
Alone	0.5	0.46	0.43	0.49
Alone	1	0.32	0.29	0.35
Alone	3	0.28	0.25	0.31
Family of origin	0.5	0.53	0.52	0.54
Family of origin	1	0.36	0.36	0.37
Family of origin	3	0.31	0.30	0.32
With couple/children	0.5	0.49	0.47	0.50
With couple/children	1	0.32	0.30	0.33
With couple/children	3	0.26	0.25	0.27

Time for this code chunk to run: 0.9 minutes

Show code

kmfit3coupchil2 <- function(data, i){
      dat <- dplyr::mutate(data[i,], yrs_to_tr_dropout= dplyr::case_when(yrs_to_tr_dropout==0~ 0.0001, T~yrs_to_tr_dropout))
      set.seed(2125)
      kmlist<-list()
      design_weights <- survey::svydesign(ids = ~1, data = dat, weights = ~weight_mult_a0)

      s<-survey::svykm(Surv(time=yrs_to_tr_dropout,event=event)~con_quien_vive_joel, se=F,
                  #con_quien_vive_joel, 
                  design=design_weights, rho=0, gamma=0) 
      
      kmlist<-cbind(s$`With couple/children`$time,s$`With couple/children`$surv) %>% as_tibble() %>% dplyr::filter(V1>2.9 & V1<3.1) %>% summarise(mean=mean(V2)) %>% unlist()
      return(kmlist)
}
set.seed(2125)
bootapp3coupchil2 <- boot(data = prueba2_imp22, statistic = kmfit3coupchil2, R = B_iter)

 mestr <- 
   geex::m_estimate(estFUN = kmfit3coupchil2,  # Function of estimating functions
   data = dat, # Data to be used (must be data frame)
   root_control = setup_root_control(start = rep(c(0),145)))   # Set starting values

Time for this code chunk to run: 0 minutes

Adjusted RMSTs

Define function

Multinomial logistic

Show code

#TABLE OF RMSTS
source("https://raw.githubusercontent.com/s-conner/akm-rmst/master/AKM_rmst.R")

cat("Multinomial logistic, stabilized weights but not truncated?")
RMST11<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a00, tau=.5)

Show code

RMST12<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a00, tau=1)

Show code

RMST13<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a00, tau=3)

Time for this code chunk to run: 0 minutes

Show code

rmst_se_ind<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11$rmst_group,RMST12$rmst_group, RMST13$rmst_group)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::mutate(RMST= round(RMST,2))%>% 
  dplyr::select(-SE) %>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Group ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.))) 
rmst_se_ind_diff<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11$rmst_diff,RMST12$rmst_diff, RMST13$rmst_diff)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::select(-SE) %>%
  dplyr::mutate_if(is.numeric, ~round(.,2)) %>% 
  tidyr::separate(rn, into = c("rn1", "rn2"), sep = "vs.") %>% 
  dplyr::mutate(diff=paste0(`Est.`," (",CIL,", ", CIU,")")) %>% 
  dplyr::select(tt, rn1, rn2, diff)%>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Groups ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+"))) %>% 
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.)))
rmst_se_ind_ratio<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11$rmst_ratio,RMST12$rmst_ratio, RMST13$rmst_ratio)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::select(-SE, -`Log Est.`) %>% 
  dplyr::mutate_if(is.numeric, ~round(.,2)) %>% 
  tidyr::separate(rn, into = c("rn1", "rn2"), sep = "vs.") %>% 
  dplyr::mutate(ratio=paste0(`Est.`," (",CIL,", ", CIU,")")) %>% 
  dplyr::select(tt, rn1, rn2, ratio)%>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Groups ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+"))) %>% 
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.)))

rmst_se_ind_diff %>% 
  dplyr::left_join(rmst_se_ind,by=c("tt"="tt", "rn1"="rn")) %>% 
  dplyr::left_join(rmst_se_ind,by=c("tt"="tt", "rn2"="rn")) %>% 
  dplyr::select(tt, rn1, `RMST.x`, rn2, `RMST.y`, diff) %>% 
  dplyr::left_join(rmst_se_ind_ratio,by=c("tt"="tt", "rn1"="rn1", "rn2"="rn2")) %>% 
#  knitr::kable(format= "markdown", caption="Restricted Mean Survival Time (in years)", col.names = )
 knitr::kable(format = 'markdown',
        escape = FALSE,
        col.names= c("Time", "Cat1","RMST","Cat2","RMST","diff","ratio"), 
     caption="Restricted Mean Survival Time (in years)") #%>% 

Table 7: Restricted Mean Survival Time (in years)

Time	Cat1	RMST	Cat2	RMST	diff	ratio
Six months	Family of origin	0.39	Alone	0.36	0.03 (0.02, 0.04)	1.07 (1.04, 1.1)
Six months	With couple/children	0.38	Alone	0.36	0.01 (0, 0.03)	1.04 (1.01, 1.07)
Six months	With couple/children	0.38	Family of origin	0.39	-0.01 (-0.01, -0.01)	0.97 (0.96, 0.98)
One year	Family of origin	0.60	Alone	0.55	0.05 (0.03, 0.07)	1.09 (1.05, 1.14)
One year	With couple/children	0.57	Alone	0.55	0.02 (-0.01, 0.04)	1.03 (0.99, 1.07)
One year	With couple/children	0.57	Family of origin	0.60	-0.03 (-0.04, -0.02)	0.94 (0.93, 0.96)
Three years	Family of origin	1.24	Alone	1.13	0.11 (0.03, 0.19)	1.1 (1.02, 1.17)
Three years	With couple/children	1.11	Alone	1.13	-0.02 (-0.1, 0.06)	0.98 (0.91, 1.05)
Three years	With couple/children	1.11	Family of origin	1.24	-0.13 (-0.17, -0.1)	0.89 (0.86, 0.92)

Show code

    # kableExtra::kable_styling(font_size = 27) %>%
    # gsub("font-size: initial !important;", 
    #      "font-size: 30pt !important;", 
    #      .)

Time for this code chunk to run: 0.1 minutes

Show code

#TABLE OF RMSTS

cat("Multinomial logistic, stabilized weights but not truncated?")
RMST11b<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a0, tau=.5)

Show code

RMST12b<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a0, tau=1)

Show code

RMST13b<-
akm_rmst_2(time=prueba2_imp22$yrs_to_tr_dropout, status=prueba2_imp22$event, group=prueba2_imp22$con_quien_vive_joel,weight=prueba2_imp22$weight_mult_a0, tau=3)

Time for this code chunk to run: 0 minutes

BART

Show code

rmst_se_ind_b<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11b$rmst_group,RMST12b$rmst_group, RMST13b$rmst_group)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::mutate(RMST= round(RMST,2))%>% 
  dplyr::select(-SE) %>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Group ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.))) 
rmst_se_ind_diff_b<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11b$rmst_diff,RMST12b$rmst_diff, RMST13b$rmst_diff)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::select(-SE) %>%
  dplyr::mutate_if(is.numeric, ~round(.,2)) %>% 
  tidyr::separate(rn, into = c("rn1", "rn2"), sep = "vs.") %>% 
  dplyr::mutate(diff=paste0(`Est.`," (",CIL,", ", CIU,")")) %>% 
  dplyr::select(tt, rn1, rn2, diff)%>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Groups ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+"))) %>% 
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.)))
rmst_se_ind_ratio_b<-
cbind(tt=c(rep("Six months",3),rep("One year",3),rep("Three years",3)),rbind(RMST11b$rmst_ratio,RMST12b$rmst_ratio, RMST13b$rmst_ratio)) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  dplyr::select(tt, rn, everything()) %>% 
  dplyr::select(-SE, -`Log Est.`) %>% 
  dplyr::mutate_if(is.numeric, ~round(.,2)) %>% 
  tidyr::separate(rn, into = c("rn1", "rn2"), sep = "vs.") %>% 
  dplyr::mutate(ratio=paste0(`Est.`," (",CIL,", ", CIU,")")) %>% 
  dplyr::select(tt, rn1, rn2, ratio)%>% 
    dplyr::mutate(across(starts_with("rn"), ~str_remove(., "Groups ")))%>%
  dplyr::mutate(across(starts_with("rn"), ~str_remove_all(., "\\d+"))) %>% 
  dplyr::mutate(across(starts_with("rn"), ~str_trim(.)))

rmst_se_ind_diff_b %>% 
  dplyr::left_join(rmst_se_ind_b,by=c("tt"="tt", "rn1"="rn")) %>% 
  dplyr::left_join(rmst_se_ind_b,by=c("tt"="tt", "rn2"="rn")) %>% 
  dplyr::select(tt, rn1, `RMST.x`, rn2, `RMST.y`, diff) %>% 
  dplyr::left_join(rmst_se_ind_ratio_b,by=c("tt"="tt", "rn1"="rn1", "rn2"="rn2")) %>% 
#  knitr::kable(format= "markdown", caption="Restricted Mean Survival Time (in years)", col.names = )
 kableExtra::kbl(format = 'markdown',
        escape = FALSE,
        col.names= c("Time", "Cat1","RMST","Cat2","RMST","diff","ratio"), 
    caption="Restricted Mean Survival Time (in years) (BART weights)")

Table 8: Restricted Mean Survival Time (in years) (BART weights)

Time	Cat1	RMST	Cat2	RMST	diff	ratio
Six months	Family of origin	0.39	Alone	0.36	0.03 (0.02, 0.04)	1.07 (1.04, 1.1)
Six months	With couple/children	0.38	Alone	0.36	0.01 (0, 0.03)	1.04 (1.01, 1.07)
Six months	With couple/children	0.38	Family of origin	0.39	-0.01 (-0.01, -0.01)	0.97 (0.96, 0.98)
One year	Family of origin	0.60	Alone	0.55	0.05 (0.03, 0.07)	1.09 (1.05, 1.14)
One year	With couple/children	0.57	Alone	0.55	0.02 (-0.01, 0.04)	1.03 (0.99, 1.07)
One year	With couple/children	0.57	Family of origin	0.60	-0.03 (-0.04, -0.02)	0.94 (0.93, 0.96)
Three years	Family of origin	1.24	Alone	1.13	0.11 (0.04, 0.19)	1.1 (1.03, 1.18)
Three years	With couple/children	1.10	Alone	1.13	-0.02 (-0.1, 0.06)	0.98 (0.91, 1.05)
Three years	With couple/children	1.10	Family of origin	1.24	-0.14 (-0.17, -0.1)	0.89 (0.86, 0.92)

Time for this code chunk to run: 0.2 minutes

Session Info

June 29, 2023

Show code

Sys.getenv("R_LIBS_USER")

[1] "C:/Users/CISS Fondecyt/OneDrive/Documentos/R/win-library/4.1"

Show code

rstudioapi::getSourceEditorContext()

Document Context: 
- id:        '2728DA37'
- path:      'C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/joel/analisis_joel22.Rmd'
- contents:  <2091 rows>
Document Selection:
- [822, 13] -- [822, 13]: ''

Show code

save.image("__analisis_joel22.RData")

unlink("*_cache", recursive = T, force = T, expand = TRUE)

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) %>% 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('Paquetes estadísticos utilizados')),
      options=list(
initComplete = htmlwidgets::JS(
      "function(settings, json) {",
      "$(this.api().tables().body()).css({'font-size': '80%'});",
      "}")))

Time for this code chunk to run: 0.1 minutes