Linear mixed models to handle missing at random data in trial-based economic evaluations.

Trial-based cost-effectiveness analyses (CEAs) are an important source of evidence in the assessment of health interventions. In these studies, cost and effectiveness outcomes are commonly measured at multiple time points, but some observations may be missing. Restricting the analysis to the participants with complete data can lead to biased and inefficient estimates. Methods, such as multiple imputation, have been recommended as they make better use of the data available and are valid under less restrictive Missing At Random (MAR) assumption. Linear mixed effects models (LMMs) offer a simple alternative to handle missing data under MAR without requiring imputations, and have not been very well explored in the CEA context. In this manuscript, we aim to familiarize readers with LMMs and demonstrate their implementation in CEA. We illustrate the approach on a randomized trial of antidepressants, and provide the implementation code in R and Stata. We hope that the more familiar statistical framework associated with LMMs, compared to other missing data approaches, will encourage their implementation and move practitioners away from inadequate methods.

INTRODUCTION

Cost‐effectiveness analyses (CEAs) conducted alongside randomized controlled trials are an important source of information for decision‐makers in the process of technology appraisal (Ramsey et al., 2015). The analysis is based on healthcare outcome data and health service use, typically collected at multiple time points and then combined into overall measures of effectiveness and cost (Drummond et al., 2005). The derivation of these measures is problematic in the presence of missing outcome data, which are a common issue in trials, with on average only 60%–75% of randomized participants having complete CEA data (Leurent et al., 2020; Noble et al., 2012). A popular approach to handle missingness is to discard the participants with incomplete observations (complete case analysis (CCA)), allowing for derivation of the overall measures based on the completers alone. We note that slightly different definitions of CCA are possible, depending on the form of the model of interest, the type of missingness and the inclusion of observed covariates. Without loss of generality, throughout the paper, we refer to CCA as to when only individuals with fully‐observed effect and cost data (i.e., the completers) are included in the analysis, thus requiring the removal of all cases with partially‐observed outcome data. This approach, although appealing by its simplicity, has well‐recognized limitations including loss of efficiency and an increased risk of bias (Carpenter & Kenward, 2012; Faria et al., 2014; Rubin, 1987). In recent years there has been an increment in the uptake of more appropriate statistical methods, such as multiple imputation (MI) (Leurent et al., 2018b) or Bayesian methods (Gabrio et al., 2019), which improve efficiency, rely on less restrictive missingness assumptions, and facilitate the task of conducting sensitivity analysis to alternative missingness assumptions (Gabrio et al., 2020; Leurent et al., 2020). Among these, Missing At Random [MAR ; Rubin (1987)] often provides a desirable starting point as it implies that valid inferences can be drawn based on the observed data. Both CCA and MI are generally valid under some form of MAR assumptions. When estimating a difference between arms adjusted for baseline covariates a CCA will be valid if being a complete‐case is only dependent on these covariates. Multiple imputation will be valid under less restrictive assumptions, as missingness could also depend on the observed outcome at the other time‐points and other variables not in the analysis model. Although MI has become more frequent and accessible in CEAs, its' implementation in routine analyses remains limited compared to CCA (Gabrio et al., 2017; Leurent et al., 2018a), perhaps because of a lack of familiarity, computational time, or the analytical difficulties of the approaches. Indeed, it requires particular care when combined with other statistical procedures, such as taking into account clustering (R. Gomes et al., 2012), or using bootstrap re‐sampling (Brand et al., 2019).

We propose the use of linear mixed effects models (LMMs) as an alternative approach under MAR. Linear mixed effects models are commonly used for the modeling of dependent data (e.g., repeated‐measures) and belong to the general class of likelihood‐based methods. Linear mixed effects models have been occasionally used in CEA under the denomination of multilevel or hierarchical models (Manca et al., 2005; Rice & Jones, 1997) to account for the dependence between observations with a hierarchical structure (e.g., cluster‐randomized trials). Repeated‐measures also follow a hierarchical structure since data within each individual are correlated, and LMM are increasingly used to analyze longitudinal outcome in trials. However, LMMs appear surprisingly uncommon for the analysis of repeated measures in trial‐based CEA, perhaps because of a lack of awareness or familiarity with fitting LMMs. To our knowledge, only (Faria et al., 2014) briefly examined the use of LMMs in trial‐based CEA as an alternative to MI, and concluded that they offer a valid approach under MAR.

