/***SAS code for the analysis of
simulated data****
Please save the accompanying
sim_data.csv (with 10,015 observations) to your local hard drive and update the
path in the PROC IMPORT below.
sim_data has four
columns
1. id: subject id
2. tte0: time to event from baseline (not
from previous event as this is total time analysis)
3. event: 0:censoring event, 1:mild event,
2: severe event
4. covariate: a simulated covariate with a
certain relation with rate and severity (see below)
*/;
ODS HTML CLOSE;
ODS HTML;
PROC IMPORT OUT= WORK.sim_data
DATAFILE= “C:\sim_data.csv” DBMS=CSV
REPLACE; *UPDATE THIS PATH
TO WHERE THE CSV FILE IS SAVED;
GETNAMES=YES;
DATAROW=2;
RUN;
PROC SORT DATA=sim_data; BY id tte0; RUN;
PROC NLMIXED DATA=sim_data gconv=0;
*Coefficients starting with br
are for the rate component;
*those starting with bs are for the severity component;
*v_r is the variance of
random-effects for rate, and v_s is the same for the
severity component;
*PARMS statement declares model parameters and
assigns their initial values;
PARMS
gamma=1
br_0=0 br_x=0
bs_0=0 bs_x=0
v_r=1 v_s=1 cov_r_s=0;
BOUNDS gamma>0,v_r>=0,v_s>=0;
***RATE component***;
*lin_rate: the linear
component of the rate function;
lin_rate = br_0+br_x*x+zr; *zr is the
random-effect term for rate;
alpha = exp(lin_rate);
*Survival and hazard functions of Weibull;
S_t=exp(-(alpha*tte0)**gamma);
h_t = gamma*alpha*((tte0*alpha)**(gamma-1));
*ll1: the log-likelihood for the rate component;
*note the total
time analysis means each event contributes its hazard function,
*and there is a
survival-type contribution for the censoring event;
ll1 =
(event>0)*log(h_t)+(event=0)*log(S_t);
***SEVERITY component***;
*lin_severity: the
linear component of the severity function;
lin_severity=bs_0+bs_x*x+zs; *zs is the
random-effect term for severity;
*p is the probability of having severe
exacerbation, conditional on having an exacerbation;
p=EXP(lin_severity)/(1+EXP(lin_severity));
*ll2: the log-likelihood for the severity
component;
*If for the censoring even (coded as event=0),
there is no contribution to the severity component;
ll2=LOG((event=0)*1+(event>0)*((event=1)*(1-p)+(event=2)*p));
***Likelihood calculations***;
ll=ll1+ll2;
*We use ‘general’ as we have constructed the
likelihood ourselves;
*It tells SAS
just to maximize ll and do not worry about its type etc;
model tte0 ~ general(ll);
random zr zs
~ normal([0,0],[v_r,cov_r_s,v_s]) subject=id;
*Estimate uses Delta method for inference about
functionals of model parameters;
ESTIMATE ‘zr_sd‘ SQRT(v_r);
ESTIMATE ‘zs_sd‘ SQRT(v_s);
ESTIMATE ‘Correlation coefficient (zr_zs_rho)’ cov_r_s/SQRT(v_r*v_s);
RUN;
*End of code;