A Bayesian Modeling Framework for Health Care Resource Use and Costs in Trial-Based Economic Evaluations

Abstract

Individual-level data are routinely used in trial-based economic evaluations to assess the effectiveness and costs of a given intervention. While effectiveness measures are often expressed via utility scores derived from health-related quality-of-life instruments (e.g., EQ-5D questionnaires), information on different types of health care resource use (HRU) measures (e.g., number and types of services) are collected to compute the costs. Partially complete HRU data, particularly for self-reported questionnaires, are handled via ad hoc methods that rely on some assumptions (fill in a zero) that are typically hard to justify. Although methods have been proposed to account for the uncertainty surrounding missing data, particularly in the form of multiple imputation or Bayesian methods, these have mostly been implemented at the level of costs at different times or over the entire study period, while little attention has been given to how missing values at the level of HRUs should be addressed and their implications on the final analysis. We present a general Bayesian framework for the analysis of partially observed HRUs in trial-based economic evaluations, which can accommodate the typical complexities of the data (e.g., excess zeros, skewness, missingness) and quantify the impact of missingness uncertainty on the results. We show the benefits of our approach with a motivating example and compare the results to those from more standard analyses fitted at the level of cost variables after adopting some ad hoc imputation. This article highlights the importance of adopting a comprehensive modeling approach to handle partially observed HRU data in economic evaluations and the strategic advantages of building these models within a Bayesian framework.

Highlights

• Missing health care service data in trial-based economic evaluations are often removed or imputed using quite restrictive assumptions (e.g., no use of service).

• We propose a flexible Bayesian approach to account for missing health care service uncertainty and compare the results with models fitted at more aggregated levels (e.g., total costs) using a real case study.

• Our results show that, depending on the (assumed) missingness assumptions and the level of data aggregation at which analyses are performed, results may be considerably changed.

• When feasible, analyses should be conducted at the most disaggregated level to ensure that all available information collected in the trial is used in the analysis without relying on (often) restrictive ad hoc imputation approaches.

Economic evaluations represent an increasingly important component in the process of health technology assessment (HTA) of new health care interventions in many countries.^1,2 In the context of trial-based economic evaluations, HTA institutions recommend cost-utility analysis as the “reference case” approach, using health-related quality-of-life (HRQol) utility scores as the health outcome to ensure direct comparability of interventions across disease areas. Alongside utilities, costs are computed by applying the appropriate unit prices, retrieved from the literature or national sources, to the amount of health care resource use (HRU) consumed by the patients over a certain period of time (e.g., number and type of visits over the past 3 mo). Self-reported questionnaires, such as the Client Service Receipt Inventory ³ and iMedical Cost Questionnaire, ⁴ are routinely used to collect HRU data but often suffer from low completion rates.

Previous reviews have advocated for the importance of dealing with missing outcome data in trial-based economic evaluations^5,6 and have highlighted an overall improvement over time in current practice toward the adoption of approaches that formally take into account the uncertainty around missing values, such as multiple imputation. ⁷ A more recent review provided an updated analysis of the missingness approaches used in economic evaluations and revealed a somewhat unclear picture about the consistency of the methods used to address missingness at different levels of aggregation. ⁸ More specifically, the authors highlighted a clear difference between the typical methods used to handle missing values at the level of cost/utility data, where multiple imputation is often used, compared with the level of questionnaire items, where some form of single imputation approach is used. While this may not pose a problem for EQ-5D questionnaires, often characterized by unit nonresponse (all answers skipped), the same cannot be said for HRU questionnaires, which are typically lengthy and subject to item nonresponse (some answers skipped). As a result, when confronted with partially missing HRU responses, analysts may decide to “fill in the gaps” based on some specific assumptions such as no use of service (i.e., imputed as 0). Although this can be justified in some cases based on external information, when this is lacking, it is essential that the uncertainty around missing HRU data is properly quantified and its impact on the analysis results reflected.

Building on the existent literature, we propose a novel Bayesian framework for trial-based economic evaluations that allows a flexible model specification and the handling of missing data at the HRU level. The choice of a Bayesian approach has practical advantages compared with a standard frequentist framework,^9,10 including 1) the use of a modular structure to increase model complexity in a relatively easy way, ¹¹ 2) natural interpretation of cost-effectiveness results in probabilistic terms, ¹² and 3) direct implementation of probabilistic sensitivity analysis, consisting in the quantification of the impact of parameters’ uncertainty on the conclusions. ¹³ We show the benefits of using our framework on a real case study, with a focus on appropriately modeling partially observed HRU values and its implications in terms of inferences and, crucially, cost-effectiveness results.

The article is structured as follows: the second section presents the case study and describes the data. The third section defines the general framework of the statistical models used to analyze costs/utilities and HRU values. The fourth section compares the results of models fitted at different levels of aggregation and shows how different assumptions about partially missing cases may affect the model estimates. The fifth section performs the economic evaluation, summarizes the inferences for each model, and compares the cost-effectiveness results. The sixth section discusses the proposed framework and suggests some improvements for future work. Finally, the Appendix includes additional material related to model assessment and results, while the software code is provided in the online supplementary material.

Case Study: The PBS Trial

The Positive Behaviour Support (PBS) ¹⁴ trial was a multicenter, randomized controlled trial involving community intellectual disability services and service users with mild to severe intellectual disabilities and challenging behaviors. Positive behavior support is a multicomponent intervention that is designed to foster prosocial actions and enhance the person’s quality of life and their integration within the local community. Participants (N = 244) were enrolled and randomly allocated to staff teams trained to deliver PBS in addition to treatment as usual (reference intervention, $n_2=108$ ) or to staff teams trained to deliver treatment as usual alone (comparator, $n_1=136$ ). Measures for quality of life in the form of EQ-5D-3L questionnaires ¹⁵ and HRU information based on Client Service Receipt Inventory (CSRI) questionnaires ³ were collected for each individual ( $i=1,\dots ,n$ ) at baseline ( $j=0$ ) and at $6$ and $12$ ( $j=1,2$ ) mo of follow-up.

