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

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

* *

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