In this letter, we aim to familiarize readers with the implementation of LMMs in trial‐based CEA using standard software and summarize the statistical and economic results from a case study. Finally, we discuss the proposed approach and provide some suggestions for future work.

METHODS

Linear mixed effects models for repeated measurements

Linear mixed effects model extends the usual linear model framework by the addition of “random effect” terms, which can take into account the dependence between observations. A simple model for a repeated measure can be written as where Y denotes the outcome repeatedly collected for each individual i = 1,…,N at multiple times j = 1,…,J. The model parameters commonly referred to as fixed effects include an intercept β and the coefficients (β ,…,β ) associated with the predictors X ,…,X , while ω and ϵ are two random terms: ϵ is the usual error term and ω is a random intercept which captures variation in outcomes between individuals. The random terms are typically assumed to be normally distributed such that and . Equation (1) treats the data as having a 2‐level structure, where and capture the variance of the responses within (level 1) and between (level 2) individuals, respectively. The models can be extended to deal with more complex structures, for example, by allowing the effect of the covariates to vary across individuals (random slope) or a different covariance structure of the errors. Linear mixed effects models can be fitted even if some outcome data are missing and provide correct inferences under MAR (Schafer & Graham, 2002). We refer interested readers to the available LMM literature for an in‐depth description of the methods (Molenberghs et al., 2004).

One particular type of LMMs commonly used in the analysis of repeated measures in clinical trials is referred to as a Mixed Model for Repeated Measurement (MMRM; Wolfinger, 1993). The model includes a categorical effect for time, an interaction between time and treatment arm, and allows errors to have different variance and correlation over time (i.e., unstructured covariance structure). Consider the case where we want to model health related quality of life data (i.e., utilities) collected at three times (baseline and two follow‐ups). The model can then be expressed as: where U is the utility measured for patient i at time j, TIME and TRT are the indicators for the time and treatment arm, and TIME TRT are the interaction terms between time and treatment arm at each follow‐up (j ≥ 2). Within this model, β , represents the mean utility at baseline, β and β represent the mean utility for the control arm (TRT  = 0) at the first and second follow‐up, while β and β capture the mean difference between the arms at the two follow‐ups. The errors at each time follow a multivariate normal distribution ϵ ∼ N(0, Σ) with an unstructured covariance matrix Σ, that is, errors are allowed to have different variance and correlation between them. Note that no treatment parameter was included at baseline, this makes sense as no treatment effect is expected at baseline, but is also important to obtain β and β parameters adjusted for baseline utility (see Section 2.2). The model can be easily extended to handle additional time‐points, or to model any other repeatedly measured continuous outcome such as cost.

Incremental (between‐group) or marginal (within‐group) estimates for aggregated outcomes over the trial period, such as quality‐adjusted life years (QALYs) or total costs can be retrieved as linear combinations of the parameter estimates from Equation (2). For example, the mean difference in total cost is obtained by summing up the estimated differences at each follow‐up point, while differences on a QALYs scale can be obtained as weighted linear combinations of the coefficient estimates of the utility model (see Appendix C).

Adjusting for baseline variables

It is standard practice to adjust the analysis for imbalances in some baseline characteristics, to control for potential imbalance and as well as to increase precision (Manca et al., 2005) while also strengthening the plausibility of the MAR assumption (Little et al., 2012). We note that, differences in baseline outcome values between arms are already adjusted for within the model specification in Equation (2). Other baseline variables can be adjusted for by adding them as predictors to the model, although LMMs require covariates to be completely observed. However, in randomized controlled trials, missing baseline data can be usually addressed by implementing single imputation techniques (e.g., mean‐imputation) to obtain complete data prior to fitting the model, without loss of validity or efficiency (White & Thompson, 2005).

Assessing cost‐effectiveness

Once estimates for the average effectiveness and total cost differences are obtained, uncertainty can be assessed through resampling methods (e.g., bootstrapping). Results are then often summarized using a Cost‐Effectiveness Plane (CEP; Black [1990]) and a Cost‐Effectiveness Acceptability Curve (CEAC; Van Hout et al. [1994]).