Individual-level utility scores at each time $u_{\mathrm{ij}}$ were computed based on the patients’ answers to the EQ-5D questionnaires and national value sets, ¹⁶ while HRU data were collected through CSRI questionnaires and covered a wide range of health care services. For the purpose of this study, we consider individual-level HRUs collected at each time ${\mathrm{HRU}}_{\mathrm{ij}}^k$ on $k=1,\dots ,K=9$ services: number of visits with a psychiatric doctor (PSYDR), psychologist (PSYCH), physiotherapist (PHYSI), dentist (DENT), social worker (SOCWORK), community worker (COMWORK), general practitioner (GP), nurse (NURSE), and private therapist (THERAP). CSRI data were characterized both by unit and item nonresponse, meaning that individuals may either provide full or only partial information in terms of the services used at each time. A statistical summary of the missingness rates of the utility and HRU variables by time point in the PBS trial is shown in Table 1.

Table 1

Number (Proportion) of Missing Cases for the Utilities and Different Health Care Resource Use Services (Self-Reported Questionnaires), Presented by Time Point of the Trial ^a

Outcome	Baseline (j = 0)	6 mo ( $\mathit{j}=1$ )	12 mo ( $\mathit{j}=2$ )	ic (Outcome)
Utilities ( $u$ )	14 (5.7%)	23 (9.4%)	16 (6.6%)	40 (16.4%)
PSYDR ( $\mathrm{hr}{\mathrm{u}}^1$ )	0	18 (7.4%)	12 (4.9%)	22 (9%)
PSYCH ( $\mathrm{hr}{\mathrm{u}}^2$ )	2 (1.4%)	12 (4.9%)	9 (3.7%)	14 (5.7%)
PHYSI ( $\mathrm{hr}{\mathrm{u}}^3$ )	1 (0.4%)	13 (5.3%)	10 (4.1%)	16 (6.6%)
DENT ( $\mathrm{hr}{\mathrm{u}}^4$ )	1 (0.4%)	14 (5.7%)	12 (4.9%)	19 (7.8%)
SOCWORK ( $\mathrm{hr}{\mathrm{u}}^5$ )	1 (0.4%)	14 (5.7%)	11 (4.5%)	18 (7.4%)
COMWORK ( $\mathrm{hr}{\mathrm{u}}^6$ )	0	13 (5.3%)	11 (4.5%)	15 (6.1%)
GP ( $\mathrm{hr}{\mathrm{u}}^7$ )	3 (2.1%)	15 (6.2%)	14 (5.7%)	22 (9%)
NURSE ( $\mathrm{hr}{\mathrm{u}}^8$ ) 2	(1.4%)	17 (7%)	13 (5.3%)	24 (9.8%)
THERAP ( $\mathrm{hr}{\mathrm{u}}^9$ )	1 (0.4%)	15 (6.1%)	12 (4.9%)	19 (7.8%)
ic (time)	10 (4.2%)	22 (9%)	17 (7%)	59 (24.2%)

a

In the table, the number (%) of individuals with partially observed data at each time across outcomes (bottom row), for each outcome across times (rightmost column), and across both times and outcomes (bottom right cell) is also reported. At each time j, only unit nonresponse characterizes EQ-5D questionnaires (and so utilities), while both unit and item nonresponse affects Client Service Receipt Inventory questionnaires (and so health care resource uses).

Overall, the proportions of missing values are moderate and never exceed $10\%$ at any time point. However, when aggregating all variables, missingness becomes more substantial, with a proportion of incomplete cases of $\approx 24\%$ (bottom right cell in Table 1).

In addition to missingness, a typical feature of HRU data is the amount of “extreme” values that remain constant over time (e.g., zero values), typically denoted as “structural values.” ¹⁷ Table 2 reports the number (proportions) of structural zeros observed for each type of HRU service in the PBS study.

Table 2

Number (Proportion) of Structural Zeros Observed for Each Type of HRU Service in the PBS Study

Outcome	Structural Zeros
PSYDR ( $\mathrm{hr}{\mathrm{u}}^1$ )	37 (15%)
PSYCH ( $\mathrm{hr}{\mathrm{u}}^2$ )	159 (65%)
PHYSI ( $\mathrm{hr}{\mathrm{u}}^3$ )	189 (77%)
DENT ( $\mathrm{hr}{\mathrm{u}}^4$ )	35 (14%)
SOCWORK ( $\mathrm{hr}{\mathrm{u}}^5$ )	47 (19%)
COMWORK ( $\mathrm{hr}{\mathrm{u}}^6$ )	196 (80%)
GP ( $\mathrm{hr}{\mathrm{u}}^7$ )	1 (0.4%)
NURSE ( $\mathrm{hr}{\mathrm{u}}^8$ )	53 (22%)
THERAP ( $\mathrm{hr}{\mathrm{u}}^9$ )	85 (35%)

HRU, health care resource use; PBS, Positive Behaviour Support.

The proportions of structural zeros in most HRUs range from $15\%\mathrm{to}80\%$ , with the only exception being the number of GP visits ( $0.4\%$ ). The presence of these zero values induces a high degree of skewness in each of these variables’ empirical distributions, which needs to be appropriately addressed in the analysis stage.

Modeling Framework

We first present our framework building upon previously introduced approaches. We start by considering models fitted at the level of total costs and quality-adjusted life-years (QALYs)¹⁸⁻²¹ or costs and utilities at different times^22,23 and then extend these to account for partially observed HRU data. The framework allows for typical features of cost-effectiveness data (e.g., skewness, correlation, structural values) while also handling missing HRUs without relying on ad hoc assumptions (e.g., assumed zeros). Although different study designs in terms of data collection could be considered (e.g., multiple times for the effectiveness but only 1 for costs), we argue that the PBS data reflect the typical design used in many trial-based economic evaluations. ²⁴ Thus, we illustrate the framework using the trial design (see the second section) as reference, where disaggregated outcomes are collected at baseline and equally spaced follow-up times. However, we note that the framework can also be modified to handle different types of data structures, including different collection times between outcomes and/or a cross-sectional design for some variables.

