// ********************************************************** //
// ASSOCIATION OF FALLS WITH LONG-TERM MORTALITY - STATA CODE //
//                                                            //
//    #0: Set PLUS directory and load CSV dataset.            //
//    #1: Define survey design and survival time data.        //
//    #2: Baseline characteristics by severity of fall.       //
//    #3: Mortality rates by severity of fall.                //
//    #4: Cox proportional hazards models.                    //
//    #5: Adjusted non-parametric survival curves.            //
//    #6: Cumulative hazard ratios over time.                 //
//    #7: Subgroup analyses.                                  //
//                                                            //
// ********************************************************** //

// *************************************************************
// #0
// Set PLUS directory and load CSV dataset.

* Set PLUS directory where ado-file "stpm" is installed 
sysdir set PLUS "c:\...\ado\plus"

* Load CSV dataset
import delimited "c:\...\Padrón-Monedero_et_al_2020_Plos_One_Falls_and_mortality_in_older_residents_Data.csv"

// ********************************************
// #1
// Define survey design and survival time data.

* Stratified cluster sampling with inverse probability weights
svyset cluster [pweight=peso], strata(resid2) vce(linearized) singleunit(missing)
svydescribe

* Follow-up restricted to 5 years
stset seg, id(id) failure(dead==1) exit(time 5)
stdescribe

// *********************************************
// #2
// Baseline characteristics by severity of fall.

* Unweighted sample counts
tabulate fallcat
foreach x of varlist agecat sex factype dementia ndiscat antidep nmedcat ui funcdep {
   tabulate `x' fallcat, missing
   }

* Weighted percentages
svy: tabulate fallcat, percent format(%6.1f)
foreach x of varlist agecat sex factype dementia ndiscat antidep nmedcat ui funcdep {
   svy: tabulate `x' fallcat, column pearson percent format(%6.1f)
   }

// ************************************
// #3
// Mortality rates by severity of fall.

* Unweighted deaths and person-years
stptime, by(fall) per(100) dd(1)
stptime, by(fallcat) per(100) dd(1)

* Weighted mortality rates
streset [pweight=peso] // Include weights in survival time setting

stptime, by(fall) per(100) dd(1)
stptime, by(fallcat) per(100) dd(1)

streset [pweight=1] // Remove weights from survival time setting

// ********************************
// #4
// Cox proportional hazards models.

* Cox proportional hazards models adjusted for age, sex, and type of facility
svy: stcox i.fall i.agecat sex i.factype

svy: stcox i.fallcat i.agecat sex i.factype
svy: stcox fallcat i.agecat sex i.factype // Test for linear trend

* Cox proportional hazards models further adjusted for dementia, number of chronic conditions,
* use of antidepressants, number of prescribed medications, urinary incontinence, and functional dependency
svy: stcox i.fall i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep

svy: stcox i.fallcat i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep
svy: stcox fallcat i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep

// ****************************************
// #5
// Adjusted non-parametric survival curves.

* Center indicators for categorical covariates around their weighted percentages
foreach x of varlist agecat factype ndiscat nmedcat ui funcdep {
   tabulate `x', gen(`x')
   }

foreach x of varlist agecat2 agecat3 agecat4 agecat5 sex factype2 factype3 dementia ndiscat2 ndiscat3 antidep nmedcat2 nmedcat3 ui2 ui3 funcdep2 funcdep3 {
   svy: mean `x'
   matrix wtdm = e(b)
   generate `x'_c = `x' - wtdm[1,1]
   }
matrix drop _all

* Fall-stratified Cox model
svy: stcox agecat2_c agecat3_c agecat4_c agecat5_c sex_c factype2_c factype3_c dementia_c ndiscat2_c ndiscat3_c antidep_c nmedcat2_c nmedcat3_c ui2_c ui3_c funcdep2_c funcdep3_c, strata(fallcat) basesurv(surv_cox)

* Export adjusted non-parametric survival curves for Figure 1
export delimited surv_cox _t fallcat using "c:\...\Figure1.csv" if surv_cox!=., nolabel replace

// *******************************************************************************
// #6
// Cumulative hazard ratios over time obtained from a spline-based survival model.

* Include weights in survival time setting
streset [pweight=peso]
stdescribe

* Create indicators for categorical covariates
foreach x of varlist fallcat agecat factype ndiscat nmedcat ui funcdep {
   tabulate `x', gen(`x')
   }

