Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
44,872 characters · 24 sections · 43 citation commands
Nonparametric regression for cost-effectiveness analyses with observational data - a tutorial
Key points
A central challenge in healthcare decision-making is determining which interventions to implement or reimburse within the limits of constrained budgets. This often necessitates choosing between two or more competing treatment options for the same medical condition, where one treatment proves to be both more effective and more costly than the other. Cost-effectiveness analyses (CEAs) are designed to inform such choices by assessing whether the additional health gains of a treatment justify its additional costs. A key component of a CEA is the estimation of the treatment effects on both patients’ health outcomes and associated costs using available data.
If cost-effectiveness data are obtained from a randomized trial, where treatment allocation does not depend on the patients' baseline characteristics, the statistical estimation of treatment effects is relatively simple. If the data are observational, however, and treatment allocation is hence related to the patients' baseline characteristics, the estimation procedure is much more complex. Despite the availability of several textbooks on the analysis of observational data imbens2015pate,aronow2019foundations,hernan2024causal, appropriate methodology is not widely used in the field of CEA. That is, many researchers rely on simple linear regression or propensity score matching, preventing them from capturing complex, nonlinear relationships or interactions within the data. While there are some proposals for more sophisticated statistical methods in CEAs kreif2013regression,chen2024tutorial, they still depend on restrictive parametric models, which may limit them from appropiately modeling complex data structures.
We see several potential reasons for the reliance on such relatively simple models. One is simply researchers' potential unawareness of useful alternatives, or of the advantages which these alternatives offer. Another is the added complexity of nonparametric models, with respect to both usage and interpretation. In this paper, we aim to enlargen the toolkit available to cost-effectiveness researchers, by providing an accessible guideline on how to use state-of-the-art nonparametric regression methods for CEAs. Our focus will be on the Bayesian Additive Regression Trees (BART) model, which has recently been extended for application in the cost-effectiveness setting.
We begin by explaining what researchers aim to learn from CEAs using observational data. Our presentation will be quite informal; a more rigorous version may be found in esser2024seemingly. We shall not consider more general issues in the analysis of observational data, such as well-defined patient populations, treatment options, follow-up periods and so on. Such aspects are important, arguably more important than the statistical aspects, but they apply to CEAs in exactly the same way as they do to other epidemiological analyses. We thus refer the reader to canonical sources on these issues, such as hernan2016using.
In a typical CEA, the researcher is faced with two treatment options, $t = 0$ and $t = 1$. Often the former is called the control, and the latter the intervention. The fundamental goal of a CEA is to estimate, given available data, the probability that the intervention is cost-effective, compared with the control. To make this possible, we consider two essential quantities:
As indicated above, we frequently encounter the situation that both treatment effects are positive, meaning that the intervention is both more effective and more expensive than the intervention. A natural question then arises: do the added health gains of the intervention justify the increase in costs? To answer this, we need to quantify how much we are willing to pay for an additional unit of health. This quantity is called the willingness-to-pay (WTP) and is customarily assigned the symbol $\lambda$. \\ \\ Example 1. Suppose that costs are measured in EUR, and $\Delta_c = 1000$; health is measured in QALYs, and $\Delta_q = 0.1$. Then, by implementing the intervention, we gain an average health benefit of $0.1$ QALYs per patient, while, having to pay an additional $1000$ EUR. We hence need to pay $10000$ EUR in exchange for one QALY gained (because $1000/0.1=10000$). The intervention is cost-effective if and only if $\lambda$ is greater than $10000$. In other words, the intervention is cost-effective if and only if $\lambda * 0.1 - 1000 > 0$. \\ \\ The decision process in the preceding example can be formalized by working with the the incremental net benefit (INB): $$ \text{INB}_{\lambda} \coloneqq \lambda \Delta_q - \Delta_c. $$
By putting a monetary value on health ($\lambda$) and checking if the benefits of the intervention outweigh its additional costs, we can decide wether the additional health we get from the intervention is worth the extra cost. More precisely: if the INB is zero or negative, we conclude that the intervention is not cost-effective and should hence not be implemented. Conversely, if the INB is positive, we consider the intervention to be cost-effective and worthy of implementation.
In most applications, the true values of $\Delta_c$ and $\Delta_q$ are unknown, and so is the true value of $\text{INB}_{\lambda}$ for any given $\lambda$. They thus must be estimated from available data, and there will inevitably be uncertainty about these estimates. In this paper, we adopt the Bayesian point of view when conceptualizing this uncertainty, which means that we treat all unknown quantities as random variables. This allows us to speak about the probability of cost-effectiveness\footnote{In the frequentist or "classical" paradigm, this is not possible, and we would instead have to argue in terms of p-values. See lothgren2000definition.}, which is given by $\Pr(\text{INB}_{\lambda} > 0)$, the probability that the INB is positive given a specific $\lambda$ value.
While there are many different methods to estimate $\Delta_c$ and $\Delta_q$, they are all ultimately based on empirical observations of patient's cost and health. These observations are often positively or negatively correlated, as illustrated by the following examples.
A proper statistical analysis should account for this correlation, and doing so will also induce correlation in the distribution of $\Delta_c$ and $\Delta_q$. The following example illustrates how this correlation can affect the results. \\ \\ Example 2. We now consider a situation where $\Delta_c$ and $\Delta_q$ are estimated from available data, and we are hence uncertain about their values. This uncertainty is expressed through a probability distribution; an example of which is illustrated in Figure (ref). The left column shows two bivariate density plots of $\Delta_q$ and $\Delta_c$, while the right column shows the probability of cost-effectiveness as a function of $\lambda$. These plots are called the cost-effectiveness plane (CE-plane) and the cost-effectiveness acceptability curve (CEAC), respectively. While the mean and variance of $\Delta_c$ and $\Delta_q$ are the same in both rows, the correlation between them differs: it is negative in the top row and positive in the bottom row. This difference in correlation leads to a noticeable change in the CEACs and, consequently, in how the results are interpreted. For example, at a willingness-to-pay of $50,000 \text{\euro}$ per QALY gained the probability of cost-effectiveness is approximately $0.7$ in the top row, but nearly $1$ in the bottom row.
When cost-effectiveness data are obtained from a randomized trial, model specification is not critical. By design, the patients' baseline characteristics are unrelated to which treatment the patient receives, and we hence need not adjust for them; simply comparing mean costs and effects in the treatment groups will produce good estimates of the treatment effects. If a researcher insists on adjusting anyway, a simple linear model will perform well lin2013agnostic. Without randomization, however, there will generally be certain baseline characteristics, called confounders, which influence both the treatment choice and the outcomes (i.e. cost and health) and hence need to be adjusted for in the analysis. This adjustment necessitates the use of a statistical model, as we need to specify how exactly the confounders influence the outcomes. Even with relevant confounders measured, a badly chosen model can lead to biased and inconsistent estimates of the treatment effects aronow2025nonparametric, making model choice a critical matter when analysing observational data.
With that in mind, we will now critically examine some of the commonly used model choices in CEAs. A frequently used model is seemingly unrelated regression (SUR, zellner1962efficient, ben2023conducting), which asserts that
With SUR, we model the patients' costs $c_i$ and health $q_i$ as linear functions of a set of covariates $\mathbf{x_i}$, consisting the patients' baseline characteristics or some transformation thereof (so it may include, say, quadratic or interaction terms), plus a bivariate-normally distributed distributed error term. In a CEA, both the linearity and the normality assumption are questionable, and we shall now scrutinize each of them in turn.
The normality assumption is rarely realistic in the CEA setting, particularly for the cost outcome, which can only take non-negative values. However, violations of this assumption are not critical: in a landmark paper, white1980heteroskedasticity showed that, as long as certain mild conditions are met, the estimates from linear regression models (such as SUR) remain unbiased and consistent\footnote{Consistency means, roughly speaking, that with a sufficiently large sample the estimate will be close to the true value.}, even when the distribution of the error term is misspecified. Nonetheless, the normality assumption can be weakened by through the use of bootstrapping, as is frequently done in CEA lothgren2000definition.
Another way to avoid the normality assumption, which has been promoted by many authors but not applied much in practice, is the use of generalized linear models (GLMs). Those allow the specification of non-normal error terms—for example, by modeling costs using a lognormal or gamma distribution baio2018statistical,gabrio2019bayesian. However, this approach has a significant drawback: the results can be very sensitive to the choice of the distribution. A valuable simulation study by briggs2005parametric showed that, if a lognormal distribution is assumed for the costs, but the true distribution is gamma, the estimated mean costs are strongly biased. The sample mean on the other hand, which is essentially a linear regression without covariates, performs well regardless of the true distribution. The same patterns will persist in the more observational setting, with adjustment for baseline covariates. The reason is that, for GLMs with non-normal error terms, there are no theoretical guarantees analogous to those in linear models — that is, if the error term in a GLM is not specified correctly, the resulting estimates may be inaccurate even with large sample sizes. In consequence, choosing the wrong distribution for the outcome can lead to both biased and inconsistent treatment effect estimates.
The more fundamental problem with the SUR model is the assumed linearity. Violations of this assumption can lead to severe bias in the estimated treatment effects dorie2019automated; we illustrate this through a simple simulation experiment in Appendix (ref). GLMs suffer from analogous problems. For example, a GLM for costs often assumes that the logarithm of the expected costs depends linearly on the covariates baio2018statistical; an assumption just as restrictive as the linearity in the SUR model. More generally, both SUR and GLMs are examples of parametric models, which make strong assumptions about how the covariates are related to the outcomes. Violations of these assumptions lead to distorted estimates of treatment effects dorie2019automated.
In order to deal with the aforementioned problem, statisticians have developed regression techniques which do not make strong parametric assumptions about the relationship between the covariates and the outcomes. These so called nonparametric regression models\footnote{The term "nonparametric" is somewhat of a misnomer, since such models do in fact have parameters. However, the number of these parameters (or more technically: the dimension of the parameter space) is not chosen in advance, unlike in, say, the linear model. Nonetheless, the terminology is widely used.} are often much more accurate than linear models or GLMs. Crucially, under similar regularity conditions, consistency guarantees analogous to those of linear regression, also exist for nonparametric regression models, even when the normal error distribution is misspecified kleijn2006misspecification. These added benefits come at a cost, however, as nonparametric regression methods can be more challenging to use and interpret. In particular, unlike in the basic linear model, the treatment effects of interest are no longer represented by a single model parameter, making their interpretation les straightforward. We thus intend to make nonparametric regression accessible to applied cost-effectiveness researchers.
We will focus on a specific nonparametric regression model, namely Bayesian additive regression trees (BART), which has been shown to perform well at estimating treatment effects dorie2019automated. We have recently extended this model to account for the correlation between costs and health outcomes found in CEA data; this extended version is called suBART esser2024seemingly
As the name suggests, BART is a Bayesian statistical model and works by aggregating a set of regression trees. We shall now explain briefly what each of these terms mean.
While there regression model come in many forms, they all share the common goal of describing the relationship between an outcome variable and one or more covariates. A regression tree does this by dividing the covariate space (the set of all possible values of $\mathbf{x}_i$) into distinct "bins". Then each of these bins is assigned a single value, representing the predicted or expected outcome for all observations whose covariates fall within that bin. We illustrate this through a simple example. \\ \\ Example 3. A representation of a simple regression tree is shown on the left of Figure (ref). We may imagine this tree giving us the predicted QALY for a particular patient, based on their baseline health state (expressed as a number between 0 and 1). Given a specific value for this baseline utility --- say $x$ --- how do we find the predicted QALY? We start at the top of the tree, and then move down the levels sequentially. At each step, we follow either the left or the right direction, depending on whether $x$ is smaller or larger than the shown threshold value. Once we reach the end of a branch (where the tree no longer splits), we take the given value as the predicted QALY. We can also use the tree to plot the expected QALY as a function of $x$. This is illustrated in the right panel of Figure (ref). Where a linear model would have given us a line, the regression tree produces a so-called step function. \\ \\
Regression trees are to some extent analogous to the decision tree models which health economists sometimes use in simulation experiments (see for example briggs2006decision, chapter 2). However, among other differences, a regression tree is not a priori specified by the researcher, but rather "learned from the data".
The basic idea of BART is to express the outcome as a sum of regression trees plus a normally distributed error term\footnote{As noted above, such a normally distributed error term does not pose issues even when the conditional distribution of the outcomes is not normal.}: \[y_i = \sum^m_{t=1} f_i\left(\mathbf{x}_i\right) + \varepsilon_i,\] where each $f_i$ is a regression tree; each of which captures different nonlinear effects of the covariates on the outcome. By constraining the trees to be small, BART effectively controls overfitting, which in turn improves generalization. This framework allows for modelling complex nonlinear relationships between covariates and outcome without relying on strong parametric assumptions. The total number of trees $m=200$ is frequently used, since this choice has been found to provide good performance while keeping the computational demands of the model fitting manageable. See chipman2010bart,tan2019bayesian for detailed treatments of the original BART model, and particularly the process by which the trees are learned from the data.
The acronym suBART stands for seemingly unrelated BART, a name deliberately reminiscent of the SUR model we described in Section (ref). SUR consists of two linear models which are coupled through correlated error terms. suBART keeps the latter, but replaces the linear models with additive regression trees. More precisely, the model is given by:
where each $f_i$ and $g_i$ is a regression tree.
As the name suggests, BART is a Bayesian statistical method, which means that the model parameters --- and functions of these parameters, such as $\Delta_c$ --- are treated as random variables, rather than unknown constants. We learn about these parameters by examining their posterior distribution, which quantifies how different values of the parameters are more or less likely, taking into account both prior beliefs and the evidence provided by the data.
Bayesian estimation has a long tradition in CEA, with some authors arguing that it is the most natural framework to perform cost-effectiveness research under baio2018statistical. The main reason is that cost-effectiveness researchers are concerned with the probability of cost-effectiveness which is nothing other than $\Pr(\lambda \Delta_q - \Delta_c > 0)$. In other words, we want to make probabilistic statements about our parameters of interest, while explicitly acknowledging uncertainty about their values, and this -- strictly speaking --- is only possible under the Bayesian paradigm.
In addition, there are pragmatic considerations that motivate our use of Bayesian estimation. It allows for the use of complex statistical models that are computationally difficult to estimate using alternative approaches. Morever, Bayesian estimation often performs well even when assessed by frequentist or other non-Bayesian criteria. BART is a prime example, as shown through the results of the American Causal Inference Conference (ACIC) competition dorie2019automated.
In practice, Bayesian estimation works by taking a large number of random draws from the posterior distribution. We then use these draws to approximate our quantities of interest, such as the aforementioned probability of cost-effectiveness or confidence intervals for the treatment effects. This approximation will be quite accurate if the number of draws is sufficiently large. One popular procedure to obtain these draws is is called Markov chain Monte Carlo (MCMC), and the collection of draws is called an MCMC sample. While superficially similar to bootstrapping, which is a commonly used tool in CEA, the underlying theory of MCMC is quite different. While bootstrap samples are resampled datasets used to approximate sampling variability in a frequentist framework, MCMC samples represent draws from the posterior distribution and quantify uncertainty in the parameters conditional on the observed data. In practice, however, the use of MCMC is quite similar to that of bootstrapping. For example, to construct a 90% credible interval\footnote{The Bayesian analogue to the frequentist confidence interval.}, we can simply take the 5th and 95th percentiles of the MCMC sample. Also analogously to bootstrapping, there are no strict guidelines for what the aforementioned "large number" of draws should be. The usual choice is a few thousand, and we will use $4000$ in the example below. For more information on MCMC, we refer to greenberg2012introduction, Chapter 7. Details on the specific MCMC method for the suBART model may be found in esser2024seemingly.
In the causal inference literature, one frequently encounters the so-called propensity score (PS), which is the probability that a patient received the treatment, given their baseline characteristics. The PS summarizes all the confounding variables into a single number, and its inclusion as an additional covariate helps control for confounding more effectively hahn2020bayesian,dorie2019automated. We therefore begin by estimating the PS for each patient in the data, using a BART model to predict treatment assignment based on baseline characteristics. This estimated PS is then included as an additional covariate when fitting the suBART model for the outcomes. This approach has been shown to perform well—and better than a suBART model without PS adjustment—in the context of observational CEAs esser2024seemingly.
In this section we guide the reader through the use of suBART. We use a simulated dataset, which is designed as an example and not meant to represent any specific real-world research endeavour. In the supplementary material\footnote{See \url{https://github.com/Jonas-Esser/nonparametric_CEA}}, we provide the full dataset, together with all of the following code. The dataset contains the treatment indicator t (which indicates which treatment a patient received), the two outcome variables (costs and health --- represented by c and q, respectively) and the confounding variables age, sex, and education. For illustration, we display a subset of the data in Table (ref).
Some observations in the outcome variables are missing. The suBART method can handle these automatically, provided we assume that the data are missing at random ---that is, the probability of a value being missing does not depend on the value itself, but may depend on other observed variables (such as confounders).
In order for the suBART code to work correctly, it is important that the data are formatted properly. Unordered categorical variables with more than two levels (i.e., variables that represent categories without any natural order, such as blood type, region, or treatment center) should be declared as factors. In the dataset, this applies to the education variable. All other variables should be declared as numbers. Furthermore, all missing observations should be properly formatted as NA values, not as text strings or numbers such as $-99$, which is occasionally observed in practice.
We put all covariates into a data frame called X, including the treatment indicator but not the outcomes c and q. The X object should contain only those variables which are supposed to enter into the regression model, and no others. Likewise, we put the outcomes c and \texttt{q} into their own data frame \texttt{Y}.
We proceed to estimate the PS using a BART model, where treatment assignment is modeled as a function of the previously defined confounders. The number of trees and MCMC iterations was specified earlier. In line with standard practice, we discard the first 1000 posterior draws (usually referred to as burn-in samples) following standard practice in the MCMC literature. This means that although we set the total number of MCMC iterations (n_mcmc) to $5000$ in the code, only the remaining $4000$ samples after burn-in are used for inference.
Once estimated, we append the PS to our covariate data frame, X. For clarity and to keep the original data intact, we store this extended data frame in a new object called X_ps.
As mentioned above, in nonlinear models such as BART, treatment effects estimates are not represented by a single model parameter, unlike with the simple linear model without interaction terms. Instead, they are derived from the model predictions as follows: for each patient, we estimate the expected outcomes under both treatment options and then calculate the average difference between these expected outcomes. The procedure is commonly known as g-computation hernan2024causal.
In algorithmic form, our strategy may be outlined as follows:
These steps will be performed for every MCMC draw. For the sake of step $3$, we create another data frame, which contains each patient twice: once with the treatment indicator set to $0$, and once with it set to $1$. This data frame is called X_test in the code.
We now fit the model, using the previously constructed data frames.
The resulting suBART_fit object contains a variety of objects; we only need the entry suBART_fit\$y_hat_test, which corresponds to the the predictions from Step (ref).
In this step, we apply the calculations from Section (ref) to all draws in the MCMC sample. We use the resulting draws of $\Delta_c$ and $\Delta_q$ to estimate $\Pr\left(\text{INB}_{\lambda} > 0\right)$ for a specified range of $\lambda$ values, which will allow us to compute and plot the cost-effectiveness acceptability curve (CEAC), as described by willan2006statistical.
Finally, we plot the cost-effectiveness plane and CEAC, using the ggplot2 package wickham2016ggplot.
In this tutorial, we demonstrated how to perform statistical estimation in CEAs with observational data. Our approach avoids strong parametric assumptions on the relationship between treatment and baseline variables with the expected costs and health state of a patient, which is very important when data is not obtained from a randomized trial. While we do not explicitly address the common issue of non-normal error terms, we argued that our approach is likely robust to violations of the normality assumption.
We used the suBART model in our presentation, both because of its strong empirical performance in this context and its methodological advantages. We stress, however, that our primary goal was not to present this model, but rather to provide a general framework for cost-effectiveness researchers analyzing observational data. suBART is merely an example, and there are alternatives which can be used in a similar way, and which we also expect to perform well. In principle, any non-Bayesian method for estimating $\Delta_c$ and $\Delta_q$ can be used --- for example, the Causal Forest of wager2018estimation or Double Machine Learning knaus2022double. The MCMC sample referenced in Section (ref) would then be replaced by a bootstrap sample, and the rest of the analysis be carried out analogously. However, such an analysis would be no longer Bayesian, and hence the interpretation of the results changes: strictly speaking, we could no longer speak of the "probability of cost-effectiveness", and must instead argue in terms of frequentist p-values; see lothgren2000definition. Additionally, such approaches may be infeasible in practice due to computational constraints.
We now discuss some limitations of our presented approach, relating to missing data, and outline possible solutions to address them.
Despite these limitations, the presented framework is sufficient for use in many CEAs with observational data and produces more credible results than commonly used parametric models. We thus hope that readers will find it useful in their research.
No funding was received for the completion of this paper.
All authors have no conflicts of interest to declare.
All data and code used are available under \url{https://github.com/Jonas-Esser/nonparametric_CEA}.
Not applicable.
The initial idea for the paper was proposed by Jonas Esser; the scope and design were then worked out in cooperation with Johanna van Dongen and Judith Bosmans. The text was written primarily by Jonas Esser, while incorporating suggestions from all other authors. The R code, figures and computer experiments were created by Jonas Esser and Mateus Maia.
\interlinepenalty=1000
\setcounter{subsection}{0}