Assume that some patient-level HRU and HRQol data from a trial are collected via self-reported questionnaires (e.g., CSRI and EQ-5D) at equally spaced times $j=1,\dots ,J$ on $i=1,\dots ,n$ individuals, who are randomly allocated to either a control ( $t=0$ ) or intervention ( $t=1$ ) group with sample sizes $n_1$ and $n_2$ , respectively. For each individual and time point, questionnaire answers on multiple health care services and HRQol domains are collected. Let ${\mathrm{HRU}}_{\mathrm{ij}}^k$ denote the health resource use information collected on service $k$ (for $k=1,\dots K$ ) for the $i$ -th individual and $j$ -th time point in the trial, and let ${\mathrm{HrQol}}_{\mathrm{ij}}^l$ denote the corresponding HRQol information collected on domain $l$ (for $l=1,\dots ,L$ ). HRU and HRQol data are then, respectively, combined with national unit prices for each type of service $p_k$ and value sets to derive service-specific costs $c_{\mathrm{ij}}^k$ and utility scores $u_{\mathrm{ij}}$ for each individual at each time in the study, with the total individual costs obtained as the sum of all service-specific costs at each time, that is, $c_{\mathrm{ij}}=\sum _{k=1}^K{c}_{\mathrm{ij}}^k$ . Next, total costs ( $t{c}_i$ ) and QALYs ( $e_i$ ) for each individual over the trial period are computed as

$$\begin{array}{ccc}t{c}_i=\sum _{j=1}^J{c}_{\mathrm{ij}}\hfill & \mathrm{and}\hfill & e_i=\sum _{j=1}^J(u_{\mathrm{ij}}+u_{\mathrm{ij}-1})\frac{{\delta }_j}{2},\hfill \end{array}$$

(1)

where ${\delta }_j=\frac{{\mathrm{Time}}_j-{\mathrm{Time}}_{j-1}}{\mathrm{Unit}\mathrm{of}\mathrm{time}}$ is the proportion of the time unit (typically 1 y) that is covered between time $j-1$ and $j$ in the trial.

The framework allows the modeling of different types of variables: HRU and HRQoL, costs and utilities at each time, or total costs and QALYs. In principle, it is arguably easier to specify a model for aggregated (e.g., total costs) compared with disaggregated variables (e.g., costs or HRU), since it allows for performing the analysis based on a fewer number of variables and a cross-sectional setting. However, in reality, the occurrence of missing values at the most disaggregated level (i.e., HRU and HRQoL) is unavoidable, with proportions of unobserved values that are usually substantial.^5,6 In the presence of partially missing questionnaire answers, modeling at aggregated levels (e.g., costs or total costs) does not allow direct incorporation of the evidence from the partially observed items, which need to be discarded unless their values are imputed prior to the analysis. While item nonresponse is often not a problem for HRQoL, ²⁵ it severely affects HRU data, with many individuals providing information about the usage of only some of the services. ²⁶ To our knowledge, guidelines on the handling of partially missing HRU data are not available, with analysts often not reporting detailed information about their proportions of missing values. ⁸ When information on HRUs is lacking, it is plausible to assume that some deterministic imputation of item nonresponses is used to improve data completeness. A typical example occurs when individuals with a missing value on a specific service are assumed to not have used it, thus suggesting the replacement of the unobserved value with a zero. However, such assumption can never be verified based on the observed data and, when repeatedly applied for different services, times, and individuals, may have a substantial effect on the analysis results.

We separately present our approach by level of data aggregation: 1) total costs and QALYs $(e_i,t{c}_i)$ , 2) costs and utilities computed at each time $(u_{\mathrm{ij}},c_{\mathrm{ij}})$ , and 3) HRUs and utilities computed at each time $(u_{\mathrm{ij}},{\mathrm{HRU}}_{\mathrm{ij}}^k)$ . For each level of aggregation, we show how the model can be specified and compared with simpler approaches to handle partially missing cases prior to model fitting. For illustrative purposes, and given that item nonresponse is unlikely to affect HRQoL questionnaires, in the reminder we will assume that item nonresponses affect only HRU data.

Modeling of Total Costs and QALYs

At the most aggregated level, data are summarized into effectiveness and total costs ( $e_i,t{c}_i$ ), which are computed over the study period using equation 1. Most of the methodological literature has focused on the modeling of these variables with an emphasis on the need to specify flexible models to handle the typical complexities of the data. These include correlation between the outcomes, skewness, and presence of structural values.^{20,21,27–30} For instance, dependence between outcomes can be captured by specification of the joint distribution $p(e_i,t{c}_i)$ as

$$p(e_i,t{c}_i)=p\left(e_i\right)p(t{c}_i\mid e_i)=p\left(t{c}_i\right)p(e_i\mid t{c}_i),$$

(2)

where, for example, $p\left(e_i\right)$ is the marginal distribution of the effectiveness and $p(t{c}_i\mid e_i)$ is the conditional distribution of the total costs given the effectiveness. Note that, while it is possible to use interchangeably either factorization in equation 2, without loss of generality, we describe our analysis in the following through a marginal distribution for the effectiveness and a conditional distribution for the costs. Choice of the distribution for $e_i$ and $t{c}_i$ can be guided by the specific characteristics of the data, for example, gamma or lognormal to handle skewness, with each distribution being indexed by a set of location $({\varphi }_{\mathrm{ie}},{\varphi }_{\mathrm{itc}})$ and nuisance $({\psi }_e,{\psi }_{\mathrm{tc}})$ parameters. A generalized linear structure based on some link functions $g(·)$ , such as log or logit, is often used to incorporate outcome-specific covariates at the location level:

$$\begin{array}{ccc}{g}_e\left({\varphi }_{\mathrm{ie}}\right)={\alpha }_0+{\alpha }_1{X}_{\mathrm{ie}}\hfill & \mathrm{and}\hfill & g_{\mathrm{tc}}\left({\varphi }_{\mathrm{itc}}\right)={\beta }_0+{\beta }_1{X}_{\mathrm{itc}},\hfill \end{array}$$

