// ************************************************************** // // ANALYSIS OF SOCIAL ENGAGEMENT AND MORTALITY - STATA CODE // // // // #0: Set PLUS directory and load CSV dataset. // // #1: Baseline characteristics by level of social engagement. // // #2: Summary follow-up statistics. // // #3: Create inverse probability weights. // // #4: Standardized nonparametric and smooth survival curves. // // #5: 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:\...\Pastor-Barriuso_et_al_2020_BMC_Geriatrics_Social_engagement_and_mortality_Data.csv" // ************************************************************** // #1 // Table 1. Baseline characteristics by level of social engagement. * Define sampling weight equal to the inverse of the probability of selection svyset [pweight=sweight] * Overall distribution by level of social engagement tabulate socpart svy: tabulate socpart, percent format(%6.1f) * Baseline characteristics by level of social engagement foreach x of varlist agecat sex educ marital facowner facsize staycat caregiver extvisit ndiscat barthelcat { tabulate `x' socpart svy: tabulate `x' socpart, column pearson percent format(%6.1f) } svyset, clear // ************************************************************** // #2 // Summary follow-up statistics. * Set survival time data stset time, id(id) failure(dth==1) stdescribe * Unweighted person-years and deaths by level of social engagement stptime, by(socpart) per(100) dd(1) * Sampling-weighted mortality rates by level of social engagement streset [pweight=sweight] stptime, by(socpart) per(100) dd(1) * Sampling-weighted cumulative mortality at 2, 5, and 10 years of follow-up by level of social engagement sts list, failure at(2 5 10) by(socpart) stset, clear // ************************************************************** // #3 // Create inverse probability weights to standardize survival curves for each level of social engagement // to the weighted distribution of baseline confounders in the overall study population. * Set sampling weights svyset [pweight=sweight] * Sampling-weighted marginal proportions in each level of social engagement svy: proportion socpart matrix mp = e(b) * Standardization 1 by age, sex, educational level, and marital status * Sampling-weighted polytomous logistic model to estimate each resident's population probability * of being in its own level of social engagement conditional on the observed confounders svy: mlogit socpart i.agecat sex i.educ i.marital predict cp0-cp2 if e(sample), pr * Inverse probability weights further stabilized by the sampling-weighted marginal proportions in each level of social engagement generate stdweight1 = . forvalues i = 0/2 { replace stdweight1 = mp[1,`i'+1]/cp`i' if socpart==`i' } drop cp0-cp2 * 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 distributions of baseline confounders for each level of social engagement svyset [pweight=cweight1] foreach x of varlist agecat sex educ marital { svy: tabulate `x' socpart, column pearson percent format(%6.1f) } * Standardization 2 by age, sex, educational level, marital status, * facility ownership, size, length of stay, assigned caregiver, and frequency of external visits svyset [pweight=sweight] svy: mlogit socpart i.agecat sex i.educ i.marital facowner i.facsize i.staycat caregiver i.extvisit predict cp0-cp2 if e(sample), pr generate stdweight2 = . forvalues i = 0/2 { replace stdweight2 = mp[1,`i'+1]/cp`i' if socpart==`i' } drop cp0-cp2 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 marital facowner facsize staycat caregiver extvisit { svy: tabulate `x' socpart, column pearson percent format(%6.1f) } * Standardization 3 by age, sex, educational level, marital status, * facility ownership, size, length of stay, assigned caregiver, frequency of external visits, * number of chronic conditions, and functional dependency svyset [pweight=sweight] svy: mlogit socpart i.agecat sex i.educ i.marital facowner i.facsize i.staycat caregiver i.extvisit i.ndiscat i.barthelcat predict cp0-cp2 if e(sample), pr generate stdweight3 = . forvalues i = 0/2 { replace stdweight3 = mp[1,`i'+1]/cp`i' if socpart==`i' } drop cp0-cp2 matrix drop _all generate cweight3 = sweight*stdweight3 label variable cweight3 "Combined weight for standardization 3" drop stdweight3 summarize cweight3, detail svyset [pweight=cweight3] foreach x of varlist agecat sex educ marital facowner facsize staycat caregiver extvisit ndiscat barthelcat { svy: tabulate `x' socpart, column pearson percent format(%6.1f) } svyset, clear // ************************************************************** // #4 // Nonparametric and smooth survival curves weighted by combined weights and stratified by level of social engagement. * Standardization 1 * Set survival time data with combined weights stset time [pweight=cweight1], id(id) failure(dth==1) stdescribe * Kaplan-Meier survival curves by level of social engagement sts generate s_km_std1 = s, by(socpart) * Spline-based survival model stratified by level of social engagement (with two internal knots at the 33th and 67th percentiles) tabulate socpart, generate(socpart) stpm socpart2 socpart3, df(3) scale(hazard) stratify(socpart2 socpart3) * Predicted log cumulative hazard functions predict lnch_spl_std1, cumhazard predict se_lnch_spl_std1, cumhazard stdp * Predicted median survival times predict t50_spl_std1, centile(50) predict se_t50_spl_std1, centile(50) stdp tab1 t50_spl_std1 se_t50_spl_std1 drop t50_spl_std1 se_t50_spl_std1 * Standardization 2 streset [pweight=cweight2] stdescribe sts generate s_km_std2 = s, by(socpart) stpm socpart2 socpart3, df(3) scale(hazard) stratify(socpart2 socpart3) predict lnch_spl_std2, cumhazard predict se_lnch_spl_std2, cumhazard stdp predict t50_spl_std2, centile(50) predict se_t50_spl_std2, centile(50) stdp tab1 t50_spl_std2 se_t50_spl_std2 drop t50_spl_std2 se_t50_spl_std2 * Standardization 3 streset [pweight=cweight3] stdescribe sts generate s_km_std3 = s, by(socpart) stpm socpart2 socpart3, df(3) scale(hazard) stratify(socpart2 socpart3) predict lnch_spl_std3, cumhazard predict se_lnch_spl_std3, cumhazard stdp predict t50_spl_std3, centile(50) predict se_t50_spl_std3, centile(50) stdp tab1 t50_spl_std3 se_t50_spl_std3 drop t50_spl_std3 se_t50_spl_std3 * Export data for further analyses in R export delimited id socpart cweight1 cweight2 cweight3 dth time s_km_std1 lnch_spl_std1 se_lnch_spl_std1 s_km_std2 lnch_spl_std2 se_lnch_spl_std2 s_km_std3 lnch_spl_std3 se_lnch_spl_std3 /// using "c:\...\survival.csv", nolabel replace drop socpart1 socpart2 socpart3 s_km_std1 lnch_spl_std1 se_lnch_spl_std1 s_km_std2 lnch_spl_std2 se_lnch_spl_std2 s_km_std3 lnch_spl_std3 se_lnch_spl_std3 stset, clear // ************************************************************** // #5 // Subgroup analyses: Spline-based survival models weighted by combined weights and stratified by social engagement and covariate subgroup. * Overall stset time [pweight=cweight3], id(id) failure(dth==1) stdescribe tabulate socpart, generate(socpart) stpm socpart2 socpart3, df(3) scale(hazard) stratify(socpart2 socpart3) predict lnch_std3, cumhazard time(5) predict se_lnch_std3, cumhazard time(5) stdp forvalues i = 0/2 { tab1 lnch_std3 se_lnch_std3 if socpart==`i' } drop lnch_std3 se_lnch_std3 stset, clear * By baseline age group (<85 or >=85 years) svyset [pweight=sweight] generate agedic = age>=85 svy: mlogit socpart agedic predict mp0-mp2 if e(sample), pr svy: mlogit socpart agedic##(sex educ marital facowner facsize staycat caregiver extvisit ndiscat barthelcat) i.agecat predict cp0-cp2 if e(sample), pr generate stdweight = . forvalues i = 0/2 { replace stdweight = mp`i'/cp`i' if socpart==`i' } drop mp0-mp2 cp0-cp2 generate cweight_age = sweight*stdweight drop stdweight centile cweight_age, centile(1 99) replace cweight_age = r(c_1) if cweight_ager(c_2) svyset, clear stset time [pweight=cweight_age], id(id) failure(dth==1) stdescribe forvalues i = 1/3 { generate age_socpart`i'= agedic*socpart`i' } stpm socpart2 socpart3 agedic age_socpart2 age_socpart3, df(3) scale(hazard) stratify(socpart2 socpart3 agedic age_socpart2 age_socpart3) predict lnch_age, cumhazard time(5) predict se_lnch_age, cumhazard time(5) stdp forvalues i = 0/1 { forvalues j = 0/2 { tab1 lnch_age se_lnch_age if agedic==`i' & socpart==`j' } } drop agedic age_socpart1-age_socpart3 lnch_age se_lnch_age stset, clear drop cweight_age * By sex (women or men) svyset [pweight=sweight] svy: mlogit socpart sex predict mp0-mp2 if e(sample), pr svy: mlogit socpart sex##(agecat educ marital facowner facsize staycat caregiver extvisit ndiscat barthelcat) predict cp0-cp2 if e(sample), pr generate stdweight = . forvalues i = 0/2 { replace stdweight = mp`i'/cp`i' if socpart==`i' } drop mp0-mp2 cp0-cp2 generate cweight_sex = sweight*stdweight drop stdweight centile cweight_sex, centile(1 99) replace cweight_sex = r(c_1) if cweight_sexr(c_2) svyset, clear stset time [pweight=cweight_sex], id(id) failure(dth==1) stdescribe forvalues i = 1/3 { generate sex_socpart`i'= sex*socpart`i' } stpm socpart2 socpart3 sex sex_socpart2 sex_socpart3, df(3) scale(hazard) stratify(socpart2 socpart3 sex sex_socpart2 sex_socpart3) predict lnch_sex, cumhazard time(5) predict se_lnch_sex, cumhazard time(5) stdp forvalues i = 0/1 { forvalues j = 0/2 { tab1 lnch_sex se_lnch_sex if sex==`i' & socpart==`j' } } drop sex_socpart1-sex_socpart3 lnch_sex se_lnch_sex stset, clear drop cweight_sex * By facility ownership (public/subsidized or private) svyset [pweight=sweight] svy: mlogit socpart facowner predict mp0-mp2 if e(sample), pr svy: mlogit socpart facowner##(agecat sex educ marital staycat caregiver extvisit ndiscat barthelcat) i.facsize predict cp0-cp2 if e(sample), pr generate stdweight = . forvalues i = 0/2 { replace stdweight = mp`i'/cp`i' if socpart==`i' } drop mp0-mp2 cp0-cp2 generate cweight_owner = sweight*stdweight drop stdweight centile cweight_owner, centile(1 99) replace cweight_owner = r(c_1) if cweight_ownerr(c_2) svyset, clear stset time [pweight=cweight_owner], id(id) failure(dth==1) stdescribe forvalues i = 1/3 { generate owner_socpart`i'= facowner*socpart`i' } stpm socpart2 socpart3 facowner owner_socpart2 owner_socpart3, df(3) scale(hazard) stratify(socpart2 socpart3 facowner owner_socpart2 owner_socpart3) predict lnch_owner, cumhazard time(5) predict se_lnch_owner, cumhazard time(5) stdp forvalues i = 0/1 { forvalues j = 0/2 { tab1 lnch_owner se_lnch_owner if facowner==`i' & socpart==`j' } } drop owner_socpart1-owner_socpart3 lnch_owner se_lnch_owner stset, clear drop cweight_owner * By facility size (<300 or >=300 beds) svyset [pweight=sweight] generate facsizedic = facsize==2 svy: mlogit socpart facsizedic predict mp0-mp2 if e(sample), pr svy: mlogit socpart facsizedic##(agecat sex educ marital staycat caregiver extvisit ndiscat barthelcat) facowner i.facsize predict cp0-cp2 if e(sample), pr generate stdweight = . forvalues i = 0/2 { replace stdweight = mp`i'/cp`i' if socpart==`i' } drop mp0-mp2 cp0-cp2 generate cweight_size = sweight*stdweight drop stdweight centile cweight_size, centile(1 99) replace cweight_size = r(c_1) if cweight_sizer(c_2) svyset, clear stset time [pweight=cweight_size], id(id) failure(dth==1) stdescribe forvalues i = 1/3 { generate size_socpart`i'= facsizedic*socpart`i' } stpm socpart2 socpart3 facsizedic size_socpart2 size_socpart3, df(3) scale(hazard) stratify(socpart2 socpart3 facsizedic size_socpart2 size_socpart3) predict lnch_size, cumhazard time(5) predict se_lnch_size, cumhazard time(5) stdp forvalues i = 0/1 { forvalues j = 0/2 { tab1 lnch_size se_lnch_size if facsizedic==`i' & socpart==`j' } } drop facsizedic size_socpart1-size_socpart3 lnch_size se_lnch_size stset, clear drop cweight_size * By functional dependency at baseline (no or mild/moderate) svyset [pweight=sweight] generate bartheldic = barthelcat>=1 svy: mlogit socpart bartheldic predict mp0-mp2 if e(sample), pr svy: mlogit socpart bartheldic##(agecat sex educ marital facowner facsize staycat caregiver extvisit ndiscat) i.barthelcat predict cp0-cp2 if e(sample), pr generate stdweight = . forvalues i = 0/2 { replace stdweight = mp`i'/cp`i' if socpart==`i' } drop mp0-mp2 cp0-cp2 generate cweight_bar = sweight*stdweight drop stdweight centile cweight_bar, centile(1 99) replace cweight_bar = r(c_1) if cweight_barr(c_2) svyset, clear stset time [pweight=cweight_bar], id(id) failure(dth==1) stdescribe forvalues i = 1/3 { generate bar_socpart`i'= bartheldic*socpart`i' } stpm socpart2 socpart3 bartheldic bar_socpart2 bar_socpart3, df(3) scale(hazard) stratify(socpart2 socpart3 bartheldic bar_socpart2 bar_socpart3) predict lnch_bar, cumhazard time(5) predict se_lnch_bar, cumhazard time(5) stdp forvalues i = 0/1 { forvalues j = 0/2 { tab1 lnch_bar se_lnch_bar if bartheldic==`i' & socpart==`j' } } drop socpart1-socpart3 bartheldic bar_socpart1-bar_socpart3 lnch_bar se_lnch_bar stset, clear drop cweight_bar