// ************************************************************* //
// ANALYSIS OF TYPE OF FACILITY AND MORTALITY - STATA CODE       //
//                                                               //
//    #0: Set PLUS directory and load CSV dataset.               //
//    #1: Baseline characteristics by type of facility.          //
//    #2: Summary follow-up statistics.                          //
//    #3: Create inverse probability weights.                    //
//    #4: Standardized nonparametric and smooth survival curves. //
//    #5: Subgroup analyses.                                     //
//    #6: Secondary analyses with profit status.                 //
//                                                               //
// ************************************************************* //

// *************************************************************
// #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:\...\Damian_et_al_2018_Plos_One_Facility_ownership_and_mortality_Data.csv"

// *************************************************************
// #1
// Table 1. Baseline characteristics by type of facility.
// Note: Unweighted counts and weighted percentages.

* Define sampling weight equal to the inverse of the probability of selection
svyset [pweight=sweight]

* Weighted mean (SD) age at baseline
svy: mean age
estat sd

* Overall distribution by type of facility
tabulate factype
svy: tabulate factype

* Baseline characteristics by type of facility
foreach x of varlist agecat sex educ staycat dementia ndiscat barthelcat {
   tabulate `x' factype, column chi2
   svy: tabulate `x' factype, column pearson
   }

svyset, clear

// *************************************************************
// #2
// Summary follow-up statistics.

* Set survival time data
stset time, id(id) failure(dth==1)
stdescribe

* Median and interquartile range of follow-up time
summarize time, detail

* Unweighted person-years and deaths by type of facility
stptime, by(factype) per(100) dd(1)

* Sampling-weighted mortality rates by type of facility
streset [pweight=sweight] // Include sampling weights in survival time setting
stptime, by(factype) per(100) dd(1)

* Sampling-weighted cumulative mortality at 2, 5, and 10 years of follow-up by type of facility
sts list, failure at(2 5 10) by(factype)

* Sampling-weighted follow-up times to 35%, 50%, and 65% cumulative mortality risks by type of facility
sts generate s_km = s, by(factype)

forvalues i = 0/3 {
   summarize _t if factype==`i' & s_km<=0.65
   summarize _t if factype==`i' & s_km<=0.50
   summarize _t if factype==`i' & s_km<=0.35
   }

drop s_km
stset, clear

// *************************************************************
// #3
// Create inverse probability weights to standardize survival curves for each type of facility
// to the weighted distribution of baseline confounders in the overall study population.

* Set sampling weights
svyset [pweight=sweight]

* Sampling-weighted marginal proportions in each type of facility
svy: proportion factype
matrix mp = e(b)

* Standardization 1 by age, sex, educational level, and length of stay

* Sampling-weighted multinomial logistic regression to estimate each resident's population probability 
* of being in its own type of facility conditional on the observed confounders
svy: mlogit factype i.agecat sex educ staycat
predict cp0-cp3 if e(sample), pr

* Inverse probability weights further stabilized by the sampling-weighted marginal proportions in each type of facility
generate stdweight1 = .
forvalues i = 0/3 {
   replace stdweight1 = mp[1,`i'+1]/cp`i' if factype==`i'
   }
drop cp0-cp3

* Define combined weights as the product of sampling weights and standardization weights
generate cweight1 = sweight*stdweight1
label variable cweight1 "Combined weight for standardization 1"
drop stdweight1

* Descriptive statistics of combined weights
summarize cweight1, detail

* Fully-weighted frequency distributions of baseline confounders for each type of facility
svyset [pweight=cweight1] // Set combined weights

foreach x of varlist agecat sex educ staycat {
   svy: tabulate `x' factype, column percent
   }

* Standardization 2 by age, sex, educational level, length of stay, dementia, number of chronic conditions other than dementia, and functional dependency
svyset [pweight=sweight]

svy: mlogit factype i.agecat sex educ staycat dementia ndiscat i.barthelcat
predict cp0-cp3 if e(sample), pr