(3)

where the sets of regression parameters $\alpha$ and $\beta$ include the intercept and covariate-specific coefficients ( $X_{\mathrm{ie}},X_{\mathrm{itc}}$ ) for each model. Posterior estimates of interest, that is, the mean incremental QALYs ( ${\mathrm{\Delta }}_e={\mu }_{e0}-{\mu }_{e1}$ ) and total costs ( ${\mathrm{\Delta }}_{\mathrm{tc}}={\mu }_{\mathrm{tc}0}-{\mu }_{\mathrm{tc}1}$ ) between treatment groups, can be derived in terms of linear combinations of the regression parameters or through simulation approaches such as Monte Carlo methods.^9,10 Figure 1 shows a graphical representation of the modeling framework for the aggregated variables, where the effectiveness and total cost distributions are represented in terms of combined “modules”—the blue and the red boxes—in which the random quantities are linked through logical relationships.

Figure 1

Joint distribution, expressed in terms of a marginal distribution for the effectiveness and a conditional distribution for the total costs, respectively, indicated with a blue and red line. The parameters indexing the corresponding distributions or “modules” are indicated with different Greek letters, while the black and magenta arrows show the dependence relationships between the parameters within and between the 2 models, respectively. Note that $i$ denotes the individual index, while the treatment index $t$ is omitted to ease notation.

An appealing feature of this strategy is the need to specify only 2 distributions. However, when HRU nonresponse occurs, analysts need to handle the partially observed disaggregated data prior to fitting the model. As a simple solution, they could either discard all cases with partially observed HRUs or impute as zero all HRU nonresponses ( ${\overline{\mathrm{HRU}}}_{\mathrm{ij}}^k=0$ ) and use these to generate costs at each time $c_{\mathrm{ij}}$ . In the second case, analysts must also decide how to handle missing cases at the cost level by either discarding all cases with partially observed costs at any time or impute as zero all partially observed costs at each time ( ${\overline{c}}_{\mathrm{ij}}=0$ ) and use these to generate total costs $t{c}_i$ .

Modeling of Costs and Utilities at Each Time

At an intermediate aggregation level, the data consist of 2 longitudinal variables, namely, the utilities and costs computed at each time ( $u_{\mathrm{ij}},c_{\mathrm{ij}}$ ). In recent years, attention has been given to the specification of longitudinal models that can deal with missing outcomes at different times while also addressing the complexities of the data.^22,23 For example, the joint outcome distribution $p(u_{\mathrm{ij}},c_{\mathrm{ij}})$ at time $j>0$ can be specified as

$$p(u_{\mathrm{ij}},c_{\mathrm{ij}})=p(u_{\mathrm{ij}}\mid u_{i,j-1},c_{\mathrm{ij}})p(c_{\mathrm{ij}}\mid u_{\mathrm{ij}-1},c_{i,j-1}),$$

(4)

where $p(u_{\mathrm{ij}}\mid u_{i,j-1},c_{\mathrm{ij}})$ is the conditional distribution of the utilities at $j>0$ given utilities at time $j-1$ and costs at $j$ , while $p(c_{\mathrm{ij}}\mid u_{\mathrm{ij}-1},c_{\mathrm{ij}-1})$ is the conditional distribution of the costs at $j>0$ given utilities and costs at $j-1$ . Similarly to the modeling of aggregated variables, distributions for $u_{\mathrm{ij}}$ and $c_{\mathrm{ij}}$ are indexed by a set of time-specific location $({\varphi }_{\mathrm{iju}},{\varphi }_{\mathrm{ijc}})$ and nuisance $({\psi }_{\mathrm{ju}},{\psi }_{\mathrm{jc}})$ parameters and should be chosen according to the specific features of the data. When a generalized linear structure based on some link functions $g_u(·)$ and $g_c(·)$ is used to incorporate outcome-specific covariates at the location level, the conditional mean at time $j$ for the utilities and costs can be expressed as

$$\begin{array}{ccc}{g}_u\left({\varphi }_{\mathrm{iju}}\right)={\alpha }_{j0}+{\alpha }_{j1}{X}_{\mathrm{iu}}\hfill & \mathrm{and}\hfill & g_c\left({\varphi }_{\mathrm{ijc}}\right)={\beta }_{j0}+{\beta }_{j1}{X}_{\mathrm{ic}},\hfill \end{array}$$

(5)

where the sets of regression parameters $\alpha$ and $\beta$ include the time-specific intercept and covariate-related coefficients for each model. Once estimates for the marginal mean utility and cost for each treatment group and time ( ${\mu }_{\mathrm{ju}},{\mu }_{\mathrm{jc}}$ ) are retrieved from the model, posterior estimates at the levels of QALYs and total costs ( ${\mathrm{\Delta }}_e,{\mathrm{\Delta }}_{\mathrm{tc}}$ ) can be derived based on ${\mu }_{\mathrm{ju}}$ and ${\mu }_{\mathrm{jc}}$ and equation 1. Figure 2 shows a graphical representation of the modeling framework for the utility and cost distributions at time $j$ in terms of combined modules denoted by the blue and red boxes, respectively.

Figure 2

Joint distribution, expressed in terms of a conditional distribution for utilities and costs at time $j$ given their values at time $j-1$ , respectively, indicated with a blue and red boxes. The parameters indexing the corresponding distributions or “modules” are indicated with different Greek letters, while the black, magenta, and green arrows show the dependence relationships between the parameters within and between the models and between times, respectively. Note that $i$ and $j$ denote the individual and time index, while the treatment index $t$ is omitted to ease notation.