TRIAL OVERVIEW

Overview

The Health Technology Assessment Study of the Use of Antidepressants for Depression in Dementia (SADD, ISRCTN88882979) was a placebo‐controlled randomized trial of participants from old‐age psychiatry services in England (Banerjee et al., 2011). A total of n = 326 participants were enrolled and randomized to receive placebo, or either one of two antidepressants (sertraline and mirtazapine) together with treatment as usual, and followed for 9 months. Outcomes were collected at baseline, 3 and 9 months, including EQ‐5D‐3L which was converted into utility (Dolan, 1997), and resource use for 0–3 months and 3–9 months which were converted into costs (Romeo et al., 2013). Details about the objective and conclusions from the trial as well as the study CEA are reported elsewhere (Banerjee et al., 2013). We restrict our attention to the 9 months cost‐utility analysis of placebo (n  = 111) versus mirtazapine (n  = 108) from the health and social care cost perspective.

Table 1 shows the missing data patterns of the utility and cost data (U ,C ) at baseline (j = 1), 3 months (j = 2) and 9 months (J = 3) follow‐up.

TABLE 1

Missingness patterns for the utility and cost variables in the SADD study

Missing data patterns	Placebo (N 1  = 111)	Mirtazapine (N 2  = 108)	Total (N = 219)
U 1 U 2 U 3 C 1 C 2 C 3	N 1 (%)	N 2 (%)	N(%)
‐ ‐ ‐ ‐ ‐ ‐	54 (48%)	47 (43%)	101 (46%)
‐ X X ‐ X X	12 (11%)	17 (16%)	29 (13%)
‐ X X ‐ ‐ ‐	15 (13%)	14 (13%)	29 (13%)
‐ ‐ X ‐ ‐ ‐	11 (10%)	13 (12%)	24 (11%)
‐ ‐ X ‐ ‐ X	9 (8%)	8 (7%)	17 (8%)
‐ X ‐ ‐ ‐ ‐	4 (4%)	4 (4%)	8 (4%)
‐ X X ‐ ‐ X	2 (2%)	2 (2%)	4 (2%)
‐ ‐ ‐ ‐ X X	2 (2%)	3 (3%)	5 (2%)
‐ ‐ ‐ ‐ ‐ X	2 (2%)	0 (0%)	2 (1%)

Note: For each pattern and treatment group the corresponding number and proportions of participants are reported. U and C indicate utility or cost variable at time j, while ‐ and X indicate an observed or missing variable, respectively.

Abbreviation: SADD, Study of the Use of Antidepressants for Depression in Dementia.

About half of the trial participants had complete cost‐effectiveness data, with utility values which tended to be more frequently missing than costs.

Cost‐effectiveness analysis of SADD

We applied the LMM model described above to the SADD data to estimate the difference in QALYs and cost, adjusting for baseline values (utility or cost). A set of 10,000 bootstrap replications were performed to derive cost‐effectiveness results. All models were fitted in R, using the packages and functions described in Appendix A.1. The R and Stata code for implementing the methods are given in Appendix A.2 and Appendix A.3, respectively, and on an online repository (Gabrio, 2022).

RESULTS

Table 2 reports the empirical means and standard deviations associated with the utility and cost variables in each treatment group over the study period.

TABLE 2

Observed means (and standard deviations) of the utility and cost variables at baseline, 3 and 9 months follow‐up for the placebo and mirtazapine group in the SADD trial

	Placebo	(N 1  = 111)	Mirtazapine	(N 2  = 108)
Utilities
Baseline	0.670 (0.268)	N = 111	0.688 (0.292)	N = 108
3 months	0.733 (0.286)	N = 78	0.763 (0.267)	N = 71
9 months	0.734 (0.274)	N = 62	0.827 (0.195)	N = 64
Costs
Baseline	1514 (3153)	N = 111	1191 (2256)	N = 108
3 months	1437 (3338)	N = 97	1130 (1868)	N = 88
9 months	2146 (4401)	N = 84	2550 (4288)	N = 78

Abbreviation: SADD, Study of the Use of Antidepressants for Depression in Dementia.

People in the mirtazapine group had higher mean utilities at all time‐points, lower mean costs at the first follow‐up and higher mean costs at the last follow‐up.

