Correlated multistate models for multiple processes: an application to renal disease progression in systemic lupus erythematosus

Summary Bidirectional changes over time in the estimated glomerular filtration rate and in urine protein content are of interest for the treatment and management of patients with lupus nephritis. Although these processes may be modelled by separate multistate models, the processes are likely to be correlated within patients. Motivated by the lupus nephritis application, we develop a new multistate modelling framework where subject-specific random effects are introduced to account for the correlations both between the processes and within patients over time. Models are fitted by using bespoke code in standard statistical software. A variety of forms for the random effects are introduced and evaluated by using the data from the Systemic Lupus International Collaborating Clinics.


R Code for the Maximisation of the Likelihood Function from Section 4
We present the R code that was used to evaluate the likelihood function for example presented in Section 5 of the manuscript. We assume that we have longitudinal data for a number of subjects across a number of ordered time points (t) such that the eGFR state (gfstate), proteinuria state (pustate) are recorded at each time point together with a unique patient identifier (ptno). Hence, the dataset is of the form given in Table 1 ptno t gfstate pustate 001 0.00 The commented R code for the calculation of the likelihood function for the correlated eGFR and proteinuria models is given below. We note that this code will fit the 'simple random effects' model from Section 5 of the paper. Slight modifications can be made to the code to fit the inverse, power inverse and separate random effects models from the manuscript.
We make the following definitions for use in the presented R code: • v= the log-transformed gamma random effect • data = longitudinal dataset containing information on eGFR state and proteinuria level state, at the patient level, over time. Similar to that shown in Table 1.
• Q GF = non-stochastic transition intensity matrix for eGFR.
• Q PU = non-stochastic transition intensity matrix for proteinuria level.
The function 'corr gf pu msm1'outputs patient-level likelihood function contribution as a function of the random effect, v, for user-supplied values of the random effect variance and transition intensity matrix parameters for eGFR and proteinuria level.
#Evaluate the patient-level contribution to the likelihood #function from the eGFR model, for a #given value of the random effect v. Here we use the msm #library to evaluate the likelihood contribution. The #outputted quantity from this step is-2*log-likelihood #at the patient level, for a given value of v.
mod_GF<-msm(gfstate~t,data=data,opt.method="optim",method="BFGS", control=list(maxit=0),qmatrix=Q1,hessian=FALSE,use.deriv=FALSE)$opt$value #Evaluate the patient-level contribution to the likelihood #function from the proteinuria level model, for a given value #of the random effect v. Here we use the 'msm' library to #evaluate the likelihood contribution. The outputted #quantity from this step is -2*log-likelihood at the patient #level, for a given value of v.
mod_PU<-msm(pustate~t,data=data,opt.method="optim",method="BFGS", control=list(maxit=0),qmatrix=Q2,hessian=FALSE,use.deriv=FALSE)$opt$value #We transform the -2*log-likelihood values from above to form #likelihood function contributions from each of the eGFR #and proteinuria level parts of the model, at the #patient level and for a given value of the random effect, v.
pdf<-(1/gamma(tau))*(1/v)*(tau^tau)*((-log(v))^(tau-1))*exp(-tau*(-log(v))) #Multiply this evaluated probability density function by #the patient-level likelihood function contributions from #each of the eGFR and proteinuria level #parts of the model. out<-lik_GF*lik_PU*pdf #Outputted below is a patient-level likelihood function #contribution as a function of the random effect, v, for #user-supplied values of the random effect variance and #transition intensity matrix parameters for eGFR #and proteinuria level.

return(out) }
The function below integrates takes the patient-level outputted likelihood function contribution from 'corr gf pu msm1' and integrates out the random effect.
• full data = a dataset containing longitudinal records of eGFR state and proteinuria state for all patients. Patients are indexed by the column 'ptno'.
• ind= patient number index. #Pass the patient-level data, eGFR #transition intensity matrix, proteinuria level transition #intensity matrix, random effect variance value and patient-level #data to 'corr_gf_pu_msm3' for the formation of the patient-level #contribution to the likelihood function and integration step. #Return this patient-level likelihood function contribution.
lik_cont<-corr_gf_pu_msm3(Q_GF,Q_PU,dat,theta) return(lik_cont) } The function 'corr gf pu msm5' evaluates each patient-level contribution to the likelihood function and outputs -2*log-likelihood function. The inputs are as follows • pts= list of patient ids (ptno) • full data= a dataset containing longitudinal records of eGFR state and proteinuria state for all patients. Patients are indexed by the column 'ptno'.
• pars= a 9x1 vector of real-valued user-defined numbers. These are the numbers that will be transformed to transition intensity matrix entries and the random effect variance value during the likelihood maximisation process.