Although this strategy can be more challenging to implement than that in the “Modelling of Total Costs and QALYs” section, it provides a more efficient way to handle partially observed cases. Indeed, by modeling longitudinal variables, missingness uncertainty at any time can be directly quantified through a probabilistic approach without relying on ad hoc imputation methods. However, if participants provide partial information at the level of the questionnaires (i.e., item nonresponse), analysts are still required to handle partially observed responses prior to model fitting. As a simple solution, they could either discard all cases with partially observed HRUs or impute as zero all partially observed HRUs at each time ( ${\overline{\mathrm{HRU}}}_{\mathrm{ij}}^k=0$ ) and use these to generate costs at each time $c_{\mathrm{ij}}$ .

Modeling of HRU Services and Utilities at Each Time

At the most disaggregated level, the data consist in the series of longitudinal variables ${\mathrm{hrqol}}_{\mathrm{ij}}^l,{\mathrm{hru}}_{\mathrm{ij}}^k$ , corresponding to the HRU and HRQoL individual responses at each time provided to the $k$ -th and $l$ -th item of the questionnaires, for $k=1,\dots ,K$ and $l=1,\dots ,L$ . Given that HRQoL questionnaires are mostly affected by unit nonresponse, no substantial gain can be obtained by focusing the modeling at the level of the items compared with the utilities, ³¹ whereas HRU questionnaires are often characterized by item nonresponse. As a result, we can specify the model at the level of $u_{\mathrm{ij}}$ and ${\mathrm{hru}}_{\mathrm{ij}}^k$ to make full use of the available information in the study. For example, the joint outcome distribution $p(u_{\mathrm{ij}},{\mathrm{hru}}_{\mathrm{ij}}^k)$ at time $j>0$ can be specified as

$$p(u_{\mathrm{ij}},{\mathrm{hru}}_{\mathrm{ij}}^k)=p(u_{\mathrm{ij}}\mid u_{i,j-1},{\mathrm{hru}}_{\mathrm{ij}}^k)p({\mathrm{hru}}_{\mathrm{ij}}^k\mid u_{\mathrm{ij}-1},{\mathrm{hru}}_{i,j-1}^k),$$

(6)

where $p({\mathrm{hru}}_{\mathrm{ij}}^k\mid u_{\mathrm{ij}-1},{\mathrm{hru}}_{\mathrm{ij}-1}^k)$ is the conditional distribution of the $k$ -th service at $j>0$ given utilities and the past HRU value at $j-1$ . Similarly to the modeling of $c_{\mathrm{ij}}$ , the distributions of ${\mathrm{hru}}_{\mathrm{ij}}^k$ are indexed by the set of time and service-specific location $({\varphi }_{\mathrm{ijhru}}^k)$ and nuisance $({\psi }_{\mathrm{jhru}}^k)$ parameters and should be chosen to handle the specific features of the data. While the generalized linear model structure for longitudinal outcomes in the “Modeling of Costs and Utilities at Each Time” section can also be applied to HRU data, model parameters should be given appropriate interpretations. For example, some link function $g_{\mathrm{hru}}(·)$ can be used to incorporate HRU-specific covariates ( $X_{\mathrm{ihru}}$ ) at the level of the conditional usage rate for service $k$ at time $j$ :

$$g_{\mathrm{hru}}\left({\varphi }_{\mathrm{ijhru}}^k\right)={\beta }_{j0}^k+{\beta }_{j1}^k{X}_{\mathrm{ihru}},$$

(7)

where the sets of regression parameters $\beta$ include the time- and service-specific intercept and covariate-related coefficients for the model. An additional complication is that service-specific HRUs are often characterized by large proportions of zero values. To overcome this problem, the use of 2-part regressions or Hurdle models has been suggested in the literature, although it has mostly been applied to cost variables. ²⁰ These consist in mixture models defined by 2 components: the first is a mass distribution at the spike, while the second is a parametric model applied to the natural range of the relevant variable. Typically, a logistic regression is used to estimate the probability of incurring a “structural” value (e.g., $0$ ); this is then used to weight the mean of the “nonstructural” values estimated in the second component.

To specify a Hurdle model for HRUs, we first define an individual indicator variable $d_{\mathrm{ihru}}^k$ taking value $1$ if the $i-$ th individual is associated with a zero value on the $k$ -th service at all times ( ${\overline{hru}}_{i\forall j}^k=0$ ) and 0 otherwise ( ${\mathrm{hru}}_{i\forall j}^k>0$ ). This is then modeled as

$$\begin{array}{c}\multicolumn{1}{c}{d_{\mathrm{ihru}}^k:=\mathbb{I}({\overline{\mathrm{hru}}}_{i\forall j}=0)~\mathrm{Bernoulli}\left({\eta }_i^k\right)}\\ \multicolumn{1}{c}{\mathrm{logit}\left({\eta }_i^k\right)={\gamma }_0^k+{\gamma }_1^k{X}_{i0},}\end{array}$$

(8)

where ${\eta }_i^k$ is the conditional probability of a zero structural value in the $k$ -th service, which is estimated on the logit scale as a function of a set of regression parameters $\mathbf{\gamma }$ and covariates $X_{i0}$ . Estimates for the marginal probability of a structural zero ${\pi }^k$ can be obtained either by linear combination of the regression parameters or by simulation based on equation 8. Estimates for the marginal usage rate of nonzero HRU values ${\mu }_{j>0}^k$ can be obtained by fitting equation 7 to only nonzero HRUs. Next, we compute the overall population average HRU usage rate for each service and time ${\mu }_{\mathrm{jhru}}^k$ as the linear combination

$${\mu }_{\mathrm{jhru}}^k=(1-{\pi }^k){\mu }_{j>0}^k.$$

These quantities can then be combined with unit prices $p^k$ associated with each service $k$ to generate corresponding estimates of marginal mean costs ${\mu }_{\mathrm{jc}}^k$ and summed across all services and times to obtain estimates of the marginal mean total costs ${\mu }_{\mathrm{tc}}$ . Figure 3 shows a graphical representation of the modeling framework for the utility, non-zero and zero cost distributions at time j in terms of modules denoted by the blue, red and green boxes, respectively.

Figure 3

