Database (step 1)

Date created: 22:33:32 25 Jul 2021.

Install commands that are unavailable or out of date.

. *<< dd_do : noout > >
. clear all

. *https://onlinelibrary.wiley.com/doi/epdf/10.1002/sim.8894
. *https://pclambert.net/pdf/Stata_Nordic2019_Lambert.pdf
. *https://slidetodoc.com/automated-reports-using-stata-chuck-huber-ph-d/
. *~Mi unidad\Alvacast\SISTRAT 2019 (github)\_supp_mstates\stata\12874_2020_1192_MOESM1_ESM.docx
. *https://opr.princeton.edu/workshops/Downloads/2015May_StataGraphicsKoffman.pdf
. *http://www.bruunisejs.dk/StataHacks/My%20commands/matprint/matprint_demo/
. *https://pure.au.dk/portal/files/140882936/ScientificWorkInStataGoneEasy.pdf
. *https://www.stata.com/meeting/nordic-and-baltic18/slides/nordic-and-baltic18_Bruun.pdf
. *https://github.com/dvorakt/TIER_exercises/blob/master/dyndoc_debt_growth/debt%20and%20growth%20stata%20dyndoc.do
. 
. cap noi which predictms
c:\ado\plus\p\predictms.ado
*! version 4.3.0 14mar2021 MJC

. if _rc==111 {
.         cap noi net install multistate, from("https://www.mjcrowther.co.uk/code/multistate") 
.         }

. cap noi which merlin
c:\ado\plus\m\merlin.ado
*! version 2.0.2 19mar2021 MJC

. if _rc==111 {
.         cap noi net install merlin, from("https://www.mjcrowther.co.uk/code/merlin/") 
.         }

. cap noi which sumat
c:\ado\plus\s\sumat.ado
*! Part of package matrixtools v. 0.28
*! Support: Niels Henrik Bruun, niels.henrik.bruun@gmail.com
*! 2021-01-03 toxl added

. if _rc==111 {
.         cap noi scc install matrixtools
.         }

. cap noi which estwrite
c:\ado\plus\e\estwrite.ado
*! version 1.2.4 04sep2009
*! version 1.0.1 15may2007 (renamed from -eststo- to -estwrite-; -append- added)
*! version 1.0.0 29apr2005 Ben Jann (ETH Zurich)

. if _rc==111 {
.         cap noi ssc install estwrite
.         }

. cap noi which winsor2
c:\ado\plus\w\winsor2.ado
*! Inspirit of -winsor-(NJ Cox) and -winsorizeJ-(J Caskey)
*! Lian Yujun, arlionn@163.com, 2013-12-25
*! 1.1 2014.12.16

. if _rc==111 {
.         cap noi ssc install winsor2
.         }       

.         

We need to obtain the file and the work folder.

. mata : st_numscalar("OK", direxists("/volumes/sdrive/data//"))

. if scalar(OK) == 1 {
.         cap noi cd "/volumes/sdrive/data//"
.         global pathdata "/volumes/sdrive/data//"
.         di "Location= ${pathdata}; Date: `c(current_date)', considering an OS `c(os)' for the user: `c(username)'"
. }

. else display "This file does not exist"
This file does not exist

