 Research
 Open Access
 Published:
Joint modelling of multivariate longitudinal clinical laboratory safety outcomes, concomitant medication and clinical adverse events: application to artemisininbased treatment during pregnancy clinical trial
BMC Medical Research Methodology volume 21, Article number: 208 (2021)
Abstract
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 aminotransferase, total bilirubin) and concomitant medication are associated with clinical AEs over time following artemisininbased antimalarial therapy.
Methods
We used data from a trial of artemisininbased 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 submodels from four outcomes: longitudinal submodel for alanine aminotransferase, longitudinal submodel for total bilirubin, Poisson submodel for concomitant medication and Poisson submodel for clinical AEs. Since the clinical AEs was our primary outcome, the longitudinal submodels and concomitant medication submodel were linked to the clinical AEs submodel 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
Background
There has been limited analysis of safety data being conducted in drug clinical trials compared to efficacy data. During drug clinical trials, scheduled and nonscheduled visits are conducted where multiple safety outcomes are collected including longitudinal clinical laboratory outcomes and clinical AEs [1, 2]. Data on concomitant medications defined as any medication/supplement taken by the trial participant other than the investigational product, are also repeatedly collected in order to support drug safety evidence. The repeatedly measured clinical laboratory outcomes are considered as objective and proxy indicators for druginduced human body organ dysfunction and tend to be strongly associated with incidence of clinical AEs. For example, classic liver function biomarkers such as total bilirubin and alanine aminotransferase are usually used in assessment of drug safety [3]. However, each of the clinical laboratory safety outcomes is usually analysed separately based on summary statistics such as means or proportion of patients with elevated values for the safety markers [4, 5]. The clinical AE data is also usually analysed separately using the descriptive nonparametric methods (e.g. cumulative incidence function and crude incidence rate) [6,7,8]. Although ignored in most safety analyses, concomitant medication taken during the trial period can also influence the magnitude of drug safety estimates since the concomitant medication can interact with the treatment under investigation. Concomitant medication is usually also analysed separately (e.g. using mean cumulative function) [9]. The standard approaches that analysed each of the clinical laboratory safety data, clinical AEs and concomitant medication separately are inefficient and may yield biased estimates of AE incidence rate (i.e. clinical AEs occurrence over followup time) since they ignore the correlation structure of the multiple outcomes over the followup time, leading to information loss. Furthermore, the separate analyses present challenges in interpreting overall drug safety profile. Integrated approach in analysis of the clinical laboratory safety data, concomitant medication and clinical AEs through the use of joint models, can offer an opportunity to efficiently harness the available information towards enriching the drug safety estimates. The joint models can efficiently quantify how the incidence rate of AEs is associated with the multivariate longitudinal process of the multiple continuous laboratorybased measurements, accounting for baseline covariates, concomitant medications and unobserved heterogeneity that is not captured in the conventional analytical approaches of safety data, all of which can potentially confound the estimates.
Recently, there have been increased proposals for improved analysis of multivariate clinical laboratory data, concomitant medications and clinical AEs. Schildrout et al. [10] and Rosenkranz [11] discuss the utility of longitudinal modelling of clinical laboratory data using estimating equations and maximum likelihood. However, such separate longitudinal models ignore correlation between the multiple safety outcomes. Southworth and Heffernan proposed joint modelling approach that captures the joint behaviour of clinical laboratory safety data elevated values indicative of potential liver damage [12]. However, this was limited as the model did not use all the available information due to focus on extreme values and, further, such modelling is mainly suited for large sample sizes. Merz et al. [3] also recommended accounting for both multivariate longitudinal clinical laboratory data and clinical AEs, although they do not explicitly discuss how this can be practically achieved. Barker [9] demonstrates how mean cumulative function can be used as an exploratory tool in analysis number of concomitant medications taken over the followup time in order to support drug safety evidence interpretation. Building on these recent developments, this paper explores whether jointly modelling multivariate longitudinal clinical laboratory, concomitant medication and clinical AEs adjusted for baseline characteristics can be efficient in estimating AE incidence rate and useful in understanding joint evolution of the potentially correlated multiple safety outcomes.
In clinical trials with frequent AE occurrence (e.g. antimalarial treatment in pregnancy trials), event count models such as Poisson or negative binomial regression are considered important in quantifying the AE rate [13,14,15,16]. However, the AEs are also highly correlated with the cumulative concomitant medications and clinical laboratory safety outcomes. Joint modelling can offer an opportunity to explicitly model the correlations and account for the followup time, ensuring efficient and unbiased estimates. Despite the advanced development and evident benefits of joint modelling of multivariate longitudinal data [17,18,19], joint modelling has received limited attention in drug safety assessment. Recent statistical software developments have provided an opportunity for extending joint models to more different types of outcomes [20, 21], beyond the traditionally implemented joint models with single continuous longitudinal data and single survival outcome. In this paper, we present a joint model for longitudinal continuous clinical laboratory safety data (alanine amino transferase and total bilirubin), individual clinical AEs and concomitant medication. We focus on demonstrating the utility of the proposed model in assessing the association between the clinical AEs in relation to these outcomes controlling for baseline covariates. Furthermore, we assessed whether the joint model improved the estimates in log incidence rate ratio (IRR) on effect of treatment on clinical AEs compared to the conventional Poisson model. The proposed joint model is useful in providing evidencebased decisions in determining the most important clinical safety measurements to focus on for drug safety assessment especially in resourceconstrained settings.
Methods
Data
Our proposed joint model was applied to data from PREGACT trial as the motivating data. PREGACT trial (ClinicalTrials.gov number, NCT00852423) was a multicentre, openlabel 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 polymerasechain reaction (PCR) adjusted cure rates at day 63.
The safety outcomes collected at baseline and during followup were AEs, total bilirubin, alanine aminotransferase, white blood cell and red blood cell counts. The patients were directly observed on days 02 (after receiving the dose). The patients were then asked to return to the clinic for followup 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 followup 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.
Ethical considerations
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.
Notation and joint model specification
Let each patient who was randomized and received at least a dose be denoted as i = 1, … …, n. Let \({V}_i=\left({v}_{i1}^{\mathrm{T}},\kern0.5em {v}_{i2}^{\mathrm{T}}\right)\) 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 (v_{i1}, v_{i2}) are of (n_{ik}x1) dimension for the observed longitudinal measurements of the kth outcome; v_{ik} = (v_{i1k}, ……., v_{ink})^{T}. We accommodate the situations where observation times, t_{ijk} may differ between individuals and outcomes (j = 1, … . ., n_{ik}). For each patient, we let \({C}_i=\left({c}_{i1}^{\mathrm{T}},\kern0.5em {c}_{i2}^{\mathrm{T}}\right)\) represent a bivariate count outcome vector where k = 1 is for total number of AEs experienced during their followup time T_{i} and k = 2 for total number concomitant medication used over the followup time T_{i}. A set of covariates that were collected at baseline for each individual are defined as X_{i} = {X_{1i}, X_{2i}, … . ., X_{pi}}.
Conventional Poisson model for clinical AEs
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 safetyrelated 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 nonnegative 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, \({\boldsymbol{X}}_{ik}^{\mathrm{T}}\), 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 timevarying covariates (with time t0, representing baseline covariates in M1), the vector of covariates adjusted in model M2 is \({\boldsymbol{X}}_{ik}^{\mathrm{T}}\left(\boldsymbol{t}\right)\) 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 submodels; longitudinal submodel and timetoevent submodel (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 submodel is replaced with count outcome submodel. Since we assumed Poisson distribution of the clinical AEs count, the count outcome submodel was specified as Poissondistributed. Therefore, the joint model presented in this paper consist of two longitudinal continuous laboratory safety biomarkers submodels (one for total bilirubin and one for alanine aminotransferase) two Poisson submodels (one for the individual clinical AEs count and one for the individual concomitant medication count). The connecting of these outcomes submodels is very flexible such that we could use random effects or expected value of outcomes [23] as detailed in the model specification below.
Joint model formulation
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 submodels that made the joint model considered in this work. The estimates in all the models are considered on a logarithm scale.
Longitudinal submodel
We modelled two continuous longitudinal clinical laboratory safety outcomes (i.e. total bilirubin and alanine aminotransferase) using a longitudinal submodel 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 v_{i1}, v_{i2} was assumed to be normally distributed with mean μ_{k} and variance \({\delta}_k^2\) (k = 1, 2) on a logarithm scale. Each continuous longitudinal clinical laboratory outcome is modelled using the submodel 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 \({\boldsymbol{X}}_{ik}^{\mathrm{T}}(t)\boldsymbol{\beta}\), specified as a linear function of the covariates at a given time for the outcome k is flexible such that it can accommodate timevarying covariates. The Z_{i}(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 patientspecific shared random effects b_{i} is assumed to have a multivariate normal distribution with mean 0 and variancecovariance matrix Σ, s.t. (b_{i}) ∼ MVN(0, ∑).
Alanine aminotransferase submodel
Based on M1, considering patientspecific random intercept as our shared parameter of interest, linking the outcomes, the longitudinal submodel for alanine aminotransferase can be written simply as M4 below;
where b_{0i1} is the patient specific random intercept corresponding to the alanine aminotransferase outcome.
Total bilirubin submodel
Similarly, the longitudinal submodel for total bilirubin, M5, can be formulated as;
where b_{0i2} is the patient specific random intercept corresponding to total bilirubin outcome.
Count submodels
We considered two event count submodels; one for concomitant medication count and the other for the AE count. The two count submodels were derived from a modified form of the generic Poisson model. The count submodel was as follows;
Concomitant medication count submodel
Assuming that the concomitant medication count had a Poisson distribution, its submodel was defined as
such that log(φ b_{03i}) represents the logarithm of the Poisson mean response, φ, conditional on the patientspecific random intercept, b_{03i} for concomitant medication count (i.e. the concomitant medication count submodel is linked with the AE count submodel in the joint model via the random intercept). The \({\boldsymbol{X}}_{ik}^{\mathrm{T}}(t)\boldsymbol{\beta}\) 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 patientspecific deviation from the mean random intercept of the reported concomitant medication use.
Clinical AE count submodel
Considering that AE count was the primary outcome of interest and assuming a Poisson distribution for the AE count, the AE count submodel was formulated as;
The log(φ b_{01i}, b_{02i}, b_{03i}) 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 submodels and concomitant medication submodel (b_{01i}, b_{02i}, b_{03i}) is linked to the linear predictor \({\boldsymbol{X}}_{ik}^{\mathrm{T}}(t)\boldsymbol{\beta}\) 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 patientspecific 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 patientspecific deviation from the mean random intercept of a respective log alanine aminotransferase.
Alternative clinically meaningful formulation of the AE count submodel 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) v_{1}, v_{2}, under the current value association structure, we can modify the AE count submodel as;
The E[v_{1}(t) μ_{1}(t)] and E[v_{2}(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, v_{3}, submodel is still linked to the AE count submodel via a random intercept b_{03i}. We maintained the random intercept (b_{03i}α_{3}) parameterization since this efficiently deals with any potential overdispersion 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.
Joint model likelihood and estimation
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(b_{i} θ) is a normal density function for the random effects conditional on θ and f(V_{i} θ) is a normal density function for the log total bilirubin and log alanine aminotransferase conditional on θ and b_{i}. The f(T_{i}, C_{i} b_{i}, θ) represents Poisson density function conditional on θ and b_{i}; patientspecific total follow up days is T_{i} and C_{i} represents total events recorded (for AE count submodel and concomitant count submodel). 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 GaussHermite 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.
Statistical data analysis
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, bednet 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 followup 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 followup 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 betweenoutcome 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.
Results
The trial recruited 870 women who were randomly assigned to received AL, ASAQ, DHAPQ; each arm recruited 290 women. Baseline characteristics were similar across and between the trial arms (Table 1). Maternal age was skewed to the right. Overall, median maternal age was 20 years (IQR: 18.024.0). Overall mean haemoglobin concentration at enrolment (g/dl) was 10.1 (SD: 1.3) and distribution was also consistently similar across all the treatment arms. The overall reported bed net use in a night before enrolment was 23.2% such and it was similar across and between the study arms (Table 1).
In total 1512 AEs were observed over the whole followup time and 108 of these were definitely not related to the treatment. Proportion of women who experienced at least one AE in the AL, ASAQ and DHAPQ treatment arms were 64.5% (187 of 300 women), 70.7% (205 of 290 women) and 59.3% (172 of 290 women), respectively. Almost all patients took at least a concomitant medication over the followup such that for each arm 99% (287 out of 290) took at least a concomitant medication (Table 2). Across all the treatment arms, the reported median concomitant medications taken was 3; for AL 3(IQR: 25), for ASAQ 3(IQR: 14) and for DHAPQ 3(IQR: 25).
We observed that iron supplement (26.3%), paracetamol (26.2%) and albendazole (10.7%) were frequently reported concomitant medications (Table 3). Almost all the patients analysed, 99.2%, took the iron supplements and approximately twothirds, 60.3%, of the patients took paracetamol apart from the antimalarial drug under investigation. During the followup period, a total of 353 patients reported having used albendazole once.
Overall median of logtransformed alanine aminotransferase was 2.7 (IQR: 2.53.0) such that the medians for AL was 2.8 (IQR: 2.53.0), for ASAQ was 2.7 (IQR: 2.53.0) and for DHAPQ was 2.8 (IQR: 2.53.0). Overall median for logtransformed total bilirubin was 2.0 (IQR: 1.72.4). Figure 1 below shows the detailed distribution of each of clinical laboratory safety outcomes over the followup days by antimalarial drug treatment arm. As shown in graph 1 within Fig. 1, the alanine aminotransferase did not vary much throughout the followup days. The general trend indicated that across the followup days, the alanine aminotransferase concentration was slightly higher among the women who received DHAPQ compared to those women who received either AL or ASAQ. After day 0 the total bilirubin concentration gradually decreased across all the treatment arms and, on average, women treated with AL had the highest bilirubin concentration across the followup days (see Graph 2 within Fig. 1).
Parameter estimates for the models
In this section, we report the estimates from the conventional Poisson model (where AEs count is an outcome) (Table 4) and the proposed joint model (Table 5). Our interest was to establish whether the joint model yielded better results than the conventional Poisson model and how the jointly modelled outcomes evolved over the followup time in relation to clinical AEs; the results from the AE count submodel in the joint model (Table 5) are of main interest in our reporting.
Parameter estimates from the conventional Poisson model for the clinical AE count
Using the conventional Poisson model, M9, the clinical AE incidence rate did not differ significantly across the treatment groups (p < 0.149). However, the estimates suggested that the AE incidence rate was higher in the ASAQ treatment group compared to the AL treatment group (log IRR: 0.1838; 95% CI: 0.0574, 0.3102; p = 0.004) (Table 4). Exponentiating the point estimate of the log IRR yields IRR of 1.2018, implying that the patients treated with ASAQ were expected to have 20.18% higher rate of clinical AEs occurrence compared to those treated with AL. The AE incidence rate did not differ between DHAPQ and AL treatment groups (log IRR: 0.0255; 95% CI: − 0.1587, 0.1076; p = 0.707). A year increase in maternal age was associated with 0.0208 increase in log IRR of the AEs (95% CI: 0.0068, 0.0347; p = 0.004). The AE incidence rate was significantly higher among the those who enrolled in the second trimester compared to those who enrolled in the third trimester (log IRR: 0.3768; 95% CI: 0.1943, 0.5192; p < 0.001).
Parameter estimates from the joint model
We report the fixed and random parameter estimates for the final fitted joint model based on mixed association structure, M9, (i.e. random effects and current value association structure) in Table 5. AE count submodel indicated that the AE incidence rate did not vary across all the treatment groups (p = 0.376). However, the estimates indicated that AE incidence rate was higher in ASAQ treatment group compared to AL (log IRR: 0.2118, 95% CI: 0.0082, 0.4154; P = 0.041) suggesting that patients treated with ASAQ were expected to have 23.59% higher rate of clinical AEs occurrence compared to those treated with AL (i.e. IRR = exp.(0.2118) = 1.2359). The AE incidence rate did not significantly differ in DHAPQ treatment compared to AL (log IRR: 0.0508; 95% CI: − 0.2559, 0.1542). Our results show that the magnitude of the treatment effect was higher in the joint model, M9, compared to the estimates in the conventional Poisson model, M1, above. Among the association parameters describing joint evolution between the outcomes, only concomitant medication was associated with clinical AEs. The association parameter, α_{3}, linking concomitant medication and clinical AEs suggested that the clinical AEs outcome was positively and strongly associated with the concomitant medications (log IRR: 1.7487; 95% CI: 1.5471, 1.9503; p < 0.001). A unit increase for a patient in deviation from the mean random intercept of concomitant medication use was associated with 5fold increase in rate of clinical AEs (i.e.IRR = exp.(1.7487) = 5.7471). We found that the clinical AEs were not associated the alanine aminotransferase nor total bilirubin over the followup time as demonstrated by their respective association parameters that were nonsignificant. For the association between alanine aminotransferase and clinical AEs over followup time the estimated log IRR was − 0.0288 (95% CI: − 0.5045, 0.4469; p = 0.906) and for the association between total bilirubin and clinical AEs the log IRR was 0.1153 (95%CI: − 0.0889, 0.3194; p = 0.269). As also observed in the conventional Poisson model above, the clinical AE submodel indicated that those women who enrolled in the second trimester had a higher AE incidence rate compared to those who enrolled in the third trimester (log IRR 0.3986; 95% CI: 0.1592, 0.6381; p = 0.001); the women who enrolled in the second trimester had 48.97% higher rate of clinical AEs (i.e. IRR = exp.(0.3986) = 1.4897) compared to those who enrolled in the third trimester.
In the submodel for log alanine aminotransferase, we observed no significant difference across all the variables included in that submodel. For the log total bilirubin submodel, we observed 0.1786 lower log total bilirubin among the women who enrolled in the second trimester compared to those who enrolled in the third trimester (log IRR: 0.1786; 95% CI: − 0.2752, − 0.0820, p < 0.001). The concomitant submodel showed that concomitant medication was associated with the treatment group. Women enrolled in ASAQ treatment group reported lower incidence rate of concomitant medication use compared to those in AL treatment group (log IRR: 0.1841; 95% CI: − 0.3054, − 0.0629; p = 0.003). Women in DHAPQ treatment group also reported lower incidence rate of concomitant medication use compared to those in AL treatment group. Those women who enrolled in the second trimester reported a higher incidence rate of concomitant medication use than those who enrolled in the trial while in the third trimester (log IRR 0.2861; 95% CI: 0.1454, 0.4268; p < 0.001).
The standard deviations for patientspecific random intercepts, across the outcomes, suggested that concomitant medication had the highest unobserved heterogeneity with standard deviation of 0.5378 (95% CI: 0.4914, 0.5887). The lowest unobserved heterogeneity was observed in log alanine aminotransferase outcome with random intercept standard deviation of 0.2547 (0.2324, 0.2792).
Sensitivity analyses
As part of sensitivity analysis, we attempted to relax the normality assumption of the random effects. We fitted a robust joint model for the mixture association structure, M9, where the random effects were assumed to have a tdistribution. The model yielded similar results as those in Table 5 and Alkaike Information Criterion (AIC) of 18,230.28 which was comparable to an AIC of 18,278.53 for the model in Table 5. Secondly, fitting a joint model with random effects association structure, M5, we obtained similar estimates as those in Table 5 and with negligibly lower AIC of 18,277.87.
Following a reviewer’s suggestion, we further investigated the effect of the concomitant medications with high frequencies on AEs. Although iron supplement was the most commonly reported concomitant medication, we did not consider it in assessing the impact of specific concomitant medication on AEs since it was taken by almost all the patients, 99.2%. Instead, we considered paracetamol since it was the secondhighest reported concomitant medication where 60.3% of the patients reported it at least once. Assessment of the effect of treatment on clinical AEs using a joint model of a similar structure to M9 but where the concomitant medication count was confined to paracetamol yielded similar results as those obtained in Table 5 (see Table 6). The paracetamol concomitant medication was positively and strongly associated with clinical AEs (log IRR: 1.2969; 95% CI: 1.1328, 1.4609; p < 0.001). However, this observed effect of paracetamol concomitant medication was slightly lower compared to the overall concomitant medication effect as reported in Table 4 (log IRR: 1.7487; 95% CI: 1.5471, 1.9503; p < 0.001). No interaction was observed between the paracetamol concomitant medication and the treatment (p = 0.078) such that we focussed on reporting results with the main effect only model as shown in Table 6. The model with interaction terms had a negligibly lower AIC of 16,556.19 than the model without the interaction terms with an AIC of 16,557.29.
The interest and structure of the proposed joint model could not permit us to assess the impact of albendazole on AEs. The patients who took albendazole took it once; therefore, this concomitant medication arose as a binary outcome that does not meet the joint model specification (i.e. where the concomitant medication is incorporated in the model as a count outcome). Since the objective of the exploratory analysis was to establish the impact of specific frequent concomitant medication on AEs, we did not undertake any additional analysis on the effect of the “other” concomitant medication on AEs.
Discussion
In clinical trials, multiple safety outcomes collected over the followup time require advanced analyses in order to develop a valid drug safety profile that makes full use of the information from the multiple safety outcomes. This paper introduced a framework for joint modelling of multivariate continuous longitudinal clinical laboratory safety outcomes, concomitant medication counts and adverse event count, exploiting the advantages of the existing extended multivariate generalised linear and nonlinear mixed effects models framework [23]. Our work focussed on investigating whether clinical AEs varied by treatment arm and how the laboratory outcomes and concomitant medication were associated with clinical AEs over followup time in artemisininbased treatment of malaria during pregnancy clinical trial. We found that the concomitant medication was strongly associated with the clinical AEs incidence, suggesting that the patients who experienced more AEs were more likely to report the use of more concomitant medications. If not wellaccounted for such increased use of concomitant medication can distort the treatment effect on the clinical AE occurrence. From methodological perspective, we discovered that the conventional Poisson model underestimated the treatment effect on the AE incidence rate ratio; the effect of treatment on AE incidence rate estimates for the joint model were higher than those from conventional Poisson model. The difference in the models can be attributed to the fact that the joint model uses more of the available information through the explicit modelling of the observed correlation, unobserved heterogeneity and any nonlinear effects [25,26,27]. The joint model efficiently uses the available information, including the unobserved heterogeneity, leading to more improved estimates compared to the conventional Poisson model involving separate analysis.
Controlling for baseline covariate to address any potential confounding in AE incidence rate ratio estimates, we found that maternal age and enrolling in the second trimester were factors associated with the AE occurrence. A year increase in maternal age was associated with 2% increased rate of AE incidence. This could be as a result that as the women become older they are more likely to report AEs out of experience than the younger women. As expected women who enrolled in the second trimesters had 49% increased rate of AE incidence compared to those who enrolled in the third trimester. This could probably be due to the fact that physiological changes that take place during pregnancy are more impactful in the second than in the third trimester. Hence, the susceptibility to AE occurrence after taking either the antimalarial drug or the concomitant medication (that was strongly associated with increase rate of AE occurrence) since the pregnancy related body physiological changes can interact with the drug.
In assessment of the impact of concomitant medication on AEs, identifying specific concomitant medication is usually of interest. In this paper, we observed that paracetamol greatly contributed to the overall association between concomitant medication and AEs. Although the difference in the effect of overall concomitant medication on AEs and the effect of paracetamol concomitant medication on AEs was not big, the difference could not be considered clinically negligible. Therefore, undertaking both analyses (i.e. effect of overall concomitant medication and effect of specific most frequent concomitant medication on AEs) should be considered useful in understanding antimalarial drug safety profile. Alternatively, an interesting future joint model development can consider quantifying the contribution of each of the specific frequent concomitant medication on AEs occurrence. This can be helpful in identifying and handling any potential noisy concomitant medications that contribute less to the AEs occurrence.
Within the joint modelling framework, our proposed model provides additional knowledge by extending the number of outcomes that can be analysed simultaneously, including multiple longitudinal outcomes and multiple count outcomes, focussing on drug safety assessment. We acknowledge that other researchers have previously proposed similar kind of modelling but applied in the efficacy context [28,29,30,31,32]. Buu et al. considered a joint model of count and binary outcome [30]. Li et al. considered jointly modelling proportion, count and continuous longitudinal outcomes [31]. Yang and Kang introduced a joint model for mixed Poisson outcome and continuous outcome [32]. Our model is considered as a special case of these previously proposed models focussing on drug safety assessment. Key unique features of our proposed joint model include incorporating mixed association structure of random effects and expected value for linking the outcomes, use of restricted cubic splines to accommodate the capturing of any potential timedependent effects, accommodating multiple safety outcomes and novel application to drug safety assessment. A case study based paper like ours can further facilitate the adoption of improved statistical analysis of safety data in clinical trials [33].
In this paper, data missingness mechanism was ignorable. Our future work will consider investigating the impact of missing data on the model performance under varying followup times and model association structures. Secondly, in this paper we worked with a large sample size. An interesting further methodological work would be to assess how the sample size impacts the performance of the model that we applied. Unfortunately, the key limitation of the joint model introduced and applied in this paper is that it is computationally demanding such that it takes longer time to converge. This directly affects the scope of the simulations that can be done to investigate the highlighted potential methodological issues associated with such complicated models. An extension can also be made to the number of outcomes included the joint model depending on the availability of the data. For example, aspartate amino transferase, highly correlated with alanine amino transference, can be added in our joint model.
Conclusion
Jointly modelling of multivariate longitudinal clinical laboratory safety data, concomitant medication and clinical AEs efficiently harnesses the safety data in order to better understand drug safety profile. This paper has demonstrated the utility of joint modelling to yield improved AE incidence rate estimates compared to the conventional Poisson model. The proposed joint model helps to quantify the association between the laboratory safety outcomes, concomitant medication, and clinical AEs over followup time. Key public health relevance of the joint modelling of the safety outcomes includes facilitating the choice of the most important clinical laboratory safety and other safety outcomes for profiling drug safety, especially in resourceconstrained settings. In this exploratory analysis, we have also established that the AEs that were observed in the PREGACT trial in Malawi were more linked to concomitant medications that the patients took. Our findings provide an assurance that artemisininbased treatments can be safely used in second and third trimester during pregnancy.
Availability of data and materials
The data that support the findings of this study are available from Institute of Tropical Medicine, Antwerp, Belgium but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the corresponding author upon reasonable request and with permission of Institute of Tropical Medicine, Antwerp, Belgium (http://www.itg.be).
Abbreviations
 AEs:

Adverse events
 AL:

Artemether–lumefantrine
 ASAQ:

Amodiaquine–artesunate
 DHAPQ:

Dihydroartemisinin–piperaquine
 IRR:

Incidence rate ratio
 CI:

Confidence interval
 PCR:

Polymerasechain reaction
 AIC:

Akaike Information Criterion
References
 1.
Zink RC, Marchenko O, SanchezKam M, Ma H, Jiang Q. Sources of safety data and statistical strategies for design and analysis: clinical trials. Ther Innov Reg Sci. 2018;52(2):141–58.
 2.
Munsaka MS. A questionbased approach to the analysis of safety data. In: Biopharmaceutical applied statistics symposium. Biostatistical analysis of clinical trials, vol. 2. Singapore: Springer; 2018. p. 193–216.
 3.
Merz M, Lee KR, KullakUblick GA, Brueckner A, Watkins PB. Methodology to assess clinical liver safety data. Drug Saf. 2014;37(1):33–45.
 4.
ChuangStein C. Laboratory data in clinical trials: a statistician's perspective. Control Clin Trials. 1998;19(2):167–77.
 5.
Klepper MJ, Cobert B. Drug safety data: how to analyze, summarize and interpret to determine risk. Sudbury: Jones & Bartlett Learning; 2010.
 6.
Siddiqui O. Statistical methods to analyze adverse events data of randomized clinical trials. J Biopharm Stat. 2009;19(5):889–99.
 7.
Unkel S, Amiri M, Benda N, Beyersmann J, Knoerzer D, Kupas K, et al. On estimands and the analysis of adverse events in the presence of varying followup times within the benefit assessment of therapies. Pharm Stat. 2019;18(2):166–83.
 8.
Patson N, Mukaka M, Otwombe KN, Kazembe L, Mathanga DP, Mwapasa V, et al. Systematic review of statistical methods for safety data in malaria chemoprevention in pregnancy trials. Malar J. 2020;19(1):1–11.
 9.
Barker C. Exploratory method for summarizing concomitant medication data–the mean cumulative function. Pharm Stat. 2010;9(4):331–6.
 10.
Schildcrout JS, Jenkins CA, Ostroff JH, Gillen DL, Harrell FE, Trost DC. Analysis of longitudinal laboratory data in the presence of common selection mechanisms: a view toward greater emphasis on premarketing pharmaceutical safety. Stat Med. 2008;27(12):2248–66.
 11.
Rosenkranz GK. Modeling laboratory data from clinical trials. Computat Stat Data Anal. 2009;53(3):812–9.
 12.
Southworth H, Heffernan JE. Multivariate extreme value modelling of laboratory safety data from clinical studies. Pharm Stat. 2012;11(5):367–72.
 13.
Cook RJ, Lawless J. The statistical analysis of recurrent events. New York: Springer; 2007.
 14.
Therneau TM, Grambsch PM. Modeling survival data: extending the cox model. Newyork: Springer; 2000.
 15.
Fu H, Luo J, Qu Y. Hypoglycemic events analysis via recurrent timetoevent (HEART) models. J Biopharm Stat. 2016;26(2):280–98.
 16.
Luo J, Qu Y. Analysis of hypoglycemic events using negative binomial models. Pharm Stat. 2013;12(4):233–42.
 17.
Hickey G, Philipson P, Jorgensen A, KolamunnageDona R. A comparison of joint models for longitudinal and competing risks data, with application to an epilepsy drug randomized controlled trial. J R Stat Soc Ser A (Stat Soc). 2018;181(4):1105–23.
 18.
Hickey GL, Philipson P, Jorgensen A, KolamunnageDona R. Joint modelling of timetoevent and multivariate longitudinal outcomes: recent developments and issues. BMC Med Res Methodol. 2016;16(1):117.
 19.
Syrjälä E, Nevalainen J, Peltonen J, Takkinen HM, Hakola L, Åkerlund M, et al. A joint modeling approach for childhood meat, fish and egg consumption and the risk of advanced islet autoimmunity. Sci Rep. 2019;9(1):1–10.
 20.
Hickey GL, Philipson P, Jorgensen A, KolamunnageDona R. joineRML: a joint model and software package for timetoevent and multivariate longitudinal outcomes. BMC Med Res Methodol. 2018;18(1):50.
 21.
Crowther MJ. merlin—a unified modeling framework for data analysis and methods development in Stata. Stata J. 2020;20(4):763–84.
 22.
Pekyi D, Ampromfi AA, Tinto H, TraoreCoulibaly M, Tahita MC, Valea I, et al. Four artemisininbased treatments in African pregnant women with malaria. New England J Med. 2016;374(10):913–27.
 23.
Crowther MJ. Extended multivariate generalised linear and nonlinear mixed effects models. 2017. https://arxiv.org/abs/1710.02223. Accessed 5 Dec 2020.
 24.
Tuerlinckx F, Rijmen F, Verbeke G, De Boeck P. Statistical inference in generalized linear mixed models: a review. Br J Math Stat Psychol. 2006;59(2):225–55.
 25.
Henderson R, Diggle P, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics. 2000;1(4):465–80.
 26.
Ibrahim JG, Chu H, Chen LM. Basic concepts and methods for joint models of longitudinal and survival data. J Clin Oncol. 2010;28(16):2796.
 27.
Long JD, Mills JA. Joint modeling of multivariate longitudinal data and survival data in several observational studies of Huntington’s disease. BMC Med Res Methodol. 2018;18(1):1–15.
 28.
Chi YY, Ibrahim JG. Joint models for multivariate longitudinal and multivariate survival data. Biometrics. 2006;62(2):432–45.
 29.
Andrinopoulou ER, Clancy JP, Szczesniak R. Multivariate joint modeling to identify markers of growth and lung function decline that predict cystic fibrosis pulmonary exacerbation onset. BMC Pulm Med. 2020;20:1–11.
 30.
Buu A, Li R, Tan X, Zucker RA. Statistical models for longitudinal zeroinflated count data with applications to the substance abuse field. Stat Med. 2012;31(29):4074–86.
 31.
Li H, Zhang Y, Carroll RJ, Keadle SK, Sampson JN, Matthews CE. A joint modeling and estimation method for multivariate longitudinal data with mixed types of responses to analyze physical activity data generated by accelerometers. Stat Med. 2017;36(25):4028–40.
 32.
Yang Y, Kang J. Joint analysis of mixed Poisson and continuous longitudinal data with nonignorable missing values. Computat Stat Data Anal. 2010;54(1):193–207.
 33.
Phillips R, Cornelius V. Understanding current practice, identifying barriers and exploring priorities for adverse event analysis in randomised controlled trials: an online, crosssectional survey of statisticians from academia and industry. BMJ Open. 2020;10(6):e036875.
Acknowledgements
We thank all the study participants and data collectors in the clinical trial.
Funding
The clinical trial was supported by the European and Developing Countries Clinical Trials Partnership, the Malaria in Pregnancy Consortium (which is funded through a grant from the Bill and Melinda Gates Foundation to the Liverpool School of Tropical Medicine), the Belgian Development Cooperation Agency, the Liverpool School of Tropical Medicine, the Medical Research Council UK, the Netherlands Organization for Scientific Research, and SanofiAventis NP was funded by the Fogarty International Center of the National Institutes of Health [Grant number D43TW010075 to MKL, TC and DM].
Author information
Affiliations
Contributions
Study idea conceptualization: NP, TC, MM. Methodology development: NP, MM, LK, TC, MKL, DM, VM, GC, UD. Analysis and first draft development: NP. Supervision, writing revised draftsreviewing and editing: MM, LK, TC, MKL, DM, VM, GC, UD. Enrolled patients and conducted clinical investigation: VM, GC, UD. Reading and approving the final manuscript: NP, MM, LK, TC, MKL, DM, VM, GC, UD.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
This study was approved by the University of the Witwatersrand’s Human Research Ethics Committee. The PREGACT trial received ethical clearance from ethics committee at the Antwerp University Hospital and College of Medicine Research Ethics Committee at the University of Malawi in Malawi [22]. Informed consent was obtained from all the participants.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Patson, N., Mukaka, M., D’Alessandro, U. et al. Joint modelling of multivariate longitudinal clinical laboratory safety outcomes, concomitant medication and clinical adverse events: application to artemisininbased treatment during pregnancy clinical trial. BMC Med Res Methodol 21, 208 (2021). https://doi.org/10.1186/s12874021014129
Received:
Accepted:
Published:
Keywords
 Adverse events
 Randomised controlled trials
 Joint model
 Drug safety
 Concomitant medication