Joint distribution, expressed in terms of conditional distributions for utilities, nonzero and structural zero HRUs at time $j$ given outcomes at time $j-1$ , respectively, indicated with a blue, red, and green boxes. The parameters indexing the corresponding distributions or “modules” are indicated with different Greek letters, while black, magenta, and green solid arrows show the dependence relationships between the parameters within and between the models and times, respectively. The dashed arrows denote deterministic relationships. Note that $i$ and $j$ denote the individual and time index, while the treatment index $t$ is omitted to ease notation.

The proposed strategy is more challenging to implement compared with those given in the “Modeling of Total Costs and QALYs” and “Modeling of Costs and Utilities at Each Time” sections due to the presence of multiple HRU variables, each consisting of a 2-part mixture for the nonzero and structural zero component. However, defining the model at the most disaggregated level provides the most efficient way to handle partially observed HRUs. We note that the dependence structure in equation 4 and equation 6 is simply one possible choice, which was selected based on the comparison of the model fit to the PBS data and with respect to alternative specifications. In general, we recommend that analysts explore different specifications and select the one that provides the best compromise between flexible dependence structures between the modeled variables and a feasible implementation to the available data.

Results

All models were fitted using JAGS, ³² a software specifically designed for the analysis of Bayesian models using Markov chain Monte Carlo (MCMC) simulation, ³³ which can be interfaced with R through the package R2jags. ³⁴ Samples from the posterior distribution of the parameters of interest generated by JAGS and saved to the R workspace were then used to produce summary statistics and plots. We ran 2 chains with $20,000$ iterations per chain, using a burn in of $10,000$ , for a total sample of $20,000$ iterations for posterior inference. For each unknown quantity in the model, we assessed convergence and autocorrelation of the MCMC simulations using diagnostic measures such as the potential scale reduction factor and the effective sample size. ³⁵ For each modeling approach, the fit of different distributions for the HRU/cost variables was compared using standard measures such as the widely applicable information criterion, or WAIC, ³⁶ and posterior predictive checks. Based on these comparisons, normal distributions for QALYs and utility variables, gamma distributions (compared with normal and lognormal distributions) for all cost variables, and normal distributions (compared with Poisson and negative binomial) for all HRU variables, using a hurdle approach to handle structural zeros was selected as the best-fitting distributions. Alternative prior distributions were considered to check that any unintended information was not incorporated into the models through the priors.

We also compared the results from our analyses to those of a traditional modeling approach based on a frequentist statistical framework. Although it is difficult to identify a standard approach for handling missing values and modeling health economic data in routine analyses, we used information collected from previous reviews^5,6,8 to identify the methods used in this analysis. A detailed description of the methods used under the traditional modeling approach is provided in Appendix A alongside a summary of the statistical and health economic results derived from it. In Appendix B, we explain in detail how the hurdle model was implemented in JAGS, while in the online supplementary material, we provide the full JAGS code for the models.

Model Estimates

Table 3 compares the posterior results of models fitted to the PBS data at the level of S1) total costs and QALYs, S2) costs and utilities at each time, and S3) HRU categories and utilities at each time. Within the first 2 strategies, alternative approaches to handle partially observed data are implemented. For strategy 1, missing cases are either included in the analysis (ALL), imputed as zero only for HRUs (IMP-H), or imputed as zero for both HRUs and costs (IMP-HC). For strategy 2, missing cases are either included in the analysis (ALL) or imputed as zero only for HRUs (IMP-H). Finally, for strategy 3, all data with no ad hoc imputation are used (ALL). Table 4 shows the posterior estimates and $95\%$ credible intervals of the marginal mean total costs and QALYs by treatment group ( ${\mu }_{\mathrm{tc}},{\mu }_e$ ) obtained from the different modeling strategies. Estimates (and credible intervals) for mean costs, utilities, and HRU rates, which can be obtained only from models fitted under strategy 2 and/or 3, are reported in Appendix A.

Table 3

Posterior Means and $95\%$ Credible Intervals of ${\mu }_{\mathrm{tc}}$ and ${\mu }_e$ in the Control ( $t=1$ ) and Intervention ( $t=2$ ) Groups Obtained under Alternative Missingness Approaches When Fitting the Model: at the Level of Total Costs and QALYs (Strategy 1), at the Level of Costs and Utilities at Each Time (Strategy 2), and at the Level of HRU and Utilities at Each Time (Strategy 3)

Approach	${\mathbf{\mu }}_{\mathit{tc}\mathbf{(}\mathit{t}\mathbf{=}1\mathbf{)}}$		${\mathbf{\mu }}_{\mathit{tc}\mathbf{(}\mathit{t}\mathbf{=}2\mathbf{)}}$		${\mathbf{\mu }}_{\mathit{e}\mathbf{(}\mathit{t}\mathbf{=}1\mathbf{)}}$		${\mathbf{\mu }}_{\mathit{e}\mathbf{(}\mathit{t}\mathbf{=}2\mathbf{)}}$
	Mean	95% CI	Mean	95% CI	Mean	95% CI	Mean	95% CI
Strategy 1: Total cost and QALY
ALL	2,543	(2,157; 2,938)	2,754	(2,249; 3,310)	0.487	(0.449; 0.575)	0.609	(0.569; 0.649)
IMP-H	2,888	(2,411; 3,393)	2,379	(1,939; 2,872)	0.488	(0.45; 0.526)	0.61	(0.571; 0.651)
IMP-HC	2,395	(1,897; 2,899)	2,237	(1,175; 2,749)	0.486	(0.449; 0.523)	0.61	(0.572; 0.649)
Strategy 2: Cost and utility at each time
ALL	2,607	(2,253; 2,971)	2,701	(2,278; 3,145)	0.494	(0.463; 0.527)	0.6	(0.566; 0.635)
IMP-H	2,453	(2,087; 2,843)	2,273	(1,874; 2,664)	0.494	(0.462; 0.526)	0.6	(0.565; 0.633)
Strategy 3: HRU category and utility at each time
ALL	2687	(2,173; 3,194)	2,587	(1,995; 3,206)	0.513	(0.475; 0.55)	0.599	(0.565; 0.634)

