Joint rate-severity modeling of (recurrent) events

This page contains the SAS code and an exemplary simulated dataset for fitting joint rate-severity models for recurrent event, as described here.

 

Link to SAS code

 

Link to exemplary data set

 

 

/***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.csvDBMS=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;

 

 

Note: parameters that have been used to generate the simulated dataset:

Follow-up time ~ Normal(1,0.1)    //average follow-up time is 1 year in each person

x~Normal(0,1) //the single covariate modeled in this example

Shape (gamma parameter of the Weibull): 1.2

br_0=0 (log-rate of event for those with covariate=0)

br_x=-1 (log-acceleration factor for the covariate. The log-hazard ratio associating one unit increase in x to the rate of outcome is b_x*gamma=-1.2)

bs_0=1 (logit of probability of severe event for those with covariate=0)

bs_x=1 (log-odds ratio for the association with one unit increase in covariate with the probability of severe v. mild event)

zr_sd (SD of random-effect term for rate)=0.5

zs_sd (SD of random-effect term for severity)=1

zr_zs_rho: Correlation coefficient between rate and severity effects=-0.5