Joint modelling of multivariate longitudinal clinical laboratory safety outcomes, concomitant medication and clinical adverse events: application to artemisinin-based treatment during pregnancy clinical trial

Background: In drug trials, clinical adverse events (AEs), concomitant medication and laboratory safety outcomes are repeatedly collected to support drug safety evidence. Despite the potential correlation of these outcomes, they are typically analysed separately, potentially leading to misinformation and inefficient estimates due to partial assessment of safety data. Using joint modelling, we investigated whether clinical AEs vary by treatment and how laboratory outcomes (alanine amino-transferase, total bilirubin) and concomitant medication are associated with clinical AEs over time following artemisinin-based antimalarial therapy. Methods: We used data from a trial of artemisinin-based treatments for malaria during pregnancy that randomized 870 women to receive artemether–lumefantrine (AL), amodiaquine–artesunate (ASAQ) and dihydroartemisinin–piperaquine (DHAPQ). We fitted a joint model containing four sub-models from four outcomes: longitudinal sub-model for alanine aminotransferase, longitudinal sub-model for total bilirubin, Poisson sub-model for concomitant medication and Poisson sub-model for clinical AEs. Since the clinical AEs was our primary outcome, the longitudinal sub-models and concomitant medication sub-model were linked to the clinical AEs sub-model via current value and random effects association structures respectively. We fitted a conventional Poisson model for clinical AEs to assess if the effect of treatment on clinical AEs (i.e. incidence rate ratio (IRR)) estimates differed between the conventional Poisson and the joint models, where AL was reference treatment. Results: Out of the 870 women, 564 (65%) experienced at least one AE. Using joint model, AEs were associated with the concomitant medication (log IRR 1.7487; 95% CI: 1.5471, 1.9503; p < 0.001) but not the total bilirubin (log IRR: -0.0288; 95% CI: − 0.5045, 0.4469; p = 0.906) and alanine aminotransferase (log IRR: 0.1153; 95% CI: − 0.0889, 0.3194; p = 0.269). The Poisson model underestimated the effects of treatment on AE incidence such that log IRR for ASAQ was 0.2118 (95% CI: 0.0082, 0.4154; p = 0.041) for joint model compared to 0.1838 (95% CI: 0.0574, 0.3102; p = 0.004) for Poisson model. Conclusion: We demonstrated that although the AEs did not vary across the treatments, the joint model yielded efficient AE incidence estimates compared to the Poisson model. The joint model showed a positive relationship between the AEs and concomitant medication but not with laboratory outcomes. Trial registration: ClinicalTrials.gov: NCT00852423.

