/***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;