Readmission risk among women in women-only versus mixed gender treatment programs for substance use disorder in Chile (4th step)
07 julio, 2021
Transition intensities
`%>%`<-magrittr::`%>%`
transition_label_4s<- c(`1`="Transition 1: Admission to Treatment completion",`2`="Transition 2: Admission to Discharge without completion",`3`="Transition 3: Admission to Readmission",`4`="Transition 4: Treatment completion to Readmission", `5`="Transition 5: Discharge without completion to Readmission")
layout(matrix(1:6, nc = 2, byrow = F))
for (i in 1:n_trans2){
plot(km.lc2[[i]], col="red", conf.int=F, xlim=c(0,12));
lines(fits_wei2[[i]], col="coral4", ci=F);
lines(fits_weiph2[[i]], col="navyblue", ci=F);
lines(fits_gomp2[[i]], col="lightpink", ci=F);
lines(fits_llogis2[[i]], col="gray25", ci=F);##A0A36D
lines(fits_gam2[[i]], col="darkorchid4", ci=F);##886894
lines(fits_ggam2[[i]], col="#496A72", ci=F);
lines(fits_logn2[[i]], col="gray70", ci=F);
lines(fits_exp2[[i]],col="#A0A36D", ci=F);
lines(fits_genf2[[i]],col="cadetblue", ci=F)
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Lognormal", "Exponential", "Generalized F"), col =
c("red","coral4","navyblue","lightpink","gray25",#"#B89673",##886894
"darkorchid4","#496A72","gray70","#A0A36D","cadetblue"),
title = "Distributions", cex = .95, bty = "n", lty=1.2,lwd=3)# lty = 1:2,
title(main=transition_label_4s[[i]])
}
path_jpg_1<-rstudioapi::getSourceEditorContext()$path
if(no_mostrar==1){
jpeg(paste0(sub("Proyecto_carla34.Rmd","",path_jpg_1),"pred_surv_4st_int_only.jpg"), height=18, width= 10, res= 320, units = "in")
layout(matrix(1:6, nc = 2, byrow = F))
for (i in 1:n_trans2){
plot(km.lc2[[i]], col="red", conf.int=F, xlim=c(0,12));
lines(fits_wei2[[i]], col="coral4", ci=F);
lines(fits_weiph2[[i]], col="navyblue", ci=F);
lines(fits_gomp2[[i]], col="lightpink", ci=F);
lines(fits_llogis2[[i]], col="gray25", ci=F);##A0A36D
lines(fits_gam2[[i]], col="darkorchid4", ci=F);##886894
lines(fits_ggam2[[i]], col="#496A72", ci=F);
lines(fits_logn2[[i]], col="gray70", ci=F);
lines(fits_exp2[[i]],col="#A0A36D", ci=F);
lines(fits_genf2[[i]],col="cadetblue", ci=F)
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Lognormal", "Exponential", "Generalized F"), col =
c("red","coral4","navyblue","lightpink","gray25",#"#B89673",##886894
"darkorchid4","#496A72","gray70","#A0A36D","cadetblue"),
title = "Distributions", cex = .95, bty = "n", lty=1.2,lwd=3)# lty = 1:2,
title(main=transition_label_4s[[i]])
}
dev.off()
}
Figure 1. Vissual Assessment of Hazards, Four-states Model (w/o covars)
haz_plot_4s_int_only<-
haz_4s %>%
dplyr::mutate(dist=factor(dist,levels=c("Kernel density",dists_wo_covs_4s$formal))) %>%
ggplot()+
geom_line(aes(time, est, color=dist),size=1)+
#geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=.5, linetype=0) +
scale_color_manual(name="Distributions", values = c("black","#f54b96","#00e9b1","#69b763",
"#166000","#b27ff9","#fa863b","#013eab","#a7aa48","#b34b40","darkred"))+
#c("black",brewer.pal(n = 9, name = 'Paired')))+
#c("#112A60","#085754","#D3A347","#4F3C91","red","#112A60","#085754","#8F630D","#251363")) +
facet_wrap(.~trans,labeller = labeller(trans = transition_label_4s), ncol=2, scales = "free_y")+
sjPlot::theme_sjplot2()+
theme(legend.position="bottom",
strip.background = element_rect(fill = "white", colour = "white"))+
scale_x_continuous(breaks = seq(1, 15, 2))+
#ylim(0,1.25)+
#theme(axis.text.x = element_blank(),
# panel.grid.major = element_blank(),
# panel.grid.minor = element_blank()) +
labs(y="Hazard",x="Time (years)")
haz_plot_int<-list()
for (i in 1:(length(transition_label_4s))){
haz_plot_int[[i]]<-
haz_4s %>%
dplyr::mutate(dist=factor(dist,levels=c("Kernel density",dists_wo_covs_4s$formal))) %>%
dplyr::filter(trans==i) %>%
ggplot()+
geom_line(aes(time, est, color=dist),size=1)+
scale_color_manual(name="Distributions", values = c("red", "coral4","navyblue","lightpink","gray25","darkorchid4","#496A72","gray70","#A0A36D","cadetblue"))+
sjPlot::theme_sjplot2()+
theme(legend.position="bottom",
strip.background = element_rect(fill = "white", colour = "white"))+
scale_x_continuous(breaks = seq(1, 15, 2))+
labs(y="",x="")+ theme(legend.position="none")+theme(plot.title = element_text(size = 11, face = "bold"))+
ggtitle(transition_label_4s[[i]])
}
legend_int_only <-
haz_4s %>%
dplyr::mutate(dist=factor(dist,levels=c("Kernel density",dists_wo_covs_4s$formal))) %>%
dplyr::filter(trans==i) %>%
ggplot()+
geom_line(aes(time, est, color=dist),size=1)+
scale_color_manual(name="Distributions", values = c("red", "coral4","navyblue","lightpink","gray25","darkorchid4","#496A72","gray70","#A0A36D","cadetblue"))+
sjPlot::theme_sjplot2()+
theme(legend.position="bottom",
strip.background = element_rect(fill = "white", colour = "white"))+
scale_x_continuous(breaks = seq(1, 15, 2))+
lims(x = c(0,0), y = c(0,0))+
theme_void()+
theme(legend.position = c(0.5,0.5),
legend.key.size = unit(.2, "cm"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"))+
guides(colour = guide_legend(ncol = 2,override.aes = list(size=7)))## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.haz_plot_int[[6]]<-legend_int_only
fig_haz_plot_int<-
ggarrange(
plotlist=haz_plot_int, ncol=2, nrow=3)## Warning: Removed 1721 row(s) containing missing values (geom_path).fig_haz_plot_int_lab<-
annotate_figure(fig_haz_plot_int, left = textGrob("Cumulative hazards", rot = 90, vjust = 1, gp = gpar(cex = 1.3)),
bottom = textGrob("Years", gp = gpar(cex = 1.3)))
fig_haz_plot_int_lab
if(no_mostrar==1){
jpeg(paste0(sub("Proyecto_carla34.Rmd","",path_jpg_1),"pred_haz_4st_int_only.jpg"), height=15, width= 10, res= 320, units = "in")
fig_haz_plot_int_lab
#PPredicted survival curves for the five transitions of the four-states multistate model
#grid.arrange(haz_plot_4s_int_only,fit_plot_int_4s_only, heights=c(2,1.5))
dev.off()
}`%>%`<-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 -> Treatment completion",
"Admission -> Readmission",
"Treatment completion -> 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 -> Treatment completion",
"Admission -> Readmission",
"Treatment completion -> 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 -> Treatment completion",
"Admission -> Discharge without completion",
"Admission -> Readmission",
"Treatment completion -> Readmission",
"Discharge without completion -> 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 -> Treatment completion",
"Admission -> Discharge without completion",
"Admission -> Readmission",
"Treatment completion -> Readmission",
"Discharge without completion -> 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"))`%>%`<-magrittr::`%>%`
trans_for_three_st<-c("Admission -> Treatment completion","Admission -> Readmission","Treatment completion -> Readmission")
trans_for_four_st<-c("Admission -> Treatment completion","Admission -> Discharge without completion","Admission -> Readmission","Treatment completion -> Readmission","Discharge without completion -> Readmission")
abc_let<-c("A","B","C","D","E")
st_plot<-list()
for (i in 1:(length(trans_for_four_st))){
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]]), 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]]), 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]]))
}
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)))
st_plot2<-list()
for (i in 1:(length(trans_for_three_st))){
st_plot2[[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]]))
}
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[[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")),
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_plot2[[1]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot2[[2]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+theme(plot.title = element_text(size = 10, face = "bold")),
st_plot2[[3]] + theme(legend.position="none")+xlim(0,5)+ylim(0,1)+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 2. 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/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/")
} 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_paper_carla.pdf"), 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")
ggsave(paste0(jpg_path,"Fig2_alt_paper_carla.pdf"), fig_cox_abcd_three_states_lab, width = 10, height = 11.5, dpi = 600, units= "in")
}plots2<- data.frame(title=rep(c("(a) Admission -> Treatment completion","(b) Admission -> Treatment w/o completion","(c) Admission -> Readmission", "(d) Treatment completion -> Readmission","(e) Treatment w/o completion -> Readmission"),1),trans=1:5)
layout(matrix(1:10, nc = 2, byrow = F))
for(i in c(1:5)){
plot(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv)), type="l",
xlab="log(Years)", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv)), lty=2)
legend(0,-4, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("LOG CUMULATIVE HAZARD VS LOG TIME PLOT \n",plots2[i,"title"]), cex.main=1.5)
}
for(i in c(1:5)){
plot(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv), type="l",
xlab="Years", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv), lty=2)
legend(7.5,.23, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("CUMULATIVE HAZARD PLOT: -LOG(KM SURVIVAL) \n",plots2[i,"title"]), cex.main=1.5)
}
Figure 3. Vissual Assessment of Hazards, Four-states Model (w/o covars)
if(no_mostrar==1){
jpeg(paste0(sub("Proyecto_carla34.Rmd","",path_jpg),"ph1.jpg"), height=15, width= 10, res= 320, units = "in")
layout(matrix(1:10, nc = 2, byrow = F))
for(i in c(1:5)){
plot(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv), type="l",
xlab="Years", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv), lty=2)
legend(7.5,.31, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("CUMULATIVE HAZARD PLOT: -LOG(KM SURVIVAL) \n",plots2[i,"title"]), cex.main=1.5)
}
for(i in c(1:5)){
plot(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv)), type="l",
xlab="log(Years)", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv)), lty=2)
legend(0,-3, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("LOG CUMULATIVE HAZARD VS LOG TIME PLOT \n",plots2[i,"title"]), cex.main=1.5)
}
dev.off()
}coxzph<-data.frame()
for(i in c(1:5)){
coxzph<-
rbind(coxzph,cox.zph(coxph(Surv(time,status)~factor(tipo_de_programa_2), data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==i)), transform="km", terms=TRUE, singledf=FALSE, global=TRUE)$table)
}
coxzph$p<-round(coxzph$p, 3)
data.table::data.table(coxzph,keep.rownames = T) %>%
dplyr::filter(grepl("GLOBAL",rn)) %>%
dplyr::mutate(chisq=round(chisq,3)) %>%
dplyr::mutate(transition=1:5) %>%
dplyr::select(transition, chisq, df, p) %>%
knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
caption = paste0("Table 1. Proportionality (Arrival, ~3 years + RP)"),
col.names = c("State","Time","Residential","Outpatient"),
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 1)Admission", 1, 5) %>%
kableExtra::pack_rows("Starting From 2)Readmission", 6, 10) %>%
kableExtra::pack_rows("Starting From 3)2nd Readmission", 11, 15) %>%
kableExtra::pack_rows("Starting From 4)3rd Readmission", 16, 20) %>%
# kableExtra::add_footnote("Note. Readmission state will not be considered because is an absorbing state", notation="none") %>%
kableExtra::scroll_box(width = "100%", height = "375px")## Error in UseMethod("nodeset_apply"): no applicable method for 'nodeset_apply' applied to an object of class "NULL"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_medians2<-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)
#newdat5a_2 newdat5b_2
set.seed(1234)
simfinal_fmsm2_4st_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_2[1,],newdat5b_2[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_fmsm2_4st_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_2[1,],newdat5b_2[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_fmsm2_4st_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_2[1,],newdat5b_2[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_fmsm2_4st_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_2[1,],newdat5b_2[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_medians2<-Sys.time()
paste0("Time in process: ");time_after_sim_medians2-time_before_sim_medians2## [1] "Time in process: "## Time difference of 19.51208 mins# probability of each final outcome
# mean and quantiles of the time to that outcome for people who experience itsimfinal_fmsm2_4st_df<-
rbind.data.frame(
simfinal_fmsm2_4st_30d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.08),
simfinal_fmsm2_4st_90d%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=.25),
simfinal_fmsm2_4st_1y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=1),
simfinal_fmsm2_4st_3y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=3),
simfinal_fmsm2_4st_5y%>%dplyr::select("tipo_de_programa_2","quantity","val","lower","upper")%>%mutate(year=5))## Error in eval(lhs, parent, parent): objeto 'simfinal_fmsm2_4st_30d' no encontradosimfinal_fmsm2_4st_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")## Error in eval(lhs, parent, parent): objeto 'simfinal_fmsm2_4st_df' no encontradoif(no_mostrar==1){
jpeg(paste0(sub("Proyecto_carla34.Rmd","",path_jpg),"ph1.jpg"), height=15, width= 10, res= 320, units = "in")
layout(matrix(1:10, nc = 2, byrow = F))
for(i in c(1:5)){
plot(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv), type="l",
xlab="Years", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5, col=1)
lines(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time,
-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv), lty=2)
legend(7.5,.31, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("CUMULATIVE HAZARD PLOT: -LOG(KM SURVIVAL) \n",plots2[i,"title"]), cex.main=1.5)
}
for(i in c(1:5)){
plot(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==0))$surv)), type="l",
xlab="log(Years)", ylab="", xaxs="i",yaxs="i",
las=1,cex.lab=1.5, cex.axis=1.5)
lines(log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$time),
log(-log(survfit(Surv(time,status)~1, data=subset(ms2_CONS_C1_SEP_2020_women_imputed, trans==plots2[i,"trans"] & tipo_de_programa_2==1))$surv)), lty=2)
legend(0,-3, c("MG", "WO"), bty="n", lty=c(2,1), cex=1.5)
title(main=paste0("LOG CUMULATIVE HAZARD VS LOG TIME PLOT \n",plots2[i,"title"]), cex.main=1.5)
}
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: 'B04D31FA'
## - path: 'C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Proyecto_carla34.Rmd'
## - contents: <838 rows>
## Document Selection:
## - [101, 15] -- [101, 15]: ''#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/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/mult_state_carla_4.RData")
} else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_carla_4.RData")
} else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla_4.RData")
} else {
save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla_4.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