generate stdweight2 = .
forvalues i = 0/3 {
   replace stdweight2 = mp[1,`i'+1]/cp`i' if factype==`i'
   }
drop cp0-cp3
matrix drop _all

generate cweight2 = sweight*stdweight2
label variable cweight2 "Combined weight for standardization 2"
drop stdweight2

summarize cweight2, detail

svyset [pweight=cweight2]

foreach x of varlist agecat sex educ staycat dementia ndiscat barthelcat { // Supplementary Table S2
   svy: tabulate `x' factype, column percent
   }

svyset, clear

// *************************************************************
// #4
// Nonparametric and smooth survival curves weighted by combined weights and stratified by type of facility.

* Standardization 1 by age, sex, educational level, and length of stay

* Set survival time data with combined weights
stset time [pweight=cweight1], id(id) failure(dth==1)
stdescribe

* Kaplan-Meier survival curves by type of facility 
sts generate s_km_std1 = s, by(factype)

* Spline-based survival model stratified by type of facility (with a single internal knot at the 50th percentile)
tabulate factype, generate(factype)
stpm factype2 factype3 factype4, df(2) scale(hazard) stratify(factype2 factype3 factype4)

* Predicted log cumulative hazard function
predict lnch_spl_std1, cumhazard
predict se_lnch_spl_std1, cumhazard stdp

* Predicted centiles of survival times
predict t35_spl_std1, centile(35)
predict t50_spl_std1, centile(50)
predict t65_spl_std1, centile(65)
tab1 t35_spl_std1 t50_spl_std1 t65_spl_std1

predict se_t35_spl_std1, centile(35) stdp
predict se_t50_spl_std1, centile(50) stdp
predict se_t65_spl_std1, centile(65) stdp
tab1 se_t35_spl_std1 se_t50_spl_std1 se_t65_spl_std1

drop t35_spl_std1 t50_spl_std1 t65_spl_std1 se_t35_spl_std1 se_t50_spl_std1 se_t65_spl_std1

* Standardization 2 by age, sex, educational level, length of stay, dementia, number of chronic conditions other than dementia, and functional dependency
streset [pweight=cweight2]
stdescribe

sts generate s_km_std2 = s, by(factype)

stpm factype2 factype3 factype4, df(2) scale(hazard) stratify(factype2 factype3 factype4)

predict lnch_spl_std2, cumhazard
predict se_lnch_spl_std2, cumhazard stdp

predict t35_spl_std2, centile(35)
predict t50_spl_std2, centile(50)
predict t65_spl_std2, centile(65)
tab1 t35_spl_std2 t50_spl_std2 t65_spl_std2
predict se_t35_spl_std2, centile(35) stdp
predict se_t50_spl_std2, centile(50) stdp
predict se_t65_spl_std2, centile(65) stdp
tab1 se_t35_spl_std2 se_t50_spl_std2 se_t65_spl_std2
drop t35_spl_std2 t50_spl_std2 t65_spl_std2 se_t35_spl_std2 se_t50_spl_std2 se_t65_spl_std2

* Export data for further analyses in R
export delimited id factype _d _t s_km_std1 lnch_spl_std1 se_lnch_spl_std1 s_km_std2 lnch_spl_std2 se_lnch_spl_std2 ///
       using "c:\...\survival.csv", nolabel replace

drop factype1 factype2 factype3 factype4 s_km_std1 lnch_spl_std1 se_lnch_spl_std1 s_km_std2 lnch_spl_std2 se_lnch_spl_std2
stset, clear

// *************************************************************
// #5
// Subgroup analyses: Spline-based survival models weighted by combined weights and stratified by type of facility and covariate subgroup.

* By baseline age group (<84 or >=85 years)

* Create stratum-specific combined weights
svyset [pweight=sweight]

svy: mlogit factype agedic
predict mp0-mp3 if e(sample), pr

svy: mlogit factype agedic##(sex educ staycat dementia ndiscat barthelcat)
predict cp0-cp3 if e(sample), pr

generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3

generate cweight_age = sweight*stdweight
label variable cweight_age "Combined weight for age strata"
drop stdweight

svyset, clear

* Set survival time data with combined weights
stset time [pweight=cweight_age], id(id) failure(dth==1)
stdescribe

* Spline-based survival model stratified by type of facility and age group
tabulate factype, generate(factype)
forvalues i = 1/4 {
   generate age_factype`i'= agedic*factype`i'
   }
stpm factype2 factype3 factype4 agedic age_factype2 age_factype3 age_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 agedic age_factype2 age_factype3 age_factype4)