HRU, health care resource use; QALY, quality-adjusted life-year.

In general, estimates of mean QALYs in both treatment groups are always very similar. This is expected in that, within each strategy, zero imputation was considered only for HRU/cost variables, thus not affecting utilities and/or QALYs. In addition, since only unit nonresponse affects HRQoL data, no substantial information gain is obtained by modeling utilities compared with QALYs, with average estimates remaining almost unchanged across all strategies.

A comparison of the estimates obtained under strategy 1 suggests that an analysis fitted to all cases without zero imputation (ALL) is associated with mean total cost estimates for the intervention that are higher compared with those of the control group ( $\approx 211$ £), whereas after zero imputing the data (IMP-H or IMP-HC), the situation is reversed. A similar trend is also observed when comparing the estimates obtained under strategy 2, with the intervention being on average more expensive than the control under ALL ( $\approx 94$ £) but less expensive under IMP-H ( $\approx -180$ £). Finally, mean estimates obtained under ALL from strategy 3 suggest that the intervention is cheaper than the control ( $\approx -100$ £) with comparable credible interval widths to those from strategy 2. Figure 4 shows the mean differences in total costs between treatments derived from each model, distinguished according to the type of strategy (S1 = blue, S2 = green, S3 = red) and the missingness approach used (ALL, IMP-H, IMP-HC).

Figure 4

Estimated mean total cost difference between treatment groups in the Positive Behaviour Support study based on different modeling strategies (S1, S2, S3) and alternative approaches to handling missing health care resource use and/or cost data.

Mean estimates derived from models fitted at the level of total costs (S1, blue color) show a considerable degree of variability and slightly positive (ALL), slightly negative (IMP-HC), or quite negative (IMP-H) values. Similar conclusions can be drawn based on the results from models fitted at the level of the costs at each time (S2, green color). However, results derived from direct modeling of HRU data (S3, red color) show average estimates in favor of the intervention (slightly cheaper than the control) and located between the results under ALL and the other zero-imputation approaches (for S1 and S2). Although no formal bias assessment is possible, due to the empirical nature of the analyses, we note that for aggregated modeling strategies (S1 and S2), estimates based on some form of zero imputation (IMP-H, IMP-HC, IMP-H) are systematically lower compared with those based on models in which no ad hoc imputation was done (ALL). This suggests that reliance on the zero imputation of HRU/costs, which distorts the original data and likely leads to an underestimation of missingness uncertainty, can substantially drive the model estimates.

Economic Evaluation

We end by assessing the cost-effectiveness of the intervention, comparing the results under each modeling strategy. We specifically rely on the examination of the cost-effectiveness plane (CEP) ³⁷ and the cost-effectiveness acceptability curve (CEAC) ³⁸ to summarize the health economic analysis.

The CEP (Figure 5a) is a graphical representation of the joint distribution of the population average effectiveness and costs increments between the arms. In the graph, we show the results only under ALL (light blue for S1, light green for S2, and light red for S3) to ease presentation. The slope of the straight line crossing the plane is the willingness-to-pay threshold (often indicated as $k$ ), with points lying below the line defining the sustainability area, where the treatment is considered more cost-effective than the control is. In the graph, we also show the incremental cost-effectiveness ratio (ICER) computed under each approach and denoted it with a darker colored dot. For all strategies, almost all samples fall in the North-East and South-East quadrants quite close to each other, although the clouds of dots are progressively shifted downward moving from S1 to S3, with the latter also displaying a slightly negative ICER. This suggests that, even though results obtained under the different modeling strategies do not largely differ, slightly more favorable conclusions for the intervention group are observed under S3.

Figure 5

(a) Cost-effectiveness planes (CEPs) and (b) cost-effectiveness acceptability curves (CEACs) associated with alternative modeling strategies and missingness approaches. In the CEPs, incremental cost-effectiveness ratios based on the results from ALL under the 3 modeling strategies are indicated with corresponding darker colored dots (blue for S1, green for S2, and red for S3), while the straight line passing through the plot (evaluated at $k=ϵ25,000$ ) denotes the acceptance threshold value. In the CEACs, in addition to the results from ALL under each modeling strategy (solid lines), the probability values for the alternative missingness approaches (IMP-H and IMP-HC) are denoted with different types of blue (within S1) and green (within S2) dashed lines.

The CEAC (Figure 5b) is obtained as the proportion of dots lying in the sustainability area upon varying the acceptance threshold $k$ . For this analysis, based on current practice, we considered a range of $k$ up to £10,000 per QALY gained. The CEAC displays the probability of cost-effectiveness, thus providing a simple summary of the uncertainty associated with the “optimal” decision making suggested by the ICER. Interestingly, results obtained under ALL (solid lines) for all strategies are associated with estimates that are systematically lower compared with zero-imputation strategies under S1 and S2, that is, either IMP-H (dashed lines) or IMP-HC (dotted line), for all $k$ values. Finally, we note that results from different strategies without zero imputation (solid lines) show quite similar cost-effectiveness conclusions for $k>\pounds 3,000$ . This suggests how, across the scenarios explored and dataset used, differences in cost-effectiveness conclusions are mainly driven by the different approaches used to handle missing cases rather than the modeling strategies.

Discussion

Trial-based economic evaluations are typically conducted on quantities that are derived from disaggregated data, such as self-reported EQ-5D and HRU questionnaires, which are almost inevitably affected by missingness. For HRU outcomes, the lack of a “gold standard” method of data collection ³⁹ and the typical occurrence of item nonresponse patterns make the task of handling missingness particularly challenging. Analysts routinely rely on some ad hoc methods to handle missing values that often involve specific assumptions, for example, zero-imputed corresponding to no use of service. Although this might be reasonable in some cases, zero imputations is often performed to ensure a higher completeness rate for more aggregated outcomes. This, however, can be a dangerous practice in that it fails to fully recognize the impact of missingness uncertainty and may even distort estimates and mislead cost-effectiveness assessments. Although this problem may be less relevant in some cases, for example, when disaggregated data are either fully observed or fully missing at each moment of collection (e.g., EQ-5D questionnaires), it becomes crucial when data are affected by item-level missingness, which is typically the case for HRU questionnaires.