Table 3 reports the LMM‐estimated marginal means of the utilities and costs at each follow‐up by group, as well as QALYs and total costs. Incremental results, incremental cost‐effectiveness ratio and probability of cost‐effectiveness (at an acceptance threshold of £25,000) between the two groups are also reported.

TABLE 3

Linear mixed model results for the SADD trial

	Placebo		Mirtazapine		Incremental
	Estimate	95% CI	Estimate	95% CI	Estimate	95% CI
U 2	0.731	(0.675; 0.786)	0.752	(0.694; 0.810)	0.021	(−0.053; 0.096)
U 3	0.727	(0.671; 0.783)	0.781	(0.721; 0.841)	0.054	(−0.024; 0.132)
QALYs	0.540	(0.510; 0.571)	0.561	(0.530; 0.593)	0.021	(−0.016; 0.059)
C 2	1369	(840; 1899)	1252	(695; 1809)	−117	(−866; 631)
C 3	2092	(1226; 2959)	2760	(1856; 3664)	668	(−544; 1880)
Total costs	3462	(2236; 4688)	4012	(2735; 5289)	550	(−1156; 2257)
ICER					57,774
Cost‐effectiveness probability					34%

Note: Point estimates and 95% confidence intervals for the marginal means and mean difference of the utilities and costs at 3 (j = 2) and 9 (j = 3) months follow‐up as well as for the quality‐adjusted life years (QALYs) and total costs over the duration of the study. Incremental results are also reported in terms of the incremental cost‐effectiveness ratio (ICER) and cost‐effectiveness probability at k = £25,000 per QALY gained. Costs are expressed in British pounds.

Abbreviation: SADD, Study of the Use of Antidepressants for Depression in Dementia.

The results suggest that the mirtazapine group is associated with higher utilities compared with the placebo group at all follow‐ups, which leads to higher mean QALYs estimates. Although mean costs at 13 weeks are on average higher in the placebo group, the mirtazapine group is associated with relatively higher mean costs at 9 months as well as higher mean total costs. Overall, the mean difference in QALYs was 0.021(95%CI:−0.016;0.059), and difference in cost £550 (95%CI:−1156;2257). We additionally compare the results in terms of mean QALYs and total costs obtained under LMM with those from two alternative analyses: CCA and MI, the latter producing estimates in line with those of LMM with standard errors reduced by approximately13% for QALYs and 25% for total cost estimates compared to CCA (see Appendix B.2).

Figure 1 shows the associated CEP and CEAC. Mirtazapine appears to be slightly more effective and costly, but with important uncertainty, resulting in a probability of being cost‐effective of around 50% at a threshold of £25,000 per QALY gained.

FIGURE 1

Cost‐Effectiveness Plane (CEP) and Cost‐Effectiveness Acceptability Curve from the Study of the Use of Antidepressants for Depression in Dementia based on a linear mixed effects model fitted on 10,000 bootstrap replications. In the CEP, an acceptance threshold of k = £25,000 per quality‐adjusted life year (QALY) gained is used and the position of the incremental cost‐effectiveness ratio (ICER) on the plot is denoted with a darker colored dot. (a) Cost‐effectiveness plane. (b) Cost‐effectiveness acceptability curve

DISCUSSION

In this article we proposed the use of LMMs as an alternative approach to conduct trial‐based CEAs Although repeatedly criticized in the literature, analysts commonly handle missingness in QALYs and total costs via CCA (Gabrio et al., 2017; Leurent et al., 2018a), which causes some loss of information and potentially biases the inferences. The use of imputation methods, especially MI, have been recommended in the literature since they can obtain valid inferences using all observed responses (Faria et al., 2014). However, there may be practical obstacles to the spread of these methods among practitioners. First, analysts may be unfamiliar with the implementation of MI, and feel more confident using simpler approaches. Second, MI can be time‐consuming, particularly when combined with bootstrapping, for which there are alternative implementation strategies possible (Brand et al., 2019).