* Predicted log cumulative hazards at 5 years of follow-up
predict lnch_age, cumhazard time(5)
predict se_lnch_age, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_age se_lnch_age if agedic==`i' & factype==`j'
      }
   }
drop age_factype1-age_factype4 lnch_age se_lnch_age

stset, clear
drop cweight_age

* By sex (female or male)

svyset [pweight=sweight]

svy: mlogit factype sex
predict mp0-mp3 if e(sample), pr
svy: mlogit factype sex##(agecat educ staycat dementia ndiscat barthelcat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_sex = sweight*stdweight
label variable cweight_sex "Combined weight for sex strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_sex], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate sex_factype`i'= sex*factype`i'
   }
stpm factype2 factype3 factype4 sex sex_factype2 sex_factype3 sex_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 sex sex_factype2 sex_factype3 sex_factype4)
predict lnch_sex, cumhazard time(5)
predict se_lnch_sex, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_sex se_lnch_sex if sex==`i' & factype==`j'
      }
   }
drop sex_factype1-sex_factype4 lnch_sex se_lnch_sex

stset, clear
drop cweight_sex

* By educational level (<primary or >=primary)

svyset [pweight=sweight]

svy: mlogit factype educ
predict mp0-mp3 if e(sample), pr
svy: mlogit factype educ##(agecat sex staycat dementia ndiscat barthelcat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_educ = sweight*stdweight
label variable cweight_educ "Combined weight for education strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_educ], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate educ_factype`i'= educ*factype`i'
   }
stpm factype2 factype3 factype4 educ educ_factype2 educ_factype3 educ_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 educ educ_factype2 educ_factype3 educ_factype4)
predict lnch_educ, cumhazard time(5)
predict se_lnch_educ, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_educ se_lnch_educ if educ==`i' & factype==`j'
      }
   }
drop educ_factype1-educ_factype4 lnch_educ se_lnch_educ

stset, clear
drop cweight_educ

* By length of stay (<3 or >=3 years)

svyset [pweight=sweight]

svy: mlogit factype staycat
predict mp0-mp3 if e(sample), pr
svy: mlogit factype staycat##(agecat sex educ dementia ndiscat barthelcat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_stay = sweight*stdweight
label variable cweight_stay "Combined weight for stay strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_stay], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate stay_factype`i'= staycat*factype`i'
   }
stpm factype2 factype3 factype4 staycat stay_factype2 stay_factype3 stay_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 staycat stay_factype2 stay_factype3 stay_factype4)
predict lnch_stay, cumhazard time(5)
predict se_lnch_stay, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_stay se_lnch_stay if staycat==`i' & factype==`j'
      }
   }
drop stay_factype1-stay_factype4 lnch_stay se_lnch_stay

stset, clear
drop cweight_stay

* By dementia (no or yes)

svyset [pweight=sweight]

svy: mlogit factype dementia
predict mp0-mp3 if e(sample), pr
svy: mlogit factype dementia##(agecat sex educ staycat ndiscat barthelcat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_dem = sweight*stdweight
label variable cweight_dem "Combined weight for dementia strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_dem], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate dem_factype`i'= dementia*factype`i'
   }
stpm factype2 factype3 factype4 dementia dem_factype2 dem_factype3 dem_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 dementia dem_factype2 dem_factype3 dem_factype4)
predict lnch_dem, cumhazard time(5)
predict se_lnch_dem, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_dem se_lnch_dem if dementia==`i' & factype==`j'
      }
   }
drop dem_factype1-dem_factype4 lnch_dem se_lnch_dem

stset, clear
drop cweight_dem