. 
. mata : st_numscalar("OK", direxists("E:\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG\"))

. if scalar(OK) == 1 {
.         cap noi cd "E:\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG"
E:\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG
.         global pathdata "E:\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG"
.         global pathdata2 "E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/"
.         di "Location= ${pathdata}; Date: `c(current_date)', considering an OS `c(os)' for the user: `c(username)'"
Location= E:\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG; Date: 25 Jul 2021, considering an OS Windows for the user: andre
. }

. else display "This file does not exist"

.         
. mata : st_numscalar("OK", direxists("C:\Users\CISS Fondecyt\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG\"))

. if scalar(OK) == 1 {
.         cap noi cd "C:\Users\CISS Fondecyt\Mi unidad\Alvacast\SISTRAT 2019 (github)"
.         global pathdata "C:\Users\CISS Fondecyt\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG"
.         global pathdata2 "C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/"
.         di "Location= ${pathdata}; Date: `c(current_date)', considering an OS `c(os)' for the user: `c(username)'"
. }

. else display "This file does not exist"
This file does not exist

. 
. mata : st_numscalar("OK", direxists("C:\Users\andre\Desktop\_WO vs MG\"))

. if scalar(OK) == 1 {
.         cap noi cd "C:\Users\andre\Desktop\_WO vs MG"
.         global pathdata "C:\Users\andre\Desktop\_mult_state_ags"
.         global pathdata2 "C:/Users/andre/Desktop/_mult_state_ags/"
.         di "Location= ${pathdata}; Date: `c(current_date)', considering an OS `c(os)' for the user: `c(username)'"
. }

. else display "This file does not exist"
This file does not exist

. 
. mata : st_numscalar("OK", direxists("C:\Users\CISS Fondecyt\OneDrive\Documentos\"))

. if scalar(OK) == 1 {
.         cap noi cd "C:\Users\CISS Fondecyt\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG"
.         global pathdata "C:\Users\CISS Fondecyt\Mi unidad\Alvacast\SISTRAT 2019 (github)\_WO vs MG"
.         global pathdata2 "C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/"
.         di "Location= ${pathdata}; Date: `c(current_date)', considering an OS `c(os)' for the user: `c(username)'"
. }

. else display "This file does not exist"
This file does not exist

. 

Path data= E:\Mi unidad\Alvacast\SISTRAT 2019 (github)_WO vs MG;

Timestamp: 25 Jul 2021, considering that is a Windows OS for the username: andre

First we open the files and drop the variables that would mistakenly amplify the sample, and define labels.

The file is located and named as: E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/four_st_msprep.dta

Then we define the transition matrix:

and transform the database in a long format, according to the specifications and the transition matrix.

Finally, the database adopt the following structure:

Set the database as a renewal time.

. *stset _stop, enter(_start) failure(_status=1) //* scale(365.25) id(id) 
. 
. *file:///G:/Mi%20unidad/Alvacast/SISTRAT%202019%20(github)/_supp_mstates/stata/crowther2017%20(1).pdf
. stset _time, failure(_status==1)

     failure event:  _status == 1
obs. time interval:  (0, _time]
 exit on or before:  failure

------------------------------------------------------------------------------
     81,243  total observations
          0  exclusions
------------------------------------------------------------------------------
     81,243  observations remaining, representing
     22,520  failures in single-record/single-failure data
   46651531  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =     4,130

. 

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

Aalen-Johanssen, Non-parametric transition probabilities

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

Generated an Aalen-Johanssen estimator to obtain the transition probabilities of the data from the time 0 (from admission). For this, we separated the transition probabilities according to the setting at baseline.

. *http://fmwww.bc.edu/repec/bocode/m/msaj.ado
. msaj, transmatrix(mat_four_states) by(tipo_de_programa_2) ci

. rename (P_AJ_*) (ajprob*)

To generate figures, we select the valid transitions

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

Parametric Models

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

Intercept-only

Generamos una lista de modelos paramétricos sin variables predictivas. Entre ellas, el Exponencial, weibull, gompertz, log-logistico, log-normal y gama generalizado. Adicionalmente, se definió una serie de modelos royston-palmar con una función de splines cúblicos restringidos, en que los knots (#gl -1) se definen en cada percentil de la distribución. En este caso, no existe una variable perdictiva más que la función de riesgo acumulado. Guardamos las estimaciones en el archivo `parmodels_m_ago_c’.

. forvalues i = 1/5 {
  2.         // Exponential
.         set seed 2125
  3.         qui cap merlin (_time if _trans == `i', family(exponential, fail(_status)))
  4.         estimates store m`i'_exp
  5.         *estimates save "${pathdata2}parmodels.ster", replace   
.         // Weibull
.         set seed 2125
  6.         qui cap merlin (_time if _trans == `i', family(weibull, fail(_status)))
  7.         estimates store m`i'_weib
  8.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Gompertz
.         set seed 2125
  9.         qui cap merlin (_time if _trans == `i', family(gompertz, fail(_status)))
 10.         estimates store m`i'_gom
 11.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Log logistic
.         set seed 2125
 12.         qui cap merlin (_time if _trans == `i', family(loglogistic, fail(_status)))
 13.         estimates store m`i'_logl
 14.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Log normal
.         set seed 2125
 15.         qui cap merlin (_time if _trans == `i', family(lognormal, fail(_status)))
 16.         estimates store m`i'_logn
 17.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Generalised gamma
.         set seed 2125
 18.         qui cap merlin (_time if _trans == `i', family(ggamma, fail(_status)))
 19.         estimates store m`i'_ggam
 20.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Royston Parmar models
.         set seed 2125
 21.         forvalues j=2/10 {
 22.                 qui cap merlin (_time if _trans == `i', family(rp, df(`j') fail(_status)))
 23.                 estimates store m`i'_rp`j'
 24.                 *estimates save "${pathdata2}parmodels.ster", append    
.         }       
 25. }

. 
. estwrite _all using "${pathdata2}parmodels_m_ago_c.sters", replace
(saving m1_exp)
(saving m1_weib)
(saving m1_gom)
(saving m1_logl)
(saving m1_logn)
(saving m1_ggam)
(saving m1_rp2)
(saving m1_rp3)
(saving m1_rp4)
(saving m1_rp5)
(saving m1_rp6)
(saving m1_rp7)
(saving m1_rp8)
(saving m1_rp9)
(saving m1_rp10)
(saving m2_exp)
(saving m2_weib)
(saving m2_gom)
(saving m2_logl)
(saving m2_logn)
(saving m2_ggam)
(saving m2_rp2)
(saving m2_rp3)
(saving m2_rp4)
(saving m2_rp5)
(saving m2_rp6)
(saving m2_rp7)
(saving m2_rp8)
(saving m2_rp9)
(saving m2_rp10)
(saving m3_exp)
(saving m3_weib)
(saving m3_gom)
(saving m3_logl)
(saving m3_logn)
(saving m3_ggam)
(saving m3_rp2)
(saving m3_rp3)
(saving m3_rp4)
(saving m3_rp5)
(saving m3_rp6)
(saving m3_rp7)
(saving m3_rp8)
(saving m3_rp9)
(saving m3_rp10)
(saving m4_exp)
(saving m4_weib)
(saving m4_gom)
(saving m4_logl)
(saving m4_logn)
(saving m4_ggam)
(saving m4_rp2)
(saving m4_rp3)
(saving m4_rp4)
(saving m4_rp5)
(saving m4_rp6)
(saving m4_rp7)
(saving m4_rp8)
(saving m4_rp9)
(saving m4_rp10)
(saving m5_exp)
(saving m5_weib)
(saving m5_gom)
(saving m5_logl)
(saving m5_logn)
(saving m5_ggam)
(saving m5_rp2)
(saving m5_rp3)
(saving m5_rp4)
(saving m5_rp5)
(saving m5_rp6)
(saving m5_rp7)
(saving m5_rp8)
(saving m5_rp9)
(saving m5_rp10)
(file E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/parmodels_m_ago_c.sters saved)

Selected the models with best AICs.

In case of the more flexible parametric models (non-standard), we selected the models that showed the best trade-off between lower complexity and better fit.

With covariables

Posteriorly, we compared AICs with covariates. We added the time of arrival as a covariate to each transition (excluding the first).

. local stname `" "tipo_de_programa_2 edad_al_ing_grupos escolaridad_rec sus_principal_mod freq_cons_sus_prin compromiso_biopsicosocial tenencia_de_la_viv
> ienda_mod num_otras_sus_mod numero_de_hijos_mod_rec tipo_de_plan_res" "tipo_de_programa_2 edad_al_ing_grupos escolaridad_rec sus_principal_mod freq_cons
> _sus_prin compromiso_biopsicosocial tenencia_de_la_vivienda_mod num_otras_sus_mod numero_de_hijos_mod_rec tipo_de_plan_res" "tipo_de_programa_2 edad_al_
> ing_grupos escolaridad_rec sus_principal_mod freq_cons_sus_prin compromiso_biopsicosocial tenencia_de_la_vivienda_mod num_otras_sus_mod numero_de_hijos_
> mod_rec tipo_de_plan_res" "tipo_de_programa_2 edad_al_ing_grupos escolaridad_rec sus_principal_mod freq_cons_sus_prin compromiso_biopsicosocial tenencia
> _de_la_vivienda_mod num_otras_sus_mod numero_de_hijos_mod_rec tipo_de_plan_res _start" "tipo_de_programa_2 edad_al_ing_grupos escolaridad_rec sus_princi
> pal_mod freq_cons_sus_prin compromiso_biopsicosocial tenencia_de_la_vivienda_mod num_otras_sus_mod numero_de_hijos_mod_rec tipo_de_plan_res _start" "' 

. 
. forvalues i = 1/5 {
  2.         gettoken state stname : stname 
  3.         // Exponential
.         set seed 2125
  4.         qui cap merlin (_time `state' if _trans == `i', family(exponential, fail(_status)))
  5.         estimates store m2_`i'_exp
  6.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Weibull
.         set seed 2125
  7.         qui cap merlin (_time `state' if _trans == `i', family(weibull, fail(_status)))
  8.         estimates store m2_`i'_weib
  9.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Gompertz
.         set seed 2125
 10.         qui cap merlin (_time `state' if _trans == `i', family(gompertz, fail(_status)))
 11.         estimates store m2_`i'_gom
 12.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Log logistic
.         set seed 2125
 13.         qui cap merlin (_time `state' if _trans == `i', family(loglogistic, fail(_status)))
 14.         estimates store m2_`i'_logl
 15.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Log normal
.         set seed 2125
 16.         qui cap merlin (_time `state' if _trans == `i', family(lognormal, fail(_status)))
 17.         estimates store m2_`i'_logn
 18.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Generalised gamma
.         set seed 2125
 19.         qui cap merlin (_time `state' if _trans == `i', family(ggamma, fail(_status)))
 20.         estimates store m2_`i'_ggam
 21.         *estimates save "${pathdata2}parmodels.ster", append    
.         // Royston Parmar models
.         set seed 2125
 22.         forvalues j=2/10 {
 23.                 qui cap merlin (_time `state' if _trans == `i', family(rp, df(`j') fail(_status)))
 24.                 estimates store m2_`i'_rp`j'
 25.                 *estimates save "${pathdata2}parmodels.ster", append    
.         }       
 26. }

. 
. estwrite _all using "${pathdata2}parmodels_m2_ago_c.sters", replace
(saving m2_1_exp)
(saving m2_1_weib)
(saving m2_1_gom)
(saving m2_1_logl)
(saving m2_1_logn)
(saving m2_1_ggam)
(saving m2_1_rp2)
(saving m2_1_rp3)
(saving m2_1_rp4)
(saving m2_1_rp5)
(saving m2_1_rp6)
(saving m2_1_rp7)
(saving m2_1_rp8)
(saving m2_1_rp9)
(saving m2_1_rp10)
(saving m2_2_exp)
(saving m2_2_weib)
(saving m2_2_gom)
(saving m2_2_logl)
(saving m2_2_logn)
(saving m2_2_ggam)
(saving m2_2_rp2)
(saving m2_2_rp3)
(saving m2_2_rp4)
(saving m2_2_rp5)
(saving m2_2_rp6)
(saving m2_2_rp7)
(saving m2_2_rp8)
(saving m2_2_rp9)
(saving m2_2_rp10)
(saving m2_3_exp)
(saving m2_3_weib)
(saving m2_3_gom)
(saving m2_3_logl)
(saving m2_3_logn)
(saving m2_3_ggam)
(saving m2_3_rp2)
(saving m2_3_rp3)
(saving m2_3_rp4)
(saving m2_3_rp5)
(saving m2_3_rp6)
(saving m2_3_rp7)
(saving m2_3_rp8)
(saving m2_3_rp9)
(saving m2_3_rp10)
(saving m2_4_exp)
(saving m2_4_weib)
(saving m2_4_gom)
(saving m2_4_logl)
(saving m2_4_logn)
(saving m2_4_ggam)
(saving m2_4_rp2)
(saving m2_4_rp3)
(saving m2_4_rp4)
(saving m2_4_rp5)
(saving m2_4_rp6)
(saving m2_4_rp7)
(saving m2_4_rp8)
(saving m2_4_rp9)
(saving m2_4_rp10)
(saving m2_5_exp)
(saving m2_5_weib)
(saving m2_5_gom)
(saving m2_5_logl)
(saving m2_5_logn)
(saving m2_5_ggam)
(saving m2_5_rp2)
(saving m2_5_rp3)
(saving m2_5_rp4)
(saving m2_5_rp5)
(saving m2_5_rp6)
(saving m2_5_rp7)
(saving m2_5_rp8)
(saving m2_5_rp9)
(saving m2_5_rp10)
(file E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/_WO vs MG/parmodels_m2_ago_c.sters saved)

Selected the models with lower fit indices.

Saved at= 22:40:12 25 Jul 2021