/* Osteosarcopenia as a risk factor for fractures and mortality – 19-year follow-up of a population-based sample M. Blomqvist, M. S. Nuotio, K. Sääksjärvi, J. Pentti, S. Koskinen, S. Stenholm All analyses were conducted using SAS software version 9.4 (SAS Institute Inc., Cary, North Carolina, United States). This file contains the code used to generate the results presented in the manuscript submitted for publication in Springer Nature Aging Clinical and Experimental Research. Dataset preparation was performed in a separate script. Further information is available upon reasonable request. Corresponding author: Matias Blomqvist, MD Email: maerbl@utu.fi. *\ /*#################################################################################*/ /*### BASELINE CHARACTERISTICS /*#################################################################################*/ proc freq data=baseline_data; table osp_group; run; proc sort data=baseline_data; by osp_group; run; /* Age */ proc univariate data=baseline_data; where osp_group > .; var age; by osp_group; run; proc anova data=baseline_data; class osp_group; model age = osp_group; means osp_group / tukey; run; /* Sex */ proc freq data=baseline_data; table osp_group*sex /chisq; run; /* Smoking */ proc freq data=baseline_data; table osp_group*MB_SMOKING2 /chisq; run; /* Education */ proc freq data=baseline_data; table osp_group*education /chisq; run; /* Physical activity */ proc freq data=baseline_data; table osp_group*MB_PHYSACT /chisq; run; /* Mobility limitation */ proc freq data=baseline_data; table osp_group*MB_MOBILITY_LIMITATION /chisq; run; /******************************************************************/ /* KAPLAN-MEIER ***************************************************/ /******************************************************************/ proc format; value spop_fmt 0 = "No PS, no osteoporosis" 1 = "PS only" 2 = "Osteoporosis only" 3 = "Osteosarcopenia"; run; /* ANY FRACTURE */ proc lifetest data=baseline_data plots=survival(cl) outsurv=km_data; where osp_group in (1, 2); format osp_group spop_fmt.; time time_to_event*MB_KM_EVENT(0); strata osp_group; run; ods graphics / reset=all; ods html style=statistical; /*ods listing gpath='C:\ ... '; */ proc sort data=km_data; by osp_group; run; proc sgplot data=km_data; title 'Any fracture'; highlow x=time_to_event low=sdf_lcl high=sdf_ucl / group=osp_group lineattrs=(thickness=0) fillattrs=(transparency=0.5) name="ci_band"; band x=time_to_event lower=sdf_lcl upper=sdf_ucl / group=osp_group transparency=0.5 name="ci_band" modelname="km"; step x=time_to_event y=survival / group=osp_group lineattrs=(pattern=solid thickness=1) name="km"; xaxis label="Time (yr)" valueattrs=(size=12) labelattrs=(size=12); yaxis label="Survival Probability" valueattrs=(size=12) labelattrs=(size=12) min=0 max=1; keylegend "ci_band" / valueattrs=(size=12) title="" location=inside position=bottomleft across=1; run; ods html close; ods listing close; /* MOF */ proc lifetest data=mof plots=survival(cl) outsurv=km_data; where osp_group in (0, 1, 2, 3); /* change accordingly */ format osp_group spop_fmt.; time time_to_event*MB_KM_EVENT(0); strata osp_group; run; ods graphics / reset=all; ods html style=statistical; /*ods listing gpath='C:\ ... '; */ proc sort data=km_data; by osp_group; run; proc sgplot data=km_data; title 'Major osteoporotic fracture'; highlow x=time_to_event low=sdf_lcl high=sdf_ucl / group=osp_group lineattrs=(thickness=0) fillattrs=(transparency=0.5) name="ci_band"; band x=time_to_event lower=sdf_lcl upper=sdf_ucl / group=osp_group transparency=0.5 name="ci_band" modelname="km"; step x=time_to_event y=survival / group=osp_group lineattrs=(pattern=solid thickness=1) name="km"; xaxis label="Time (yr)" valueattrs=(size=12) labelattrs=(size=12); yaxis label="Survival Probability" valueattrs=(size=12) labelattrs=(size=12) min=0 max=1; keylegend "ci_band" / valueattrs=(size=12) title="" location=inside position=bottomleft across=1; run; ods html close; ods listing close; /* HIP FRACTURE */ proc lifetest data=lonkkamurtumat plots=survival(cl) outsurv=km_data; where osp_group in (3, 2); /* change accordingly */ format osp_group spop_fmt.; time time_to_event*MB_KM_EVENT(0); strata osp_group; run; ods graphics / reset=all; ods html style=statistical; /*ods listing gpath='C:\ ... '; */ proc sort data=km_data; by osp_group; run; proc sgplot data=km_data; title 'Hip fracture'; highlow x=time_to_event low=sdf_lcl high=sdf_ucl / group=osp_group lineattrs=(thickness=0) fillattrs=(transparency=0.5) name="ci_band"; band x=time_to_event lower=sdf_lcl upper=sdf_ucl / group=osp_group transparency=0.5 name="ci_band" modelname="km"; step x=time_to_event y=survival / group=osp_group lineattrs=(pattern=solid thickness=1) name="km"; xaxis label="Time (yr)" valueattrs=(size=12) labelattrs=(size=12); yaxis label="Survival Probability" valueattrs=(size=12) labelattrs=(size=12) min=0 max=1; keylegend "ci_band" / valueattrs=(size=12) title="" location=inside position=bottomleft across=1; run; ods html close; ods listing close; /* DEATH */ proc lifetest data=mbdata_death plots=survival(cl) outsurv=km_data; where osp_group in (2, 1); /* change accordingly */ format osp_group spop_fmt.; time time_to_event*death(0); strata osp_group; run; ods graphics / reset=all; ods html style=statistical; /*ods listing gpath='C:\ ... '; */ proc sort data=km_data; by osp_group; run; proc sgplot data=km_data; title 'Death'; highlow x=time_to_event low=sdf_lcl high=sdf_ucl / group=osp_group lineattrs=(thickness=0) fillattrs=(transparency=0.5) name="ci_band"; band x=time_to_event lower=sdf_lcl upper=sdf_ucl / group=osp_group transparency=0.5 name="ci_band" modelname="km"; step x=time_to_event y=survival / group=osp_group lineattrs=(pattern=solid thickness=1) name="km"; xaxis label="Time (yr)" valueattrs=(size=12) labelattrs=(size=12); yaxis label="Survival Probability" valueattrs=(size=12) labelattrs=(size=12) min=0 max=1; keylegend "ci_band" / valueattrs=(size=12) title="" location=inside position=bottomleft across=1; run; ods html close; ods listing close; /******************************************************************/ /* COX MACRO ******************************************************/ /******************************************************************/ %macro cox(dataset, reference, covariates); title 'Cox proportional hazards, dataset = ' &dataset.' , reference group = ' &reference. ' , covariates = '&covariates.; proc freq data=&dataset.; table osp_group*event; run; proc phreg data=&dataset.; where osp_group > .; class osp_group(ref=&reference.); model time_to_event*event(0,1) = osp_group &covariates. / rl; run; title 'Follow-up time'; proc means data=&dataset.; var time_to_event; run; %mend; /******************************************************************/ /* FINE-GRAY MACRO ************************************************/ /******************************************************************/ %macro finegray(dataset, reference, covariates); title 'Cox proportional hazards, Fine-Gray, dataset = ' &dataset.' , reference group = ' &reference. ' , covariates = '&covariates. ; proc freq data=&dataset.; table osp_group*event; run; proc phreg data=&dataset.; where osp_group > .; class osp_group(ref=&reference.); model time_to_event*event(0) = osp_group &covariates. / eventcode=2 rl; run; proc means data=&dataset.; var time_to_event; run; %mend; /******************************************************************/ /* DATA ANALYSES **************************************************/ /******************************************************************/ /* FRACTURE OUTCOMES */ /* model 1 */ %cox(mof, '0', age sex); %cox(baseline_data, '0', age sex); %cox(hip_fractures, '0', age sex); /* model 2 */ %cox(mof, '0', age sex education MB_SMOKING2 MB_PHYSACT); %cox(baseline_data, '0', age sex education MB_SMOKING2 MB_PHYSACT); %cox(hip_fractures, '0', age sex education MB_SMOKING2 MB_PHYSACT); /*model 3 */ %cox(mof, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %cox(baseline_data, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %cox(hip_fractures, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); /*model 3 + Fine-Gray */ %finegray(mof, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(baseline_data, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(hip_fractures, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); /*model 3 + Fine-Gray - OSP as reference*/ %finegray(mof, '3', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(baseline_data, '3', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(hip_fractures, '3', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); /* DEATH */ proc freq data=mbdata_death; table osp_group*event; run; /* model 1 */ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='0'); model time_to_event*death(0) = osp_group age sex / rl; run; /* model 2 */ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='0'); model time_to_event*death(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT / eventcode=1 rl; run; /* model 3 */ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='0'); model time_to_event*death(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION/ eventcode=1 rl; run; /* model 3 - OSP as reference*/ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='3'); model time_to_event*death(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION/ eventcode=1 rl; run; /* model 1 - OSP as reference*/ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='3'); model time_to_event*death(0) = osp_group age sex / rl; run; /* model 2 - OSP as reference*/ proc phreg data=mbdata_death; where osp_group > .; class osp_group(ref='2'); model time_to_event*death(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT / eventcode=1 rl; run; /******************************************************************/ /* SENSITIVITY ANALYSES *******************************************/ /******************************************************************/ /* excluded all those with previous fracture of same type*/ /* model 3 + Fine-Grey */ %finegray(prev_mof_excl, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(prev_any_frac_excl, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(prev_hip_frac_excl, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); /* include only low-energy fractures */ %finegray(any_lowenergyonly, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(mof_lowenergyonly, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(hip_lowenergyonly, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); /* restrict follow-up period to 10 yr */ %finegray(any_fracture_10y, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(mof_10y, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); %finegray(hip_fracture_10y, '0', age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION); proc freq data=mbdata_death10v; table osp_group*event; run; proc phreg data=mbdata_death10v; where osp_group > .; class osp_group(ref='0'); model time_to_event*event(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION/ eventcode=1 rl; run; /* Schoenfeld residuals: proportional hazards assumption */ proc phreg data=baseline_data; where osp_group > .; class osp_group(ref='0'); model time_to_event*event(0) = osp_group age sex education MB_SMOKING2 MB_PHYSACT MB_MOBILITY_LIMITATION / rl; output out=cox_output ressch=osp_group1_res osp_group2_res osp_group3_res age_res sex_res education_res MB_SMOKING2_res MB_PHYSACT_res MB_MOBILITY_LIMITATION_res; run; proc corr data=cox_output pearson; var osp_group1_res osp_group2_res osp_group3_res age_res sex_res education_res MB_SMOKING2_res MB_PHYSACT_res MB_MOBILITY_LIMITATION_res; with time_to_event; run;