* By number of chronic conditions other than dementia (0–2 or >=3)

svyset [pweight=sweight]

svy: mlogit factype ndiscat
predict mp0-mp3 if e(sample), pr
svy: mlogit factype ndiscat##(agecat sex educ staycat dementia barthelcat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_ndis = sweight*stdweight
label variable cweight_ndis "Combined weight for number of disease strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_ndis], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate ndis_factype`i'= ndiscat*factype`i'
   }
stpm factype2 factype3 factype4 ndiscat ndis_factype2 ndis_factype3 ndis_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 ndiscat ndis_factype2 ndis_factype3 ndis_factype4)
predict lnch_ndis, cumhazard time(5)
predict se_lnch_ndis, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_ndis se_lnch_ndis if ndiscat==`i' & factype==`j'
      }
   }
drop ndis_factype1-ndis_factype4 lnch_ndis se_lnch_ndis

stset, clear
drop cweight_ndis

* By functional dependency (no/mild or moderate/severe/total)

svyset [pweight=sweight]

svy: mlogit factype bartheldic
predict mp0-mp3 if e(sample), pr
svy: mlogit factype bartheldic##(agecat sex educ staycat dementia ndiscat)
predict cp0-cp3 if e(sample), pr
generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp`i'/cp`i' if factype==`i'
   }
drop mp0-mp3 cp0-cp3
generate cweight_bar = sweight*stdweight
label variable cweight_bar "Combined weight for functional dependency strata"
drop stdweight

svyset, clear

stset time [pweight=cweight_bar], id(id) failure(dth==1)
stdescribe

forvalues i = 1/4 {
   generate bar_factype`i'= bartheldic*factype`i'
   }
stpm factype2 factype3 factype4 bartheldic bar_factype2 bar_factype3 bar_factype4, df(2) scale(hazard) ///
     stratify(factype2 factype3 factype4 bartheldic bar_factype2 bar_factype3 bar_factype4)
predict lnch_bar, cumhazard time(5)
predict se_lnch_bar, cumhazard time(5) stdp
forvalues i = 0/1 {
   forvalues j = 0/3 {
      tab1 lnch_bar se_lnch_bar if bartheldic==`i' & factype==`j'
      }
   }
drop factype1-factype4 bar_factype1-bar_factype4 lnch_bar se_lnch_bar

stset, clear
drop cweight_bar

// *************************************************************
// #6
// Secondary analyses classifying facility according to profit status.

* Create combined weights for standardization 2
svyset [pweight=sweight]

svy: proportion factypeb
matrix mp = e(b)

svy: mlogit factypeb i.agecat sex educ staycat dementia ndiscat i.barthelcat
predict cp0-cp3 if e(sample), pr

generate stdweight = .
forvalues i = 0/3 {
   replace stdweight = mp[1,`i'+1]/cp`i' if factypeb==`i'
   }
drop cp0-cp3
matrix drop _all

generate cweight_profit = sweight*stdweight
label variable cweight_profit "Combined weight for profit status"
drop stdweight

summarize cweight_profit, detail

svyset [pweight=cweight_profit]

foreach x of varlist agecat sex educ staycat dementia ndiscat barthelcat {
   svy: tabulate `x' factypeb, column percent
   }

svyset, clear

* Spline-based survival model weighted by combined weights and stratified by type of facility
stset time [pweight=cweight_profit], id(id) failure(dth==1)
stdescribe

tabulate factypeb, generate(factypeb)
stpm factypeb2 factypeb3 factypeb4, df(2) scale(hazard) stratify(factypeb2 factypeb3 factypeb4)

predict lnch_profit, cumhazard time(5)
predict se_lnch_profit, cumhazard time(5) stdp
tab1 lnch_profit se_lnch_profit

predict t50_profit, centile(50)
predict se_t50_profit, centile(50) stdp
tab1 t50_profit se_t50_profit

drop factypeb1 factypeb2 factypeb3 factypeb4 lnch_profit se_lnch_profit t50_profit se_t50_profit

stset, clear