* Spline-based parametric survival model (with a single internal knot at the 50th percentile) stratified by severity of fall
stpm fallcat2 fallcat3 ///
     agecat2 agecat3 agecat4 agecat5 sex factype2 factype3 ///
     dementia ndiscat2 ndiscat3 antidep nmedcat2 nmedcat3 ui2 ui3 funcdep2 funcdep3, ///
     df(2) scale(hazard) cluster(cluster) stratify(fallcat2 fallcat3)

* Tests for proportional hazards
test [s0]fallcat2 [s1]fallcat2
test [s0]fallcat3 [s1]fallcat3

* Spline coefficients and their covariance matrices
matrix b = e(b)
matrix V = e(V)

matrix b1 = b[1,7],b[1,1],b[1,4]
matrix V1 = V[7,7],V[7,1],V[7,4]\ ///
            V[1,7],V[1,1],V[1,4]\ ///
            V[4,7],V[4,1],V[4,4]
matrix list b1
matrix list V1, nohalf

matrix b2 = b[1,8],b[1,2],b[1,5]
matrix V2 = V[8,8],V[8,2],V[8,5]\ ///
            V[2,8],V[2,2],V[2,5]\ ///
            V[5,8],V[5,2],V[5,5]
matrix list b2
matrix list V2, nohalf
matrix drop _all

* Extract natural cubic splines of log time
generate double lnt = I__b_0
label variable lnt "Log time"
generate double lnt_s1 = I__b_1
label variable lnt_s1 "Orthogonalized natural cubic spline term 1 for log time"

* Export data for Figure 2
export delimited _t lnt lnt_s1 using "c:\...\Figure2.csv" if e(sample), nolabel replace

streset [pweight=1] // Remove weights from survival time setting

// ***********************************************************
// #7
// Cox proportional hazards models in pre-specified subgroups.

* Overall
svy: stcox i.fallcat i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep, nohr

* By age
recode agecat (0/2 = 0) (3/4 = 1), generate(agedic)
svy: stcox fallcat##agedic sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#1.agedic
lincom 2.fallcat + 2.fallcat#1.agedic
test 1.fallcat#1.agedic 2.fallcat#1.agedic // Test for interaction

* By sex
svy: stcox fallcat##sex i.agecat i.factype dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#1.sex
lincom 2.fallcat + 2.fallcat#1.sex
test 1.fallcat#1.sex 2.fallcat#1.sex

* By type of facility
svy: stcox fallcat##2.factype i.agecat sex dementia i.ndiscat antidep i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#2.factype
lincom 2.fallcat + 2.fallcat#2.factype
test 1.fallcat#2.factype 2.fallcat#2.factype

* By dementia
svy: stcox fallcat##dementia i.agecat sex i.factype i.ndiscat antidep i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#1.dementia
lincom 2.fallcat + 2.fallcat#1.dementia
test 1.fallcat#1.dementia 2.fallcat#1.dementia

* By no. of medical conditions
svy: stcox fallcat##2.ndiscat i.agecat sex i.factype dementia antidep i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#2.ndiscat
lincom 2.fallcat + 2.fallcat#2.ndiscat
test 1.fallcat#2.ndiscat 2.fallcat#2.ndiscat

* By use of antidepressants
svy: stcox fallcat##antidep i.agecat sex i.factype dementia i.ndiscat i.nmedcat i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#1.antidep
lincom 2.fallcat + 2.fallcat#1.antidep
test 1.fallcat#1.antidep 2.fallcat#1.antidep

* By no. of prescribed medications
svy: stcox fallcat##2.nmedcat i.agecat sex i.factype dementia i.ndiscat antidep i.ui i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#2.nmedcat
lincom 2.fallcat + 2.fallcat#2.nmedcat
test 1.fallcat#2.nmedcat 2.fallcat#2.nmedcat

* By urinary incontinence
svy: stcox fallcat##0.ui i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.funcdep, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#0.ui
lincom 2.fallcat + 2.fallcat#0.ui
test 1.fallcat#0.ui 2.fallcat#0.ui

* By functional dependency
svy: stcox fallcat##0.funcdep i.agecat sex i.factype dementia i.ndiscat antidep i.nmedcat i.ui, nohr
lincom 1.fallcat
lincom 2.fallcat
lincom 1.fallcat + 1.fallcat#0.funcdep
lincom 2.fallcat + 2.fallcat#0.funcdep
test 1.fallcat#0.funcdep 2.fallcat#0.funcdep