Decision whether to use Markov or Semi-Markov

#state arrival extended (semi-)Markov to mean that the i → j transition hazard depends on thetime of arrival at state i. 

#Build a Cox proportional hazard model including treatment and time in previous state as covariates

tab_cox_markov<- data.frame()
for (i in c(2:max(trans_matrix,na.rm=T))){
coxph(Surv(Tstart,Tstop,status)~factor(tipo_de_plan_res_1)+Tstart_yr,
                  data=subset(ms_d_match_surv_res, trans==i) %>% dplyr::mutate(Tstart_yr=Tstart/365.25),method = "breslow") %>% 
    assign(paste0("CoxMarkov",i),.,envir=.GlobalEnv)
  round(exp(coef(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("HR",i),.,envir=.GlobalEnv)
  round(exp(confint(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  round(coef(summary(get(paste0("CoxMarkov",i))))[,5],4)%>% assign(paste0("P",i),.,envir=.GlobalEnv)
  data.frame(get(paste0("CI",i))) %>% dplyr::rename("Lower 95% CI"="X2.5..","Upper 95% CI"="X97.5..")%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  tab_cox_markov_add<- cbind.data.frame(plots[i, "title"],get(paste0("HR",i)),get(paste0("CI",i)),get(paste0("P",i)))
  tab_cox_markov<-rbind.data.frame(tab_cox_markov,tab_cox_markov_add)
}
invisible(formatC(c(exp(coef(get(paste0("CoxMarkov",2))))[2], exp(confint(get(paste0("CoxMarkov",2))))[2,]), format = "f", digits = 8))
invisible(scales::scientific(c(0.000416, 0.000240,  0.000592), digits = 3))
invisible(scales::scientific(c(0.000329, 0.00000404,  0.000653), digits = 3))

tab_cox_markov2<- data.frame()
for (i in c(2:max(trans_matrix,na.rm=T))){
coxph(Surv(time,status)~factor(tipo_de_plan_res_1)+Tstart_yr,
                  data=subset(ms_d_match_surv_res, trans==i) %>% dplyr::mutate(Tstart_yr=Tstart/365.25),method = "breslow") %>% 
    assign(paste0("CoxMarkov",i),.,envir=.GlobalEnv)
  round(exp(coef(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("HR",i),.,envir=.GlobalEnv)
  round(exp(confint(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  round(coef(summary(get(paste0("CoxMarkov",i))))[,5],4)%>% assign(paste0("P",i),.,envir=.GlobalEnv)
  data.frame(get(paste0("CI",i))) %>% dplyr::rename("Lower 95% CI"="X2.5..","Upper 95% CI"="X97.5..")%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  tab_cox_markov_add<- cbind.data.frame(plots[i, "title"],get(paste0("HR",i)),get(paste0("CI",i)),get(paste0("P",i)))
  tab_cox_markov2<-rbind.data.frame(tab_cox_markov2,tab_cox_markov_add)
}
#2021-08-21: con renewal da lo mismo

tab_cox_markov %>% 
  data.table(keep.rownames=T) %>% 
  dplyr::rename("Terms"="rn", "Transition"="plots[i, \"title\"]",
                "HR"="get(paste0(\"HR\", i))","P"="get(paste0(\"P\", i))") %>% 
  dplyr::mutate(Terms=dplyr::case_when(grepl("tipo_de_", Terms)~"Setting (Residential)",
                                    grepl("Tstart",Terms)~"Time in previous state(in years)")) %>% 
  dplyr::mutate(P=ifelse(P<.001,"<.001",sprintf("%1.3f",P))) %>% 
  dplyr::rename("Sig."="P") %>% 
  dplyr::mutate(`95% CIs`=paste0(sprintf("%2.2f",`Lower 95% CI`),", ",sprintf("%2.2f",`Upper 95% CI`))) %>% 
  dplyr::select(-`Lower 95% CI`,-`Upper 95% CI`) %>% 
  dplyr::select(Transition, Terms, HR, `95% CIs`, Sig.) %>% 
      knitr::kable(format= "html", format.args= list(decimal.mark= ".", big.mark= ","),
               caption="Table 11. PH Model incluiding time in previous state & Setting as a covariate",
               align= c("c",rep('c', 5)))%>%
  #kableExtra::pack_rows("Three-states", 1, 2) %>% 
  #kableExtra::pack_rows("Four-states", 3, 4) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size= 11)%>% 
  kableExtra::kable_classic() %>% 
  kableExtra::scroll_box(width = "100%", height = "350px")

Table 11. PH Model incluiding time in previous state & Setting as a covariate
Transition	Terms	HR	95% CIs	Sig.
Readmission to Second Readmission	Setting (Residential)	1.04	0.95, 1.13	0.451
Readmission to Second Readmission	Time in previous state(in years)	1.18	1.13, 1.22	<.001
Second to Third Readmission	Setting (Residential)	1.06	0.90, 1.25	0.459
Second to Third Readmission	Time in previous state(in years)	1.16	1.09, 1.24	<.001
Third to Fourth Readmission	Setting (Residential)	1.09	0.82, 1.46	0.551
Third to Fourth Readmission	Time in previous state(in years)	1.09	0.97, 1.23	0.148

#a variable appears on both the left and right sides of the formula
#this warning should be normal, since we are dealing with time to arrival at a determined state.

Time for this code chunk to run: 0 minutes

#state arrival extended (semi-)Markov to mean that the i → j transition hazard depends on thetime of arrival at state i. 

#Build a Cox proportional hazard model including treatment and time in previous state as covariates

tab_cox_markov_9s<- data.frame()
for (i in c(3:max(trans_matrix2,na.rm=T))){
coxph(Surv(Tstart,Tstop,status)~factor(tipo_de_plan_res_1)+Tstart_yr,
                  data=subset(ms_d_match_surv_oct_2022, trans==i) %>% dplyr::mutate(Tstart_yr=Tstart/365.25),method = "breslow") %>% 
    assign(paste0("CoxMarkov",i),.,envir=.GlobalEnv)
  round(exp(coef(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("HR",i),.,envir=.GlobalEnv)
  round(exp(confint(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  round(coef(summary(get(paste0("CoxMarkov",i))))[,5],4)%>% assign(paste0("P",i),.,envir=.GlobalEnv)
  data.frame(get(paste0("CI",i))) %>% dplyr::rename("Lower 95% CI"="X2.5..","Upper 95% CI"="X97.5..")%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  tab_cox_markov_add2<- cbind.data.frame(plots2[i, "title"],get(paste0("HR",i)),get(paste0("CI",i)),get(paste0("P",i)))
  tab_cox_markov_9s<-rbind.data.frame(tab_cox_markov_9s,tab_cox_markov_add2)
}
invisible(formatC(c(exp(coef(get(paste0("CoxMarkov",2))))[2], exp(confint(get(paste0("CoxMarkov",2))))[2,]), format = "f", digits = 8))
invisible(scales::scientific(c(0.000416, 0.000240,  0.000592), digits = 3))
invisible(scales::scientific(c(0.000329, 0.00000404,  0.000653), digits = 3))

tab_cox_markov2_9s<- data.frame()
for (i in c(3:max(trans_matrix2,na.rm=T))){
coxph(Surv(time,status)~factor(tipo_de_plan_res_1)+Tstart_yr,
                  data=subset(ms_d_match_surv_oct_2022, trans==i) %>% dplyr::mutate(Tstart_yr=Tstart/365.25),method = "breslow") %>% 
    assign(paste0("CoxMarkov",i),.,envir=.GlobalEnv)
  round(exp(coef(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("HR",i),.,envir=.GlobalEnv)
  round(exp(confint(get(paste0("CoxMarkov",i)))),2)%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  round(coef(summary(get(paste0("CoxMarkov",i))))[,5],4)%>% assign(paste0("P",i),.,envir=.GlobalEnv)
  data.frame(get(paste0("CI",i))) %>% dplyr::rename("Lower 95% CI"="X2.5..","Upper 95% CI"="X97.5..")%>% assign(paste0("CI",i),.,envir=.GlobalEnv)
  tab_cox_markov_add22<- cbind.data.frame(plots2[i, "title"],get(paste0("HR",i)),get(paste0("CI",i)),get(paste0("P",i)))
  tab_cox_markov2_9s<-rbind.data.frame(tab_cox_markov2_9s,tab_cox_markov_add22)
}
#2021-08-21: con renewal da lo mismo

tab_cox_markov_9s %>% 
  data.table(keep.rownames=T) %>% 
  dplyr::rename("Terms"="rn", "Transition"="plots2[i, \"title\"]",
                "HR"="get(paste0(\"HR\", i))","P"="get(paste0(\"P\", i))") %>% 
  dplyr::mutate(Terms=dplyr::case_when(grepl("tipo_de_", Terms)~"Setting (Residential)",
                                    grepl("Tstart",Terms)~"Time in previous state(in years)")) %>% 
  dplyr::mutate(P=ifelse(P<.001,"<.001",sprintf("%1.3f",P))) %>% 
  dplyr::rename("Sig."="P") %>% 
  dplyr::mutate(`95% CIs`=paste0(sprintf("%2.2f",`Lower 95% CI`),", ",sprintf("%2.2f",`Upper 95% CI`))) %>% 
  dplyr::select(-`Lower 95% CI`,-`Upper 95% CI`) %>% 
  dplyr::select(Transition, Terms, HR, `95% CIs`, Sig.) %>% 
      knitr::kable(format= "html", format.args= list(decimal.mark= ".", big.mark= ","),
               caption="Table 11. PH Model incluiding time in previous state & Setting as a covariate",
               align= c("c",rep('c', 5)))%>%
  #kableExtra::pack_rows("Three-states", 1, 2) %>% 
  #kableExtra::pack_rows("Four-states", 3, 4) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size= 11)%>% 
  kableExtra::kable_classic() %>% 
  kableExtra::scroll_box(width = "100%", height = "350px")

Table 11. PH Model incluiding time in previous state & Setting as a covariate
Transition	Terms	HR	95% CIs	Sig.
Completion to Readmission	Setting (Residential)	1.69	1.51, 1.90	<.001
Completion to Readmission	Time in previous state(in years)	1.87	1.67, 2.10	<.001
Non-completion to Readmission	Setting (Residential)	1.55	1.45, 1.65	<.001
Non-completion to Readmission	Time in previous state(in years)	2.23	2.06, 2.41	<.001
Readmission to Second Readmission (TC)	Setting (Residential)	1.16	0.93, 1.44	0.199
Readmission to Second Readmission (TC)	Time in previous state(in years)	1.17	1.08, 1.27	<.001
Readmission to Second Readmission (TNC)	Setting (Residential)	1.00	0.90, 1.12	0.938
Readmission to Second Readmission (TNC)	Time in previous state(in years)	1.19	1.13, 1.25	<.001
Second to Third Readmission (TC)	Setting (Residential)	0.91	0.61, 1.37	0.661
Second to Third Readmission (TC)	Time in previous state(in years)	1.34	1.16, 1.56	<.001
Second to Third Readmission (TNC)	Setting (Residential)	1.14	0.94, 1.38	0.174
Second to Third Readmission (TNC)	Time in previous state(in years)	1.12	1.03, 1.21	0.005

#a variable appears on both the left and right sides of the formula
#this warning should be normal, since we are dealing with time to arrival at a determined state.

Time for this code chunk to run: 0 minutes

overall_test<-function(x){
  x<-deparse(substitute(x))
  paste0("Chi2(", get(x)[["Nsub"]],")=",
  round(get(x)[["orig_ch_stat"]],2),", p=",
  sprintf("%1.4f",get(x)[["p_ch_stat_wb"]]))
} 

#_#_#_#__#_#_#_#_#_#__###### markov test ##### _#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

ms_d_match_surv_res_tm<- ms_d_match_surv_res
ms_d_match_surv_oct_2022_tm<- ms_d_match_surv_oct_2022


attr(ms_d_match_surv_res_tm, "trans")<-trans_matrix
attr(ms_d_match_surv_oct_2022_tm, "trans")<-trans_matrix2


#Time grid
tseq <- seq(1,1837,by=30) 

#Three approaches to testing are considered; i) A simple method based on including 
#time of entry into the state as a covariate in a Cox model for each transition 
#intensity ii) Use of the stratified version of the Commenges-Andersen test 2 
#for a univariate frailty, and iii) A novel class of tests based on families of 
#log-rank statistics, where patients are grouped by their state occupancy at landmark times.

invisible("Pruebas")

start_time <- Sys.time()

#ms_d_match_surv_oct_2022_tm[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),], 
cox_markov_test_prueba12<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm,
  formula="tipo_de_plan_res_1", id = "id", transition=1, grid=tseq, B=10)

overall_test(cox_markov_test_prueba12)

## [1] "Chi2(22452)=0, p=0.0000"

Time for this code chunk to run: 0 minutes

cox_markov_test_prueba13<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, formula="tipo_de_plan_res_1", id = "id", transition=2, grid=tseq, B=200)

overall_test(cox_markov_test_prueba13)

## [1] "Chi2(22452)=0, p=0.0000"

cox_markov_test_prueba24<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, formula="tipo_de_plan_res_1", id = "id", transition=3, grid=tseq, B=200)

overall_test(cox_markov_test_prueba24)

## [1] "Chi2(5356)=35.3, p=0.0100"

cox_markov_test_prueba35<-
mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, formula="tipo_de_plan_res_1", id = "id", transition=4, grid=tseq, B=200)

overall_test(cox_markov_test_prueba35)

## [1] "Chi2(13067)=68.71, p=0.0000"

end_time <- Sys.time()

print("Time taken in process")

## [1] "Time taken in process"

end_time - start_time

## Time difference of 5.317755 hours

#The model considered the transition from intermediate states to our absorbing state (being readmitted at the fourth time) is explained by the time spent in the previous health state. This covariate (time in the previous state) was shown to be statistically significant (p<.001); results indicated a longer duration spent in the first treatment is associated with increased risk of readmission. Therefore, a semi-Markov (called a Markov renewal model) or clock reset approach should be undertaken for both models.

Time for this code chunk to run: 0.1 minutes

#ms_d_match_surv_oct_2022_tm[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),], 
cox_markov_test_prueba46<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm,#[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),], 
  formula="tipo_de_plan_res_1", id = "id", transition=5, grid=tseq, B=200)

overall_test(cox_markov_test_prueba46)

## [1] "Chi2(1518)=8.42, p=0.0050"

cox_markov_test_prueba57<-
mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, #[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),], 
  formula="tipo_de_plan_res_1", id = "id", transition=6, grid=tseq, B=200)

overall_test(cox_markov_test_prueba57)

## [1] "Chi2(4046)=22.96, p=0.0000"

cox_markov_test_prueba68<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, #[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),], 
  formula="tipo_de_plan_res_1", id = "id", transition=7, grid=tseq, B=200)

overall_test(cox_markov_test_prueba68)

## [1] "Chi2(447)=8.85, p=0.0100"

cox_markov_test_prueba79<-
  mstate::MarkovTest(ms_d_match_surv_oct_2022_tm, #[which(ms_d_match_surv_oct_2022_tm$id %in% 1:20),],
  formula="tipo_de_plan_res_1", id = "id", transition=8, grid=tseq, B=200)

overall_test(cox_markov_test_prueba79)

## [1] "Chi2(1305)=9.21, p=0.0100"

end_time2 <- Sys.time()

Time for this code chunk to run: 0.6 minutes

Session Info

Sys.getenv("R_LIBS_USER")

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

rstudioapi::getSourceEditorContext()

## Document Context: 
## - id:        '3D1B67AC'
## - path:      'C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Matching_Process1_25_NOV_22.Rmd'
## - contents:  <403 rows>
## Document Selection:
## - [151, 1] -- [151, 1]: ''

if (grepl("CISS Fondecyt",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_125_apr22.RData")
  } else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_125_apr22.RData")
  } else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_125_apr22.RData")
  } else if (grepl("G:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_125_apr22.RData")
  } else {
    save.image("~/mult_state_125_apr22.RData")
    path.expand("~/mult_state_125_apr22.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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] etm_1.1.1         data.table_1.14.2 forcats_0.5.1     stringr_1.4.0    
##  [5] dplyr_1.0.7       purrr_0.3.4       readr_2.0.1       tidyr_1.1.3      
##  [9] tibble_3.0.3      ggplot2_3.3.5     tidyverse_1.3.1   mstate_0.3.1     
## [13] survival_3.1-12  
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.0          lubridate_1.7.10  webshot_0.5.2     httr_1.4.2       
##  [5] tools_4.0.2       backports_1.1.7   bslib_0.2.5.1     utf8_1.1.4       
##  [9] R6_2.5.1          DBI_1.1.1         colorspace_1.4-1  withr_2.5.0      
## [13] tidyselect_1.1.1  gridExtra_2.3     emmeans_1.6.3     curl_4.3         
## [17] compiler_4.0.2    cli_3.0.1         rvest_1.0.1       xml2_1.3.2       
## [21] sandwich_3.0-1    sass_0.4.0        scales_1.1.1      survMisc_0.5.5   
## [25] mvtnorm_1.1-2     systemfonts_1.0.2 digest_0.6.25     foreign_0.8-81   
## [29] svglite_2.0.0     rmarkdown_2.11    rio_0.5.27        pkgconfig_2.0.3  
## [33] htmltools_0.5.1.1 highr_0.9         dbplyr_2.1.1      rlang_0.4.11     
## [37] readxl_1.3.1      rstudioapi_0.13   jquerylib_0.1.4   generics_0.1.3   
## [41] zoo_1.8-9         jsonlite_1.7.2    zip_2.2.0         car_3.0-11       
## [45] magrittr_2.0.1    kableExtra_1.3.4  Matrix_1.2-18     Rcpp_1.0.7       
## [49] munsell_0.5.0     fansi_0.4.1       abind_1.4-5       lifecycle_1.0.1  
## [53] stringi_1.4.6     multcomp_1.4-17   yaml_2.2.1        carData_3.0-4    
## [57] MASS_7.3-51.6     grid_4.0.2        parallel_4.0.2    crayon_1.4.1     
## [61] survminer_0.4.9   lattice_0.20-41   haven_2.4.3       splines_4.0.2    
## [65] hms_1.1.0         knitr_1.33        pillar_1.7.0      ggpubr_0.4.0     
## [69] estimability_1.3  ggsignif_0.6.3    codetools_0.2-16  reprex_2.0.1     
## [73] glue_1.4.1        evaluate_0.14     modelr_0.1.8      vctrs_0.3.8      
## [77] tzdb_0.1.2        cellranger_1.1.0  gtable_0.3.1      km.ci_0.5-2      
## [81] assertthat_0.2.1  xfun_0.25         openxlsx_4.2.4    xtable_1.8-4     
## [85] broom_0.7.9       coda_0.19-4       rstatix_0.7.0     ggcorrplot_0.1.3 
## [89] viridisLite_0.4.1 KMsurv_0.1-5      TH.data_1.0-10    ellipsis_0.3.2

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

Time for this code chunk to run: 0 minutes

Ambulatory or residential? a multi-state analysis of treatments for substance use disorders (Step 1)

December 01, 2022

Decision whether to use Markov or Semi-Markov

Session Info