Readmission risk among women in women-only versus mixed gender treatment programs for substance use disorder in Chile (3rd step)
25 mayo, 2021
Transition intensities
#note_flexsurv_1d_dwca <-
# paste0("N= ",format(as.numeric(flexsurv_1d_dwca$N), big.mark=","),"; Events= ",format(as.numeric(flexsurv_1d_dwca$events), big.mark=","),"; Censored= ",format(as.numeric(flexsurv_1d_dwca$N - flexsurv_1d_dwca$events), big.mark=","),"; Time at risk= ",format(as.numeric(flexsurv_1d_dwca$trisk), big.mark=","),"; Df= ",format(as.numeric(flexsurv_1d_dwca$npars), big.mark=","), "; AIC= ",format(as.numeric(flexsurv_1d_dwca$AIC), big.mark=","))
dt_coefs_simple_surv<-
data.frame(
real_vars=c("tipo_de_programa_2Women specific", "edad_al_ing_grupos.L", "edad_al_ing_grupos.Q", "edad_al_ing_grupos.C","escolaridad_rec.L", "escolaridad_rec.Q", "sus_principal_modCocaine hydrochloride", "sus_principal_modCocaine paste", "sus_principal_modMarijuana", "sus_principal_modOther", "freq_cons_sus_prin.L", "freq_cons_sus_prin.Q", "freq_cons_sus_prin.C", "freq_cons_sus_prin^4", "compromiso_biopsicosocial.L", "compromiso_biopsicosocial.Q", "tenencia_de_la_vivienda_modOthers", "tenencia_de_la_vivienda_modOwner/Transferred dwellings/Pays Dividends", "tenencia_de_la_vivienda_modRenting", "tenencia_de_la_vivienda_modStays temporarily with a relative", "num_otras_sus_mod.L", "num_otras_sus_mod.Q", "numero_de_hijos_mod_recYes", "tipo_de_plan_resResidential", "comp_statusDischarge without clinical advice", "tipo_de_programa_2Women specific:comp_statusDischarge without clinical advice","time_to_outcome","factor(tipo_de_programa_2)1"),
formal_vars= c('Type of program-Women specific', 'Age at admission to treatment, grouped- 30-39', 'Age at admission to treatment, grouped- 40-49', 'Age at admission to treatment, grouped- 50+', 'Ed. Attainment- Completed high school or less', 'Ed. Attainment- More than high school', 'Primary or main substance- Cocaine hydrochloride', 'Primary or main substance- Cocaine paste', 'Primary or main substance- Marijuana', 'Primary or main substance- Other', 'Consumption frequency of primary or main substance- 2 to 3 days a week', 'Consumption frequency of primary or main substance- 4 to 6 days a week', 'Consumption frequency of primary or main substance- 1 day a week or more', 'Consumption frequency of primary or main substance- Daily', 'Biopsychosocial involvement- 2-Moderate', 'Biopsychosocial involvement- 3-Severe', 'Tenure status of households- Others', 'Tenure status of households- Owner/Transferred dwellings/Pays Dividends', 'Tenure status of households- Renting', 'Tenure status of households- Stays temporarily with a relative', 'Co-occurring SUD- One additional substance', 'Co-occurring SUD- More than one additional substance', 'Have children (Dichotomized)- Yes', 'Setting of Treatment- Residential', 'Cause of Discharge- W/o clinical advice', 'Type of Program x Cause of Discharge', 'Days in baseline treatment','Women-specific program')#'Educational Attainment', tipo_de_plan_resResidential
)
trans_1_3s<-
data.table::data.table(fits_c_gomp[[1]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_2_3s<-
data.table::data.table(fits_c_genf[[2]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_3_3s<-
data.table::data.table(fits_c_ggam[[3]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
coefs_psa_3s<-
trans_1_3s%>%
dplyr::full_join(trans_2_3s,by="formal_vars")%>%
dplyr::full_join(trans_3_3s,by="formal_vars") %>%
dplyr::mutate(formal_vars=factor(formal_vars,levels=c("shape","scale","rate","mu","sigma","Q","P","Women-specific program"))) %>% #, "meanlog", "sdlog"
dplyr::arrange(formal_vars)
trans_1_4s<-
data.table::data.table(fits_c_gomp2[[1]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_2_4s<-
data.table::data.table(fits_c_gomp2[[2]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_3_4s<-
data.table::data.table(fits_c_logn2[[3]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_4_4s<-
data.table::data.table(fits_c_ggam2[[4]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
trans_5_4s<-
data.table::data.table(fits_c_logn2[[5]]$res[,1:3],keep.rownames = T) %>%
dplyr::left_join(dt_coefs_simple_surv, by=c("rn"="real_vars")) %>%
#dplyr::mutate(p.value=ifelse(p.value<.001,"<0.001",as.character(sprintf("%1.3f",p.value)))) %>%
dplyr::mutate(formal_vars=ifelse(is.na(formal_vars),rn, formal_vars)) %>%
dplyr::mutate_at(vars(est, `L95%`,`U95%`),~dplyr::case_when(grepl("tipo_de_programa_2",rn)~exp(.),T~.)) %>%
dplyr::mutate(conf.low=sprintf("%1.2f",`L95%`),
conf.high=sprintf("%1.2f",`U95%`),
statistic=sprintf("%1.2f",est),
`CI 95%`=paste0(conf.low,", ",conf.high)) %>%
dplyr::filter(rn %in% c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","factor(tipo_de_programa_2)1")) %>%
dplyr::select(formal_vars, statistic, `CI 95%`)
coefs_psa_4s<-
dplyr::full_join(trans_1_4s,trans_2_4s,by="formal_vars")%>% #
dplyr::full_join(trans_3_4s,by="formal_vars")%>%
dplyr::full_join(trans_4_4s,by="formal_vars")%>%
dplyr::full_join(trans_5_4s,by="formal_vars")%>%
dplyr::mutate(formal_vars=factor(formal_vars,levels=c("shape","scale","rate","mu","sigma","Q","P", "meanlog", "sdlog","Women-specific program"))) %>% #
dplyr::arrange(formal_vars)
options(knitr.kable.NA = '')
coefs_psa_3s %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 0a. Partitioned Survival Results, Three-states model"),
col.names = c("Actual state", rep(c("Estimate","IC 95%"),3)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Adm -> Ther.Disch" = 2,"Adm -> Readm" = 2, "Ther.Disch -> Readm" = 2))| Actual state | Estimate | IC 95% | Estimate | IC 95% | Estimate | IC 95% |
|---|---|---|---|---|---|---|
| shape | -0.60 | -0.63, -0.57 | ||||
| rate | 0.32 | 0.25, 0.42 | ||||
| mu | 1.76 | 1.41, 2.10 | 5.35 | 4.10, 6.60 | ||
| sigma | 1.67 | 1.54, 1.82 | 2.63 | 2.30, 3.01 | ||
| Q | -2.43 | -2.81, -2.06 | 0.39 | 0.16, 0.62 | ||
| P | 1.67 | 1.20, 2.33 | ||||
| Women-specific program | 1.26 | 1.18, 1.34 | 0.91 | 0.84, 0.98 | 0.90 | 0.68, 1.20 |
| Note. We removed Readmission because it was an absorbing state |
coefs_psa_4s%>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 0b. Partitioned Survival Results, Four-states model"),
col.names = c("Actual state", rep(c("Estimate","IC 95%"),5)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Adm -> Ther.Disch" = 2,"Adm -> Disch.w/oClin.Adv."=2, "Adm -> Readm" = 2, "Ther.Disch -> Readm" = 2, "Disch.w/oClin.Adv.-> Readm"=2))| Actual state | Estimate | IC 95% | Estimate | IC 95% | Estimate | IC 95% | Estimate | IC 95% | Estimate | IC 95% |
|---|---|---|---|---|---|---|---|---|---|---|
| shape | -0.20 | -0.22, -0.18 | -0.89 | -0.92, -0.85 | ||||||
| rate | 0.39 | 0.30, 0.51 | 0.84 | 0.72, 0.99 | ||||||
| mu | 5.35 | 4.10, 6.60 | ||||||||
| sigma | 2.63 | 2.30, 3.01 | ||||||||
| Q | 0.39 | 0.16, 0.62 | ||||||||
| meanlog | 3.16 | 2.60, 3.72 | 4.27 | 3.65, 4.88 | ||||||
| sdlog | 1.31 | 1.26, 1.37 | 3.05 | 2.97, 3.13 | ||||||
| Women-specific program | 1.29 | 1.21, 1.37 | 0.96 | 0.92, 1.00 | 0.93 | 0.83, 1.05 | 0.90 | 0.68, 1.20 | 0.82 | 0.69, 0.96 |
| Note. We removed Readmission because it was an absorbing state |
Cumulative Baseline Hazards of Joint Models
`%>%`<-magrittr::`%>%`
# 3 states -Women-specific
cumhaz_data_a <- rbind(data.frame(cox_cumhaz_0_a_we$Haz, #2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_a$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_a$Haz,
model = "Parametric"))
cumhaz_data_a$trans <- factor(cumhaz_data_a$trans,
levels = seq(1, 3),
labels = c("Adm -> Ther.Disch",
"Adm -> Readm",
"Ther.Disch -> Readm"))
# 3 states -General-population
cumhaz_data_b <- rbind(data.frame(cox_cumhaz_0_a_gp$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_b$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_b$Haz,
model = "Parametric"))
cumhaz_data_b$trans <- factor(cumhaz_data_b$trans,
levels = seq(1, 3),
labels = c("Adm -> Ther.Disch",
"Adm -> Readm",
"Ther.Disch -> Readm"))
## 4 states -Women specific
cumhaz_data_c <- rbind(data.frame(cox_cumhaz_0_c_we$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_c$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_c$Haz,
model = "Parametric"))
cumhaz_data_c$trans <- factor(cumhaz_data_c$trans,
levels = seq(1, 5),
labels = c("Adm -> Ther.Disch",
"Adm -> Disch.w/oClin.Adv.",
"Adm -> Readm",
"Ther.Disch -> Readm",
"Disch.w/oClin.Adv.-> Readm"
))
## 4 states -General population
cumhaz_data_d <- rbind(data.frame(cox_cumhaz_0_c_gp$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_d$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_d$Haz,
model = "Parametric"))
cumhaz_data_d$trans <- factor(cumhaz_data_d$trans,
levels = seq(1, 5),
labels = c("Adm -> Ther.Disch",
"Adm -> Disch.w/oClin.Adv.",
"Adm -> Readm",
"Ther.Disch -> Readm",
"Disch.w/oClin.Adv.-> Readm"
))
fig_cox_ab2<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_a,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_b,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="solid"),size=1, alpha=.65) +
facet_wrap(trans~., ncol=1) +
scale_color_manual(name ="Model",values=c("#0E53A7","#200772","#BF8830"),labels= c("Cox","NP Cox","Par"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women-specific","General Population"))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("Years") +
ylab("Cumulative hazards") +
theme_minimal()+
theme(legend.position=c(.9,.59),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
fig_cox_cd2<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_c,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_d,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="solid"),size=1, alpha=.65) +
facet_wrap(trans~., ncol=1) +
scale_color_manual(name ="Model",values=c("#0E53A7","#200772","#BF8830"),labels= c("Cox","NP Cox","Par"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women-specific","General Population"))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("Years") +
ylab("Cumulative") +
theme_minimal()+
theme(legend.position=c(.9,.6),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))fig_cox_ab2 + theme(legend.position=c(.9,.6))+ ylim(0,1)
Figure 1a. Estimate of Cumulative Hazards, Three States Model
fig_cox_cd2+ theme(legend.position=c(.9,.585))+ ylim(0,2)
Figure 1b. Estimate of Cumulative Hazards, Four States Model
cumhaz_plot1<-
ggplot()+
geom_step(data=data.frame(flexsurv_cumhaz_a$Haz) %>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=data.frame(flexsurv_cumhaz_b$Haz)%>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="solid"),size=1, alpha=.65) +
scale_color_manual(name ="Transition",values=c("#0E53A7","#200772","#BF8830"),labels= c("Admission -> Ther. Discharge","Admission -> Readmission","Therapeutic Discharge -> Readmission"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women-specific","General Population"))+
scale_x_continuous(breaks=seq(0, 11, by = 1))+
ylim(0,2)+
ylab("Cumulative hazards")+
theme_minimal()+
theme(axis.title.x = element_blank(),
legend.position=c(.65,.6),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
cumhaz_plot2<-
ggplot()+
geom_step(data=data.frame(flexsurv_cumhaz_c$Haz) %>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=data.frame(flexsurv_cumhaz_d$Haz)%>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="solid"),size=1, alpha=.65) +
scale_color_manual(name ="Transitions",values=c("gray40","#1B0773","#6E0069","#025167","#A68600"),labels= c("Admission-> Therapeutic Discharge","Admission-> Discharge Without Clinical Advice","Admission->Readmission","Therapeutic Discharge->Readmission","Discharge Without Clinical Advice->Readmission"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women-specific","General Population"))+
#values=c("#BF8F30","#1D7373","#876ED7","#1240AB")
scale_x_continuous(breaks=seq(0, 11, by = 1))+
theme_minimal()+
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position=c(.65,.28),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))+
ylim(0,2)
ggarrange(cumhaz_plot1,cumhaz_plot2)
Figure 1d. Estimate of Cumulative Hazards of the Parametric Model
if(no_mostrar==1){
jpeg("eso14.jpg", height=15, width= 10, res= 320, units = "in")
ggarrange(cumhaz_plot1,cumhaz_plot2)
dev.off()
}Transition Probabilities
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#This functionality is implemented in pmatrix.simfs,
#but the tcovs argument actually has no impact on the transition probabilities, as evidenced below.
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
pmatrix_t_a_lo<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_t_a_",t)),"lower")
}))
pmatrix_t_b_lo<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_t_b_",t)),"lower")
}))
pmatrix_4s_t_a_lo<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_4s_t_a_",t)),"lower")
}))
pmatrix_4s_t_b_lo<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_4s_t_b_",t)),"lower")
}))
pmatrix_t_a_up<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_t_a_",t)),"upper")
}))
pmatrix_t_b_up<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_t_b_",t)),"upper")
}))
pmatrix_4s_t_a_up<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_4s_t_a_",t)),"upper")
}))
pmatrix_4s_t_b_up<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
attr(get(paste0("pmatrix_4s_t_b_",t)),"upper")
}))
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
pmatrix_t_a_df_final<-
cbind.data.frame(t=rep(seq(0, 11, by = 1/12),each=dim(mat_3_states)[2]),
trans=rep(1:dim(mat_3_states)[2]),
pmatrix_t_a,pmatrix_t_a_lo,pmatrix_t_a_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]"))
pmatrix_t_b_df_final<-
cbind.data.frame(t=rep(seq(0, 11, by = 1/12),each=dim(mat_3_states)[2]),
trans=rep(1:dim(mat_3_states)[2]),
pmatrix_t_b,pmatrix_t_b_lo,pmatrix_t_b_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]"))
pmatrix_4s_t_a_df_final<-
cbind.data.frame(t=rep(seq(0, 11, by = 1/12),each=dim(mat_4_states)[2]),
trans=rep(1:dim(mat_4_states)[2]),
pmatrix_4s_t_a,pmatrix_4s_t_a_lo,pmatrix_4s_t_a_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3","Trans_4"="X4",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1","Trans_4_lo"="X4.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2","Trans_4_up"="X4.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]")) %>%
dplyr::mutate(Trans_4_cis=paste0(sprintf("%1.2f",Trans_4)," [",
sprintf("%1.2f",Trans_4_lo),"-",
sprintf("%1.2f",Trans_4_up),"]"))
pmatrix_4s_t_b_df_final<-
cbind.data.frame(t=rep(seq(0, 11, by = 1/12),each=dim(mat_4_states)[2]),
trans=rep(1:dim(mat_4_states)[2]),
pmatrix_4s_t_b,pmatrix_4s_t_b_lo,pmatrix_4s_t_b_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3","Trans_4"="X4",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1","Trans_4_lo"="X4.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2","Trans_4_up"="X4.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]")) %>%
dplyr::mutate(Trans_4_cis=paste0(sprintf("%1.2f",Trans_4)," [",
sprintf("%1.2f",Trans_4_lo),"-",
sprintf("%1.2f",Trans_4_up),"]"))#A 2021-04-23, se corrigieron los modelos no-paramétricos para que calzaran. Hubo una confusión con cuál era de cuál programa en función de la letra. Para homologar, la letra A es con programa mujeres
set.seed(1234)
pmatrix_cox_a_0<-
mstate::mssample(msf1$Haz, trans = mat_3_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:16:38 2021
## Replication 1000 finished at Tue May 25 13:16:41 2021
## Replication 1500 finished at Tue May 25 13:16:45 2021
## Replication 2000 finished at Tue May 25 13:16:48 2021
## Replication 2500 finished at Tue May 25 13:16:51 2021
## Replication 3000 finished at Tue May 25 13:16:55 2021
## Replication 3500 finished at Tue May 25 13:16:58 2021
## Replication 4000 finished at Tue May 25 13:17:01 2021
## Replication 4500 finished at Tue May 25 13:17:05 2021
## Replication 5000 finished at Tue May 25 13:17:09 2021set.seed(1234)
pmatrix_cox_b_0<-
mstate::mssample(msf0$Haz, trans = mat_3_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:17:13 2021
## Replication 1000 finished at Tue May 25 13:17:17 2021
## Replication 1500 finished at Tue May 25 13:17:21 2021
## Replication 2000 finished at Tue May 25 13:17:25 2021
## Replication 2500 finished at Tue May 25 13:17:29 2021
## Replication 3000 finished at Tue May 25 13:17:35 2021
## Replication 3500 finished at Tue May 25 13:17:38 2021
## Replication 4000 finished at Tue May 25 13:17:42 2021
## Replication 4500 finished at Tue May 25 13:17:46 2021
## Replication 5000 finished at Tue May 25 13:17:50 2021#A 2021-04-23, se corrigieron los modelos no-paramétricos para que calzaran. Hubo una confusión con cuál era de cuál programa en función de la letra. Para homologar, la letra A es con programa mujeres
set.seed(1234)
pmatrix_cox_a_4s_0<-
mstate::mssample(msf1_4s$Haz, trans = mat_4_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:17:55 2021
## Replication 1000 finished at Tue May 25 13:18:00 2021
## Replication 1500 finished at Tue May 25 13:18:04 2021
## Replication 2000 finished at Tue May 25 13:18:09 2021
## Replication 2500 finished at Tue May 25 13:18:14 2021
## Replication 3000 finished at Tue May 25 13:18:19 2021
## Replication 3500 finished at Tue May 25 13:18:25 2021
## Replication 4000 finished at Tue May 25 13:18:30 2021
## Replication 4500 finished at Tue May 25 13:18:35 2021
## Replication 5000 finished at Tue May 25 13:18:40 2021set.seed(1234)
pmatrix_cox_b_4s_0<-
mstate::mssample(msf0_4s$Haz, trans = mat_4_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:18:46 2021
## Replication 1000 finished at Tue May 25 13:18:51 2021
## Replication 1500 finished at Tue May 25 13:18:56 2021
## Replication 2000 finished at Tue May 25 13:19:01 2021
## Replication 2500 finished at Tue May 25 13:19:06 2021
## Replication 3000 finished at Tue May 25 13:19:11 2021
## Replication 3500 finished at Tue May 25 13:19:16 2021
## Replication 4000 finished at Tue May 25 13:19:22 2021
## Replication 4500 finished at Tue May 25 13:19:27 2021
## Replication 5000 finished at Tue May 25 13:19:32 2021#flexsurv_cumhaz_c$Haz
set.seed(1234)
pmatrix_flexsurv_a<-
mstate::mssample(flexsurv_cumhaz_a$Haz, trans = mat_3_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:19:35 2021
## Replication 1000 finished at Tue May 25 13:19:38 2021
## Replication 1500 finished at Tue May 25 13:19:41 2021
## Replication 2000 finished at Tue May 25 13:19:45 2021
## Replication 2500 finished at Tue May 25 13:19:48 2021
## Replication 3000 finished at Tue May 25 13:19:52 2021
## Replication 3500 finished at Tue May 25 13:19:55 2021
## Replication 4000 finished at Tue May 25 13:19:59 2021
## Replication 4500 finished at Tue May 25 13:20:03 2021
## Replication 5000 finished at Tue May 25 13:20:06 2021set.seed(1234)
pmatrix_flexsurv_b<-
mstate::mssample(flexsurv_cumhaz_b$Haz, trans = mat_3_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:20:09 2021
## Replication 1000 finished at Tue May 25 13:20:12 2021
## Replication 1500 finished at Tue May 25 13:20:15 2021
## Replication 2000 finished at Tue May 25 13:20:19 2021
## Replication 2500 finished at Tue May 25 13:20:22 2021
## Replication 3000 finished at Tue May 25 13:20:25 2021
## Replication 3500 finished at Tue May 25 13:20:29 2021
## Replication 4000 finished at Tue May 25 13:20:32 2021
## Replication 4500 finished at Tue May 25 13:20:35 2021
## Replication 5000 finished at Tue May 25 13:20:38 2021#pstate1–pstate3 are close to the first rows of the matrices returned by pmatrix.fs
set.seed(1234)
pmatrix_flexsurv_c<-
mstate::mssample(flexsurv_cumhaz_c$Haz, trans = mat_4_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:20:43 2021
## Replication 1000 finished at Tue May 25 13:20:47 2021
## Replication 1500 finished at Tue May 25 13:20:51 2021
## Replication 2000 finished at Tue May 25 13:20:56 2021
## Replication 2500 finished at Tue May 25 13:21:00 2021
## Replication 3000 finished at Tue May 25 13:21:05 2021
## Replication 3500 finished at Tue May 25 13:21:10 2021
## Replication 4000 finished at Tue May 25 13:21:14 2021
## Replication 4500 finished at Tue May 25 13:21:18 2021
## Replication 5000 finished at Tue May 25 13:21:23 2021set.seed(1234)
pmatrix_flexsurv_d<-
mstate::mssample(flexsurv_cumhaz_d$Haz, trans = mat_4_states, clock = "reset", M = n_iter/2, tvec = seq(0, 11, by = 1/12),do.trace=500)## Replication 500 finished at Tue May 25 13:21:27 2021
## Replication 1000 finished at Tue May 25 13:21:31 2021
## Replication 1500 finished at Tue May 25 13:21:36 2021
## Replication 2000 finished at Tue May 25 13:21:40 2021
## Replication 2500 finished at Tue May 25 13:21:44 2021
## Replication 3000 finished at Tue May 25 13:21:49 2021
## Replication 3500 finished at Tue May 25 13:21:53 2021
## Replication 4000 finished at Tue May 25 13:21:57 2021
## Replication 4500 finished at Tue May 25 13:22:02 2021
## Replication 5000 finished at Tue May 25 13:22:07 2021We simulated the trajectories from Semi-markov models by y generating the time to the next transition until the individual reaches an absorbing state or a specified censoring time.
invisible(c("For semi-Markov models, mstate provides the function mssample to produce both simulated trajectories and transition probability matrices from semi-Markov models, given the estimated piecewise-constant cumulative hazards (Fiocco et al. 2008), produced by msfit or msfit.flexsurvreg, though this is generally less efficient than pmatrix.simfs. In this example, 1000 samples from mssample give estimates of transition probabilities that are accurate to within around 0.02. However with pmatrix.simfs, greater precision is achieved by simulating 100 times as many trajectories in a shorter time"))
#pmatrix_t_a_df_final pmatrix_t_b_df_final pmatrix_4s_t_a_df_final pmatrix_4s_t_b_df_final
fig_ab_pmatrix<-
pmatrix_cox_a %>%
dplyr::left_join(pmatrix_cox_b, by="time", suffix=c("",".wcovs_cox.gp")) %>%
dplyr::left_join(pmatrix_cox_a_0, by="time", suffix=c("",".wo_covs.ws")) %>%
dplyr::left_join(pmatrix_cox_b_0, by="time", suffix=c("",".wo_covs.gp")) %>%
#pstate1–pstate3 are close to the first rows of the matrices returned by pmatrix.fs
dplyr::left_join(pmatrix_flexsurv_a, by="time", suffix=c("",".wcovs_par.ws")) %>%
dplyr::left_join(pmatrix_flexsurv_b, by="time", suffix=c("",".wcovs_par.gp")) %>%
dplyr::left_join(pmatrix_t_a_df_final[which(pmatrix_t_a_df_final$trans==1),c("t","Trans_1","Trans_2","Trans_3")], by=c("time"="t"))%>%
dplyr::rename("pstate1.wcovs_2par.ws"="Trans_1","pstate2.wcovs_2par.ws"="Trans_2", "pstate3.wcovs_2par.ws"="Trans_3") %>%
dplyr::left_join(pmatrix_t_b_df_final[which(pmatrix_t_b_df_final$trans==1),c("t","Trans_1","Trans_2","Trans_3")], by=c("time"="t"))%>%
dplyr::rename("pstate1.wcovs_2par.gp"="Trans_1","pstate2.wcovs_2par.gp"="Trans_2", "pstate3.wcovs_2par.gp"="Trans_3") %>%
dplyr::rename("pstate1.wcovs_cox.ws"="pstate1","pstate2.wcovs_cox.ws"="pstate2", "pstate3.wcovs_cox.ws"="pstate3") %>%
reshape2::melt(id.vars="time") %>% #janitor::tabyl(variable)
#tidyr::separate(variable, into=c("state","covs","program")) %>%
dplyr::mutate(state= dplyr::case_when(grepl("pstate1",variable)~"1)Admission",
grepl("pstate2",variable)~"2)Therapeutic Discharge",
grepl("pstate3",variable)~"3)Readmission",
T~NA_character_)) %>%
dplyr::mutate(covs= dplyr::case_when(grepl("wo_covs",variable)~"NP Cox",
grepl("wcovs_2par",variable)~"Parametric (flexsurv)",
grepl("wcovs_par",variable)~"Parametric (mssample)",
grepl("wcovs_cox",variable)~"Cox",#debe ir después, menos restrictivo
T~NA_character_)) %>%
dplyr::mutate(type_of_program= dplyr::case_when(grepl("ws",variable)~"Women-specific",
grepl("gp",variable)~"General Population",
T~NA_character_)) %>%
ggplot(aes(time, value, color= covs, linetype =type_of_program))+
geom_step(size=1, alpha=.65) +
facet_wrap(state~., ncol=1, scales="free_y") + #type_of_program
scale_color_manual(name ="Model",values=c("#A66F00","#FFD300","#6A48D7","#06266F"),labels= c("Cox","NP Cox","Par (flexsurv)", "Par (mssample)"))+
# scale_linetype_manual(name ="Type of program",values=c(1,4), labels= c("General Population","Women-specific"))+
#scale_linetype_manual(name= "Covariates",values=c(1,4), labels=c("w/covs","w/o covs"))+
scale_x_continuous(breaks=seq(0, 11,.5))+
xlab("Years") +
ylab("") +
theme_minimal()+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
theme(legend.position=c(.9,.6),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
fig_ab_pmatrix
Figure 2a. Estimate of State Occupancies at Baseline, Three-states model
#pmatrix_t_a_df_final pmatrix_t_b_df_final pmatrix_4s_t_a_df_final pmatrix_4s_t_b_df_final
fig_cd_pmatrix<-
pmatrix_cox_c %>%
dplyr::left_join(pmatrix_cox_d, by="time", suffix=c("",".wcovs_cox.gp")) %>%
dplyr::left_join(pmatrix_cox_a_4s_0, by="time", suffix=c("",".wo_covs.ws")) %>%
dplyr::left_join(pmatrix_cox_b_4s_0, by="time", suffix=c("",".wo_covs.gp")) %>%
#pstate1–pstate3 are close to the first rows of the matrices returned by pmatrix.fs
dplyr::left_join(pmatrix_flexsurv_c, by="time", suffix=c("",".wcovs_par.ws")) %>%
dplyr::left_join(pmatrix_flexsurv_d, by="time", suffix=c("",".wcovs_par.gp")) %>%
dplyr::left_join(pmatrix_4s_t_a_df_final[which(pmatrix_4s_t_a_df_final$trans==1),c("t","Trans_1","Trans_2","Trans_3","Trans_4")], by=c("time"="t"))%>%
dplyr::rename("pstate1.wcovs_2par.ws"="Trans_1","pstate2.wcovs_2par.ws"="Trans_2", "pstate3.wcovs_2par.ws"="Trans_3",
"pstate4.wcovs_2par.ws"="Trans_4") %>%
dplyr::left_join(pmatrix_4s_t_b_df_final[which(pmatrix_4s_t_b_df_final$trans==1),c("t","Trans_1","Trans_2","Trans_3","Trans_4")], by=c("time"="t"))%>%
dplyr::rename("pstate1.wcovs_2par.gp"="Trans_1","pstate2.wcovs_2par.gp"="Trans_2", "pstate3.wcovs_2par.gp"="Trans_3",
"pstate4.wcovs_2par.gp"="Trans_4") %>%
dplyr::rename("pstate1.wcovs_cox.ws"="pstate1","pstate2.wcovs_cox.ws"="pstate2", "pstate3.wcovs_cox.ws"="pstate3",
"pstate4.wcovs_cox.ws"="pstate4") %>%
reshape2::melt(id.vars="time") %>% #janitor::tabyl(variable)
#tidyr::separate(variable, into=c("state","covs","program")) %>%
dplyr::mutate(state= dplyr::case_when(grepl("pstate1",variable)~"1)Admission",
grepl("pstate2",variable)~"2)Therapeutic Discharge",
grepl("pstate3",variable)~"3)Discharge without Clinical Advice",
grepl("pstate4",variable)~"4)Readmission",
T~NA_character_)) %>%
dplyr::mutate(covs= dplyr::case_when(grepl("wo_covs",variable)~"NP Cox",
grepl("wcovs_2par",variable)~"Parametric (flexsurv)",
grepl("wcovs_par",variable)~"Parametric (mssample)",
grepl("wcovs_cox",variable)~"Cox",#debe ir después, menos restrictivo
T~NA_character_)) %>%
dplyr::mutate(type_of_program= dplyr::case_when(grepl("ws",variable)~"Women-specific",
grepl("gp",variable)~"General Population",
T~NA_character_)) %>%
ggplot(aes(time, value, color= covs, linetype =type_of_program))+
geom_step(size=1, alpha=.65) +
facet_wrap(state~., ncol=1, scales="free_y") + #type_of_program
scale_color_manual(name ="Model",values=c("#A66F00","#FFD300","#6A48D7","#06266F"),labels= c("Cox","NP Cox","Par (flexsurv)", "Par (mssample)"))+
scale_linetype_manual(name ="Type of program",values=c(1,4), labels= c("General Population","Women-specific"))+
#scale_linetype_manual(name= "Covariates",values=c(1,4), labels=c("w/covs","w/o covs"))+
scale_x_continuous(breaks=seq(0, 11,.5))+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("") +
theme_minimal()+
theme(legend.position=c(.9,.9),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
fig_cd_pmatrix
Figure 2b. Estimate of State Occupancies at Baseline, Four-states model
#_#_#_#_#_#_#_#_#_#_#__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#Transition probability matrix from a fully-parametric, semi-Markov multi-state model
#_#_#_#_#_#_#_#_#_#_#__#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_Transition probabilities of parametric models by treatment arm
pmatrix_t_a_df_final_mod<-
data.table::data.table(pmatrix_t_a_df_final)[,-c("Trans_1_cis","Trans_2_cis","Trans_3_cis")]%>%
dplyr::rename("Trans_1.lo"="Trans_1_lo", "Trans_2.lo"="Trans_2_lo",
"Trans_3.lo"="Trans_3_lo", "Trans_1.up"="Trans_1_up",
"Trans_2.up"="Trans_2_up", "Trans_3.up"="Trans_3_up") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::rename("trans_from"="trans") %>%
dplyr::mutate(variable=dplyr::case_when(!grepl("up|lo",variable)~paste0(variable,".est"),
T~as.character(variable))) %>%
tidyr::separate(variable, c("trans_to", "cis"), "\\.") %>% #View()
dplyr::mutate(cis=ifelse(cis=="est","estimate",cis)) %>%
dplyr::ungroup() %>%
#dplyr::group_by(t, trans_from, trans_to, cis) %>% summarise(n=n())%>% dplyr::filter(n>1)
tidyr::pivot_wider(names_from="cis", values_from="value") %>%
dplyr::mutate(trans_to=factor(str_sub(trans_to, start= -1)),trans_from=factor(trans_from))
pmatrix_t_b_df_final_mod<-
data.table::data.table(pmatrix_t_b_df_final)[,-c("Trans_1_cis","Trans_2_cis","Trans_3_cis")]%>%
dplyr::rename("Trans_1.lo"="Trans_1_lo", "Trans_2.lo"="Trans_2_lo",
"Trans_3.lo"="Trans_3_lo", "Trans_1.up"="Trans_1_up",
"Trans_2.up"="Trans_2_up", "Trans_3.up"="Trans_3_up") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::rename("trans_from"="trans") %>%
dplyr::mutate(variable=dplyr::case_when(!grepl("up|lo",variable)~paste0(variable,".est"),
T~as.character(variable))) %>%
tidyr::separate(variable, c("trans_to", "cis"), "\\.") %>% #View()
dplyr::mutate(cis=ifelse(cis=="est","estimate",cis)) %>%
dplyr::ungroup() %>%
#dplyr::group_by(t, trans_from, trans_to, cis) %>% summarise(n=n())%>% dplyr::filter(n>1)
tidyr::pivot_wider(names_from="cis", values_from="value") %>%
dplyr::mutate(trans_to=factor(str_sub(trans_to, start= -1)),trans_from=factor(trans_from))
rbind.data.frame(
cbind.data.frame(pmatrix_t_a_df_final_mod,program="Women-specific"),
cbind.data.frame(pmatrix_t_b_df_final_mod,program="General population")) %>%
dplyr::mutate(combi=paste0(trans_from,"-",trans_to)) %>%
dplyr::filter(!combi %in% c("1-1","2-1","2-2","3-1","3-2","3-3")) %>%
dplyr::mutate(trans_from=dplyr::case_when(trans_from==1~"1)Admission",
trans_from==2~"2)Therapeutic Discharge",
trans_from==3~"3)Readmission")) %>%
dplyr::mutate(trans_to=dplyr::case_when(trans_to==1~"to\n1)Admission",
trans_to==2~"to\n2)Therapeutic Discharge",
trans_to==3~"to\n3)Readmission")) %>%
ggplot(aes(x=t, y=estimate, fill=program, color=program))+
geom_step(size=1)+
geom_ribbon(aes(ymin=lo, ymax=up),alpha=.3)+
facet_wrap(.~trans_from+trans_to, ncol=3, scales="free_y") +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("") +
theme_minimal()+
theme(legend.position=c(.86,.8),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
Figure 3a. Estimate of State Occupancies (w/ confidence intervals), Three-states model
rbind.data.frame(
cbind.data.frame(pmatrix_t_a_df_final_mod,program="Women-specific"),
cbind.data.frame(pmatrix_t_b_df_final_mod,program="General population")) %>%
dplyr::mutate(combi=paste0(trans_from,"-",trans_to)) %>%
dplyr::filter(trans_from!=3) %>%
dplyr::mutate(trans_from=dplyr::case_when(trans_from==1~"1)From Admission",
trans_from==2~"2)From Therapeutic Discharge",
trans_from==3~"3)From Readmission")) %>%
dplyr::mutate(trans_to=dplyr::case_when(trans_to==1~"1)Admission",
trans_to==2~"2)Therapeutic Discharge",
trans_to==3~"3)Readmission")) %>%
ggplot(aes(x=t, y=estimate, fill=trans_to, color=trans_to))+
geom_area()+
facet_wrap(program~trans_from)+
sjPlot::theme_sjplot()+
theme(legend.position="bottom")+
xlab("Years")+
ylab("Probabilities")+
scale_fill_manual(name ="Transition to", values=c("#BF8F30","#1D7373","#876ED7"))+
scale_color_manual(name ="Transition to", values=c("#BF8F30","#1D7373","#876ED7"))+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_x_continuous(breaks=seq(0, 11,2))+
theme(
panel.background = element_blank(), # bg of the panel
plot.background = element_blank(), # bg of the plot
panel.grid.major = element_blank(), # get rid of major grid
panel.grid.minor = element_blank(), # get rid of minor grid
legend.background = element_blank(), # get rid of legend bg
legend.box.background = element_blank() # get rid of legend panel bg
)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
Figure 3b. Estimate of State Occupancies (Area plot), Three-states model
pmatrix_4s_t_a_df_final_mod<-
data.table::data.table(pmatrix_4s_t_a_df_final)[,-c("Trans_1_cis","Trans_2_cis","Trans_3_cis", "Trans_4_cis")]%>%
dplyr::rename("Trans_1.lo"="Trans_1_lo", "Trans_2.lo"="Trans_2_lo",
"Trans_3.lo"="Trans_3_lo", "Trans_4.lo"="Trans_4_lo", "Trans_1.up"="Trans_1_up",
"Trans_2.up"="Trans_2_up", "Trans_3.up"="Trans_3_up","Trans_4.up"="Trans_4_up") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::rename("trans_from"="trans") %>%
dplyr::mutate(variable=dplyr::case_when(!grepl("up|lo",variable)~paste0(variable,".est"),
T~as.character(variable))) %>%
tidyr::separate(variable, c("trans_to", "cis"), "\\.") %>% #View()
dplyr::mutate(cis=ifelse(cis=="est","estimate",cis)) %>%
dplyr::ungroup() %>%
#dplyr::group_by(t, trans_from, trans_to, cis) %>% summarise(n=n())%>% dplyr::filter(n>1)
tidyr::pivot_wider(names_from="cis", values_from="value") %>%
dplyr::mutate(trans_to=factor(str_sub(trans_to, start= -1)),trans_from=factor(trans_from))
pmatrix_4s_t_b_df_final_mod<-
data.table::data.table(pmatrix_4s_t_b_df_final)[,-c("Trans_1_cis","Trans_2_cis","Trans_3_cis", "Trans_4_cis")]%>%
dplyr::rename("Trans_1.lo"="Trans_1_lo", "Trans_2.lo"="Trans_2_lo",
"Trans_3.lo"="Trans_3_lo", "Trans_4.lo"="Trans_4_lo", "Trans_1.up"="Trans_1_up",
"Trans_2.up"="Trans_2_up", "Trans_3.up"="Trans_3_up","Trans_4.up"="Trans_4_up") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::rename("trans_from"="trans") %>%
dplyr::mutate(variable=dplyr::case_when(!grepl("up|lo",variable)~paste0(variable,".est"),
T~as.character(variable))) %>%
tidyr::separate(variable, c("trans_to", "cis"), "\\.") %>% #View()
dplyr::mutate(cis=ifelse(cis=="est","estimate",cis)) %>%
dplyr::ungroup() %>%
#dplyr::group_by(t, trans_from, trans_to, cis) %>% summarise(n=n())%>% dplyr::filter(n>1)
tidyr::pivot_wider(names_from="cis", values_from="value") %>%
dplyr::mutate(trans_to=factor(str_sub(trans_to, start= -1)),trans_from=factor(trans_from))
rbind.data.frame(
cbind.data.frame(pmatrix_4s_t_a_df_final_mod,program="Women-specific"),
cbind.data.frame(pmatrix_4s_t_b_df_final_mod,program="General population")) %>%
dplyr::mutate(combi=paste0(trans_from,"-",trans_to)) %>%
dplyr::filter(!combi %in% c("1-1","2-2","3-3","2-1","3-1","2-3","3-2","3-1","3-2","4-1","4-2","4-3","4-4")) %>%
dplyr::mutate(trans_from=dplyr::case_when(trans_from==1~"1)Admission",
trans_from==2~"2)Therapeutic Discharge",
trans_from==3~"3)Discharge Without Clinical Advice",
trans_from==4~"4)Readmission")) %>%
dplyr::mutate(trans_to=dplyr::case_when(trans_to==1~"to\n1)Admission",
trans_to==2~"to\n2)Therapeutic Discharge",
trans_to==3~"to\n3)Discharge Without Clinical Advice",
trans_to==4~"to\n4)Readmission")) %>%
ggplot(aes(x=t, y=estimate, fill=program, color=program))+
geom_step(size=1)+
geom_ribbon(aes(ymin=lo, ymax=up),alpha=.3)+
facet_wrap(.~trans_from+trans_to, ncol=3, scales="free_y") +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("") +
theme_minimal()+
theme(legend.position=c(.9,.2),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
Figure 3c. Estimate of State Occupancies (w/ confidence intervals), Four-states model
fig3d_plot<-
rbind.data.frame(
cbind.data.frame(pmatrix_4s_t_a_df_final_mod,program="Women-specific"),
cbind.data.frame(pmatrix_4s_t_b_df_final_mod,program="General population")) %>%
dplyr::mutate(combi=paste0(trans_from,"-",trans_to)) %>%
dplyr::filter(trans_from!=4) %>%
dplyr::mutate(trans_from=dplyr::case_when(trans_from==1~"1)From Admission",
trans_from==2~"2)From Therapeutic Discharge",
trans_from==3~"3)From Discharge Without Clinical Advice",
trans_from==4~"4)From Readmission")) %>%
dplyr::mutate(trans_to=dplyr::case_when(trans_to==1~"1)Admission",
trans_to==2~"2)Therapeutic Discharge",
trans_to==3~"3)Discharge Without Clinical Advice",
trans_to==4~"4)Readmission")) %>%
ggplot(aes(x=t, y=estimate, fill=trans_to, color=trans_to))+
geom_area()+
facet_wrap(program~trans_from)+
sjPlot::theme_sjplot()+
theme(legend.position="bottom")+
xlab("Years")+
ylab("Probabilities")+
scale_fill_manual(name ="Transition to", values=c("#BF8F30","#1D7373","#876ED7","#1240AB"))+
scale_color_manual(name ="Transition to", values=c("#BF8F30","#1D7373","#876ED7","#1240AB"))+
#scale_color_manual(name ="Model",values=c("#A66F00","#FFD300","#6A48D7","#06266F"),labels= c("Cox","NP Cox","Par (flexsurv)", "Par (mssample)"))+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_x_continuous(breaks=seq(0, 11,2))+
theme(
panel.background = element_blank(), # bg of the panel
plot.background = element_blank(), # bg of the plot
panel.grid.major = element_blank(), # get rid of major grid
panel.grid.minor = element_blank(), # get rid of minor grid
legend.background = element_blank(), # get rid of legend bg
legend.box.background = element_blank() # get rid of legend panel bg
)+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
fig3d_plot
Figure 3d. Estimate of State Occupancies (Area Plot), Four-states model
if(no_mostrar==1){
ggsave("eso13b.jpg", fig3d_plot,width = 10, height = 10, dpi = 300, units= "in")
dev.off()
}pmatrix_t_b_df_final_join<-
pmatrix_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5)) %>%
dplyr::mutate(t=factor(t, levels=c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5),
labels=c("0 days", "~30days/1 month", "~91 days/3 months", "~122days/4months", "~183 days/6 months", "~274days/9 months", "~1 year", "~3 years", "~5 years")))
pmatrix_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5)) %>%
dplyr::mutate(t=factor(t, levels=c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5),
labels=c("0 days", "~30days/1 month", "~91 days/3 months", "~122days/4months", "~183 days/6 months", "~274days/9 months", "~1 year", "~3 years", "~5 years"))) %>%
dplyr::left_join(pmatrix_t_b_df_final_join,by=c("t","trans")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
trans==3~"3)Readmission")) %>%
dplyr::select(-t) %>%
dplyr::select(`trans`,`Trans_1_cis.x`,`Trans_2_cis.x`,`Trans_3_cis.x`,`Trans_1_cis.y`,`Trans_2_cis.y`,`Trans_3_cis.y`) %>%
#2021-04-27: drop from/actual state, readmission
dplyr::filter(!grepl("Readmission",trans)) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 1. Estimated transition probabilities, Three-states model"),
col.names = c("Actual state", rep(c("Admission","Therapeutic Discharge","Readmission"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Women-specific: Transition to" = 3,"General Population: Transition to" = 3)) %>%
kableExtra::pack_rows("At 0 days", 1, 2) %>%
kableExtra::pack_rows("At ~30 days/1 month", 3, 4) %>%
kableExtra::pack_rows("At ~91 days/3 months", 5, 6) %>%
kableExtra::pack_rows("At ~122 days/4 months", 7, 8) %>%
kableExtra::pack_rows("At ~183 days/6 months", 9, 10) %>%
kableExtra::pack_rows("At ~274days/9 months", 11, 12) %>%
kableExtra::pack_rows("At ~1 year", 13, 14) %>%
kableExtra::pack_rows("At ~3 years", 15, 16) %>%
kableExtra::pack_rows("At ~5 years", 17, 18) %>%
kableExtra::scroll_box(width = "100%", height = "375px") | Actual state | Admission | Therapeutic Discharge | Readmission | Admission | Therapeutic Discharge | Readmission |
|---|---|---|---|---|---|---|
| At 0 days | ||||||
| 1)Admission | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] |
| At ~30 days/1 month | ||||||
| 1)Admission | 0.97 [0.97-0.98] | 0.03 [0.02-0.03] | 0.00 [0.00-0.00] | 0.98 [0.97-0.98] | 0.02 [0.02-0.02] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.99 [0.98-0.99] | 0.01 [0.01-0.02] | 0.00 [0.00-0.00] | 0.99 [0.98-1.00] | 0.01 [0.00-0.02] |
| At ~91 days/3 months | ||||||
| 1)Admission | 0.92 [0.91-0.93] | 0.07 [0.06-0.09] | 0.00 [0.00-0.01] | 0.94 [0.93-0.95] | 0.06 [0.05-0.07] | 0.00 [0.00-0.01] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.98 [0.96-0.99] | 0.02 [0.01-0.04] | 0.00 [0.00-0.00] | 0.98 [0.96-0.99] | 0.02 [0.01-0.04] |
| At ~122 days/4 months | ||||||
| 1)Admission | 0.90 [0.88-0.91] | 0.09 [0.08-0.11] | 0.01 [0.00-0.01] | 0.92 [0.90-0.93] | 0.08 [0.06-0.09] | 0.01 [0.00-0.01] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.97 [0.95-0.99] | 0.03 [0.01-0.05] | 0.00 [0.00-0.00] | 0.97 [0.96-0.99] | 0.03 [0.01-0.04] |
| At ~183 days/6 months | ||||||
| 1)Admission | 0.85 [0.83-0.87] | 0.13 [0.11-0.16] | 0.01 [0.01-0.02] | 0.88 [0.86-0.90] | 0.11 [0.09-0.12] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.97 [0.94-0.98] | 0.03 [0.02-0.06] | 0.00 [0.00-0.00] | 0.97 [0.94-0.98] | 0.03 [0.02-0.06] |
| At ~274days/9 months | ||||||
| 1)Admission | 0.79 [0.76-0.82] | 0.18 [0.15-0.21] | 0.03 [0.02-0.04] | 0.83 [0.81-0.86] | 0.14 [0.12-0.17] | 0.02 [0.02-0.03] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.96 [0.92-0.97] | 0.04 [0.03-0.08] | 0.00 [0.00-0.00] | 0.96 [0.93-0.97] | 0.04 [0.03-0.07] |
| At ~1 year | ||||||
| 1)Admission | 0.74 [0.71-0.78] | 0.21 [0.18-0.24] | 0.04 [0.03-0.06] | 0.79 [0.76-0.82] | 0.17 [0.15-0.20] | 0.04 [0.03-0.05] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.95 [0.91-0.97] | 0.05 [0.03-0.09] | 0.00 [0.00-0.00] | 0.95 [0.91-0.97] | 0.05 [0.03-0.09] |
| At ~3 years | ||||||
| 1)Admission | 0.53 [0.48-0.58] | 0.33 [0.28-0.38] | 0.15 [0.12-0.17] | 0.59 [0.55-0.63] | 0.27 [0.24-0.32] | 0.14 [0.11-0.16] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.90 [0.84-0.94] | 0.10 [0.06-0.16] | 0.00 [0.00-0.00] | 0.90 [0.85-0.94] | 0.10 [0.06-0.15] |
| At ~5 years | ||||||
| 1)Admission | 0.45 [0.40-0.50] | 0.34 [0.29-0.40] | 0.21 [0.17-0.24] | 0.51 [0.47-0.55] | 0.29 [0.25-0.33] | 0.20 [0.17-0.23] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.87 [0.79-0.92] | 0.13 [0.08-0.21] | 0.00 [0.00-0.00] | 0.88 [0.81-0.92] | 0.12 [0.08-0.19] |
| Note. We removed Readmission because it was an absorbing state | ||||||
#dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",estimate)," [",
# sprintf("%1.2f",lo),"-",
# sprintf("%1.2f",up),"]")) %>%
#dplyr::arrange(t,trans_from, trans_to)
#.083 = 30 days
#.25 = 91 days 3 months
#0.333 = 122 days, 4 monthspmatrix_4s_t_b_df_final_join<-
pmatrix_4s_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5)) %>%
dplyr::mutate(t=factor(t, levels=c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5),
labels=c("0 days", "~30days/1 month", "~91 days/3 months", "~122days/4months", "~183 days/6 months", "~274days/9 months", "~1 year", "~3 years", "~5 years")))
pmatrix_4s_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5)) %>%
dplyr::mutate(t=factor(t, levels=c(0.000, 0.083, .25, 0.333, 0.50, .75, 1, 3, 5),
labels=c("0 days", "~30days/1 month", "~91 days/3 months", "~122days/4months", "~183 days/6 months", "~274days/9 months", "~1 year", "~3 years", "~5 years"))) %>%
dplyr::left_join(pmatrix_4s_t_b_df_final_join,by=c("t","trans")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
trans==3~"3)From Discharge Without\nClinical Advice",
trans==4~"4)Readmission")) %>%
dplyr::select(-t) %>%
dplyr::select(trans,Trans_1_cis.x,Trans_2_cis.x,Trans_3_cis.x,Trans_4_cis.x,Trans_1_cis.y,Trans_2_cis.y,Trans_3_cis.y,Trans_4_cis.y) %>%
#2021-04-27: drop from/actual state, readmission
dplyr::filter(!grepl("Readmission",trans)) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 2. Estimated transition probabilities, Four-states model"),
col.names = c("Actual state", rep(c("Admission","Therapeutic Discharge","Discharge Without\nClinical Advice","Readmission"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Women-specific: Transition to" = 4,"General Population: Transition to" = 4)) %>%
kableExtra::pack_rows("At 0 days", 1, 3) %>%
kableExtra::pack_rows("At ~30 days/1 month", 4, 6) %>%
kableExtra::pack_rows("At ~91 days/3 months", 7, 9) %>%
kableExtra::pack_rows("At ~122 days/4 months", 10, 12) %>%
kableExtra::pack_rows("At ~183 days/6 months", 13, 15) %>%
kableExtra::pack_rows("At ~274days/9 months", 16, 18) %>%
kableExtra::pack_rows("At ~1 year", 19, 21) %>%
kableExtra::pack_rows("At ~3 years", 22, 24) %>%
kableExtra::pack_rows("At ~5 years", 25, 27) %>%
kableExtra::scroll_box(width = "100%", height = "375px") | Actual state | Admission | Therapeutic Discharge | Discharge Without Clinical Advice | Readmission | Admission | Therapeutic Discharge | Discharge Without Clinical Advice | Readmission |
|---|---|---|---|---|---|---|---|---|
| At 0 days | ||||||||
| 1)Admission | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 1.00 [1.00-1.00] | 0.00 [0.00-0.00] |
| At ~30 days/1 month | ||||||||
| 1)Admission | 0.91 [0.90-0.92] | 0.03 [0.03-0.04] | 0.06 [0.05-0.07] | 0.00 [0.00-0.00] | 0.91 [0.90-0.92] | 0.02 [0.02-0.03] | 0.07 [0.06-0.07] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.99 [0.98-0.99] | 0.00 [0.00-0.00] | 0.01 [0.01-0.02] | 0.00 [0.00-0.00] | 0.99 [0.98-1.00] | 0.00 [0.00-0.00] | 0.01 [0.00-0.02] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.98 [0.97-0.99] | 0.02 [0.01-0.03] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.99 [0.98-0.99] | 0.01 [0.01-0.02] |
| At ~91 days/3 months | ||||||||
| 1)Admission | 0.75 [0.73-0.78] | 0.08 [0.07-0.10] | 0.16 [0.14-0.18] | 0.01 [0.00-0.01] | 0.76 [0.74-0.79] | 0.06 [0.05-0.08] | 0.17 [0.15-0.19] | 0.00 [0.00-0.01] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.98 [0.96-0.99] | 0.00 [0.00-0.00] | 0.02 [0.01-0.04] | 0.00 [0.00-0.00] | 0.98 [0.96-0.99] | 0.00 [0.00-0.00] | 0.02 [0.01-0.04] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.96 [0.94-0.97] | 0.04 [0.03-0.06] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.97 [0.95-0.98] | 0.03 [0.02-0.05] |
| At ~122 days/4 months | ||||||||
| 1)Admission | 0.69 [0.67-0.72] | 0.10 [0.09-0.12] | 0.19 [0.17-0.21] | 0.01 [0.01-0.01] | 0.70 [0.68-0.73] | 0.08 [0.07-0.10] | 0.21 [0.18-0.23] | 0.01 [0.01-0.01] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.97 [0.95-0.99] | 0.00 [0.00-0.00] | 0.03 [0.01-0.05] | 0.00 [0.00-0.00] | 0.97 [0.96-0.99] | 0.00 [0.00-0.00] | 0.03 [0.01-0.04] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.95 [0.93-0.97] | 0.05 [0.03-0.07] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.96 [0.94-0.97] | 0.04 [0.03-0.06] |
| At ~183 days/6 months | ||||||||
| 1)Admission | 0.59 [0.56-0.62] | 0.14 [0.12-0.17] | 0.25 [0.22-0.28] | 0.02 [0.01-0.02] | 0.61 [0.58-0.64] | 0.11 [0.10-0.13] | 0.26 [0.24-0.29] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.96 [0.94-0.98] | 0.00 [0.00-0.00] | 0.04 [0.02-0.06] | 0.00 [0.00-0.00] | 0.97 [0.94-0.98] | 0.00 [0.00-0.00] | 0.03 [0.02-0.06] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.94 [0.91-0.95] | 0.06 [0.05-0.09] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.94 [0.93-0.96] | 0.06 [0.04-0.07] |
| At ~274days/9 months | ||||||||
| 1)Admission | 0.48 [0.44-0.51] | 0.19 [0.16-0.22] | 0.30 [0.27-0.34] | 0.03 [0.02-0.04] | 0.50 [0.47-0.53] | 0.15 [0.13-0.17] | 0.32 [0.30-0.36] | 0.03 [0.02-0.03] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.95 [0.92-0.97] | 0.00 [0.00-0.00] | 0.05 [0.03-0.08] | 0.00 [0.00-0.00] | 0.96 [0.93-0.98] | 0.00 [0.00-0.00] | 0.04 [0.02-0.07] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.92 [0.89-0.94] | 0.08 [0.06-0.11] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.93 [0.90-0.95] | 0.07 [0.05-0.10] |
| At ~1 year | ||||||||
| 1)Admission | 0.40 [0.36-0.43] | 0.22 [0.19-0.26] | 0.34 [0.31-0.37] | 0.04 [0.03-0.05] | 0.42 [0.39-0.46] | 0.18 [0.15-0.21] | 0.36 [0.33-0.40] | 0.04 [0.03-0.05] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.95 [0.91-0.97] | 0.00 [0.00-0.00] | 0.05 [0.03-0.09] | 0.00 [0.00-0.00] | 0.95 [0.91-0.97] | 0.00 [0.00-0.00] | 0.05 [0.03-0.09] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.90 [0.87-0.93] | 0.10 [0.07-0.13] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.92 [0.89-0.94] | 0.08 [0.06-0.11] |
| At ~3 years | ||||||||
| 1)Admission | 0.16 [0.13-0.19] | 0.34 [0.29-0.38] | 0.38 [0.35-0.42] | 0.12 [0.10-0.15] | 0.19 [0.16-0.22] | 0.28 [0.24-0.32] | 0.42 [0.38-0.46] | 0.11 [0.09-0.14] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.90 [0.84-0.94] | 0.00 [0.00-0.00] | 0.10 [0.06-0.16] | 0.00 [0.00-0.00] | 0.90 [0.85-0.94] | 0.00 [0.00-0.00] | 0.10 [0.06-0.15] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.83 [0.78-0.86] | 0.17 [0.14-0.22] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.85 [0.80-0.88] | 0.15 [0.12-0.20] |
| At ~5 years | ||||||||
| 1)Admission | 0.10 [0.07-0.12] | 0.36 [0.31-0.41] | 0.37 [0.33-0.41] | 0.17 [0.14-0.21] | 0.13 [0.10-0.15] | 0.31 [0.26-0.35] | 0.41 [0.36-0.45] | 0.16 [0.14-0.20] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.87 [0.80-0.92] | 0.00 [0.00-0.00] | 0.13 [0.08-0.20] | 0.00 [0.00-0.00] | 0.88 [0.81-0.92] | 0.00 [0.00-0.00] | 0.12 [0.08-0.19] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.78 [0.73-0.82] | 0.22 [0.18-0.27] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.80 [0.75-0.84] | 0.20 [0.16-0.25] |
| Note. We removed Readmission because it was an absorbing state | ||||||||
Length of Stay
Partitioned survival analyses
We first looked over the restricted mean survival time, which may be the equivalent measure of the length of stay for partitioned survival analyses of each transition (as stated by Crowther in this link). We calculated its 95% confidence intervals by 10,000 resamples.
rmst_a_3s_1<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=1,summary(sel_param_fits_b[[1]],newdata=newdat3a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_3s_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=2,summary(sel_param_fits_b[[2]],newdata=newdat3a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_3s_3<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=3,summary(sel_param_fits_b[[3]],newdata=newdat3a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_1<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=1,summary(sel_param_fits_b[[1]],newdata=newdat3b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=2,summary(sel_param_fits_b[[2]],newdata=newdat3b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_3<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=3,summary(sel_param_fits_b[[3]],newdata=newdat3b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_4s_1<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=1,summary(sel_param_fits_4s_b[[1]],newdata=newdat5a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=2,summary(sel_param_fits_4s_b[[2]],newdata=newdat5a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_3<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=3,summary(sel_param_fits_4s_b[[3]],newdata=newdat5a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_4<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=4,summary(sel_param_fits_4s_b[[4]],newdata=newdat5a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_5<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=5,summary(sel_param_fits_4s_b[[5]],newdata=newdat5a[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_1<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=1,summary(sel_param_fits_4s_b[[1]],newdata=newdat5b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=2,summary(sel_param_fits_4s_b[[2]],newdata=newdat5b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_3<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=3,summary(sel_param_fits_4s_b[[3]],newdata=newdat5b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_4<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=4,summary(sel_param_fits_4s_b[[4]],newdata=newdat5b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_5<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=5,summary(sel_param_fits_4s_b[[5]],newdata=newdat5b[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_df_psa_3s<-
rbind.data.frame(rmst_a_3s_1, rmst_a_3s_2, rmst_a_3s_3, rmst_b_3s_1, rmst_b_3s_2, rmst_b_3s_3)
rmst_df_psa_4s<-
rbind.data.frame(rmst_a_4s_1, rmst_a_4s_2, rmst_a_4s_3, rmst_a_4s_4, rmst_a_4s_5, rmst_b_4s_1, rmst_b_4s_2, rmst_b_4s_3, rmst_b_4s_4, rmst_b_4s_5)transition_label<- c(`1`="Transition 1: Admission to Ther.Discharge",`2`="Transition 2: Admission to Readmission",`3`="Transition 3: Ther.Discharge to Readmission")
ggplot(rmst_df_psa_3s, aes(x=t, y=est, color=pr, fill=pr))+
geom_step(size=1)+
geom_ribbon(aes(ymin=lcl, ymax=ucl),alpha=.3)+
facet_wrap(.~tr, labeller = labeller(tr = transition_label), ncol=1, scales = "free_y")+
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position=c(.9,.1),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
Figure 4a. Restricted Mean Survival Times (w/ confidence intervals), Three-states model
ggplot(rmst_df_psa_4s, aes(x=t, y=est, color=pr, fill=pr))+
geom_step(size=1)+
geom_ribbon(aes(ymin=lcl, ymax=ucl),alpha=.3)+
facet_wrap(.~tr, labeller = labeller(tr = transition_label_4s), ncol=1, scales = "free_y")+
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position=c(.9,.1),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
Figure 4b. Restricted Mean Survival Times (w/ confidence intervals), Four-states model
#In order to estimate performances of average LOS (ALOS) estimation methods according to the accumulating data that become available over time, ALOS with their 95% confidence intervals (CI)
#rmst_df_psa_4s
rmst_df_psa_3s %>%
dplyr::mutate(pr=dplyr::case_when(pr=="ws"~"Women-specific",
T~"General population")) %>%
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",lcl),"-",
sprintf("%1.2f",ucl),"]")) %>%
dplyr::select(pr, t, tr, st, est_ci) %>%
tidyr::pivot_wider(names_from="pr", values_from=c("est_ci")) %>%
dplyr::select(-st) %>%
dplyr::mutate_at(1,~sprintf("%1.4f",.)) %>%
dplyr::filter(t %in% c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000")) %>%
dplyr::mutate(t=factor(t,levels=c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000"), labels=c("0 days","30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::select(-tr) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 3a. Restricted Mean Survival Time by Partitioned Survival Analyses, Three-states model"),
col.names = c("Time",rep(c("Estimate [95% CI]"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Transition 1: Admission to Ther.Discharge", 1, 8) %>%
kableExtra::pack_rows("Transition 2: Admission to Readmission", 9, 16) %>%
kableExtra::pack_rows("Transition 3: Ther.Discharge to Readmission", 17, 24) %>%
kableExtra::scroll_box(width = "100%", height = "375px")| Time | Estimate [95% CI] | Estimate [95% CI] |
|---|---|---|
| Transition 1: Admission to Ther.Discharge | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.24 [0.24-0.24] | 0.24 [0.24-0.24] |
| 4 months | 0.32 [0.31-0.32] | 0.32 [0.32-0.32] |
| 6 months | 0.46 [0.46-0.47] | 0.47 [0.47-0.48] |
| 1 year | 0.87 [0.85-0.89] | 0.90 [0.88-0.91] |
| 3 years | 2.24 [2.12-2.35] | 2.38 [2.28-2.47] |
| 5 years | 3.45 [3.22-3.66] | 3.72 [3.53-3.90] |
| Transition 2: Admission to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.25-0.25] | 0.25 [0.25-0.25] |
| 4 months | 0.33 [0.33-0.33] | 0.33 [0.33-0.33] |
| 6 months | 0.50 [0.50-0.50] | 0.50 [0.50-0.50] |
| 1 year | 0.98 [0.98-0.99] | 0.99 [0.98-0.99] |
| 3 years | 2.78 [2.71-2.83] | 2.80 [2.74-2.85] |
| 5 years | 4.38 [4.24-4.50] | 4.43 [4.30-4.55] |
| Transition 3: Ther.Discharge to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.24-0.25] | 0.25 [0.24-0.25] |
| 4 months | 0.33 [0.32-0.33] | 0.33 [0.32-0.33] |
| 6 months | 0.49 [0.48-0.49] | 0.49 [0.48-0.49] |
| 1 year | 0.97 [0.94-0.98] | 0.97 [0.95-0.98] |
| 3 years | 2.81 [2.68-2.89] | 2.82 [2.70-2.89] |
| 5 years | 4.57 [4.31-4.74] | 4.60 [4.36-4.75] |
#kableExtra::add_footnote("Note. NA= Null values", notation="none") #In order to estimate performances of average LOS (ALOS) estimation methods according to the accumulating data that become available over time, ALOS with their 95% confidence intervals (CI)
#rmst_df_psa_4s
rmst_df_psa_4s %>%
dplyr::mutate(pr=dplyr::case_when(pr=="ws"~"Women-specific",
T~"General population")) %>%
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",lcl),"-",
sprintf("%1.2f",ucl),"]")) %>%
dplyr::select(pr, t, tr, st, est_ci) %>%
tidyr::pivot_wider(names_from="pr", values_from=c("est_ci")) %>%
dplyr::select(-st) %>%
dplyr::mutate_at(1,~sprintf("%1.4f",.)) %>%
dplyr::filter(t %in% c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000")) %>%
dplyr::mutate(t=factor(t,levels=c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000"), labels=c("0 days", "30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::select(-tr) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 3b. Restricted Mean Survival Time by Partitioned Survival Analyses, Four-states model"),
col.names = c("Time",rep(c("Estimate [95% CI]"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Transition 1: Admission to Ther.Discharge", 1, 8) %>%
kableExtra::pack_rows("Transition 2: Admission to Discharge w/o Clinical Advice", 9, 16) %>%
kableExtra::pack_rows("Transition 3: Admission to Readmission", 17, 24) %>%
kableExtra::pack_rows("Transition 4: Ther.Discharge to Readmission", 25, 32) %>%
kableExtra::pack_rows("Transition 5: Discharge w/o Clinical Advice to Readmission", 33, 40) %>%
kableExtra::scroll_box(width = "100%", height = "375px")| Time | Estimate [95% CI] | Estimate [95% CI] |
|---|---|---|
| Transition 1: Admission to Ther.Discharge | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.24 [0.24-0.24] | 0.24 [0.24-0.24] |
| 4 months | 0.31 [0.31-0.32] | 0.32 [0.31-0.32] |
| 6 months | 0.45 [0.45-0.46] | 0.46 [0.46-0.47] |
| 1 year | 0.83 [0.81-0.86] | 0.87 [0.85-0.89] |
| 3 years | 1.89 [1.74-2.03] | 2.09 [1.96-2.20] |
| 5 years | 2.57 [2.29-2.83] | 2.94 [2.69-3.19] |
| Transition 2: Admission to Discharge w/o Clinical Advice | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.23 [0.22-0.23] | 0.23 [0.22-0.23] |
| 4 months | 0.30 [0.29-0.30] | 0.29 [0.29-0.30] |
| 6 months | 0.42 [0.41-0.43] | 0.42 [0.41-0.43] |
| 1 year | 0.74 [0.71-0.77] | 0.73 [0.71-0.76] |
| 3 years | 1.69 [1.58-1.80] | 1.65 [1.54-1.77] |
| 5 years | 2.51 [2.30-2.71] | 2.44 [2.23-2.65] |
| Transition 3: Admission to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.25-0.25] | 0.25 [0.25-0.25] |
| 4 months | 0.33 [0.33-0.33] | 0.33 [0.33-0.33] |
| 6 months | 0.50 [0.50-0.50] | 0.50 [0.50-0.50] |
| 1 year | 1.00 [0.99-1.00] | 1.00 [0.99-1.00] |
| 3 years | 2.91 [2.85-2.95] | 2.92 [2.86-2.96] |
| 5 years | 4.68 [4.49-4.81] | 4.71 [4.54-4.83] |
| Transition 4: Ther.Discharge to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.24-0.25] | 0.25 [0.24-0.25] |
| 4 months | 0.33 [0.32-0.33] | 0.33 [0.32-0.33] |
| 6 months | 0.49 [0.48-0.49] | 0.49 [0.48-0.49] |
| 1 year | 0.97 [0.94-0.98] | 0.97 [0.95-0.98] |
| 3 years | 2.81 [2.68-2.89] | 2.82 [2.71-2.89] |
| 5 years | 4.57 [4.32-4.74] | 4.60 [4.35-4.75] |
| Transition 5: Discharge w/o Clinical Advice to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.24 [0.24-0.25] | 0.25 [0.24-0.25] |
| 4 months | 0.32 [0.32-0.33] | 0.33 [0.32-0.33] |
| 6 months | 0.48 [0.47-0.49] | 0.48 [0.48-0.49] |
| 1 year | 0.94 [0.92-0.96] | 0.95 [0.93-0.96] |
| 3 years | 2.66 [2.56-2.75] | 2.70 [2.61-2.77] |
| 5 years | 4.27 [4.07-4.44] | 4.34 [4.16-4.49] |
#kableExtra::add_footnote("Note. NA= Null values", notation="none") Multistate
tolos_t_a_str2_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_3_states)[2],get(paste0("tolos_t_a_str2_",t)))
}))
tolos_t_b_str2_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_3_states)[2],get(paste0("tolos_t_b_str2_",t)))
}))
tolos_4s_t_a_str2_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_a_str2_",t)))
}))
tolos_4s_t_b_str2_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_b_str2_",t)))
}))
tolos_4s_t_a_str3_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_a_str3_",t)))
}))
tolos_4s_t_b_str3_df<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_b_str3_",t)))
}))
tolos_t_a_df2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_3_states)[2],get(paste0("tolos_t_a_",t)))
}))
tolos_t_b_df2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_3_states)[2],get(paste0("tolos_t_b_",t)))
}))
tolos_4s_t_a_df2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_a_",t)))
}))
tolos_4s_t_b_df2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_4_states)[2],get(paste0("tolos_4s_t_b_",t)))
}))#tolos_t_a_df2 tolos_t_b_df2 tolos_4s_t_a_df2 %>% dplyr::filter(t==10) tolos_4s_t_b_df2 tolos_4s_t_a_10
los_3s<-
rbind.data.frame(
cbind.data.frame(tolos_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos_t_a_str2_df,program="Women-specific", start=2)%>% dplyr::filter(state==2),
cbind.data.frame(tolos_t_b_str2_df,program="General Population", start=2)%>% dplyr::filter(state==2)
) %>%
dplyr::filter(state<3) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Readmission"))%>%
ggplot(aes(x=t, y=est, color=program, fill=program))+
geom_step(size=1)+
geom_ribbon(aes(ymin=L, ymax=U),alpha=.3)+
facet_wrap(.~start+state, ncol=3, scales="free_y") +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position=c(.87,.1),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
los_3s
Figure 5a. Estimated Length of Stay (w/ confidence intervals), Three-states model
if(no_mostrar==1){
jpeg("eso17.jpg", height=15, width= 10, res= 96, units = "in")
los_3s
dev.off()
}rbind.data.frame(
cbind.data.frame(tolos_4s_t_a_df2,program="Women-specific",start=1),
cbind.data.frame(tolos_4s_t_b_df2,program="General Population",start=1),
cbind.data.frame(tolos_4s_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos_4s_t_b_str2_df,program="General Population", start=2),
cbind.data.frame(tolos_4s_t_a_str3_df,program="Women-specific", start=3),
cbind.data.frame(tolos_4s_t_b_str3_df,program="General Population", start=3))%>%
dplyr::filter(state<4) %>%
dplyr::mutate(comb=dplyr::case_when(
start==2 & state==1~"a",
start==2 & state==3~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(is.na(comb)) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Discharge Without Clinical Advice",
T~"to\n4)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Discharge Without Clinical Advice",
T~"Starting from Readmission"))%>%
dplyr::arrange(start, state) %>%
ggplot(aes(x=t, y=est, color=program, fill=program))+
geom_step(size=1)+
geom_ribbon(aes(ymin=L, ymax=U),alpha=.3)+
facet_wrap(.~start+state, ncol=3, scales="free_y") +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("General Population","Women-specific"))+
scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position=c(.9,.1),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
Figure 5b. Estimated Length of Stay (w/ confidence intervals), Four-states model
if(no_mostrar==1){
jpeg("C:/Users/andre/Desktop/SUD_CL/eso17.jpg", height=15, width= 10, res= 96, units = "in")
plot_comp_mssurvfit2
dev.off()
}rbind.data.frame(
cbind.data.frame(tolos_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos_t_b_str2_df,program="General Population", start=2))%>%
dplyr::mutate(comb=dplyr::case_when(start==2 & state==1~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(state<max(mat_3_states,na.rm=T)) %>%
dplyr::mutate(start=dplyr::case_when(start==1~"1)Admission", start==2~"2)Therapeutic\nDischarge",start==3~"3)Readmission",T~NA_character_))%>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic\nDischarge",
state==3~"3)Readmission"))%>%
dplyr::filter(is.na(comb)) %>% # janitor::tabyl(start, state)
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",L),"-",
sprintf("%1.2f",U),"]")) %>%
dplyr::select(program, t, start, state, est_ci) %>%
tidyr::pivot_wider(names_from="program", values_from=c("est_ci")) %>%
dplyr::mutate(comb=paste0(start, "_", state)) %>% #janitor::tabyl(start, state)
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000"),
labels=c("0 days","30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::arrange(start, state, t) %>%
dplyr::select(-start,-comb) %>%
dplyr::select(state, t, `Women-specific`, `General Population`) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 4a. Length of Stay, Three-states model"),
col.names = c("State","Time","Women-specific","General population"),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::pack_rows("Starting From Admission", 1, 16) %>%
kableExtra::pack_rows("Starting From Therapeutic Discharge", 17, 24) %>%
kableExtra::add_footnote("Note. Readmission state will not be considered because is an absorbing state", notation="none") %>%
kableExtra::scroll_box(width = "100%", height = "375px")| State | Time | Women-specific | General population |
|---|---|---|---|
| Starting From Admission | |||
| 1)Admission | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 1)Admission | 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 1)Admission | 90 days | 0.24 [0.24-0.24] | 0.24 [0.24-0.24] |
| 1)Admission | 4 months | 0.32 [0.31-0.32] | 0.32 [0.32-0.32] |
| 1)Admission | 6 months | 0.46 [0.45-0.47] | 0.47 [0.46-0.47] |
| 1)Admission | 1 year | 0.86 [0.84-0.88] | 0.89 [0.87-0.90] |
| 1)Admission | 3 years | 2.09 [1.97-2.19] | 2.23 [2.14-2.32] |
| 1)Admission | 5 years | 3.06 [2.84-3.26] | 3.32 [3.13-3.51] |
| 2)Therapeutic Discharge | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 30 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 90 days | 0.01 [0.01-0.01] | 0.01 [0.01-0.01] |
| 2)Therapeutic Discharge | 4 months | 0.02 [0.01-0.02] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 6 months | 0.04 [0.03-0.04] | 0.03 [0.02-0.03] |
| 2)Therapeutic Discharge | 1 year | 0.12 [0.10-0.14] | 0.10 [0.08-0.12] |
| 2)Therapeutic Discharge | 3 years | 0.69 [0.60-0.81] | 0.57 [0.49-0.67] |
| 2)Therapeutic Discharge | 5 years | 1.37 [1.17-1.62] | 1.14 [0.96-1.33] |
| Starting From Therapeutic Discharge | |||
| 2)Therapeutic Discharge | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 2)Therapeutic Discharge | 90 days | 0.25 [0.24-0.25] | 0.25 [0.24-0.25] |
| 2)Therapeutic Discharge | 4 months | 0.33 [0.32-0.33] | 0.33 [0.32-0.33] |
| 2)Therapeutic Discharge | 6 months | 0.49 [0.48-0.49] | 0.49 [0.48-0.49] |
| 2)Therapeutic Discharge | 1 year | 0.97 [0.94-0.98] | 0.97 [0.95-0.98] |
| 2)Therapeutic Discharge | 3 years | 2.81 [2.69-2.89] | 2.82 [2.70-2.89] |
| 2)Therapeutic Discharge | 5 years | 4.57 [4.32-4.75] | 4.60 [4.35-4.75] |
| Note. Readmission state will not be considered because is an absorbing state | |||
rbind.data.frame(
cbind.data.frame(tolos_4s_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos_4s_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos_4s_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos_4s_t_b_str2_df,program="General Population", start=2),
cbind.data.frame(tolos_4s_t_a_str3_df,program="Women-specific", start=3),
cbind.data.frame(tolos_4s_t_b_str3_df,program="General Population", start=3))%>%
dplyr::mutate(comb=dplyr::case_when(start==2 & state==1~"a",
start==2 & state==3~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(state<max(dim(mat_4_states)[2],na.rm=T)) %>%
dplyr::mutate(start=dplyr::case_when(start==1~"1)Admission", start==2~"2)Therapeutic Discharge",start==3~"3)Discharge Without Clinical Advice",start==4~"4)Readmission",T~NA_character_))%>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Discharge Without Clinical Advice",
T~"4)Readmission"))%>%
dplyr::filter(is.na(comb)) %>% #janitor::tabyl(start, state)
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",L),"-",
sprintf("%1.2f",U),"]")) %>%
dplyr::select(program,t,start,state,est_ci) %>%
tidyr::pivot_wider(names_from="program",values_from=c("est_ci")) %>%
dplyr::mutate(comb=paste0(start,"_",state)) %>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.0000","0.0833","0.2500","0.3333","0.5000","1.0000","3.0000","5.0000")) %>%
dplyr::mutate(t=factor(t,levels=c("0.0000","0.0833","0.2500","0.3333","0.5000","1.0000","3.0000","5.0000"),labels=c("0 days","30 days","90 days","4 months","6 months","1 year","3 years","5 years"))) %>%
dplyr::arrange(start, state, t) %>%
dplyr::select(-start,-comb) %>%
dplyr::select(state, t, `Women-specific`, `General Population`) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 4b. Length of Stay, Four-states model"),
col.names = c("State","Time","Women-specific","General population"),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
#kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Starting From Admission", 1, 24) %>%
kableExtra::pack_rows("Starting From Therapeutic Discharge", 25, 32) %>%
kableExtra::pack_rows("Starting From Discharge Without Clinical Advice", 33, 40) %>%
kableExtra::add_footnote("Note. Readmission state will not be considered because is an absorbing state", notation="none") %>%
kableExtra::scroll_box(width = "100%", height = "375px")| State | Time | Women-specific | General population |
|---|---|---|---|
| Starting From Admission | |||
| 1)Admission | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 1)Admission | 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 1)Admission | 90 days | 0.22 [0.21-0.22] | 0.22 [0.22-0.22] |
| 1)Admission | 4 months | 0.28 [0.27-0.28] | 0.28 [0.27-0.28] |
| 1)Admission | 6 months | 0.38 [0.37-0.39] | 0.39 [0.38-0.40] |
| 1)Admission | 1 year | 0.63 [0.60-0.66] | 0.64 [0.61-0.67] |
| 1)Admission | 3 years | 1.12 [1.02-1.21] | 1.19 [1.10-1.29] |
| 1)Admission | 5 years | 1.37 [1.21-1.52] | 1.50 [1.33-1.65] |
| 2)Therapeutic Discharge | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 30 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 90 days | 0.01 [0.01-0.01] | 0.01 [0.01-0.01] |
| 2)Therapeutic Discharge | 4 months | 0.02 [0.02-0.02] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 6 months | 0.04 [0.03-0.05] | 0.03 [0.03-0.04] |
| 2)Therapeutic Discharge | 1 year | 0.13 [0.11-0.16] | 0.10 [0.09-0.12] |
| 2)Therapeutic Discharge | 3 years | 0.72 [0.62-0.83] | 0.58 [0.49-0.68] |
| 2)Therapeutic Discharge | 5 years | 1.42 [1.23-1.64] | 1.17 [1.01-1.34] |
| 3)Discharge Without Clinical Advice | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 3)Discharge Without Clinical Advice | 30 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 3)Discharge Without Clinical Advice | 90 days | 0.02 [0.02-0.02] | 0.02 [0.02-0.03] |
| 3)Discharge Without Clinical Advice | 4 months | 0.04 [0.03-0.04] | 0.04 [0.03-0.04] |
| 3)Discharge Without Clinical Advice | 6 months | 0.07 [0.07-0.08] | 0.08 [0.07-0.09] |
| 3)Discharge Without Clinical Advice | 1 year | 0.22 [0.20-0.25] | 0.24 [0.21-0.26] |
| 3)Discharge Without Clinical Advice | 3 years | 0.98 [0.88-1.08] | 1.06 [0.96-1.15] |
| 3)Discharge Without Clinical Advice | 5 years | 1.74 [1.55-1.91] | 1.89 [1.69-2.08] |
| Starting From Therapeutic Discharge | |||
| 2)Therapeutic Discharge | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 2)Therapeutic Discharge | 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 2)Therapeutic Discharge | 90 days | 0.25 [0.24-0.25] | 0.25 [0.24-0.25] |
| 2)Therapeutic Discharge | 4 months | 0.33 [0.32-0.33] | 0.33 [0.32-0.33] |
| 2)Therapeutic Discharge | 6 months | 0.49 [0.48-0.49] | 0.49 [0.48-0.49] |
| 2)Therapeutic Discharge | 1 year | 0.97 [0.94-0.98] | 0.97 [0.95-0.98] |
| 2)Therapeutic Discharge | 3 years | 2.81 [2.68-2.89] | 2.82 [2.70-2.89] |
| 2)Therapeutic Discharge | 5 years | 4.57 [4.31-4.75] | 4.60 [4.35-4.75] |
| Starting From Discharge Without Clinical Advice | |||
| 3)Discharge Without Clinical Advice | 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 3)Discharge Without Clinical Advice | 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 3)Discharge Without Clinical Advice | 90 days | 0.24 [0.24-0.25] | 0.25 [0.24-0.25] |
| 3)Discharge Without Clinical Advice | 4 months | 0.32 [0.32-0.33] | 0.33 [0.32-0.33] |
| 3)Discharge Without Clinical Advice | 6 months | 0.48 [0.47-0.49] | 0.48 [0.48-0.49] |
| 3)Discharge Without Clinical Advice | 1 year | 0.94 [0.92-0.96] | 0.95 [0.93-0.96] |
| 3)Discharge Without Clinical Advice | 3 years | 2.66 [2.57-2.75] | 2.70 [2.61-2.78] |
| 3)Discharge Without Clinical Advice | 5 years | 4.27 [4.08-4.45] | 4.34 [4.15-4.50] |
| Note. Readmission state will not be considered because is an absorbing state | |||
Simulate Final State
We estimated the median of stay, percentiles 25, 75 and confidence intervals based on 500 resamples, by simulating a large sample of individuals from the model depending on the type of program they were admitted at baseline.
time_before_sim_medians<-Sys.time()
#Estimates the probability of each final outcome ("absorbing" state), and the mean and quantiles of the time to that outcome for people who experience it, by simulating a large sample of individuals from the model. This can be used for both Markov and semi-Markov models.C
#n_iter=12
fmsm_sel_param_fits_4s_b<-
fmsm(fits_c_gomp2[[1]],fits_c_gomp2[[2]],fits_c_logn2[[3]],fits_c_ggam2[[4]],fits_c_logn2[[5]], trans=mat_4_states)
fmsm_sel_param_fits_b<-
fmsm(fits_c_gomp[[1]],fits_c_genf[[2]],fits_c_ggam[[3]], trans=mat_3_states)
set.seed(1234)
simfinal_fmsm_3_st_30d<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1/12, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_90d<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1/4.0583, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_4m<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= (1/12)*4, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_6m<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= (1/2), #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_1y<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_3y<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 3, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_3_st_5y<-
simfinal_fmsm(fmsm_sel_param_fits_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat3a[1,],newdat3b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 5, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=10000000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_30d<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1/12, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_90d<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1/4.0583, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_4m<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= (1/12)*4, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_6m<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1/2, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_1y<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 1, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_3y<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 3, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
set.seed(1234)
simfinal_fmsm_4_st_5y<-
simfinal_fmsm(fmsm_sel_param_fits_4s_b, #Object returned by fmsm, representing a multi-state model formed from transition-specific time-to-event models fitted by flexsurvreg.
newdata = rbind.data.frame(newdat5a[1,],newdat5b[1,])[,-"strata"], #Data frame of covariate values, with one column per covariate, and one row per alternative value.
t= 5, #itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
probs = c(0.025,0.25, 0.5, 0.75,0.975),
#M=100000, #necesita grandes números
B=n_iter/20, #PondrÃa 1,000
cores=8
)
time_after_sim_medians<-Sys.time()
paste0("Time in process: ");time_after_sim_medians-time_before_sim_medians## [1] "Time in process: "## Time difference of 1.541556 hours# probability of each final outcome
# mean and quantiles of the time to that outcome for people who experience itsimfinal_fmsm_3_st_df<-
rbind.data.frame(
simfinal_fmsm_3_st_30d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.08),
simfinal_fmsm_3_st_90d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.25),
simfinal_fmsm_3_st_4m%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.333),
simfinal_fmsm_3_st_6m%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.5),
simfinal_fmsm_3_st_1y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=1),
simfinal_fmsm_3_st_3y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=3),
simfinal_fmsm_3_st_5y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=5))
simfinal_fmsm_4_st_df<-
rbind.data.frame(
simfinal_fmsm_4_st_30d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.08),
simfinal_fmsm_4_st_90d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.25),
simfinal_fmsm_4_st_4m%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.333),
simfinal_fmsm_4_st_6m%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.5),
simfinal_fmsm_4_st_1y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=1),
simfinal_fmsm_4_st_3y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=3),
simfinal_fmsm_4_st_5y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=5))#time to that outcome for people who experience
#itempo máximo al que simular, de manera que el resumen se toman de la muestra de individuos en los datos simulados quienes están en el estado absorvente a ese tiempo.
simfinal_fmsm_3_st_df %>%
dplyr::mutate(year=factor(year,levels=c(.08,.25, .333, .5, 1, 3, 5), labels=c("30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::group_by(year, tipo_de_programa_2) %>%
summarise(mean=sum(val[quantity=="mean"],na.rm=T),
lo=sum(lower[quantity=="2.5%"],na.rm=T),
up=sum(upper[quantity=="97.5%"],na.rm=T)) %>%
ggplot(aes(x=year, y=mean, color=tipo_de_programa_2))+
geom_bar(aes(x=year, y=mean, color=tipo_de_programa_2, fill=tipo_de_programa_2), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=year, ymin=lo, ymax=up, color=tipo_de_programa_2), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
ylim(0,11)+
ylab("Time to event")+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women-specific","General Population"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women-specific","General Population"))+
xlab("Maximum time of simualted individuals in the absorbing state")## `summarise()` regrouping output by 'year' (override with `.groups` argument)
Figure 6a. Time to Readmission for Users who Experience it (Mean and 95%CIs), Three-states model
simfinal_fmsm_4_st_df %>%
dplyr::mutate(year=factor(year,levels=c(.08,.25, .333, .5, 1, 3, 5), labels=c("30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::group_by(year, tipo_de_programa_2) %>%
summarise(mean=sum(val[quantity=="mean"],na.rm=T),
lo=sum(lower[quantity=="2.5%"],na.rm=T),
up=sum(upper[quantity=="97.5%"],na.rm=T)) %>%
ggplot(aes(x=year, y=mean, color=tipo_de_programa_2))+
geom_bar(aes(x=year, y=mean, color=tipo_de_programa_2, fill=tipo_de_programa_2), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=year, ymin=lo, ymax=up, color=tipo_de_programa_2), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1)+
ylim(0,11)+
ylab("Time to event")+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women-specific","General Population"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women-specific","General Population"))+
xlab("Maximum time of simualted individuals in the absorbing state")## `summarise()` regrouping output by 'year' (override with `.groups` argument)
Figure 6b. Time to Readmission for Users who Experience it (Mean and 95%CIs), Four-states model
Other patient
simulate_other_patient_before<-Sys.time()
#30-39, 2-Completed high school or less, Alcohol, Daily, 2-Moderate, Stays temporarily with a relative, One additional substance, Have children (Dichotomized) = Yes (%)
# Database to contrast adjustments
newdat3a_2 <- data.table::data.table(tipo_de_programa_2= factor(c(rep(1,1*n_trans))),
#comp_status= factor(rep(c("Therapeutic discharge","Discharge without clinical advice"),2)),
edad_al_ing_grupos= factor(rep("30-39",1*n_trans)),
escolaridad_rec= factor(rep("2-Completed high school or less",1*n_trans)),
sus_principal_mod= factor(rep("Alcohol",1*n_trans)),
freq_cons_sus_prin= factor(rep("Daily",1*n_trans)),
compromiso_biopsicosocial= factor(rep("2-Moderate",1*n_trans)),
tenencia_de_la_vivienda_mod= factor(rep("Stays temporarily with a relative",1*n_trans)),
num_otras_sus_mod= factor(rep("One additional substance",1*n_trans)),
numero_de_hijos_mod_rec= factor(rep("Yes",1*n_trans)),
tipo_de_plan_res= factor(rep("Residential",1*n_trans)),
strata= rep(1:n_trans,1),
arrival=rep(0,)
)
newdat3b_2 <- data.table::data.table(tipo_de_programa_2= factor(c(rep(0,1*n_trans))),
#comp_status= factor(rep(c("Therapeutic discharge","Discharge without clinical advice"),2)),
edad_al_ing_grupos= factor(rep("30-39",1*n_trans)),
escolaridad_rec= factor(rep("2-Completed high school or less",1*n_trans)),
sus_principal_mod= factor(rep("Alcohol",1*n_trans)),
freq_cons_sus_prin= factor(rep("Daily",1*n_trans)),
compromiso_biopsicosocial= factor(rep("2-Moderate",1*n_trans)),
tenencia_de_la_vivienda_mod= factor(rep("Stays temporarily with a relative",1*n_trans)),
num_otras_sus_mod= factor(rep("One additional substance",1*n_trans)),
numero_de_hijos_mod_rec= factor(rep("Yes",1*n_trans)),
tipo_de_plan_res= factor(rep("Residential",1*n_trans)),
strata= rep(1:n_trans,1),
arrival=rep(0,)
)
newdat5a_2 <- data.table::data.table(tipo_de_programa_2= factor(c(rep(1,1*n_trans2))),
#comp_status= factor(rep(c("Therapeutic discharge","Discharge without clinical advice"),2)),
edad_al_ing_grupos= factor(rep("30-39",1*n_trans2)),
escolaridad_rec= factor(rep("2-Completed high school or less",1*n_trans2)),
sus_principal_mod= factor(rep("Alcohol",1*n_trans2)),
freq_cons_sus_prin= factor(rep("Daily",1*n_trans2)),
compromiso_biopsicosocial= factor(rep("2-Moderate",1*n_trans2)),
tenencia_de_la_vivienda_mod= factor(rep("Stays temporarily with a relative",1*n_trans2)),
num_otras_sus_mod= factor(rep("One additional substance",1*n_trans2)),
numero_de_hijos_mod_rec= factor(rep("Yes",1*n_trans2)),
tipo_de_plan_res= factor(rep("Residential",1*n_trans2)),
strata= rep(1:n_trans2,1),
arrival=rep(0,)
)
newdat5b_2 <- data.table::data.table(tipo_de_programa_2= factor(c(rep(0,1*n_trans2))),
#comp_status= factor(rep(c("Therapeutic discharge","Discharge without clinical advice"),2)),
edad_al_ing_grupos= factor(rep("30-39",1*n_trans2)),
escolaridad_rec= factor(rep("2-Completed high school or less",1*n_trans2)),
sus_principal_mod= factor(rep("Alcohol",1*n_trans2)),
freq_cons_sus_prin= factor(rep("Daily",1*n_trans2)),
compromiso_biopsicosocial= factor(rep("2-Moderate",1*n_trans2)),
tenencia_de_la_vivienda_mod= factor(rep("Stays temporarily with a relative",1*n_trans2)),
num_otras_sus_mod= factor(rep("One additional substance",1*n_trans2)),
numero_de_hijos_mod_rec= factor(rep("Yes",1*n_trans2)),
tipo_de_plan_res= factor(rep("Residential",1*n_trans2)),
strata= rep(1:n_trans2,1),
arrival=rep(0,)
)
flexsurv_cumhaz_a2 <- flexsurv::msfit.flexsurvreg(sel_param_fits_b,#sel_param_fits,
newdata = newdat3a_2,
t = seq(.001, max_time_a, by = .01),
B = n_iter,
trans = mat_3_states, variance = FALSE)
flexsurv_cumhaz_b2 <- flexsurv::msfit.flexsurvreg(sel_param_fits_b,#sel_param_fits,
newdata = newdat3b_2,
t = seq(.001, max_time_b, by = .01),
B = n_iter,
trans = mat_3_states, variance = FALSE)
flexsurv_cumhaz_c2 <- flexsurv::msfit.flexsurvreg(sel_param_fits_4s_b, # sel_param_fits_4s_b,#sel_param_fits_4s,
newdata = newdat5a_2,
t = seq(.001, max_time_c, by = .01),
B = n_iter,
trans = mat_4_states, variance = FALSE)
flexsurv_cumhaz_d2 <- flexsurv::msfit.flexsurvreg(sel_param_fits_4s_b,#sel_param_fits_4s_b,#sel_param_fits_4s, #Alternatively, if the parameters (including covariate effects) are assumed to be different between different transitions, then a list of transition-specific models can be formed. This list has one component for each permitted transition in the multi-state model. This is more computationally efficient, particularly for larger models and datasets. See the example below, and the vignette.
newdata = newdat5b_2, #A data frame specifying the values of covariates in the fitted model, other than the transition number. This must be specified if there are other covariates. The variable names should be the same as those in the fitted model formula. There must be either one value per covariate (the typical situation) or n values per covariate, a different one for each of the n allowed transitions.
t = seq(.001, max_time_d, by = .01), #Vector of times. These do not need to be the same as the observed event times, and since the model is parametric, they can be outside the range of the data. A grid of more frequent times will provide a better approximation to the cumulative hazard trajectory for prediction with probtrans or mssample, at the cost of greater computational expense.
B = n_iter,#Number of simulations from the normal asymptotic distribution used to calculate variances. Decrease for greater speed at the expense of accuracy.
trans = mat_4_states,
variance = FALSE)#Calculate the variances and covariances of the transition cumulative hazards (TRUE or FALSE). This is based on simulation from the normal asymptotic distribution of the estimates, which is computationally-expensive.
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#_#_#_#_#_#_#_#_#_#_PMATRIX
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::pmatrix.simfs(x = sel_param_fits_b,
t = i, #Must be a single number unlike pmatrix.fs
ci = T,
#B= n_iter/5000, #Increasing B is usually more expensive than increasing M. NUmber of simulations from the noral asymptotic distribution. MOre, more accuracy
newdat= newdat3a_2,
tvar="trans", # variable in the data representing the transition type. NOt required if x is a list of models
#tcovs= "arrival",# Predicatable time-dependent covariates such as age
cores= 8,
#M = 1,#Number of individual trajectories to simulate.
trans = mat_3_states) %>%
assign(paste0("pmatrix2_t_a_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::pmatrix.simfs(x = sel_param_fits_b,
t = i,
ci = T,
#B= n_iter/5000,
newdat= newdat3b_2,
tvar="trans",
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8,#N of processor cores
trans = mat_3_states) %>%
assign(paste0("pmatrix2_t_b_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::pmatrix.simfs(x = sel_param_fits_4s_b,
t = i, #Must be a single number unlike pmatrix.fs
ci = T,
#B= n_iter/5000, #Increasing B is usually more expensive than increasing M. NUmber of simulations from the noral asymptotic distribution. MOre, more accuracy
newdat= newdat5a_2,
tvar="trans",
#tcovs= "arrival",# Predicatable time-dependent covariates such as age
cores= 8,
#M = 1,#Number of individual trajectories to simulate.
trans = mat_4_states) %>%
assign(paste0("pmatrix2_4s_t_a_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::pmatrix.simfs(x = sel_param_fits_4s_b,
t = i,
ci = T,
# B= n_iter/5000,
newdat= newdat5b_2,
tvar="trans",
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8,#N of processor cores
trans = mat_4_states) %>%
assign(paste0("pmatrix2_4s_t_b_",i),.,envir=.GlobalEnv)
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#_#_#_#_#_#_#_#_#_#_LOS
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(sel_param_fits_b,
t = i, #Maximum time to predict to.
# start = 1, #Starting state.
newdata = newdat3a_2,
trans= mat_3_states,
ci = T, #Return a confidence interval calculated by simulating from the asymptotic normal distribution of the maximum likelihood estimates. This is turned off by default, since two levels of simulation are required. If turned on, users should adjust B and/or M until the results reach the desired precision. The simulation over M is generally vectorised, therefore increasing B is usually more expensive than increasing M.
tvar = "trans", #Variable in the data representing the transition type. Not required if x is a list of models.
#tcovs = "arrival", #Predictable time-dependent covariates such as age, see sim.fmsm.
group = NULL, #Optional grouping for the states. For example, if there are four states, and group=c(1,1,2,2), then totlos.simfs returns the expected total time in states 1 and 2 combined, and states 3 and 4 combined
#M = 1, #Number of individuals to simulate in order to approximate the transition probabilities. Users should adjust this to obtain the required precision.
#B = 1, #Number of simulations from the normal asymptotic distribution used to calculate variances. Decrease for greater speed at the expense of accuracy.
cl = 0.95)%>%
assign(paste0("tolos2_t_a_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_b,
t = i,
ci = T,
#B= n_iter/5000,
newdat= newdat3b_2,
trans = mat_3_states,
tvar="trans",
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8#N of processor cores
) %>%
assign(paste0("tolos2_t_b_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i, #Must be a single number unlike pmatrix.fs
ci = T,
#B= n_iter/5000, #Increasing B is usually more expensive than increasing M. NUmber of simulations from the noral asymptotic distribution. MOre, more accuracy
newdat= newdat5a_2,
tvar="trans",
#tcovs= "arrival",# Predicatable time-dependent covariates such as age
cores= 8,
#M = 1,#Number of individual trajectories to simulate.
trans = mat_4_states) %>%
assign(paste0("tolos2_4s_t_a_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i,
ci = T,
# B= n_iter/5000,
newdat= newdat5b_2,
trans = mat_4_states,
tvar="trans",
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8#N of processor cores
) %>%
assign(paste0("tolos2_4s_t_b_",i),.,envir=.GlobalEnv)
}
## For intermediate states
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(sel_param_fits_b,
t = i, #Maximum time to predict to.
# start = 1, #Starting state.
newdata = newdat3a_2,
trans= mat_3_states,
start= 2,
ci = T, #Return a confidence interval calculated by simulating from the asymptotic normal distribution of the maximum likelihood estimates. This is turned off by default, since two levels of simulation are required. If turned on, users should adjust B and/or M until the results reach the desired precision. The simulation over M is generally vectorised, therefore increasing B is usually more expensive than increasing M.
tvar = "trans", #Variable in the data representing the transition type. Not required if x is a list of models.
#tcovs = "arrival", #Predictable time-dependent covariates such as age, see sim.fmsm.
group = NULL, #Optional grouping for the states. For example, if there are four states, and group=c(1,1,2,2), then totlos.simfs returns the expected total time in states 1 and 2 combined, and states 3 and 4 combined
#M = 1, #Number of individuals to simulate in order to approximate the transition probabilities. Users should adjust this to obtain the required precision.
#B = 1, #Number of simulations from the normal asymptotic distribution used to calculate variances. Decrease for greater speed at the expense of accuracy.
cl = 0.95)%>%
assign(paste0("tolos2_t_a_str2_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_b,
t = i,
ci = T,
start= 2,
#B= n_iter/5000,
newdat= newdat3b_2,
trans = mat_3_states,
tvar="trans",
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8#N of processor cores
) %>%
assign(paste0("tolos2_t_b_str2_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i, #Must be a single number unlike pmatrix.fs
ci = T,
#B= n_iter/5000, #Increasing B is usually more expensive than increasing M. NUmber of simulations from the noral asymptotic distribution. MOre, more accuracy
newdat= newdat5a_2,
tvar="trans",
start= 2,
#tcovs= "arrival",# Predicatable time-dependent covariates such as age
cores= 8,
#M = 1,#Number of individual trajectories to simulate.
trans = mat_4_states) %>%
assign(paste0("tolos2_4s_t_a_str2_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i,
ci = T,
# B= n_iter/5000,
newdat= newdat5b_2,
trans = mat_4_states,
tvar="trans",
start= 2,
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8#N of processor cores
) %>%
assign(paste0("tolos2_4s_t_b_str2_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i, #Must be a single number unlike pmatrix.fs
ci = T,
#B= n_iter/5000, #Increasing B is usually more expensive than increasing M. NUmber of simulations from the noral asymptotic distribution. MOre, more accuracy
newdat= newdat5a_2,
tvar="trans",
start= 3,
#tcovs= "arrival",# Predicatable time-dependent covariates such as age
cores= 8,
#M = 1,#Number of individual trajectories to simulate.
trans = mat_4_states) %>%
assign(paste0("tolos2_4s_t_a_str3_",i),.,envir=.GlobalEnv)
}
for(i in c(.25,1,3)){
set.seed(1234)
flexsurv::totlos.simfs(x = sel_param_fits_4s_b,
t = i,
ci = T,
# B= n_iter/5000,
newdat= newdat5b_2,
trans = mat_4_states,
tvar="trans",
start= 3,
#tcovs= "arrival",
#M = 1,#Number of individual trajectories to simulate.
cores= 8#N of processor cores
) %>%
assign(paste0("tolos2_4s_t_b_str3_",i),.,envir=.GlobalEnv)
}Cumulative Baseline Hazards
#Example by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5868723/
#library(mstate)
# Semi-parametric
mrcox_a2 <- mstate::msfit(object = crcox,
newdata= newdat3a_2,
trans = mat_3_states)
mrcox_b2 <- mstate::msfit(object = crcox,
newdata= newdat3b_2,
trans = mat_3_states)mrcox_c2 <- mstate::msfit(object = crcox_4s,
newdata= newdat5a_2,
trans = mat_4_states)
mrcox_d2 <- mstate::msfit(object = crcox_4s,
newdata= newdat5b_2,
trans = mat_4_states)`%>%`<-magrittr::`%>%`
# 3 states -Women-specific
cumhaz_data_a2 <- rbind(data.frame(cox_cumhaz_0_a_we$Haz, #2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_a2$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_a2$Haz,
model = "Parametric"))
cumhaz_data_a2$trans <- factor(cumhaz_data_a2$trans,
levels = seq(1, 3),
labels = c("Admission -> Therapeutic Discharge",
"Admission -> Readmission",
"Therapeutic Discharge -> Readmission"))
# 3 states -General-population
cumhaz_data_b2 <- rbind(data.frame(cox_cumhaz_0_a_gp$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_b2$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_b2$Haz,
model = "Parametric"))
cumhaz_data_b2$trans <- factor(cumhaz_data_b2$trans,
levels = seq(1, 3),
labels = c("Admission -> Therapeutic Discharge",
"Admission -> Readmission",
"Therapeutic Discharge -> Readmission"))
## 4 states -Women specific
cumhaz_data_c2 <- rbind(data.frame(cox_cumhaz_0_c_we$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_c2$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_c2$Haz,
model = "Parametric"))
cumhaz_data_c2$trans <- factor(cumhaz_data_c2$trans,
levels = seq(1, 5),
labels = c("Admission -> Therapeutic Discharge",
"Admission -> Discharge w/o Clinical Advice",
"Admission -> Readmission",
"Therapeutic Discharge -> Readmission",
"Discharge w/o Clinical Advice -> Readmission"
))
## 4 states -General population
cumhaz_data_d2 <- rbind(data.frame(cox_cumhaz_0_c_gp$Haz,#2021-04-23, cambié el modelo NP, porque estaba ws donde debiese estar el general
model = "NP Cox"),
data.frame(mrcox_d2$Haz,
model = "Cox"),
data.frame(flexsurv_cumhaz_d2$Haz,
model = "Parametric"))
cumhaz_data_d2$trans <- factor(cumhaz_data_d2$trans,
levels = seq(1, 5),
labels = c("Admission -> Therapeutic Discharge",
"Admission -> Discharge w/o Clinical Advice",
"Admission -> Readmission",
"Therapeutic Discharge -> Readmission",
"Discharge w/o Clinical Advice -> Readmission"
))
fig_cox_ab22<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_a2,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_b2,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="solid"),size=1, alpha=.65) +
facet_wrap(trans~., ncol=1, scales = "free_y") +
scale_color_manual(name ="Model",values=c("#0E53A7","#200772","#BF8830"),labels= c("Cox","NP Cox","Par"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only ","Mixed gender "))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("Years") +
ylab("Cumulative hazards") +
theme_minimal()+
theme(legend.position=c(.9,.59),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
fig_cox_cd22<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_c2,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_d2,id.vars=c("time","trans","model")), aes(time, value, color= model,linetype="solid"),size=1, alpha=.65) +
facet_wrap(trans~., ncol=1, scales = "free_y") +
scale_color_manual(name ="Model",values=c("#0E53A7","#200772","#BF8830"),labels= c("Cox","NP Cox","Par"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only ","Mixed gender "))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("Years") +
ylab("Cumulative") +
theme_minimal()+
theme(legend.position=c(.9,.6),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))trans_for_three_st<-c("Admission -> Therapeutic Discharge","Admission -> Readmission","Therapeutic Discharge -> Readmission")
trans_for_four_st<-c("Admission -> Therapeutic Discharge","Admission -> Discharge w/o Clinical Advice","Admission -> Readmission","Therapeutic Discharge -> Readmission","Discharge w/o Clinical Advice -> Readmission")
abc_let<-c("A","B","C","A","B","C","D","E")
st_plot<-list()
for (i in 1:length(trans_for_three_st)){
st_plot[[i]]<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_a2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_three_st[[i]]), aes(time, value, linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_b2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_three_st[[i]]), aes(time, value, linetype="solid"),size=1, alpha=.65) +
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("") +
ylab("") +
theme_minimal()+
theme(axis.title = element_blank())+
theme(legend.position="none",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))+
ggtitle(paste0(abc_let[[i]],") ",trans_for_three_st[[i]]))
}
for (i in 4:(length(trans_for_four_st)+3)){
st_plot[[i]]<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_c2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[i-3]]), aes(time, value, linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_d2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[i-3]]), aes(time, value, linetype="solid"),size=1, alpha=.65) +
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("") +
ylab("") +
theme_minimal()+
theme(axis.title = element_blank())+
theme(legend.position="none",
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))+
ggtitle(paste0(abc_let[[i]],") ",trans_for_four_st[[i-3]]))
}
i<-8
st_plot[[8]]<-
ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_c2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[5]]), aes(time, value, linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_d2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[5]]), aes(time, value, linetype="solid"),size=1, alpha=.65) +
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
scale_x_continuous(breaks=seq(0, 11, by = 1/2))+
xlab("") +
ylab("") +
theme_minimal()+
theme(axis.title = element_blank())+
ggtitle(paste0("E) ",trans_for_four_st[[5]]))
p3 <- ggplot()+
geom_step(data=reshape2::melt(cumhaz_data_c2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[4]]), aes(time, value, linetype="dashed"),size=1, alpha=.65) +
geom_step(data=reshape2::melt(cumhaz_data_d2,id.vars=c("time","trans","model")) %>% dplyr::filter(model=="Parametric") %>% dplyr::filter(trans==trans_for_four_st[[4]]), aes(time, value, linetype="solid"),size=1, alpha=.65) +
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
lims(x = c(0,0), y = c(0,0))+
theme_void()+
theme(legend.position = c(0.5,0.5),
legend.key.size = unit(1, "cm"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 15, face = "bold"))+
guides(colour = guide_legend(override.aes = list(size=9)))
fig_cox_abcd<-
ggarrange(
# st_plot[[1]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[2]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[3]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[4]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[5]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[6]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[7]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[8]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
p3
)## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 631 row(s) containing missing values (geom_path).## Warning: Removed 600 row(s) containing missing values (geom_path).## Warning: Removed 631 row(s) containing missing values (geom_path).## Warning: Removed 600 row(s) containing missing values (geom_path).## Warning: Removed 631 row(s) containing missing values (geom_path).## Warning: Removed 600 row(s) containing missing values (geom_path).## Warning: Removed 631 row(s) containing missing values (geom_path).## Warning: Removed 600 row(s) containing missing values (geom_path).## Warning: Removed 631 row(s) containing missing values (geom_path).## Warning: Removed 600 row(s) containing missing values (geom_path).## Warning: Removed 1131 row(s) containing missing values (geom_path).## Warning: Removed 1100 row(s) containing missing values (geom_path).fig_cox_abcd_three_states<-
ggarrange(
st_plot[[1]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[2]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot[[3]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[4]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[5]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[6]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[7]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
# st_plot[[8]] + theme(legend.position="none")+xlim(0,5)+ylim(0,2)+theme(plot.title = element_text(size = 10, face = "bold")),
p3
)## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 706 row(s) containing missing values (geom_path).## Warning: Removed 623 row(s) containing missing values (geom_path).## Warning: Removed 706 row(s) containing missing values (geom_path).## Warning: Removed 623 row(s) containing missing values (geom_path).## Warning: Removed 706 row(s) containing missing values (geom_path).## Warning: Removed 623 row(s) containing missing values (geom_path).## Warning: Removed 1131 row(s) containing missing values (geom_path).## Warning: Removed 1100 row(s) containing missing values (geom_path).fig_cox_abcd_lab<-
annotate_figure(fig_cox_abcd, left = textGrob("Cumulative hazards", rot = 90, vjust = 1, gp = gpar(cex = 1.3)),
bottom = textGrob("Years", gp = gpar(cex = 1.3)))
fig_cox_abcd_three_states_lab<-
annotate_figure(fig_cox_abcd_three_states, left = textGrob("Cumulative hazards", rot = 90, vjust = 1, gp = gpar(cex = 1.3)),
bottom = textGrob("Years", gp = gpar(cex = 1.3)))
fig_cox_abcd_lab
Figure 7a. Estimate of Cumulative Hazards, Three & Four-States Model
path_jpg<-rstudioapi::getSourceEditorContext()$path
if (grepl("CISS Fondecyt",path_jpg)==T){
jpg_path<-paste0("C:/Users/CISS Fondecyt/OneDrive/Escritorio/SUD_CL/")
} else if (grepl("andre",path_jpg)==T){
jpg_path<-paste0('C:/Users/andre/Desktop/SUD_CL/')
} else if (grepl("E:",path_jpg)==T){
jpg_path<-paste0("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/")
} else {
jpg_path<-paste0("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/")
}
if(no_mostrar==1){
ggsave(paste0(jpg_path,"Fig2_paper_carla.jpg"), fig_cox_abcd_lab, width = 10, height = 11.5, dpi = 600, units= "in")
ggsave(paste0(jpg_path,"Fig2_paper_carla.eps"), fig_cox_abcd_lab, width = 10, height = 11.5, dpi = 600, units= "in")
ggsave(paste0(jpg_path,"Fig2_paper_carla.ps"), fig_cox_abcd_lab, width = 10, height = 11.5, dpi = 600, units= "in")
ggsave(paste0(jpg_path,"Fig2_alt_paper_carla.jpg"), fig_cox_abcd_three_states_lab, width = 10, height = 11.5, dpi = 600, units= "in")
ggsave(paste0(jpg_path,"Fig2_alt_paper_carla.eps"), fig_cox_abcd_three_states_lab, width = 10, height = 11.5, dpi = 600, units= "in")
ggsave(paste0(jpg_path,"Fig2_alt_paper_carla.ps"), fig_cox_abcd_three_states_lab, width = 10, height = 11.5, dpi = 600, units= "in")
}ggarrange(fig_cox_ab22+xlim(0,4),fig_cox_cd22+xlim(0,4), labels = "AUTO")## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 1564 row(s) containing missing values (geom_path).## Warning: Removed 1627 row(s) containing missing values (geom_path).## Warning: Removed 1397 row(s) containing missing values (geom_path).## Warning: Removed 1490 row(s) containing missing values (geom_path).
Figure 7b. Estimate of Cumulative Hazards, Three & Four-States Model
if(no_mostrar==1){
ggsave("eso13ab2.jpg", ggarrange(cumhaz_plot13+xlim(0,4),cumhaz_plot23+xlim(0,4), labels = "AUTO"), width = 10, height = 11.5, dpi = 600, units= "in")
dev.off()
}Transition probabilities
simulate_other_patient_after<-Sys.time()
paste0("Time taken in process:");simulate_other_patient_after-simulate_other_patient_before## [1] "Time taken in process:"## Time difference of 3.94604 hours#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#This functionality is implemented in pmatrix.simfs,
#but the tcovs argument actually has no impact on the transition probabilities, as evidenced below.
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
pmatrix2_t_a<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
get(paste0("pmatrix2_t_a_",t))
}))
pmatrix2_t_b<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
get(paste0("pmatrix2_t_b_",t))
}))
pmatrix2_4s_t_a<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
get(paste0("pmatrix2_4s_t_a_",t))
}))
pmatrix2_4s_t_b<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
get(paste0("pmatrix2_4s_t_b_",t))
}))
pmatrix2_t_a_lo<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_t_a_",t)),"lower")
}))
pmatrix2_t_b_lo<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_t_b_",t)),"lower")
}))
pmatrix2_4s_t_a_lo<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_4s_t_a_",t)),"lower")
}))
pmatrix2_4s_t_b_lo<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_4s_t_b_",t)),"lower")
}))
pmatrix2_t_a_up<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_t_a_",t)),"upper")
}))
pmatrix2_t_b_up<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_t_b_",t)),"upper")
}))
pmatrix2_4s_t_a_up<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_4s_t_a_",t)),"upper")
}))
pmatrix2_4s_t_b_up<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
attr(get(paste0("pmatrix2_4s_t_b_",t)),"upper")
}))
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
pmatrix2_t_a_df_final<-
cbind.data.frame(t=rep(c(.25,1,3),each=dim(mat_3_states)[2]),
trans=rep(1:dim(mat_3_states)[2]),
pmatrix2_t_a,pmatrix2_t_a_lo,pmatrix2_t_a_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]"))
pmatrix2_t_b_df_final<-
cbind.data.frame(t=rep(c(.25,1,3),each=dim(mat_3_states)[2]),
trans=rep(1:dim(mat_3_states)[2]),
pmatrix2_t_b,pmatrix2_t_b_lo,pmatrix2_t_b_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]"))
pmatrix2_4s_t_a_df_final<-
cbind.data.frame(t=rep(c(.25,1,3),each=dim(mat_4_states)[2]),
trans=rep(1:dim(mat_4_states)[2]),
pmatrix2_4s_t_a,pmatrix2_4s_t_a_lo,pmatrix2_4s_t_a_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3","Trans_4"="X4",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1","Trans_4_lo"="X4.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2","Trans_4_up"="X4.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]")) %>%
dplyr::mutate(Trans_4_cis=paste0(sprintf("%1.2f",Trans_4)," [",
sprintf("%1.2f",Trans_4_lo),"-",
sprintf("%1.2f",Trans_4_up),"]"))
pmatrix2_4s_t_b_df_final<-
cbind.data.frame(t=rep(c(.25,1,3),each=dim(mat_4_states)[2]),
trans=rep(1:dim(mat_4_states)[2]),
pmatrix2_4s_t_b,pmatrix2_4s_t_b_lo,pmatrix2_4s_t_b_up) %>%
data.frame() %>%
dplyr::rename("Trans_1"="X1","Trans_2"="X2","Trans_3"="X3","Trans_4"="X4",
"Trans_1_lo"="X1.1","Trans_2_lo"="X2.1","Trans_3_lo"="X3.1","Trans_4_lo"="X4.1",
"Trans_1_up"="X1.2","Trans_2_up"="X2.2","Trans_3_up"="X3.2","Trans_4_up"="X4.2") %>%
dplyr::mutate(Trans_1_cis=paste0(sprintf("%1.2f",Trans_1)," [",
sprintf("%1.2f",Trans_1_lo),"-",
sprintf("%1.2f",Trans_1_up),"]")) %>%
dplyr::mutate(Trans_2_cis=paste0(sprintf("%1.2f",Trans_2)," [",
sprintf("%1.2f",Trans_2_lo),"-",
sprintf("%1.2f",Trans_2_up),"]")) %>%
dplyr::mutate(Trans_3_cis=paste0(sprintf("%1.2f",Trans_3)," [",
sprintf("%1.2f",Trans_3_lo),"-",
sprintf("%1.2f",Trans_3_up),"]")) %>%
dplyr::mutate(Trans_4_cis=paste0(sprintf("%1.2f",Trans_4)," [",
sprintf("%1.2f",Trans_4_lo),"-",
sprintf("%1.2f",Trans_4_up),"]"))pmatrix2_t_b_df_final_join<-
pmatrix2_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years")))
pmatrix2_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years")))%>%
dplyr::left_join(pmatrix2_t_b_df_final_join,by=c("t","trans")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
trans==3~"3)Readmission")) %>%
dplyr::select(-t) %>%
dplyr::select(`trans`,`Trans_1_cis.x`,`Trans_2_cis.x`,`Trans_3_cis.x`,`Trans_1_cis.y`,`Trans_2_cis.y`,`Trans_3_cis.y`) %>%
#2021-04-27: drop from/actual state, readmission
dplyr::filter(!grepl("Readmission",trans)) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 5a. Estimated transition probabilities, Three-states model"),
col.names = c("Actual state", rep(c("Admission","Therapeutic Discharge","Readmission"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Women-specific: Transition to" = 3,"General Population: Transition to" = 3)) %>%
kableExtra::pack_rows("~91 days/3 months", 1, 2) %>%
kableExtra::pack_rows("At ~1 year", 3, 4) %>%
kableExtra::pack_rows("At ~3 years", 5, 6) %>%
kableExtra::scroll_box(width = "100%", height = "375px")| Actual state | Admission | Therapeutic Discharge | Readmission | Admission | Therapeutic Discharge | Readmission |
|---|---|---|---|---|---|---|
| ~91 days/3 months | ||||||
| 1)Admission | 0.87 [0.84-0.89] | 0.12 [0.10-0.14] | 0.02 [0.01-0.02] | 0.89 [0.87-0.91] | 0.10 [0.08-0.12] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.87 [0.80-0.91] | 0.13 [0.09-0.20] | 0.00 [0.00-0.00] | 0.87 [0.81-0.92] | 0.13 [0.08-0.19] |
| At ~1 year | ||||||
| 1)Admission | 0.60 [0.55-0.65] | 0.29 [0.25-0.34] | 0.11 [0.08-0.14] | 0.66 [0.62-0.71] | 0.24 [0.21-0.29] | 0.09 [0.07-0.12] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.74 [0.65-0.82] | 0.26 [0.18-0.35] | 0.00 [0.00-0.00] | 0.75 [0.66-0.83] | 0.25 [0.17-0.34] |
| At ~3 years | ||||||
| 1)Admission | 0.35 [0.30-0.40] | 0.35 [0.29-0.42] | 0.29 [0.24-0.35] | 0.43 [0.37-0.48] | 0.31 [0.25-0.38] | 0.26 [0.21-0.32] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.61 [0.50-0.70] | 0.39 [0.30-0.50] | 0.00 [0.00-0.00] | 0.62 [0.49-0.72] | 0.38 [0.28-0.51] |
| Note. We removed Readmission because it was an absorbing state | ||||||
pmatrix2_4s_t_b_df_final_join<-
pmatrix2_4s_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years")))
pmatrix2_4s_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, ends_with("cis")) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::left_join(pmatrix2_4s_t_b_df_final_join,by=c("t","trans")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
trans==3~"3)From Discharge Without\nClinical Advice",
trans==4~"4)Readmission")) %>%
dplyr::select(-t) %>%
dplyr::select(trans,Trans_1_cis.x,Trans_2_cis.x,Trans_3_cis.x,Trans_4_cis.x,Trans_1_cis.y,Trans_2_cis.y,Trans_3_cis.y,Trans_4_cis.y) %>%
#2021-04-27: drop from/actual state, readmission
dplyr::filter(!grepl("Readmission",trans)) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 5b. Estimated transition probabilities, Four-states model; patient 2"),
col.names = c("Actual state", rep(c("Admission","Therapeutic Discharge","Discharge Without\nClinical Advice","Readmission"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_footnote("Note. We removed Readmission because it was an absorbing state", notation="none") %>%
kableExtra::add_header_above(c(" ", "Women-specific: Transition to" = 4,"General Population: Transition to" = 4)) %>%
kableExtra::pack_rows("At ~91 days/3 months", 1, 3) %>%
kableExtra::pack_rows("At ~1 year", 4, 6) %>%
kableExtra::pack_rows("At ~3 years", 7, 9) %>%
kableExtra::scroll_box(width = "100%", height = "375px") | Actual state | Admission | Therapeutic Discharge | Discharge Without Clinical Advice | Readmission | Admission | Therapeutic Discharge | Discharge Without Clinical Advice | Readmission |
|---|---|---|---|---|---|---|---|---|
| At ~91 days/3 months | ||||||||
| 1)Admission | 0.79 [0.77-0.81] | 0.10 [0.09-0.12] | 0.09 [0.08-0.10] | 0.01 [0.01-0.02] | 0.81 [0.79-0.83] | 0.08 [0.07-0.10] | 0.10 [0.09-0.11] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.86 [0.81-0.91] | 0.00 [0.00-0.00] | 0.14 [0.09-0.19] | 0.00 [0.00-0.00] | 0.87 [0.81-0.92] | 0.00 [0.00-0.00] | 0.13 [0.08-0.19] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.95 [0.93-0.96] | 0.05 [0.04-0.07] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.95 [0.93-0.97] | 0.05 [0.03-0.07] |
| At ~1 year | ||||||||
| 1)Admission | 0.45 [0.40-0.48] | 0.26 [0.22-0.30] | 0.20 [0.18-0.23] | 0.09 [0.07-0.12] | 0.49 [0.45-0.53] | 0.21 [0.18-0.25] | 0.22 [0.20-0.25] | 0.08 [0.06-0.10] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.74 [0.64-0.82] | 0.00 [0.00-0.00] | 0.26 [0.18-0.36] | 0.00 [0.00-0.00] | 0.75 [0.66-0.83] | 0.00 [0.00-0.00] | 0.25 [0.17-0.34] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.88 [0.84-0.91] | 0.12 [0.09-0.16] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.89 [0.86-0.92] | 0.11 [0.08-0.14] |
| At ~3 years | ||||||||
| 1)Admission | 0.17 [0.13-0.20] | 0.34 [0.28-0.40] | 0.23 [0.20-0.26] | 0.26 [0.21-0.32] | 0.22 [0.18-0.26] | 0.30 [0.25-0.36] | 0.26 [0.23-0.29] | 0.23 [0.18-0.28] |
| 2)Therapeutic Discharge | 0.00 [0.00-0.00] | 0.60 [0.50-0.70] | 0.00 [0.00-0.00] | 0.40 [0.30-0.50] | 0.00 [0.00-0.00] | 0.62 [0.50-0.73] | 0.00 [0.00-0.00] | 0.38 [0.27-0.50] |
| 3)From Discharge Without Clinical Advice | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.79 [0.74-0.84] | 0.21 [0.16-0.26] | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] | 0.81 [0.76-0.85] | 0.19 [0.15-0.24] |
| Note. We removed Readmission because it was an absorbing state | ||||||||
pmatrix2_t_b_df_final_join2<-
pmatrix2_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_1_up, Trans_2_up, Trans_3_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3") %>%
melt(id.vars=c("t","trans"))
transprob_pat2_a<-
pmatrix2_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_1_up, Trans_2_up, Trans_3_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::left_join(pmatrix2_t_b_df_final_join2,by=c("t","trans","variable")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
trans==3~"3)Readmission")) %>%
tidyr::separate(variable, c("a", "b","c")) %>%
dplyr::mutate(transition_to=paste0(a,"_",b)) %>%
dplyr::mutate(transition_to=dplyr::case_when(transition_to=="Trans_1"~"1)Admission",
transition_to=="Trans_2"~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
transition_to=="Trans_3"~"3)Readmission")) %>%
dplyr::select(-a,-b) %>%
dplyr::rename("trans_from"="trans") %>%
tidyr::pivot_longer(col=c("value.x","value.y")) %>%
tidyr::pivot_wider(names_from="c", values_from=c("value")) %>%
dplyr::mutate(sel=paste0(trans_from,"_",transition_to),
sel=dplyr::case_when(sel %in% c("1)Admission_1)Admission",
"2)Therapeutic\nDischarge_1)Admission",
"3)Readmission_1)Admission",
"2)Therapeutic\nDischarge_2)Therapeutic\nDischarge",
"3)Readmission_2)Therapeutic\nDischarge",
"3)Readmission_3)Readmission")~0,T~1)) %>%
dplyr::filter(sel==1) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=name, fill=name), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=lo, ymax=up, color=name), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(trans_from~transition_to)+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
theme(legend.position = "bottom")+
ylab("")+
xlab("")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
transprob_pat2_a
Figure 8a. Transition Probabilities (and 95% confidence intervals), 2nd patient, Three-states model
if(no_mostrar==1){
ggsave("eso13a.jpg", transprob_pat2_a,width = 12, height = 10, dpi = 300, units= "in")
dev.off()
}pmatrix2_4s_t_b_df_final_join2<-
pmatrix2_4s_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_4, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_4_lo, Trans_1_up, Trans_2_up, Trans_3_up, Trans_4_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3", "Trans_4_est"="Trans_4") %>%
melt(id.vars=c("t","trans"))
transprob_pat2_b<-
pmatrix2_4s_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_4, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_4_lo, Trans_1_up, Trans_2_up, Trans_3_up, Trans_4_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3", "Trans_4_est"="Trans_4") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::left_join(pmatrix2_4s_t_b_df_final_join2,by=c("t","trans","variable")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
trans==3~"3)Discharge Without\nClinical Advice",
trans==4~"4)Readmission")) %>%
tidyr::separate(variable, c("a", "b","c")) %>%
dplyr::mutate(transition_to=paste0(a,"_",b)) %>%
dplyr::mutate(transition_to=dplyr::case_when(transition_to=="Trans_1"~"1)Admission",
transition_to=="Trans_2"~"2)Therapeutic\nDischarge",
transition_to=="Trans_3"~"3)Discharge Without\nClinical Advice",
transition_to=="Trans_4"~"4)Readmission")) %>%
dplyr::select(-a,-b) %>%
dplyr::rename("trans_from"="trans") %>%
tidyr::pivot_longer(col=c("value.x","value.y")) %>%
tidyr::pivot_wider(names_from="c", values_from=c("value")) %>%
dplyr::mutate(sel=paste0(trans_from,"_",transition_to),
sel=dplyr::case_when(sel %in% c("1)Admission_1)Admission",
"2)Therapeutic\nDischarge_1)Admission",
"4)Readmission_1)Admission",
"2)Therapeutic\nDischarge_2)Therapeutic\nDischarge",
"2)Therapeutic\nDischarge_3)Discharge Without\nClinical Advice",
"3)Discharge Without\nClinical Advice_1)Admission",
"3)Discharge Without\nClinical Advice_2)Therapeutic\nDischarge",
"3)Discharge Without\nClinical Advice_3)Discharge Without\nClinical Advice",
"4)Readmission_2)Therapeutic\nDischarge",
"4)Readmission_3)Discharge Without\nClinical Advice",
"4)Readmission_4)Readmission")~0,T~1)) %>%
dplyr::filter(sel==1) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=name, fill=name), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=lo, ymax=up, color=name), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(trans_from~transition_to, ncol=3)+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
theme(legend.position=c(.8,.2))+
ylab("")+
xlab("")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
transprob_pat2_b
Figure 8b. Transition Probabilities (and 95% confidence intervals), 2nd patient, Four-states model
if(no_mostrar==1){
ggsave("eso13b.jpg", transprob_pat2_b,width = 12, height = 10, dpi = 300, units= "in")
dev.off()
}pmatrix_t_b_df_final_join2<-
pmatrix_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_1_up, Trans_2_up, Trans_3_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3") %>%
melt(id.vars=c("t","trans"))
transprob_pat_a<-
pmatrix_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_1_up, Trans_2_up, Trans_3_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::left_join(pmatrix_t_b_df_final_join2,by=c("t","trans","variable")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
trans==3~"3)Readmission")) %>%
tidyr::separate(variable, c("a", "b","c")) %>%
dplyr::mutate(transition_to=paste0(a,"_",b)) %>%
dplyr::mutate(transition_to=dplyr::case_when(transition_to=="Trans_1"~"1)Admission",
transition_to=="Trans_2"~"2)Therapeutic\nDischarge",
#trans==3~"3)From Discharge Without\nClinical Advice",
transition_to=="Trans_3"~"3)Readmission")) %>%
dplyr::select(-a,-b) %>%
dplyr::rename("trans_from"="trans") %>%
tidyr::pivot_longer(col=c("value.x","value.y")) %>%
tidyr::pivot_wider(names_from="c", values_from=c("value")) %>%
dplyr::mutate(sel=paste0(trans_from,"_",transition_to),
sel=dplyr::case_when(sel %in% c("1)Admission_1)Admission",
"2)Therapeutic\nDischarge_1)Admission",
"3)Readmission_1)Admission",
"2)Therapeutic\nDischarge_2)Therapeutic\nDischarge",
"3)Readmission_2)Therapeutic\nDischarge",
"3)Readmission_3)Readmission")~0,T~1)) %>%
dplyr::filter(sel==1) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=name, fill=name), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=lo, ymax=up, color=name), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(trans_from~transition_to)+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
theme(legend.position="bottom")+
ylab("")+
xlab("")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
transprob_pat_a
Figure 9a. Transition Probabilities (and 95% confidence intervals), 1st patient, Three-states model
if(no_mostrar==1){
ggsave("eso13b.jpg", transprob_pat_a,width = 12, height = 10, dpi = 300, units= "in")
dev.off()
}pmatrix_4s_t_b_df_final_join2<-
pmatrix_4s_t_b_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_4, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_4_lo, Trans_1_up, Trans_2_up, Trans_3_up, Trans_4_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3", "Trans_4_est"="Trans_4") %>%
melt(id.vars=c("t","trans"))
transprob_pat_b<-
pmatrix_4s_t_a_df_final %>%
dplyr::mutate(t=round(t,3)) %>%
dplyr::select(t, trans, Trans_1, Trans_2, Trans_3, Trans_4, Trans_1_lo, Trans_2_lo, Trans_3_lo, Trans_4_lo, Trans_1_up, Trans_2_up, Trans_3_up, Trans_4_up) %>%
dplyr::filter(t%in% c(.25, 1, 3)) %>%
dplyr::mutate(t=factor(t, levels=c(.25, 1, 3),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::rename("Trans_1_est"="Trans_1","Trans_2_est"="Trans_2","Trans_3_est"="Trans_3", "Trans_4_est"="Trans_4") %>%
melt(id.vars=c("t","trans")) %>%
dplyr::left_join(pmatrix_4s_t_b_df_final_join2,by=c("t","trans","variable")) %>%
dplyr::mutate(trans=dplyr::case_when(trans==1~"1)Admission",
trans==2~"2)Therapeutic\nDischarge",
trans==3~"3)Discharge Without\nClinical Advice",
trans==4~"4)Readmission")) %>%
tidyr::separate(variable, c("a", "b","c")) %>%
dplyr::mutate(transition_to=paste0(a,"_",b)) %>%
dplyr::mutate(transition_to=dplyr::case_when(transition_to=="Trans_1"~"1)Admission",
transition_to=="Trans_2"~"2)Therapeutic\nDischarge",
transition_to=="Trans_3"~"3)Discharge Without\nClinical Advice",
transition_to=="Trans_4"~"4)Readmission")) %>%
dplyr::select(-a,-b) %>%
dplyr::rename("trans_from"="trans") %>%
tidyr::pivot_longer(col=c("value.x","value.y")) %>%
tidyr::pivot_wider(names_from="c", values_from=c("value")) %>%
dplyr::mutate(sel=paste0(trans_from,"_",transition_to),
sel=dplyr::case_when(sel %in% c("1)Admission_1)Admission",
"2)Therapeutic\nDischarge_1)Admission",
"4)Readmission_1)Admission",
"2)Therapeutic\nDischarge_2)Therapeutic\nDischarge",
"2)Therapeutic\nDischarge_3)Discharge Without\nClinical Advice",
"3)Discharge Without\nClinical Advice_1)Admission",
"3)Discharge Without\nClinical Advice_2)Therapeutic\nDischarge",
"3)Discharge Without\nClinical Advice_3)Discharge Without\nClinical Advice",
"4)Readmission_2)Therapeutic\nDischarge",
"4)Readmission_3)Discharge Without\nClinical Advice",
"4)Readmission_4)Readmission")~0,T~1)) %>%
dplyr::filter(sel==1) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=name, fill=name), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=lo, ymax=up, color=name), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(trans_from~transition_to, ncol=3)+
scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
scale_color_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
scale_fill_manual(name ="Type of program",values=c("#0E53A7","#FFB540"), labels= c("Women only","Mixed gender"))+
ylab("")+
xlab("")+
theme(legend.position=c(.8,.2))+
ylab("")+
xlab("")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
transprob_pat_b
Figure 9b. Transition Probabilities (and 95% confidence intervals), 1st patient, Four-states model
if(no_mostrar==1){
ggsave("eso13b.jpg", transprob_pat_b,width = 12, height = 10, dpi = 300, units= "in")
dev.off()
}Length of Stay
PSA
We first looked over the restricted mean survival time, which may be the equivalent measure of the length of stay for partitioned survival analyses of each transition (as stated by Crowther in this link). We calculated its 95% confidence intervals by 10,000 resamples.
rmst_a_3s_1_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=1,summary(sel_param_fits_b[[1]],newdata=newdat3a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_3s_2_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=2,summary(sel_param_fits_b[[2]],newdata=newdat3a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_3s_3_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=3,summary(sel_param_fits_b[[3]],newdata=newdat3a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_1_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=1,summary(sel_param_fits_b[[1]],newdata=newdat3a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_2_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=2,summary(sel_param_fits_b[[2]],newdata=newdat3b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_b_3s_3_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=3,summary(sel_param_fits_b[[3]],newdata=newdat3b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=3)
}))
rmst_a_4s_1_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=1,summary(sel_param_fits_4s_b[[1]],newdata=newdat5a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_2_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=2,summary(sel_param_fits_4s_b[[2]],newdata=newdat5a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_3_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=3,summary(sel_param_fits_4s_b[[3]],newdata=newdat5a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_4_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=4,summary(sel_param_fits_4s_b[[4]],newdata=newdat5a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_a_4s_5_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="ws",t=t,tr=5,summary(sel_param_fits_4s_b[[5]],newdata=newdat5a_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_1_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=1,summary(sel_param_fits_4s_b[[1]],newdata=newdat5b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_2_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=2,summary(sel_param_fits_4s_b[[2]],newdata=newdat5b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_3_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=3,summary(sel_param_fits_4s_b[[3]],newdata=newdat5b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_4_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=4,summary(sel_param_fits_4s_b[[4]],newdata=newdat5b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_b_4s_5_2<-
do.call('rbind', lapply(seq(0, 11, by = 1/12), function(t) {
cbind(pr="gp",t=t,tr=5,summary(sel_param_fits_4s_b[[5]],newdata=newdat5b_2[1,-"strata"],type="rmst",B=n_iter,tidy=T, start=0,t=t)[,c("est","lcl","ucl")],st=4)
}))
rmst_df_psa_3s_2<-
rbind.data.frame(rmst_a_3s_1_2, rmst_a_3s_2_2, rmst_a_3s_3_2, rmst_b_3s_1_2, rmst_b_3s_2_2, rmst_b_3s_3_2)
rmst_df_psa_4s_2<-
rbind.data.frame(rmst_a_4s_1_2, rmst_a_4s_2_2, rmst_a_4s_3_2, rmst_a_4s_4_2, rmst_a_4s_5_2, rmst_b_4s_1_2, rmst_b_4s_2_2, rmst_b_4s_3_2, rmst_b_4s_4_2, rmst_b_4s_5_2)#In order to estimate performances of average LOS (ALOS) estimation methods according to the accumulating data that become available over time, ALOS with their 95% confidence intervals (CI)
#rmst_df_psa_4s
rmst_df_psa_3s_2 %>%
dplyr::mutate(pr=dplyr::case_when(pr=="ws"~"Women-specific",
T~"General population")) %>%
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",lcl),"-",
sprintf("%1.2f",ucl),"]")) %>%
dplyr::select(pr, t, tr, st, est_ci) %>%
tidyr::pivot_wider(names_from="pr", values_from=c("est_ci")) %>%
dplyr::select(-st) %>%
dplyr::mutate_at(1,~sprintf("%1.4f",.)) %>%
dplyr::filter(t %in% c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000")) %>%
dplyr::mutate(t=factor(t,levels=c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000"), labels=c("0 days","30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::select(-tr) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 6a. Restricted Mean Survival Time by Partitioned Survival Analyses, Three-states model, 2nd patient"),
col.names = c("Time",rep(c("Estimate [95% CI]"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Transition 1: Admission to Ther.Discharge", 1, 8) %>%
kableExtra::pack_rows("Transition 2: Admission to Readmission", 9, 16) %>%
kableExtra::pack_rows("Transition 3: Ther.Discharge to Readmission", 17, 24) %>%
kableExtra::scroll_box(width = "100%", height = "375px")| Time | Estimate [95% CI] | Estimate [95% CI] |
|---|---|---|
| Transition 1: Admission to Ther.Discharge | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.23 [0.23-0.24] | 0.23 [0.23-0.24] |
| 4 months | 0.30 [0.30-0.31] | 0.30 [0.30-0.31] |
| 6 months | 0.44 [0.43-0.45] | 0.44 [0.43-0.45] |
| 1 year | 0.79 [0.75-0.82] | 0.79 [0.75-0.82] |
| 3 years | 1.80 [1.66-1.94] | 1.80 [1.65-1.95] |
| 5 years | 2.61 [2.33-2.88] | 2.61 [2.34-2.88] |
| Transition 2: Admission to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.25-0.25] | 0.25 [0.25-0.25] |
| 4 months | 0.33 [0.33-0.33] | 0.33 [0.33-0.33] |
| 6 months | 0.50 [0.50-0.50] | 0.50 [0.50-0.50] |
| 1 year | 0.98 [0.97-0.99] | 0.98 [0.97-0.99] |
| 3 years | 2.73 [2.66-2.79] | 2.76 [2.69-2.81] |
| 5 years | 4.28 [4.14-4.41] | 4.34 [4.19-4.47] |
| Transition 3: Ther.Discharge to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.23 [0.22-0.24] | 0.23 [0.22-0.24] |
| 4 months | 0.30 [0.28-0.31] | 0.30 [0.28-0.31] |
| 6 months | 0.44 [0.41-0.46] | 0.44 [0.41-0.46] |
| 1 year | 0.82 [0.75-0.88] | 0.83 [0.76-0.89] |
| 3 years | 2.15 [1.87-2.39] | 2.18 [1.89-2.43] |
| 5 years | 3.28 [2.77-3.72] | 3.34 [2.82-3.80] |
#kableExtra::add_footnote("Note. NA= Null values", notation="none") #In order to estimate performances of average LOS (ALOS) estimation methods according to the accumulating data that become available over time, ALOS with their 95% confidence intervals (CI)
#rmst_df_psa_4s
rmst_df_psa_4s_2 %>%
dplyr::mutate(pr=dplyr::case_when(pr=="ws"~"Women-specific",
T~"General population")) %>%
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",lcl),"-",
sprintf("%1.2f",ucl),"]")) %>%
dplyr::select(pr, t, tr, st, est_ci) %>%
tidyr::pivot_wider(names_from="pr", values_from=c("est_ci")) %>%
dplyr::select(-st) %>%
dplyr::mutate_at(1,~sprintf("%1.4f",.)) %>%
dplyr::filter(t %in% c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000")) %>%
dplyr::mutate(t=factor(t,levels=c("0.0000","0.0833", "0.2500", "0.3333", "0.5000", "1.0000", "3.0000", "5.0000"), labels=c("0 days", "30 days", "90 days", "4 months", "6 months", "1 year", "3 years", "5 years"))) %>%
dplyr::select(-tr) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 6b. Restricted Mean Survival Time by Partitioned Survival Analyses, Four-states model, 2nd patient"),
col.names = c("Time",rep(c("Estimate [95% CI]"),2)),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Transition 1: Admission to Ther.Discharge", 1, 8) %>%
kableExtra::pack_rows("Transition 2: Admission to Discharge w/o Clinical Advice", 9, 16) %>%
kableExtra::pack_rows("Transition 3: Admission to Readmission", 17, 24) %>%
kableExtra::pack_rows("Transition 4: Ther.Discharge to Readmission", 25, 32) %>%
kableExtra::pack_rows("Transition 5: Discharge w/o Clinical Advice to Readmission", 33, 40) %>%
kableExtra::scroll_box(width = "100%", height = "375px")| Time | Estimate [95% CI] | Estimate [95% CI] |
|---|---|---|
| Transition 1: Admission to Ther.Discharge | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.23 [0.23-0.24] | 0.24 [0.24-0.24] |
| 4 months | 0.31 [0.30-0.31] | 0.31 [0.31-0.32] |
| 6 months | 0.44 [0.43-0.45] | 0.45 [0.45-0.46] |
| 1 year | 0.79 [0.76-0.82] | 0.83 [0.80-0.86] |
| 3 years | 1.66 [1.51-1.81] | 1.88 [1.73-2.02] |
| 5 years | 2.15 [1.89-2.41] | 2.54 [2.27-2.81] |
| Transition 2: Admission to Discharge w/o Clinical Advice | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.24 [0.23-0.24] | 0.24 [0.23-0.24] |
| 4 months | 0.31 [0.31-0.31] | 0.31 [0.31-0.31] |
| 6 months | 0.45 [0.45-0.46] | 0.45 [0.44-0.46] |
| 1 year | 0.84 [0.82-0.85] | 0.83 [0.81-0.85] |
| 3 years | 2.13 [2.04-2.22] | 2.10 [2.01-2.20] |
| 5 years | 3.32 [3.15-3.48] | 3.26 [3.08-3.43] |
| Transition 3: Admission to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.25 [0.25-0.25] | 0.25 [0.25-0.25] |
| 4 months | 0.33 [0.33-0.33] | 0.33 [0.33-0.33] |
| 6 months | 0.50 [0.50-0.50] | 0.50 [0.50-0.50] |
| 1 year | 0.99 [0.99-1.00] | 0.99 [0.99-1.00] |
| 3 years | 2.87 [2.79-2.92] | 2.88 [2.81-2.93] |
| 5 years | 4.56 [4.35-4.72] | 4.60 [4.38-4.75] |
| Transition 4: Ther.Discharge to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.23 [0.22-0.24] | 0.23 [0.22-0.24] |
| 4 months | 0.30 [0.28-0.31] | 0.30 [0.28-0.31] |
| 6 months | 0.44 [0.41-0.46] | 0.44 [0.41-0.46] |
| 1 year | 0.82 [0.75-0.88] | 0.83 [0.76-0.89] |
| 3 years | 2.15 [1.87-2.39] | 2.18 [1.90-2.43] |
| 5 years | 3.28 [2.77-3.72] | 3.34 [2.81-3.80] |
| Transition 5: Discharge w/o Clinical Advice to Readmission | ||
| 0 days | 0.00 [0.00-0.00] | 0.00 [0.00-0.00] |
| 30 days | 0.08 [0.08-0.08] | 0.08 [0.08-0.08] |
| 90 days | 0.24 [0.24-0.24] | 0.24 [0.24-0.25] |
| 4 months | 0.32 [0.32-0.32] | 0.32 [0.32-0.33] |
| 6 months | 0.48 [0.47-0.48] | 0.48 [0.47-0.49] |
| 1 year | 0.92 [0.90-0.94] | 0.93 [0.91-0.95] |
| 3 years | 2.58 [2.46-2.68] | 2.62 [2.51-2.71] |
| 5 years | 4.10 [3.87-4.29] | 4.19 [3.97-4.38] |
#kableExtra::add_footnote("Note. NA= Null values", notation="none") MSM
tolos2_t_a_str2_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_3_states)[2],get(paste0("tolos2_t_a_str2_",t)))
}))
tolos2_t_b_str2_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_3_states)[2],get(paste0("tolos2_t_b_str2_",t)))
}))
tolos2_4s_t_a_str2_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_a_str2_",t)))
}))
tolos2_4s_t_b_str2_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_b_str2_",t)))
}))
tolos2_4s_t_a_str3_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_a_str3_",t)))
}))
tolos2_4s_t_b_str3_df<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t=t,state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_b_str3_",t)))
}))
tolos2_t_a_df2<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_3_states)[2],get(paste0("tolos2_t_a_",t)))
}))
tolos2_t_b_df2<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_3_states)[2],get(paste0("tolos2_t_b_",t)))
}))
tolos2_4s_t_a_df2<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_a_",t)))
}))
tolos2_4s_t_b_df2<-
do.call('rbind', lapply(c(.25,1,3), function(t) {
cbind.data.frame(t= t, state=1:dim(mat_4_states)[2],get(paste0("tolos2_4s_t_b_",t)))
}))rbind.data.frame(
cbind.data.frame(tolos2_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos2_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos2_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos2_t_b_str2_df,program="General Population", start=2))%>%
dplyr::mutate(comb=dplyr::case_when(start==2 & state==1~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(state<max(mat_3_states,na.rm=T)) %>%
dplyr::mutate(start=dplyr::case_when(start==1~"1)Admission", start==2~"2)Therapeutic\nDischarge",start==3~"3)Readmission",T~NA_character_))%>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic\nDischarge",
state==3~"3)Readmission"))%>%
dplyr::filter(is.na(comb)) %>% # janitor::tabyl(start, state)
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",L),"-",
sprintf("%1.2f",U),"]")) %>%
dplyr::select(program, t, start, state, est_ci) %>%
tidyr::pivot_wider(names_from="program", values_from=c("est_ci")) %>%
dplyr::mutate(comb=paste0(start, "_", state)) %>% #janitor::tabyl(start, state)
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500", "1.0000", "3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::arrange(start, state, t) %>%
dplyr::select(-start,-comb) %>%
dplyr::select(state, t, `Women-specific`, `General Population`) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 7a. Length of Stay, Three-states model; patient 2"),
col.names = c("State","Time","Women only","Mixed gender"),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
kableExtra::pack_rows("Starting From Admission", 1, 6) %>%
kableExtra::pack_rows("Starting From Therapeutic Discharge", 7, 9) %>%
kableExtra::add_footnote("Note. Readmission state will not be considered because is an absorbing state", notation="none") %>%
kableExtra::scroll_box(width = "100%", height = "375px")| State | Time | Women only | Mixed gender |
|---|---|---|---|
| Starting From Admission | |||
| 1)Admission | ~91 days/3 months | 0.23 [0.23-0.24] | 0.24 [0.23-0.24] |
| 1)Admission | ~1 year | 0.77 [0.74-0.80] | 0.81 [0.78-0.84] |
| 1)Admission | ~3 years | 1.66 [1.53-1.79] | 1.85 [1.72-1.98] |
| 2)Therapeutic Discharge | ~91 days/3 months | 0.02 [0.01-0.02] | 0.01 [0.01-0.02] |
| 2)Therapeutic Discharge | ~1 year | 0.18 [0.15-0.21] | 0.15 [0.12-0.17] |
| 2)Therapeutic Discharge | ~3 years | 0.86 [0.72-1.00] | 0.73 [0.61-0.87] |
| Starting From Therapeutic Discharge | |||
| 2)Therapeutic Discharge | ~91 days/3 months | 0.23 [0.22-0.24] | 0.23 [0.22-0.24] |
| 2)Therapeutic Discharge | ~1 year | 0.82 [0.75-0.88] | 0.83 [0.76-0.89] |
| 2)Therapeutic Discharge | ~3 years | 2.15 [1.86-2.39] | 2.18 [1.90-2.43] |
| Note. Readmission state will not be considered because is an absorbing state | |||
rbind.data.frame(
cbind.data.frame(tolos2_4s_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos2_4s_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos2_4s_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos2_4s_t_b_str2_df,program="General Population", start=2),
cbind.data.frame(tolos2_4s_t_a_str3_df,program="Women-specific", start=3),
cbind.data.frame(tolos2_4s_t_b_str3_df,program="General Population", start=3))%>%
dplyr::mutate(comb=dplyr::case_when(start==2 & state==1~"a",
start==2 & state==3~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(state<max(dim(mat_4_states)[2],na.rm=T)) %>%
dplyr::mutate(start=dplyr::case_when(start==1~"1)Admission", start==2~"2)Therapeutic Discharge",start==3~"3)Discharge Without Clinical Advice",start==4~"4)Readmission",T~NA_character_))%>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Discharge Without Clinical Advice",
T~"4)Readmission"))%>%
dplyr::filter(is.na(comb)) %>% #janitor::tabyl(start, state)
dplyr::mutate(est_ci=paste0(sprintf("%1.2f",est)," [",
sprintf("%1.2f",L),"-",
sprintf("%1.2f",U),"]")) %>%
dplyr::select(program,t,start,state,est_ci) %>%
tidyr::pivot_wider(names_from="program",values_from=c("est_ci")) %>%
dplyr::mutate(comb=paste0(start,"_",state)) %>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500", "1.0000", "3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
dplyr::arrange(start, state, t) %>%
dplyr::select(-start,-comb) %>%
dplyr::select(state, t, `Women-specific`, `General Population`) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 7b. Length of Stay, Four-states model; patient 2"),
col.names = c("State","Time","Women only","Mixed gender"),
align =rep('c', 101)) %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
#kableExtra::add_header_above(c(rep("",1),"Women-specific"=1, "General population"=1)) %>%
kableExtra::pack_rows("Starting From Admission", 1, 9) %>%
kableExtra::pack_rows("Starting From Therapeutic Discharge", 10, 12) %>%
kableExtra::pack_rows("Starting From Discharge Without Clinical Advice", 13, 15) %>%
kableExtra::add_footnote("Note. Readmission state will not be considered because is an absorbing state", notation="none") %>%
kableExtra::scroll_box(width = "100%", height = "375px")| State | Time | Women only | Mixed gender |
|---|---|---|---|
| Starting From Admission | |||
| 1)Admission | ~91 days/3 months | 0.22 [0.22-0.22] | 0.22 [0.22-0.23] |
| 1)Admission | ~1 year | 0.67 [0.64-0.70] | 0.69 [0.67-0.72] |
| 1)Admission | ~3 years | 1.21 [1.09-1.32] | 1.34 [1.23-1.44] |
| 2)Therapeutic Discharge | ~91 days/3 months | 0.01 [0.01-0.02] | 0.01 [0.01-0.01] |
| 2)Therapeutic Discharge | ~1 year | 0.16 [0.13-0.19] | 0.13 [0.11-0.15] |
| 2)Therapeutic Discharge | ~3 years | 0.79 [0.67-0.93] | 0.67 [0.56-0.79] |
| 3)Discharge Without Clinical Advice | ~91 days/3 months | 0.01 [0.01-0.01] | 0.01 [0.01-0.02] |
| 3)Discharge Without Clinical Advice | ~1 year | 0.13 [0.12-0.15] | 0.14 [0.13-0.16] |
| 3)Discharge Without Clinical Advice | ~3 years | 0.59 [0.52-0.66] | 0.65 [0.57-0.73] |
| Starting From Therapeutic Discharge | |||
| 2)Therapeutic Discharge | ~91 days/3 months | 0.23 [0.22-0.24] | 0.23 [0.22-0.24] |
| 2)Therapeutic Discharge | ~1 year | 0.82 [0.75-0.88] | 0.83 [0.75-0.89] |
| 2)Therapeutic Discharge | ~3 years | 2.15 [1.87-2.40] | 2.18 [1.89-2.43] |
| Starting From Discharge Without Clinical Advice | |||
| 3)Discharge Without Clinical Advice | ~91 days/3 months | 0.24 [0.24-0.24] | 0.24 [0.24-0.25] |
| 3)Discharge Without Clinical Advice | ~1 year | 0.92 [0.90-0.95] | 0.93 [0.91-0.95] |
| 3)Discharge Without Clinical Advice | ~3 years | 2.58 [2.46-2.68] | 2.62 [2.50-2.72] |
| Note. Readmission state will not be considered because is an absorbing state | |||
#tolos_t_a_df2 tolos_t_b_df2 tolos_4s_t_a_df2 %>% dplyr::filter(t==10) tolos_4s_t_b_df2 tolos_4s_t_a_10
los2_3s_bar<-
rbind.data.frame(
cbind.data.frame(tolos2_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos2_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos2_t_a_str2_df,program="Women-specific", start=2)%>% dplyr::filter(state==2),
cbind.data.frame(tolos2_t_b_str2_df,program="General Population", start=2)%>% dplyr::filter(state==2)
) %>%
dplyr::filter(state<3) %>%
filter(case_when(start==1 ~ state==1, #When y == "", x > 3
T ~ state>0) #Otherwise, x < 3
) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Readmission"))%>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500","1.0000","3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=program, fill=program), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(.~start+state, ncol=3) +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
#scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position="bottom")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
los2_3s_bar
Figure 10a. Estimated Length of Stay (w/ confidence intervals), 2nd patient. Three-states model
if(no_mostrar==1){
jpeg("eso17.jpg", height=10, width= 12, res= 96, units = "in")
los2_3s_bar
dev.off()
}los2_4s_bar<-
rbind.data.frame(
cbind.data.frame(tolos2_4s_t_a_df2,program="Women-specific",start=1),
cbind.data.frame(tolos2_4s_t_b_df2,program="General Population",start=1),
cbind.data.frame(tolos2_4s_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos2_4s_t_b_str2_df,program="General Population", start=2),
cbind.data.frame(tolos2_4s_t_a_str3_df,program="Women-specific", start=3),
cbind.data.frame(tolos2_4s_t_b_str3_df,program="General Population", start=3))%>%
dplyr::filter(state<4) %>%
dplyr::mutate(comb=dplyr::case_when(
start==1 & state==2~"a",
start==1 & state==3~"a",
start==2 & state==1~"a",
start==2 & state==3~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(is.na(comb)) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Discharge Without Clinical Advice",
T~"to\n4)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Discharge Without Clinical Advice",
T~"Starting from Readmission"))%>%
dplyr::arrange(start, state) %>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500","1.0000","3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=program, fill=program), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(.~start+state, ncol=3) +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
#scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position="bottom")+#c(.9,.1))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
los2_4s_bar
Figure 10b. Estimated Length of Stay (w/ confidence intervals), 2nd patient. Four-states model
if(no_mostrar==1){
jpeg("eso17.jpg", height=10, width= 12, res= 96, units = "in")
los2_4s_bar
dev.off()
}#tolos_t_a_df2 tolos_t_b_df2 tolos_4s_t_a_df2 %>% dplyr::filter(t==10) tolos_4s_t_b_df2 tolos_4s_t_a_10
los_3s_bar<-
rbind.data.frame(
cbind.data.frame(tolos_t_a_df2,program="Women-specific", start=1),
cbind.data.frame(tolos_t_b_df2,program="General Population", start=1),
cbind.data.frame(tolos_t_a_str2_df,program="Women-specific", start=2)%>% dplyr::filter(state==2),
cbind.data.frame(tolos_t_b_str2_df,program="General Population", start=2)%>% dplyr::filter(state==2)
) %>%
dplyr::filter(state<3) %>%
filter(case_when(start==1 ~ state==1, #When y == "", x > 3
T ~ state>0) #Otherwise, x < 3
) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Readmission"))%>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500","1.0000","3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=program, fill=program), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(.~start+state, ncol=3) +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
#scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position="bottom")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
los_3s_bar
Figure 11a. Estimated Length of Stay (w/ confidence intervals), 1st patient. Three-states model
if(no_mostrar==1){
jpeg("eso17.jpg", height=10, width= 12, res= 96, units = "in")
los_3s_bar
dev.off()
}los_4s_bar<-
rbind.data.frame(
cbind.data.frame(tolos_4s_t_a_df2,program="Women-specific",start=1),
cbind.data.frame(tolos_4s_t_b_df2,program="General Population",start=1),
cbind.data.frame(tolos_4s_t_a_str2_df,program="Women-specific", start=2),
cbind.data.frame(tolos_4s_t_b_str2_df,program="General Population", start=2),
cbind.data.frame(tolos_4s_t_a_str3_df,program="Women-specific", start=3),
cbind.data.frame(tolos_4s_t_b_str3_df,program="General Population", start=3))%>%
dplyr::filter(state<4) %>%
dplyr::mutate(comb=dplyr::case_when(
start==1 & state==2~"a",
start==1 & state==3~"a",
start==2 & state==1~"a",
start==2 & state==3~"a",
start==3 & state==1~"a",
start==3 & state==2~"a",
T~NA_character_)) %>%
dplyr::filter(is.na(comb)) %>%
dplyr::mutate(state=dplyr::case_when(state==1~"1)Admission",
state==2~"2)Therapeutic Discharge",
state==3~"3)Discharge Without Clinical Advice",
T~"to\n4)Readmission"))%>%
dplyr::mutate(start=dplyr::case_when(start==1~"Starting from Admission",
start==2~"Starting from Therapeutic Discharge",
start==3~"Starting from Discharge Without Clinical Advice",
T~"Starting from Readmission"))%>%
dplyr::arrange(start, state) %>%
dplyr::mutate(t=sprintf("%1.4f",t)) %>%
dplyr::filter(t %in% c("0.2500","1.0000","3.0000")) %>%
dplyr::mutate(t=factor(t,
levels=c("0.2500", "1.0000", "3.0000"),
labels=c("~91 days/3 months", "~1 year", "~3 years"))) %>%
ggplot()+
geom_bar(aes(x=t, y=est, color=program, fill=program), position="dodge", stat="identity", alpha=0.4)+
sjPlot::theme_sjplot2()+
geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4, alpha=0.9, size=1.3)+
facet_wrap(.~start+state, ncol=3) +
scale_fill_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
scale_color_manual(name ="Type of program",values=c("#0D56A6","#FF9A00"), labels= c("Mixed gender","Women only"))+
#scale_x_continuous(breaks=seq(0, 11))+
#scale_y_continuous(labels = scales::percent_format(accuracy = 2), limits=c(0,1))+
xlab("Years") +
ylab("")+
#ylab("State occupancy probabilities") +
theme_minimal()+
theme(legend.position="bottom")+#c(.8,.85))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
panel.border = element_blank(),
legend.box.background = element_blank())
los_4s_bar
Figure 11b. Estimated Length of Stay (w/ confidence intervals), 1st patient. Four-states model
if(no_mostrar==1){
jpeg("eso17.jpg", height=10, width= 12, res= 96, units = "in")
los_4s_bar
dev.off()
}cumhaz_plot12<-
ggplot()+
geom_step(data=data.frame(flexsurv_cumhaz_a2$Haz) %>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=data.frame(flexsurv_cumhaz_b2$Haz)%>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="solid"),size=1, alpha=.65) +
scale_color_manual(name ="Transition",values=c("#0E53A7","#200772","#BF8830"),labels= c("Admission -> Ther. Discharge","Admission -> Readmission","Therapeutic Discharge -> Readmission"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
scale_x_continuous(breaks=seq(0, 11, by = 1))+
#ylim(0,2)+
ylab("Cumulative hazards")+
theme_minimal()+
theme(axis.title.x = element_blank(),
legend.position=c(.65,.6),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
cumhaz_plot22<-
ggplot()+
geom_step(data=data.frame(flexsurv_cumhaz_c2$Haz) %>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="dashed"),size=1, alpha=.65) +
geom_step(data=data.frame(flexsurv_cumhaz_d2$Haz)%>% dplyr::mutate(trans=factor(trans)), aes(time, Haz, color= trans,linetype="solid"),size=1, alpha=.65) +
scale_color_manual(name ="Transitions",values=c("gray40","#1B0773","#6E0069","#025167","#A68600"),labels= c("Admission-> Therapeutic Discharge","Admission-> Discharge Without Clinical Advice","Admission->Readmission","Therapeutic Discharge->Readmission","Discharge Without Clinical Advice->Readmission"))+
scale_linetype_manual(name= "Type of Plan",values=c("dashed","solid"), labels=c("Women only","Mixed gender"))+
#values=c("#BF8F30","#1D7373","#876ED7","#1240AB")
scale_x_continuous(breaks=seq(0, 11, by = 1))+
theme_minimal()+
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position=c(.38,.85),
legend.background = element_rect(fill="white",
size=0.5, linetype="solid",
colour ="black"))
#ylim(0,2)
ggarrange(cumhaz_plot12,cumhaz_plot22)
Figure 12. Estimate of Cumulative Hazards of the Parametric Model
if(no_mostrar==1){
jpeg("eso14.jpg", height=15, width= 10, res= 320, units = "in")
ggarrange(cumhaz_plot1,cumhaz_plot2)
dev.off()
}Session Info
Sys.getenv("R_LIBS_USER")## [1] "C:/Users/CISS Fondecyt/OneDrive/Documentos/R/win-library/4.0"rstudioapi::getSourceEditorContext()## Document Context:
## - id: 'BC4A16B9'
## - path: 'G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Proyecto_carla33.Rmd'
## - contents: <3538 rows>
## Document Selection:
## - [2021, 2] -- [2021, 2]: ''#save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla.RData")
if (grepl("CISS Fondecyt",rstudioapi::getSourceEditorContext()$path)==T){
save.image("C:/Users/CISS Fondecyt/OneDrive/Escritorio/SUD_CL/mult_state_carla_3.RData")
} else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_carla_3.RData")
} else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla_3.RData")
} else {
save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla_3.RData")
}
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Chile.1252 LC_CTYPE=Spanish_Chile.1252
## [3] LC_MONETARY=Spanish_Chile.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Chile.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] radiant.update_1.4.1 Metrics_0.1.4 muhaz_1.2.6.1
## [4] flexsurv_2.0 mstate_0.3.1 DiagrammeR_1.0.6.1.9000
## [7] Amelia_1.7.6 Rcpp_1.0.5 igraph_1.2.6
## [10] eha_2.8.1 cobalt_4.2.3 MatchIt_3.0.2
## [13] tableone_0.12.0 stargazer_5.2.2 reshape2_1.4.4
## [16] gridExtra_2.3 foreign_0.8-80 survMisc_0.5.5
## [19] ggfortify_0.4.10 survminer_0.4.8 ggpubr_0.4.0
## [22] epiR_1.0-15 forcats_0.5.0 purrr_0.3.4
## [25] readr_1.3.1 tibble_3.0.3 tidyverse_1.3.0
## [28] dplyr_1.0.1 treemapify_2.5.3 sf_0.9-3
## [31] ggiraph_0.7.0 finalfit_1.0.1 lsmeans_2.30-0
## [34] emmeans_1.4.8 RColorBrewer_1.1-2 panelr_0.7.3
## [37] lme4_1.1-23 Matrix_1.2-18 data.table_1.13.0
## [40] codebook_0.9.2 devtools_2.3.0 usethis_1.6.1
## [43] sqldf_0.4-11 RSQLite_2.2.0 gsubfn_0.7
## [46] proto_1.0.0 broom_0.7.0 zoo_1.8-8
## [49] rbokeh_0.5.1 janitor_2.0.1 plotly_4.9.2.1
## [52] kableExtra_1.1.0 Hmisc_4.4-0 Formula_1.2-3
## [55] survival_3.1-12 lattice_0.20-41 ggplot2_3.3.2
## [58] stringr_1.4.0 stringi_1.4.6 tidyr_1.1.1
## [61] knitr_1.29 matrixStats_0.56.0 boot_1.3-25
##
## loaded via a namespace (and not attached):
## [1] estimability_1.3 coda_0.19-3 acepack_1.4.1
## [4] bit64_0.9-7 multcomp_1.4-13 rpart_4.1-15
## [7] generics_0.0.2 callr_3.4.3 cowplot_1.0.0
## [10] TH.data_1.0-10 mice_3.11.0 ggfittext_0.9.0
## [13] chron_2.3-55 bit_1.1-15.2 webshot_0.5.2
## [16] xml2_1.3.2 lubridate_1.7.9 assertthat_0.2.1
## [19] xfun_0.16 hms_0.5.3 evaluate_0.14
## [22] fansi_0.4.1 dbplyr_1.4.4 readxl_1.3.1
## [25] km.ci_0.5-2 DBI_1.1.0 htmlwidgets_1.5.1
## [28] ellipsis_0.3.1 backports_1.1.7 insight_0.9.0
## [31] survey_4.0 vctrs_0.3.2 remotes_2.2.0
## [34] sjlabelled_1.1.6 abind_1.4-5 withr_2.2.0
## [37] pryr_0.1.4 checkmate_2.0.0 prettyunits_1.1.1
## [40] cluster_2.1.0 lazyeval_0.2.2 crayon_1.3.4
## [43] pkgconfig_2.0.3 labeling_0.3 units_0.6-6
## [46] nlme_3.1-148 pkgload_1.1.0 nnet_7.3-14
## [49] rlang_0.4.7 lifecycle_0.2.0 sandwich_2.5-1
## [52] modelr_0.1.8 cellranger_1.1.0 tcltk_4.0.2
## [55] rprojroot_1.3-2 KMsurv_0.1-5 carData_3.0-4
## [58] reprex_0.3.0 base64enc_0.1-3 processx_3.4.3
## [61] png_0.1-7 viridisLite_0.3.0 parameters_0.8.2
## [64] KernSmooth_2.23-17 visNetwork_2.0.9 pander_0.6.3
## [67] blob_1.2.1 classInt_0.4-3 jpeg_0.1-8.1
## [70] rstatix_0.6.0 ggeffects_0.15.1 ggsignif_0.6.0
## [73] scales_1.1.1 memoise_1.1.0 magrittr_1.5
## [76] plyr_1.8.6 hexbin_1.28.1 compiler_4.0.2
## [79] snakecase_0.11.0 cli_2.0.2 ps_1.3.3
## [82] htmlTable_2.0.1 MASS_7.3-51.6 tidyselect_1.1.0
## [85] highr_0.8 mitools_2.4 jtools_2.0.5
## [88] yaml_2.2.1 latticeExtra_0.6-29 tools_4.0.2
## [91] parallel_4.0.2 rio_0.5.16 rstudioapi_0.11
## [94] uuid_0.1-4 gistr_0.5.0 farver_2.0.3
## [97] sjPlot_2.8.4 digest_0.6.25 quadprog_1.5-8
## [100] car_3.0-8 performance_0.4.8 httr_1.4.2
## [103] gdtools_0.2.2 effectsize_0.3.2 sjstats_0.18.0
## [106] colorspace_1.4-1 rvest_0.3.6 fs_1.5.0
## [109] splines_4.0.2 statmod_1.4.34 sessioninfo_1.1.1
## [112] systemfonts_0.2.3 xtable_1.8-4 jsonlite_1.7.0
## [115] nloptr_1.2.2.2 testthat_2.3.2 R6_2.4.1
## [118] pillar_1.4.6 htmltools_0.5.0 glue_1.4.1
## [121] minqa_1.2.4 deSolve_1.28 class_7.3-17
## [124] codetools_0.2-16 maps_3.3.0 pkgbuild_1.1.0
## [127] mvtnorm_1.1-1 numDeriv_2016.8-1.1 curl_4.3
## [130] BiasedUrn_1.07 zip_2.1.1 openxlsx_4.1.5
## [133] rmarkdown_2.6 desc_1.2.0 munsell_0.5.0
## [136] e1071_1.7-3 labelled_2.5.0 sjmisc_2.8.5
## [139] haven_2.3.1 gtable_0.3.0 bayestestR_0.7.2