Our proposed joint model was applied to data from PREGACT trial as the motivating data. PREGACT trial (ClinicalTrials.gov number, {"type":"clinical-trial","attrs":{"text":"NCT00852423","term_id":"NCT00852423"}}NCT00852423) was a multicentre, open-label randomized trial carried out in 4 African countries (Burkina Faso, Ghana, Malawi, and Zambia). The current work uses the data from Malawi. In Malawi, the trial was implemented between June 2010 and August 2013 and enrolled 870 pregnant women during their second or third trimester with falciparum malaria. The women were randomly allocated in a 1:1:1 ratio to be treated with artemether–lumefantrine (AL), amodiaquine–artesunate (ASAQ) or dihydroartemisinin–diperaquine (DHAPQ). The primary outcome was polymerase-chain reaction (PCR) adjusted cure rates at day 63. The safety outcomes collected at baseline and during follow-up were AEs, total bilirubin, alanine aminotransferase, white blood cell and red blood cell counts. The patients were directly observed on days 0-2 (after receiving the dose). The patients were then asked to return to the clinic for follow-up visits on days 3 and 7 and then weekly thereafter until day 63. For this analysis, we focus on four outcomes; two clinical laboratory safety biomarkers (total bilirubin and alanine aminotransferase), clinical AEs and concomitant medication. The total bilirubin and alanine aminotransferase data was collected at enrolment, day 7, 14, 28 and 63. We focused on bilirubin and alanine aminotransferase biochemical parameters since they are key biochemical parameters in antimalarial drug safety assessment during pregnancy and were also measured in the PREGACT trial. The clinical AE count, defined as cumulate number of AEs experienced by the end of follow-up time is the primary outcome of interest in the current analysis. In all the subsequent discussions, the clinical AEs that were defined as definitely not related (by the study physician) to the antimalarial drug treatment are not considered in developing the joint model to avoid spurious results. The concomitant medication outcome is also defined as cumulative number of reported concomitant medication use by the end of the follow up time for each patient. In the current study, we considered the reported concomitant medication use regardless of its intended use. The PREGACT trial was conducted in accordance with the Declaration of Helsinki and Good Clinical Practice guidelines. The trial obtained ethical clearance from ethics committee at the Antwerp University Hospital in Belgium and College of Medicine Research Ethics Committee at the University of Malawi [22]. Prior to enrolment, informed consent was also sought from the mother. Ethical approval was also obtained from University of the Witwatersrand Human Research Ethics Committee, prior to access and utilization of the data for the current analysis. Let each patient who was randomized and received at least a dose be denoted as i = 1, … …, n. Let Vi=vi1T,vi2T be a bivariate continuous clinical laboratory outcome vector. Specifically, in our context we consider a bivariate scenario; k = 1 is alanine amino transferase and k = 2 is total bilirubin. Each of the two continuous outcome vectors (vi1,  vi2) are of (nikx1) dimension for the observed longitudinal measurements of the k-th outcome; vik = (vi1k, …….,  vink)T. We accommodate the situations where observation times, tijk may differ between individuals and outcomes (j = 1,  … . ., nik). For each patient, we let Ci=ci1T,ci2T represent a bivariate count outcome vector where k = 1 is for total number of AEs experienced during their follow-up time Ti and k = 2 for total number concomitant medication used over the follow-up time Ti. A set of covariates that were collected at baseline for each individual are defined as Xi = {X1i, X2i, … . ., Xpi}. In drug trials, where multiple safety outcomes are repeatedly collected to support safety evidence, clinical AEs are usually of primary interest. The occurrence of the clinical AEs is typically correlated with other safety-related outcomes such as alanine aminotransferase, total bilirubin and concomitant medication. Modelling these outcomes separately is inefficient since it insufficiently accounts for the correlations. Traditionally, modelling of clinical AEs count is done using the conventional Poisson model, where the Poisson mean response (of clinical AEs in this context), φ, is assumed to be equal to the variance Since the clinical AEs count is non-negative the logarithm of φ can naturally be linked with the baseline variables such that the conventional Poisson model for the clinical AE count can be presented as; Where φ represents Poisson mean response and β is a vector of coefficients (i.e. log IRR) for the corresponding vector of baseline covariates, XikT, that included maternal age, gravidity, trimester at enrolment and treatment arm.. The model M1 has log link of the mean response of clinical AEs count and assumes a Poisson distribution of the clinical AEs. The coefficients of the covariates are interpreted on logarithm scale as expected change in the log of mean clinical AEs count per unit change in the covariate. Exponentiating the coefficients (i.e. log IRR) yields IRR. However, the estimated effect of treatment on clinical AEs from model M1 may be insufficient since it does not account for other factors that can confound the estimates e.g. concomitant medication and clinical laboratory outcomes. Alternatively, one can consider incorporating this additional information (of concomitant medication and clinical laboratory outcomes) as part of covariates in model M1. This can yield model M2 below. Since the concomitant medication and clinical laboratory outcomes (alanine aminotransferase, total bilirubin) are added in the model as time-varying covariates (with time t0, representing baseline covariates in M1), the vector of covariates adjusted in model M2 is XikTt denoting the vector of covariates value at time t. However, both model M1 and M2 do not account for correlation overtime across the longitudinal alanine aminotransferase, longitudinal total bilirubin, concomitant medication and clinical AEs outcomes. In this context of multiple outcomes that are correlated overtime, joint modelling offers an opportunity to obtain improved and efficient estimates since it can efficiently account for the potential correlations and baseline covariates. The joint model also enables formal quantification of the relationship between correlated outcomes by estimating the strength of the association between the outcomes. Most previously proposed joint models consist of two sub-models; longitudinal sub-model and time-to-event sub-model (where the time to event outcome is of primary interest). Here, we propose a joint model where the primary outcome of interest is count outcome (i.e. clinical AEs count) such that time to event outcome sub-model is replaced with count outcome sub-model. Since we assumed Poisson distribution of the clinical AEs count, the count outcome sub-model was specified as Poisson-distributed. Therefore, the joint model presented in this paper consist of two longitudinal continuous laboratory safety biomarkers sub-models (one for total bilirubin and one for alanine aminotransferase) two Poisson sub-models (one for the individual clinical AEs count and one for the individual concomitant medication count). The connecting of these outcomes sub-models is very flexible such that we could use random effects or expected value of outcomes [23] as detailed in the model specification below. As highlighted above, the joint model is formulated in such a way that clinical AE count is a primary outcome. In this section we describe the structures of the sub-models that made the joint model considered in this work. The estimates in all the models are considered on a logarithm scale. We modelled two continuous longitudinal clinical laboratory safety outcomes (i.e. total bilirubin and alanine aminotransferase) using a longitudinal sub-model based on a flexible linear mixed effects model that could accommodate nonlinear changes of the laboratory safety outcomes. The flexibility was achieved through the use of restricted cubic splines. Each continuous longitudinal clinical laboratory outcome vi1, vi2 was assumed to be normally distributed with mean μk and variance δk2 (k = 1, 2) on a logarithm scale. Each continuous longitudinal clinical laboratory outcome is modelled using the sub-model below where εik(t) is the error term observed at time t for the kth outcome from patient i. The error terms for the model are assumed to be independent, identical and normally distributed with mean 0. The mean response model XikT(t)β, specified as a linear function of the covariates at a given time for the outcome k is flexible such that it can accommodate time-varying covariates. The Zi(t) is an indicator vector of random effects for kth outcome for patient i at time t such that it takes the value of 1 when there is a random effect and 0 otherwise. The vector of the patient-specific shared random effects bi is assumed to have a multivariate normal distribution with mean 0 and variance-covariance matrix Σ, s.t. (bi) ∼ MVN(0, ∑). Based on M1, considering patient-specific random intercept as our shared parameter of interest, linking the outcomes, the longitudinal sub-model for alanine aminotransferase can be written simply as M4 below; where b0i1 is the patient specific random intercept corresponding to the alanine aminotransferase outcome. Similarly, the longitudinal sub-model for total bilirubin, M5, can be formulated as; where b0i2 is the patient specific random intercept corresponding to total bilirubin outcome. We considered two event count sub-models; one for concomitant medication count and the other for the AE count. The two count sub-models were derived from a modified form of the generic Poisson model. The count sub-model was as follows; Assuming that the concomitant medication count had a Poisson distribution, its sub-model was defined as such that log(φ| b03i) represents the logarithm of the Poisson mean response, φ, conditional on the patient-specific random intercept, b03i for concomitant medication count (i.e. the concomitant medication count sub-model is linked with the AE count sub-model in the joint model via the random intercept). The XikT(t)β is the linear predictor for the Poisson model has a logarithm link function and contains a vector of fixed effects β corresponding to the baseline covariates. The α3 represents log incidence rate ratio per unit increase in the patient-specific deviation from the mean random intercept of the reported concomitant medication use. Considering that AE count was the primary outcome of interest and assuming a Poisson distribution for the AE count, the AE count sub-model was formulated as; The log(φ| b01i, b02i, b03i) represents the logarithm of Poisson mean response, φ, conditional on the shared random intercepts for the two longitudinal clinical laboratory outcomes and the concomitant medication count. Each random intercept for the longitudinal sub-models and concomitant medication sub-model (b01i, b02i, b03i) is linked to the linear predictor XikT(t)β for the Poisson model. This yields random effects association structure where α1, α2, α3 estimate the strength of the association between the respective continuous longitudinal clinical laboratory outcome, concomitant medication count and the AE count; the α1, α2, α3 represent log incidence rate per unit increase in the patient-specific deviation from the mean random intercept of a respective continuous longitudinal clinical laboratory outcome or the concomitant medication count outcome. For example, α1 represents log incidence rate ratio per unit increase in the patient-specific deviation from the mean random intercept of a respective log alanine aminotransferase. Alternative clinically meaningful formulation of the AE count sub-model could be linked to the expected value of the respective continuous longitudinal clinical laboratory outcome with the linear predictor of the Poisson model. This yields the current value association structure where the α1, α2, α3 represent log incidence rate ratio per unit increase of the respective continuous longitudinal clinical laboratory outcome, at time t. Current value association structure is very important and clinically plausible when linking continuous outcomes. In our cases, given the two continuous longitudinal clinical laboratory outcomes (i.e. log alanine aminotransferase and log total bilirubin) v1, v2, under the current value association structure, we can modify the AE count sub-model as; The E[v1(t)| μ1(t)] and E[v2(t)| μ2(t)]α2 are the expected current values of log alanine aminotransferase and total bilirubin respectively. It can be noted from M6 that the concomitant medication count, v3, sub-model is still linked to the AE count sub-model via a random intercept b03i. We maintained the random intercept (b03iα3) parameterization since this efficiently deals with any potential over-dispersion problem encountered in Poisson models. Therefore, model M8 contains a mixture of both the expected current value and random effects association structures. Since both continuous outcomes (the log alanine aminotransferase and log total bilirubin) are assumed to be normally distributed (i.e. with identity link), the model can further be simplified as; In this paper, we focus on reporting results from the M9 association structure since it is more clinically meaningful since it has association parameters that are directly interpretable in clinical practice. The estimates for the joint model are obtained by maximizing the marginal likelihood of the joint distribution of the observed data and the random effects [23]. The likelihood function for the observed data can be specified as; where θ represents the set of parameters of interest to be estimated including both the fixed and random effects (e.g. association parameters α1, α2, α3, and coefficients for the baseline covariates, β). Component f(bi| θ) is a normal density function for the random effects conditional on θ and f(Vi| θ) is a normal density function for the log total bilirubin and log alanine aminotransferase conditional on θ and bi. The f(Ti, Ci| bi, θ) represents Poisson density function conditional on θ and bi; patient-specific total follow up days is Ti and Ci represents total events recorded (for AE count sub-model and concomitant count sub-model). Since obtaining the log likelihood requires integrating out the random effects, it becomes challenging as number of random effects increases such that the traditionally used adaptive Gauss-Hermite quadrature would be limited (to handle the large number of random effects in our complex joint model). To mitigate this problem, we employed Monte Carlo integration technique since the number of draws it makes from the random effects do not need to change with increase in number of random effects [23, 24] (hence reducing the computation burden). Estimation of our models was done using merlin in Stata [21]. Only 7 patients had missing data points and the data missing mechanism was considered completely at random. Firstly, we computed the summary of the baseline characteristics for the women at enrolment. We presented means with respective standard deviation for the continuous variables (e.g. maternal age, gestational age) and frequencies with respective percentages for the categorical variables (e.g. trimester at enrolment, bed-net use). These summaries were also computed for the clinical AE occurrence and concomitant medication counts. The concomitant medication was also summarized based frequency of the specific patients and the number of patients who took the medication. In order to aid visualization of the trend of the data, we plotted box plots for the two continuous longitudinal clinical laboratory data (total bilirubin and alanine aminotransferase) outcomes by study visit day. We considered two different models to compare how they efficiently estimated effect of treatment on clinical AEs over the follow-up time. The models compared were conventional Poisson model M1 and joint model of M9 formulation above. Both the models adjusted for the same set of baseline covariates (maternal age, gravidity and trimester at enrolment). The traditional Poisson model was the first model to be fitted and we used logarithm of the total follow-up time for each patient as an offset. Then we fitted a joint model with a mixture of expected current value of outcome and random effects association structure as shown in M9. We considered the joint model to assess the joint evolution (/association) between the clinical AEs and other three outcomes (alanine aminotransferase, total bilirubin and concomitant medication) where α1, α2, α3 as defined in M9 are the parameters quantifying the between-outcome association. Secondly, we were also interested investigate the efficiency of the joint model in improving AE incidence rate estimate of the treatment effect obtained from the conventional Poisson model M1. As part of sensitivity analysis we also investigated how the association structure of the joint model affects the magnitude of the treatment effect estimates. Comparing the association structures was achieved by fitting a joint model with random effects association structure only as specified in model M7 above and joint model with mixture of current value and random effects association structure of joint model M9. For both joint models, in order to flexibly model nonlinear changes of the continuous longitudinal clinical laboratory outcomes we used restricted cubic splines with 3 degrees of freedom considered sufficient to enhance flexibility of the model. We compared the fit of the two joint models (M9 versus M7) using the Akaike Information Criterion (AIC), with small values of AIC suggesting better model. Additionally, we did an exploratory analysis to identify how the frequently reported concomitant medications are associated with the AEs. Such exploratory analysis was helpful in understanding how specific frequently used concomitant medication influence the overall impact of concomitant medication on AEs. Concomitant medication was defined as frequent if it constituted at least 10% of the reported medications. A joint model of a similar structure to M9 but where the concomitant medication count was confined to the specific frequently reported concomitant medication (i.e. paracetamol) was used during the exploratory analysis.

N/A

• Dima Health