In this article, we have presented a general modeling framework to handle item-level HRU missingness without requiring any ad hoc imputations. The framework takes advantage of the Bayesian setting to handle different features of the data while also directly quantifying the impact of missingness uncertainty on cost-effectiveness results. Our approach represents an improvement with respect to the current practice and can be implemented in a relatively easy way using freely available software.

In the PBS study, analyses implemented without zero imputation resulted in estimates of cost-effectiveness that were considerably lower compared with approaches based on zero-imputed data, especially for low values of the acceptance threshold (between $30\%\mathrm{and}50\%$ lower for $k<\pounds 4,000$ ). These results suggest that model estimates can be highly affected by the approach used, with models fitted at more aggregated levels possibly requiring more restrictive assumptions about missing HRUs (e.g., assumed zeros) or the discard of partially observed cases, which may lead to underestimation of the impact of missingness uncertainty and a substantial loss of information.

According to these considerations and previous guidelines, ⁴⁰ we formulate the following recommendations for analysts. First, at the design stage, appropriate strategies should be used to minimize the number of missing values in the data collected during the study ⁴¹ such as reducing the length of the follow-up period and/or the number of items within the questionnaires. Second, once the data are collected, the practice of ad hoc imputing HRU data (e.g., zero) should be avoided, unless clearly motivated, since it can distort the data and lead to incorrect inferences. Third, at the analysis stage, the choice of the modeling approach should be informed based on the observed missingness patterns. When individuals are associated with fully missing responses in questionnaires and over time points, no advantage can be obtained by focusing on disaggregated data, and modeling at the level of QALYs and total costs is likely to be appropriate. When individuals are associated with fully missing responses in questionnaires but partially observed responses across time points, a model fitted at the level of longitudinal utilities and costs should be preferred. When individuals are associated with partially missing responses within questionnaires and across time points, only a model fitted at the level of questionnaire responses (e.g., HRU categories) allows for full use of all the available evidence collected. Note that, in practice, the implementation of models at more disaggregated levels becomes more challenging as the number of outcomes, time points, and missing values increases. Thus, analysts should consider the feasibility of the ideal approach in relation to the available data and, when necessary, implement strategies to facilitate its implementation. Examples include simplification of the dependence structure among the modeled variables to reduce the number of parameters to estimate, the use of more informative priors to handle sparse data, and aggregation of some types of disaggregated data (e.g., different types of health care categories) to reduce the number of modeled variables.

Results obtained from a traditional approach, based on multiple imputation by chained equations (MICE) and bootstrapping (shown in Appendix A), lead to estimates and cost-effectiveness conclusions that were in line with those from Bayesian models. However, we encountered some practical issues that limited the implementation of the traditional approach. First, all attempts to apply the method at the level of HRU categories led to substantial convergence problems in the MICE algorithm. Second, on a standard computer with a $16$ GB RAM and $4$ cores, bootstrap methods substantially increased the computational time needed to run the analysis to about $4$  h, compared with a maximum of $1$ h required to fit any of the Bayesian models. We therefore conclude that, although traditional approaches can be used to perform the analysis, in practice their implementation is likely to become more challenging and computationally demanding compared with a Bayesian approach, especially when the complexity of the analysis model is increased to account for the features of the data (e.g., nonnormality, correlation) and to quantify the impact of different sources of uncertainty on the results (e.g., missing data, decision making).

Our results are obtained with specific reference to the case study considered. However, the PBS trial is very much representative of the “typical” dataset used in health economic evaluation alongside studies. Thus, it is highly likely that the same features (and potentially the same contradictions in the results, upon varying the complexity of the modeling assumptions) apply to many real cases. We note that the decision-theoretic framework granted by the Bayesian approach provides a natural setting to quantify the impact of the uncertainty on the results while also allowing an intuitive probabilistic interpretation of standard cost-effectiveness outputs (e.g., CEACs). For example, if results are not robust to a set of departures from the benchmark missingness assumption, further analyses can be sought by means of more advanced methods to explicitly investigate the variability in the unobserved values based on external information (e.g., selection models or pattern mixture models).^42,43

There are 2 main limitations of the proposed framework. First, the choice of treating partially observed data for each variable under a common missingness assumption may be practically convenient but not realistic. It is plausible that the presence of different missingness patterns suggests the existence of different missingness processes, for example, people who drop out of the study are likely associated with different reasons with respect to those who missed a few visits, and such processes should be specified separately. Second, the model may become computationally challenging when the number of variables (i.e., time points or HRU categories) increases. Alternative approaches could be used to overcome these limitations. For example, the sparsity of the data could be handled by either aggregating some variables based on some plausible justifications (e.g., aggregating different types of GP visits as a single HRU category and computing a weighted price) or by using shared priors to use the information from the observed data across all or some of the variables to facilitate the identification of some parameters.

In conclusion, in this work we have presented a flexible Bayesian analytic framework that can 1) jointly model HRUs and effectiveness, 2) account for the features of the data, and 3) make full use of the available evidence to quantify the impact of missingness uncertainty without the need to rely on ad hoc imputations prior to model fitting. Unless clearly justified based on some external information, analyses that rely on these ad hoc imputations prior to model fitting will produce results that are likely driven by the specific approach used (and therefore implicit assumptions made) and should therefore be avoided.

Supplemental Material

sj-pdf-1-mdm-10.1177_0272989X251376026 – Supplemental material for A Bayesian Modeling Framework for Health Care Resource Use and Costs in Trial-Based Economic Evaluations

Supplemental material, sj-pdf-1-mdm-10.1177_0272989X251376026 for A Bayesian Modeling Framework for Health Care Resource Use and Costs in Trial-Based Economic Evaluations by Andrea Gabrio in Medical Decision Making