***************************************************************************************************************************************************; * This code is to fit a parametric mixture model to characterize rate, duration and severity of exacerbations -- By Abdollah Safari, Jan. 2018 *; ***************************************************************************************************************************************************; /**** SAS Macro for the recurrent episode analysis of simulated data*** Please save the accompanying sim_data.csv (with 1109 observations) to your local hard drive and update the path in *FileDir* below. sim_data has 6 columns 1. id: subject id 2. period: exacerbation index number (0 means no exacerbation) 3. time: time to the next event since time of previous event (in year) 4. event: -1: out of exacerbation censoring event, -2: in exacerbation censoring event, 1: mild exacerbation onset event, 2: severe exacerbation onset event, 3: exacerbation termination event 5. var1: a simulated continuous covariate with a certain relation with submodels 6. var2: a simulated binary covariate with a certain relation with submodels The *RecrEpis* macro has 4 inputs: 1. FileDir: directory of the sim_data in your local computer 2. Dur_SubMod: whether the duration submodel should be included (YES or NO) 3. NoExac_SubMod: whether the non-susceptible submodel should be included (YES or NO) 4. Sev_SubMod: whether the severity submodel should be included (YES or NO) */ ODS HTML CLOSE; ODS HTML; %LET FileDir="C:\Users\asafari\Desktop\Asthma\R code\sim_data.csv"; *UPDATE THIS PATH TO WHERE THE CSV FILE IS SAVED; %MACRO RecrEpis(FileDir=&FileDir, Dur_SubMod="YES", NoExac_SubMod="YES", Sev_SubMod="YES"); PROC IMPORT OUT=WORK.SIM_DATA DATAFILE=&FileDir DBMS=CSV REPLACE; GETNAMES=YES; RUN; * Building the required parameters for the model and their associated likelihood function according to the selected submodels specified in the macro; * Coefficients starting with br are for the rate submodel, bp for the non-exacerbators component, bd for the duration submodel, and bs for the severity submodel; * v_r, v_d, and v_s are the variances of the random-effect of rate, duration, and severity submodels respectively; * v_rd, v_rs, and v_ds are the associated covariances of the random effects; %IF (&Dur_SubMod="YES" AND &NoExac_SubMod="YES" AND &Sev_SubMod="YES") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 gamma2=1 br_0=0 br_x1=0 br_x2=0 bp_0=0 bp_x1=0 bp_x2=0 bd_0=0 bd_x1=0 bd_x2=0 bs_0=0 bs_x1=0 bs_x2=0 v_r=1 v_d=1 v_s=1 v_rd=0 v_rs=0 v_ds=0); %LET BOUNDS_LIST= %STR(gamma1>0, gamma2>0, v_r>0, v_d>0, v_s>0); %LET RANDOM_LIST= %STR(Z_r Z_d Z_s ~ NORMAL([0 , 0 , 0] , [v_r , v_rd , v_d , v_rs , v_ds , v_s]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=P_NoEx+(1-P_NoEx)*S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r*(1-P_sev); IF (EVENT=2) THEN L=h_t_r*S_t_r*P_sev; IF (EVENT=3) THEN L=h_t_d*S_t_d; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=S_t_d; IF (PERIOD=1) THEN L=L*(1-P_NoEx); END;); %END; %IF (&Dur_SubMod="NO" AND &NoExac_SubMod="YES" AND &Sev_SubMod="YES") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 br_0=0 br_x1=0 br_x2=0 bp_0=0 bp_x1=0 bp_x2=0 bs_0=0 bs_x1=0 bs_x2=0 v_r=1 v_s=1 v_rs=0); %LET BOUNDS_LIST= %STR(v_r>0, v_s>0); %LET RANDOM_LIST= %STR(Z_r Z_s ~ NORMAL([0 , 0] , [v_r , v_rs , v_s]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=P_NoEx+(1-P_NoEx)*S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r*(1-P_sev); IF (EVENT=2) THEN L=h_t_r*S_t_r*P_sev; IF (EVENT=3) THEN L=1; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=1; IF (PERIOD=1) THEN L=L*(1-P_NoEx); END;); %END; %IF (&Dur_SubMod="YES" AND &NoExac_SubMod="NO" AND &Sev_SubMod="YES") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 gamma2=1 br_0=0 br_x1=0 br_x2=0 bd_0=0 bd_x1=0 bd_x2=0 bs_0=0 bs_x1=0 bs_x2=0 v_r=1 v_d=1 v_s=1 v_rd=0 v_rs=0 v_ds=0); %LET BOUNDS_LIST= %STR(v_r>0, v_d>0, v_s>0); %LET RANDOM_LIST= %STR(Z_r Z_d Z_s ~ NORMAL([0 , 0 , 0] , [v_r , v_rd , v_d , v_rs , v_ds , v_s]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r*(1-P_sev); IF (EVENT=2) THEN L=h_t_r*S_t_r*P_sev; IF (EVENT=3) THEN L=h_t_d*S_t_d; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=S_t_d; END;); %END; %IF (&Dur_SubMod="YES" AND &NoExac_SubMod="YES" AND &Sev_SubMod="NO") %THEN %DO; %LET PARMS_LIST= %STR( gamma1=1 gamma2=1 br_0=0 br_x1=0 br_x2=0 bp_0=0 bp_x1=0 bp_x2=0 bd_0=0 bd_x1=0 bd_x2=0 v_r=1 v_d=1 v_rd=0); %LET BOUNDS_LIST= %STR(v_r>0, v_d>0); %LET RANDOM_LIST= %STR(Z_r Z_d ~ NORMAL([0 , 0] , [v_r , v_rd , v_d]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=P_NoEx+(1-P_NoEx)*S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r; IF (EVENT=2) THEN L=h_t_r*S_t_r; IF (EVENT=3) THEN L=h_t_d*S_t_d; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=S_t_d; IF (PERIOD=1) THEN L=L*(1-P_NoEx); END;); %END; %IF (&Dur_SubMod="NO" AND &NoExac_SubMod="NO" AND &Sev_SubMod="YES") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 br_0=0 br_x1=0 br_x2=0 bs_0=0 bs_x1=0 bs_x2=0 v_r=1 v_s=1 v_rs=0); %LET BOUNDS_LIST= %STR(v_r>0, v_s>0); %LET RANDOM_LIST= %STR(Z_r Z_s ~ NORMAL([0 , 0] , [v_r , v_rs , v_s]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r*(1-P_sev); IF (EVENT=2) THEN L=h_t_r*S_t_r*P_sev; IF (EVENT=3) THEN L=1; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=1; END;); %END; %IF (&Dur_SubMod="NO" AND &NoExac_SubMod="YES" AND &Sev_SubMod="NO") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 br_0=0 br_x1=0 br_x2=0 bp_0=0 bp_x1=0 bp_x2=0 v_r=1); %LET BOUNDS_LIST= %STR(v_r>0); %LET RANDOM_LIST= %STR(Z_r ~ NORMAL([0] , [v_r]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=P_NoEx+(1-P_NoEx)*S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r; IF (EVENT=2) THEN L=h_t_r*S_t_r; IF (EVENT=3) THEN L=1; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=1; IF (PERIOD=1) THEN L=L*(1-P_NoEx); END;); %END; %IF (&Dur_SubMod="YES" AND &NoExac_SubMod="NO" AND &Sev_SubMod="NO") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 gamma2=1 br_0=0 br_x1=0 br_x2=0 bd_0=0 bd_x1=0 bd_x2=0 v_r=1 v_d=1 v_rd=0); %LET BOUNDS_LIST= %STR(v_r>0, v_d>0); %LET RANDOM_LIST= %STR(Z_r Z_d ~ NORMAL([0 , 0] , [v_r , v_rd , v_d]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r; IF (EVENT=2) THEN L=h_t_r*S_t_r; IF (EVENT=3) THEN L=h_t_d*S_t_d; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=S_t_d; END;); %END; %IF (&Dur_SubMod="NO" AND &NoExac_SubMod="NO" AND &Sev_SubMod="NO") %THEN %DO; %LET PARMS_LIST=%STR( gamma1=1 br_0=0 br_x1=0 br_x2=0 v_r=1); %LET BOUNDS_LIST= %STR(v_r>0); %LET RANDOM_LIST= %STR(Z_r ~ NORMAL([0] , [v_r]) SUBJECT=id); %LET LIKELIHOOD= %STR( IF (PERIOD=0) THEN L=S_t_r; ELSE DO; IF (EVENT=1) THEN L=h_t_r*S_t_r; IF (EVENT=2) THEN L=h_t_r*S_t_r; IF (EVENT=3) THEN L=1; IF (EVENT=-1) THEN L=S_t_r; IF (EVENT=-2) THEN L=1; END;); %END; PROC NLMIXED DATA=SIM_DATA; * PARMS statement declares model parameters and assign their initial values; PARM &PARMS_LIST; BOUNDS &BOUNDS_LIST; *** RATE SUBMODEL ***; * lin_rate: linear function of the log-normal distribution scale parameter for the rate submodel; lin_rate = br_0 + br_x1*var1 + br_x2*var2 + Z_r; * Z_r is the random effect term for the rate submodel; * Survival and hazard function of log-normal distribution for the rate submodel; S_t_r = 1 - CDF('NORMAL', (LOG(time) - lin_rate)/gamma1); h_t_r = EXP(-0.5*((LOG(time) - lin_rate)/gamma1)**2)/SQRT(2*CONSTANT('pi'))/gamma1/(time)/S_t_r; *** NON-EXACERBATOR COMPONENT ***; * lin_NoEx: logit function of the logistic regression for the non-exacerbator component; %IF (&NoExac_SubMod="YES") %THEN %DO; lin_NoEx = bp_0 + bp_x1*var1 + bp_x2*var2; P_NoEx = EXP(lin_NoEx)/(1 + EXP(lin_NoEx)); * probability of being non-exacerbator; %END; *** DURATION SUBMODEL ***; * lin_dur: linear function of the log-normal distribution scale parameter for the duration submodel; %IF (&Dur_SubMod="YES") %THEN %DO; lin_dur = bd_0 + bd_x1*var1 + bd_x2*var2 + Z_d; * Z_d is the random effect term for the duration submodel; * Survival and hazard function of log-normal distribution for the duration submodel; S_t_d = 1 - CDF('NORMAL', (LOG(time) - lin_dur)/gamma2); h_t_d = EXP(-0.5*((LOG(time) - lin_dur)/gamma2)**2)/SQRT(2*CONSTANT('pi'))/gamma2/(time)/S_t_d; %END; *** SEVERITY SUBMODEL ***; * lin_sev: logit function of the logistic regression for the severity submodel; %IF (&Sev_SubMod="YES") %THEN %DO; lin_sev = bs_0 + bs_x1*var1 + bs_x2*var2 + Z_s; * Z_s is the random effect term for the severity submodel; P_sev = EXP(lin_sev)/(1 + EXP(lin_sev)); * probability of having a severe exacerbation, conditional on having an exacerbation; %END; *** LIKELIHOOD CALCULATION ***; &LIKELIHOOD; * We use "GENERAL" as we have constructed the likelihood ourselves; * It tells SAS to maximize ll and do not worry about its type ets; ll = LOG(L); MODEL event ~ GENERAL(ll); RANDOM &RANDOM_LIST; *Estimate uses Delta method for inference about functionals of model parameters; ESTIMATE "Z_r_SD" SQRT(v_r); %IF (&Dur_SubMod="YES") %THEN %DO; ESTIMATE "Z_d_SD" SQRT(v_d); ESTIMATE "Correlation coefficient (Zr_Zd_RHO)" v_rd/SQRT(v_r*v_d); %END; %IF (&Sev_SubMod="YES") %THEN %DO; ESTIMATE "Z_s_SD" SQRT(v_s); ESTIMATE "Correlation coefficient (Zr_Zs_RHO)" v_rs/SQRT(v_r*v_s); %END; %IF (&Sev_SubMod="YES" AND &Dur_SubMod="YES") %THEN %DO; ESTIMATE "Correlation coefficient (Zd_Zs_RHO)" v_ds/SQRT(v_d*v_s); %END; RUN; %MEND; %RecrEpis(FileDir=&FileDir, Dur_SubMod="YES", NoExac_SubMod="YES", Sev_SubMod="YES")