* Journal of Statistical Software; * %JM: A SAS Macro to Fit Jointly Generalized Mixed Models for Longitudinal Data and Time-to-Event Responses; * This program produces the SAS output shown in Sections 4.1 & 4.2 of the manuscript as well as the results provided in Table 1 and Figures 1 to 4; * Authors: Alberto Garcia-Hernandez & Dimitris Rizopoulow; libname mydata 'JM/Data'; * Change location of the input AIDS dataset accordingly; * You may need to replace / by \ depending on your environment (PC or Unix); options mautolocdisplay mautosource sasautos = ("JM/macros" sasautos); * Change location of %JM macros accordingly; * You may need to replace / by \ depending on your environment (PC or Unix); options mprint nodate center; %let folder = JM/Output; * Change desired location of the output accordingly; * Do NOT quote the folder; * You may need to replace / by \ depending on your environment (PC or Unix); ********************************************************************************; * The code below fits the model of the example in Section 4.1 of the manuscript ; * NOTE: The execution of this model will take about 5 minutes ; ********************************************************************************; %macro model1(n=, Inp=disjoint, nl=); %JM(Data = mydata.aids, SubjectVar = patient, LongiType = normal, LongiTimeModel = linear, LongiVar = cd4, LongiTimevar = obstime, LongiTimeInter = ddi, LongiCovariates = female, LongiGMatrix = un, EventTimeVar = time, EventVar = death, EventVal = 1, EventModel = rpnaturalcubic, EventCovariates = ddi female, EventNknots = 3, InitialParameters = &inp., NLMIXEDOptions = &nl., OutputParameters = mydata.Example1_param_&n., OutputPredictions = mi, OutputPredictData = mydata.Example1_mi_&n., SharedParam = current_value slope, ListingFile = "&folder./Example_1_&n..lst", AdditionalOptions = calculateexectime); %mend model1; %model1(n=0, nl=gconv=0); * This code below is needed to store the dataset that contains the knots calculated for the event model in a permanent location; * We will need this data later to produce the lower quadrants of Figure 3; data mydata.Example1_CALC_KNOTS_EVENT_MODEL; set _JMCALCULATED_KNOTS_EVENT_MODEL; run; ********************************************************************************; * The code below fits the model of the example in Section 4.2 of the manuscript ; * WARNING: The execution of this model will take about 90 minutes ; ********************************************************************************; %macro model2(n=, nl=, inp=, gnl=); %JM(Data = mydata.aids, SubjectVar = patient, LongiType = binary, LongiTimeModel = linear, LongiVar = cd4_a30, LongiTimevar = obstime, LongiTimeInter = ddi, LongiCovariates = female, LongiGMatrix = un, LongiModelOptions = &gnl., LongiGLINLOptions = , EventTimeVar = time, EventVar = death, EventVal = 1, EventModel = rpnaturalcubic, EventCovariates = ddi female, EventNknots = 3, NLMIXEDOptions = &nl., OutputParameters = mydata.Example2_param_&n., OutputPredictions = mi, OutputPredictData = mydata.Example2_mi_&n., InitialParameters = &inp., SharedParam = coefficients, SharedCoefficients = bi0 bi1, ListingFile = "&folder./Example_2_&n..lst", AdditionalOptions = calculateexectime); %mend model2; %model2(n=0, gnl= method=quad, nl=gconv=0, inp = disjoint); /****************************************************************************************************************************/ /* */ /* NOTE: THE CODE BELOW PRODUCES MANUSCRIPT'S TABLE 1 AND FIGURES 1, 2, 3 AND 4 */ /* WARNING: EXECUTION OF ALL CODE BELOW MAY TAKE 10-15 HOURS */ /* */ /****************************************************************************************************************************/ *********************************************************************************************************************************; * These macro calls are to fit the same model of Section 4.1 with additional quadrature points (feeds manuscript's Figure 1) ; *********************************************************************************************************************************; %model1(n=1, nl=gconv=0 qpoints=1); %model1(n=2, nl=gconv=0 qpoints=5, Inp = mydata.Example1_param_1); %model1(n=3, nl=gconv=0 qpoints=9, Inp = mydata.Example1_param_2); %model1(n=4, nl=gconv=0 qpoints=13, Inp = mydata.Example1_param_3); %model1(n=5, nl=gconv=0 qpoints=15, Inp = mydata.Example1_param_4); %model1(n=6, nl=gconv=0 qpoints=17, Inp = mydata.Example1_param_5); *******************************************************************************************************; * The code below produces Figure 1:Example 1, joint model estimates versus number of quadrature points.; *******************************************************************************************************; proc format; value $parf 'L_INTERCEPT' = 'Intercept' 'L_TIME_B1' = 'Slope' 'L_TIME_B1_BY_DDI' = 'Slope by DDI' 'L_FEMALE' = 'Female' 'L_LOGSD_INTERCEPT' = 'Log SD bi0' 'L_COV_B0_B1' = 'Covar(bi0, bi1)' 'L_LOGSD_TIME_B1' = 'Log SD bi1' 'L_LOGSD_RESIDUAL' = 'Log SD Residual' 'E_DDI' = 'DDI (Event)' 'E_FEMALE' = 'Female (Event)' 'E_TIME_B0' = 'Time Param 0' 'E_TIME_B1' = 'Time Param 1' 'E_TIME_B2' = 'Time Param 2' 'E_TIME_B3' = 'Time Param 3' 'E_TIME_B4' = 'Time Param 4' 'ASSOCT_TD' = 'Current-value Dependent Association' 'ASSOCT_SL' = 'Slope Dependent Association' ; run; data parametersE1; set mydata.Example1_param_1 (keep=parameter estimate in=a) mydata.Example1_param_2 (keep=parameter estimate in=b) mydata.Example1_param_3 (keep=parameter estimate in=c) mydata.Example1_param_4 (keep=parameter estimate in=d) mydata.Example1_param_5 (keep=parameter estimate in=e) mydata.Example1_param_6 (keep=parameter estimate in=f) ; if a then qpoints=1; else if b then qpoints=5; else if c then qpoints=9; else if d then qpoints=13; else if e then qpoints=15; else if f then qpoints=17; format parameter $parf.; abs=abs(estimate); run; ods _all_ close; options nonumber orientation=landscape papersize=(6in 4in); goptions device=ACTXIMG; ods pdf file="&folder./Example1EstimatesVSQPoints.pdf" style=journal nopdfnote; ods graphics / border=off width=1200px height=800px; proc sgplot data=parametersE1; xaxis label="Number of quadrature points" values=(1 to 17 by 4); yaxis label="Abs(estimate) (log-scale)" type=log; series x=qpoints y=abs/group=parameter ; where parameter in ('ASSOCT_TD', 'ASSOCT_SL','L_LOGSD_RESIDUAL','L_LOGSD_INTERCEPT','L_COV_B0_B1','L_LOGSD_TIME_B1', 'E_DDI', 'E_FEMALE'); run; ods pdf close; *************************************************************************************************************************************; * These additional calls are to fit the same model of Section 4.1 with different optimization techniques (feeds manuscript's Table 1); *************************************************************************************************************************************; %model1(n=7, nl=gconv=0 technique=dbldog); %model1(n=8, nl=gconv=0 technique=quanew); %model1(n=9, nl=gconv=0 technique=newrap); %model1(n=10, nl=gconv=0 technique=nrridg); %model1(n=11, nl=gconv=0 technique=trureg); ****************************************************************************************; * The following code creates Figure 2 (manuscript's Section 4.1) ; * Example 1, joint model empirical Bayes estimates versus time of 20 subjects ; ****************************************************************************************; ods _all_ close; options nonumber orientation=landscape papersize=(6in 4in); goptions device=ACTXIMG; ods pdf file="&folder./EBEstimates.pdf" style=journal nopdfnote; ods graphics / border=off width=1200px height=800px; proc sgpanel data=mydata.Example1A_mi_0 noautolegend; panelby patient / rows=4 columns=5; rowaxis label="CD4 cell count" values=(0 to 25 by 5); colaxis label="Time (months)" values=(0 to 18 by 6); scatter x=obstime y=CD4; series x=obstime y=pred; where JMDist = "LONGI" and patient<=20; run; ods pdf close; ****************************************************************************************; * The following code creates Figure 3 (manuscript's Section 4.1) ; ****************************************************************************************; * Calculation of residuals (to be used in upper left quadrant of Figure 3); data _jmsbjmean2; set mydata.Example1A_mi_0; resid=cd4-pred; where JMDist ="LONGI"; run; * Standardization of residuals (to be used in upper right quadrant of Figure 3); proc sql noprint; select exp(estimate) into: SD_RESIDUAL from mydata.Example1A_param_0 where parameter='L_LOGSD_RESIDUAL'; quit; data _jmsbjmean3; set _jmsbjmean2; resid=resid/&SD_RESIDUAL.; where JMDist ="LONGI"; run; * Calculation of the normal quintiles (to be used in upper right quadrant of Figure 3); proc rank data=_jmsbjmean3 normal=blom out=Quant; var resid; ranks resid_quant; run; *************************************************************************; * The code below creates the plots in the lower quadrants ; * Estimation of population survival is not straight-forward using ; * models with random effects. ; * ; * The conditional probability of a censored TTE observation T without ; * longitudinal observations is exp{-H(T, theta, bi)}= S(T, theta, bi) ; * The marginal prob. = Integral exp{S(T, theta, bi)} p(bi) dbi ; * ; * In the code below we will apply the fitted model to a single ; * censored observation T with a given set of baseline covariates. ; * And we will repeat this process for T = 0.5 to 21 and for each ; * combination of treatment and gender. ; * ; * The JM macro provides X = -2LL for each model (= datapoint). ; * ; * Y = exp(X/-2) is thus the marginal probability of this observation T ; * that, since it is censored, represents the expected population ; * survival at T for the given set of covariates. ; *************************************************************************; * Creation of the matrix of baseline covariates for both models; proc sort data=mydata.aids (keep=ddI female) out=_xmatrix nodupkey; by ddI female; run; data _xmatrix; set _xmatrix; cov_n=_n_; if ddi=0 and female=0 then cov_l="Pcb/Males"; else if ddi=0 and female=1 then cov_l="Pcb/Females"; else if ddi=1 and female=0 then cov_l="ddI/Males"; else if ddi=1 and female=1 then cov_l="ddI/Females"; label cov_l='Covariates'; run; * Creating a set of times to use in the x axis of the survival plot; data _times; do time=0.5 to 21 by 0.5; output; end; run; * Creating the Cartesian product of times by covariates; proc sql noprint ; create table _xmatrix2 as select . as cd4, . as obstime, 0 as death, * from _times, _xmatrix; quit; data _xmatrix3; set _xmatrix2; _record=_n_; run; proc sql noprint; select max(_record) into : _numrecords from _xmatrix3; quit; * We obtain max and min event time from the data and store them as two macro variables; proc sql noprint; select min(time), max(time) into :lowb, :uppb from mydata.aids where death=1; quit; * We obtain the knots used in the time to event model from a dataset we had saved above, and store them as macro variables; proc sql noprint; select * into :kn1,:kn2,:kn3 from mydata.Example1_CALC_KNOTS_EVENT_MODEL; quit; * For each record of _matrix3 (each time T) we fit the model over this single observation; * Note the use of tech=none since this calll has no intention to maximize the likelihood function; * We just want to obtain the marginal -2LL of that single datapoint; %macro loop; %do index = 1 %to &_numrecords.; ods _all_ close; %JM(Data = _xmatrix3, Where= _record =&index., SubjectVar = _record, LongiTimeModel = linear, LongiVar = cd4, LongiTimevar = obstime, LongiTimeInter = ddi, LongiCovariates = female, LongiGMatrix = un, EventTimeVar = time, EventVar = death, EventVal = 1, EventModel = rpnaturalcubic, EventCovariates = ddi female, EventNknots = 3, Eventknot1 = &kn1., Eventknot2 = &kn2., Eventknot3 = &kn3., EventLowerknot = &lowb., EventUpperknot = &uppb., InitialParameters = mydata.Example1A_param_0, NLMIXEDOptions = tech=none qpoints=1, OutputParameters= _JMSurvivalEstimate, SharedParam = current_value slope, AdditionalOptions = untouchedlisting fitstatistics skipmacroheader); * As explained before the expected population survival at this time is obtained as exp(-2LL/-2); data _JMSurvivalEstimate_fit2 (keep= _record cum_haz surv ); set _JMSurvivalEstimate_fit; _record=&index.; if descr="-2 Log Likelihood"; cum_haz=value/2; surv=exp(-cum_haz); run; %if &index.=1 %then %do; data _JMSurvivalEstimate_fit3; set _JMSurvivalEstimate_fit2; run; %end; %else %if &index.>1 %then %do; data _JMSurvivalEstimate_fit3; set _JMSurvivalEstimate_fit3 _JMSurvivalEstimate_fit2; run; %end; %end; %mend loop; %loop; data _JMSurvivalEstimate_fit4; merge _xmatrix3 _JMSurvivalEstimate_fit3; by _record; run; *************************************************************************; * The code below produces Figure 3: Figure 3: Example 1, diagnostic plots; *************************************************************************; ods _all_ close; options nonumber orientation=landscape papersize=(6in 4in); goptions device=ACTXIMG; ods graphics / border=off width=1200px height=800px; ods pdf file="&folder./Diagnostic.pdf" style=printer nopdfnote ; ods layout start ; ods region x=0 y=0 width=50% height=50%; proc sgplot data= _jmsbjmean2 noautolegend ; scatter x=pred y=resid / markerattrs=(size=4); yaxis integer label="Residuals"; xaxis integer label="Empirical Bayes estimates"; run; ods region x=50% y=0 width=50% height=50%; proc sgplot data=Quant noautolegend; scatter x=resid_quant y=resid / markerattrs=(size=4); xaxis label="Theoretical normal quantiles" values=(-3.5 to 3.5 by 1); yaxis label="Standardized residuals" values=(-4 to 4 by 1); series x=resid y=resid ; run; ods region x=0 y=50% width=50% height=50%; proc sgplot data=_JMSurvivalEstimate_fit4; xaxis label="Time (months)"; yaxis label="Marginal survival" values=(0 to 1 by 0.2); series x=time y=surv /group=cov_l lineattrs=(color=black); run; ods region x=50% y=50% width=50% height=50%; proc sgplot data=_JMSurvivalEstimate_fit4; xaxis label="Time (months)"; yaxis label="Marginal cumulative hazard" values=(0 to 1 by 0.2); series x=time y=cum_haz /group=cov_l lineattrs=(color=black); run; ods layout end; ods pdf close; *********************************************************************************************************************************; * These additional calls are to fit the same model of Section 4.2 with additional quadrature points (feeds manuscript's figure 4); *********************************************************************************************************************************; %model2(n=1, gnl= method=quad(qpoints=3), nl=qpoints=3 gconv=0, inp = disjoint); %model2(n=2, nl=qpoints=5 gconv=0, inp = mydata.Example2_param_1); %model2(n=3, nl=qpoints=7 gconv=0, inp = mydata.Example2_param_2); %model2(n=4, nl=qpoints=9 gconv=0, inp = mydata.Example2_param_3); %model2(n=5, nl=qpoints=11 gconv=0, inp = mydata.Example2_param_4); %model2(n=6, nl=qpoints=13 gconv=0, inp = mydata.Example2_param_5); %model2(n=7, nl=qpoints=15 gconv=0, inp = mydata.Example2_param_6); %model2(n=8, nl=qpoints=17 gconv=0, inp = mydata.Example2_param_7); %model2(n=9, nl=qpoints=19 gconv=0, inp = mydata.Example2_param_8); %model2(n=10, nl=qpoints=21 gconv=0, inp = mydata.Example2_param_9); ********************************************************************************************************; * The code below produces Figure 4: Example 2, joint model estimates versus number of quadrature points ; ********************************************************************************************************; proc format; value $parf 'L_INTERCEPT' = 'Intercept' 'L_TIME_B1' = 'Slope' 'L_TIME_B1_BY_DDI' = 'Slope by DDI' 'L_FEMALE' = 'Female' 'L_LOGSD_INTERCEPT' = 'Log SD bi0' 'L_COV_B0_B1' = 'Covar(bi0, bi1)' 'L_LOGSD_TIME_B1' = 'Log SD bi1' 'L_LOGSD_RESIDUAL' = 'Log SD Residual' 'E_DDI' = 'DDI (Event)' 'E_FEMALE' = 'Female (Event)' 'E_TIME_B0' = 'Time Param 0' 'E_TIME_B1' = 'Time Param 1' 'E_TIME_B2' = 'Time Param 2' 'E_TIME_B3' = 'Time Param 3' 'E_TIME_B4' = 'Time Param 4' 'ASSOCT_RE_BI0' = 'Association (bi0)' 'ASSOCT_RE_BI1' = 'Association (bi1)' ; run; data parametersE2; set mydata.Example2_param_1 (keep=parameter estimate in=a) mydata.Example2_param_2 (keep=parameter estimate in=b) mydata.Example2_param_3 (keep=parameter estimate in=c) mydata.Example2_param_4 (keep=parameter estimate in=d) mydata.Example2_param_5 (keep=parameter estimate in=e) mydata.Example2_param_6 (keep=parameter estimate in=f) mydata.Example2_param_7 (keep=parameter estimate in=g) mydata.Example2_param_8 (keep=parameter estimate in=h) mydata.Example2_param_9 (keep=parameter estimate in=i) mydata.Example2_param_10 (keep=parameter estimate in=j); if a then qpoints=3 ; else if b then qpoints=5; else if c then qpoints=7; else if d then qpoints=9; else if e then qpoints=11; else if f then qpoints=13; else if g then qpoints=15; else if h then qpoints=17; else if i then qpoints=19; else if j then qpoints=21; format parameter $parf.; abs=abs(estimate); run; proc sort data=parametersE2; by parameter qpoints; run; ods _all_ close; options nonumber orientation=landscape papersize=(6in 4in); goptions device=ACTXIMG; ods pdf file="&folder./Example2EstimatesVSQPoints.pdf" style=journal nopdfnote; ods graphics / border=off width=1200px height=800px; proc sgplot data=parametersE2 ; xaxis label="Number of quadrature points" values=(3 to 21 by 2); yaxis label="Abs(estimate) (log-scale)" type=log; series x=qpoints y=abs/group=parameter ; where parameter in ('ASSOCT_RE_BI0', 'ASSOCT_RE_BI1','L_LOGSD_RESIDUAL','L_LOGSD_INTERCEPT','L_COV_B0_B1','L_LOGSD_TIME_B1','E_DDI', 'E_FEMALE'); run; ods pdf close;