We believe LMMs represent an alternative approach which can overcome some of these limitations. First, practitioners may be more comfortable with the standard regression framework. Second, LMMs can be tailored to address other data features (e.g., cluster‐randomized trials or non‐normal distribution) while also easily combined with bootstrapping. Third, LMMs do not rely on imputation, and results are therefore deterministic and easily reproducible, whereas the Monte Carlo error associated with MI may cause results to vary from one imputation to another, unless the number of imputations is sufficiently large.

There are also limitations associated with the use of the proposed methods. First, LMM estimates are valid under MAR conditional on the observed outcome and the baseline variables in the analysis model. Multiple‐imputation offers a more flexible framework allowing for inclusion of auxiliary variables not adjusted for in the analysis. Second, while MAR is often a reasonable starting point, sensitivity analysis under non‐at‐random assumptions should be considered (Gabrio et al., 2020; Leurent et al., 2020). While these can be conducted within a LMM framework (Daniels & Hogan, 2008), a MI or Bayesian framework are particularly well suited for this type of analyses (M. Gomes et al., 2019; Mason et al., 2018). Third, the LMM approach works well for repeatedly‐measure outcomes, but addressing missingness within each component of the outcomes (e.g., for disaggregated cost components) increases the complexity of the model. Depending on the amount of observed data and time points, the suggested model may not achieve convergence, and simpler covariance structures could be considered.

Further work can be done to extend the proposed methods to tackle some of the typical features of CEA data. An interesting extension would be to fit a bivariate LMM model to model simultaneously the costs and utilities. Furthermore, while LMMs are robust to non‐normally distributed data in large sample, generalized linear model specifications could be applied to improve the fit to the data using non‐normal distributions (Nixon & Thompson, 2005). A possible solution is to explore the use of multivariate generalized LMMs to account for the correlation and characteristics of the data while also providing a coherent modeling framework based on the same number of individuals for both outcomes.

To conclude, we have shown how LMMs can be used for the analysis of trial‐based CEAs Although the methodology illustrated is already known, particularly in the area of statistical analyses, to our knowledge LMMs have rarely been applied to health economic data collected alongside randomized trials. We believe the proposed methods is preferable to a complete‐case analysis when CEA data are incomplete, and that it can offer an interesting alternative to imputation methods.

CONFLICT OF INTEREST

Sube Banerjee reports grants from NIHR HTA Program, during the conduct of the study. Personal fees from Otsuka, Lundbeck, Lilly, Hamad Medical Services, Nutricia, Axovant, and grants from ESRC, NIHR, UKRI, outside the submitted work. Baptiste Leurent reports personal fees from GAIA AG, outside the submitted work. Andrea Gabrio and Catrin Plumpton declare that they have no conflict of interest.

TABLE A1

Parameters from the linear mixed models (LMM) fitted to the SADD utility and cost data

Fixed effects	Utilities model		Costs model
	Estimate	95% CI	Estimate	95% CI
TIME 1	0.674	(0.637; 0.711)	1355	(991; 1720)
TIME 2	0.731	(0.675; 0.786)	1369	(840; 1899)
TIME 3	0.727	(0.671; 0.783)	2092	(1226; 2959)
TIME 2 TRT	0.021	(−0.053; 0.095)	−117	(−866; 631)
TIME 3 TRT	0.054	(−0.024; 0.132)	668	(−544; 1880)

Note: Point estimates and 95% confidence intervals (CI).

Abbreviation: SADD, Study of the Use of Antidepressants for Depression in Dementia.

TABLE A2

Estimates and standard errors (SE) for the mean QALYs and total cost in the SADD trial under complete case analysis (CCA), multiple imputation (MI), and linear mixed model (LMM)

	Placebo		Mirtazapine		Incremental
	Estimate	SE	Estimate	SE	Estimate	SE
CCA
QALYs	0.557	(0.016)	0.579	(0.017)	0.022	(0.023)
Total costs	3375	(774)	3695	(831)	320	(1148)
MI
QALYs	0.541	(0.013)	0.561	(0.014)	0.021	(0.020)
Total costs	3472	(593)	4023	(623)	552	(862)
LMM
QALYs	0.540	(0.016)	0.561	(0.016)	0.021	(0.019)
Total costs	3462	(623)	4012	(649)	550	(868)

Abbreviation: QALYs, quality‐adjusted life years; SADD, Study of the Use of Antidepressants for Depression in Dementia.