# Methods for estimating complier-average causal effects for cost-
effectiveness analysis

K. DiazOrdaz, A. J. Franchini, R. Grieve
Department of Medical Statistics,
London School of Hygiene and Tropical Medicine,
<EMAIL_ADDRESS>
###### Abstract

In Randomised Controlled Trials (RCT) with treatment non-compliance,
instrumental variable approaches are used to estimate complier average causal
effects. We extend these approaches to cost-effectiveness analyses, where
methods need to recognise the correlation between cost and health outcomes.

We propose a Bayesian full likelihood (BFL) approach, which jointly models the
effects of random assignment on treatment received and the outcomes, and a
three-stage least squares (3sls) method, which acknowledges the correlation
between the endpoints, and the endogeneity of the treatment received.

This investigation is motivated by the REFLUX study, which exemplifies the
setting where compliance differs between the RCT and routine practice. A
simulation is used to compare the methods performance.

We find that failure to model the correlation between the outcomes and
treatment received correctly can result in poor CI coverage and biased
estimates. By contrast, BFL and 3sls methods provide unbiased estimates with
good coverage. Keywords: non-compliance, instrumental variables, bivariate
outcomes, cost-effectiveness

##  1 Introduction

Non-compliance is a common problem in Randomised Controlled Trials (RCTs), as
some participants depart from their randomised treatment, by for example
switching from the experimental to the control regimen. An unbiased estimate
of the effectiveness of treatment assignment can be obtained by reporting the
intention-to-treat (ITT) estimand. In the presence of non-compliance, a
complimentary estimand of interest is the causal effect of treatment received.
Instrumental variable (IV) methods can be used to obtain the complier average
causal effect (CACE), as long as random assignment meets the IV criteria for
identification (Angrist et al., 1996). An established approach to IV
estimation is two-stage least squares (2sls), which provides consistent
estimates of the CACE when the outcome measure is continuous, and non-
compliance is binary (Baiocchi et al., 2014).

Cost-effectiveness analyses (CEA) are an important source of evidence for
informing clinical decision-making and health policy. CEA commonly report an
ITT estimand, i.e. the relative cost-effectiveness of the intention to receive
the intervention (NICE, 2013). However, policy-makers may require additional
estimands, such as the relative cost-effectiveness for compliers. For example,
CEAs of new therapies for end-stage cancer, are required to estimate the cost-
effectiveness of treatment receipt, recognising that patients may switch from
their randomised allocation following disease progression. Alternative
estimates such as the CACE may also be useful when levels of compliance in the
RCT differ to those in the target population, or where intervention receipt,
rather than the intention to receive the intervention, is the principal cost
driver. Methods for obtaining the CACE for univariate survival outcomes have
been exemplified before (Latimer et al., 2014), but approaches for obtaining
estimates that adequately adjust for non-adherence in CEA more generally, have
received little attention. This has been recently identified as a key area
where methodological development is needed (Hughes et al., 2016).

The context of trial-based CEA highlights an important complexity that arises
with multivariate outcomes more widely, in that, to provide accurate measures
of the uncertainty surrounding a composite measure of interest, for example
the incremental net monetary benefit (INB), it is necessary to recognise the
correlation between the endpoints, in this case, cost and health outcomes
(Willan et al., 2003; Willan, 2006). Indeed, when faced with non-compliance,
and the requirement to estimate a causal effect of treatment on cost-
effectiveness endpoints, some CEA resort to per protocol (PP) analyses
(Brilleman et al., 2015), which exclude participants who deviate from
treatment. As non-compliance is likely to be associated with prognostic
variables, only some of which are observed, PP analyses are liable to provide
biased estimates of the causal effect of the treatment received.

This paper develops novel methods for estimating CACE in CEA that use data
from RCTs with non-compliance. First, we propose using the three stage least
squares (3sls) method (Zellner and Theil, 1962), which allows the estimation
of a system of simultaneous equations with endogenous regressors. Next, we
consider a bivariate version of the ‘unadjusted Bayesian’ models previously
proposed for Mendelian randomisation (Burgess and Thompson, 2012), which
simultaneously estimate the expected treatment received as a function of
random allocation, and the mean outcomes as a linear function of the expected
treatment received. Finally, we develop a Bayesian full likelihood approach
(BFL), whereby the outcome variables and the treatment received are jointly
modelled as dependent on the random assignment. This is an extension to the
multivariate case of what is known in the econometrics literature as the IV
unrestricted reduced form (Kleibergen and Zivot, 2003).

The aim of this paper is to present and compare these alternative approaches.
The problem is illustrated in Section 2 with the REFLUX study, a multicentre
RCT and CEA that contrasts laparoscopic surgery with medical management for
patients with Gastro-Oesophageal Reflux Disease (GORD). Section 3 introduces
the assumptions and methods for estimating CACEs. Section 4 presents a
simulation study used to assess the performance of the alternative approaches,
which are then applied to the case study in Section 5. We conclude with a
Discussion (Section 6), where we consider the findings from this study in the
context of related research.

##  2 Motivating example: Cost-effectiveness analysis of the REFLUX study

The REFLUX study was a UK multicentre RCT with a parallel design, in which
patients with moderately severe GORD, were randomly assigned to medical
management or laparoscopic surgery (Grant et al., 2008, 2013).

The RCT randomised 357 participants (178 surgical, 179 medical) from 21 UK
centres. An observational preference based study was conducted alongside it,
which involved 453 preference participants (261 surgical, 192 medical).

For the cost-effectiveness analysis within the trial, individual resource use
(costs in $\pounds$ sterling) and health-related quality of life (HRQoL),
measured using EQ5D (3 levels), were recorded annually for up to 5 years. The
HRQoL data were used to adjust life years and present quality-adjusted life
years (QALYs) over the follow-up period (Grant et al., 2013).111There was no
administrative censoring. As is typical, the costs were right-skewed. Table 1
reports the main characteristics of the data set.

The original CEA estimated the linear additive treatment effect on mean costs
and health outcomes (QALYs). The primary analysis used a system of seemingly
unrelated regression equations (SURs) (Zellner, 1962; Willan et al., 2004),
adjusting for baseline HRQoL EQ5D summary score (denoted by
$\mbox{EQ5D}_{0}$). The SURs can be written for cost $Y_{1i}$ and QALYs
$Y_{2i}$, as follows

$\begin{array}[]{c}Y_{1i}=\beta_{0,1}+\beta_{1,1}\mbox{treat}_{i}+\beta_{1,2}\mbox{$\mbox{EQ5D}_{0}$}_{i}+\epsilon_{1i}\\\
Y_{2i}=\beta_{0,2}+\beta_{1,2}\mbox{treat}_{i}+\beta_{2,2}\mbox{$\mbox{EQ5D}_{0}$}_{i}+\epsilon_{2i}\end{array}$
(1)

where $\beta_{1,1}$ and $\beta_{1,2}$ represent the incremental costs and
QALYs respectively. The error terms are required to satisfy
$E[\epsilon_{1i}]=E[\epsilon_{2i}]=0,\
E[\epsilon_{ki}\epsilon_{k^{\prime}i}]=\sigma_{kk^{\prime}},\
E[\epsilon_{ki}\epsilon_{k^{\prime}j}]=0,\ \mbox{\ for\ }k,\
k^{\prime}\in\\{1,2\\},\mbox{\ and for\ }i\neq j.$ Rather than assuming that
the errors are drawn from a bivariate normal distribution, estimation is
usually done by the feasible generalized least squares (FGLS) method222If we
are prepared to assume the errors are bivariate normal, estimation can proceed
by maximum likelihood.. This is a two-step method where, in the first step, we
run ordinary least squares estimation for equation (1). In the second step,
residuals from the first step are used as estimates of the elements
$\sigma_{kk^{\prime}}$ of the covariance matrix, and this estimated covariance
structure is then used to re-estimate the coefficients in equation (1)
(Zellner, 1962; Zellner and Huang, 1962).

In addition to reporting incremental costs and QALYs, CEA often report the
incremental cost-effectiveness ratio (ICER), which is defined as the ratio of
the incremental costs per incremental QALY, and the incremental net benefit
(INB), defined as $\mbox{INB}(\lambda)=\lambda\beta_{1,2}-\beta_{1,1}$, where
$\lambda$ represents the decision-makers’ willingness to pay for a one unit
gain in health outcome. Thus the new treatment is cost-effective if
$\mbox{INB}>0$. For a given $\lambda$, the standard error of INB can be
estimated from the estimated increments $\hat{\beta}_{1,1}$ and
$\hat{\beta}_{1,2}$, together with their standard errors and their correlation
following the usual rules for the variance of a linear combination of two
random variables. The willingness to pay $\lambda$ generally lies in a range,
so it is common to compute the estimated value of INB for various values of
$\lambda$. In REFLUX, the reported INB was calculated using $\lambda=\pounds
30000$, which is within the range of cost-effectiveness thresholds used by the
UK National Institute for Health and Care Excellence (NICE, 2013).

The original ITT analysis concluded that, compared to medical management, the
arm assigned to surgery had a large gain in average QALYs, at a small
additional cost and was relatively cost-effective with a positive mean INB,
albeit with 95% confidence intervals (CI) that included zero. However, these
ITT results cannot be interpreted as a causal effect of the treatment, since
within one year of randomisation, 47 of those randomised to surgery switched
and received medical management, while in the medical treatment arm, 10
received surgery. The reported reasons for not having the allocated surgery
were that in the opinion of the surgeon or the patient, the symptons were not
“sufficiently severe” or the patient was judged unfit for surgery (e.g.
overweight). The preference-based observational study conducted alongside the
RCT reported that in routine clinical practice, the corresponding proportion
who switched from an intention to have surgery and received medical
management, was relatively low (4%), with a further 2% switching from control
to intervention (Grant et al., 2013). Since the percentage of patients who
switched in the RCT was higher than in the target population and the costs of
the receipt of surgery are relatively large, there was interest in reporting a
causal estimate of the intervention. Thus, the original study also reported a
PP analysis on complete-cases, adjusted for baseline $\mbox{EQ5D}_{0}$, which
resulted in an ICER of $\pounds 7263$ per additional QALY (Grant et al.,
2013). This is not an unbiased estimate of the causal treatment effect, so in
Section 5, we re-analyse the REFLUX dataset to obtain a CACE of the cost-
effectiveness outcomes, recognising the joint distribution of costs and QALYs,
using the methods described in the next section.

##  3 Complier Average Causal effects with bivariate outcomes

We begin by defining more formally our estimands and assumptions. Let $Y_{1i}$
and $Y_{2i}$ be the continuous bivariate outcomes, and $Z_{i}$ and $D_{i}$ the
binary random treatment allocation and treatment received respectively,
corresponding to the $i$-th individual. The bivariate endpoints $Y_{1i}$ and
$Y_{2i}$ belong to the same individual $i$, and thus are correlated. We assume
that there is an unobserved confounder $U$, which is associated with the
treatment received and either or both of the outcomes. From now on, we will
assume that the (i) Stable Unit Treatment Value Assumption (SUTVA) holds: the
potential outcomes of the $i$-th individual are unrelated to the treatment
status of all other individuals (known as no interference), and that for those
who actually received treatment level $z$, their observed outcome is the
potential outcome corresponding to that level of treatment.

Under SUTVA, we can write the potential treatment received by the $i$-th
subject under the random assignment at level $z_{i}\in\\{0,1\\}$ as
$D_{i}\left(z_{i}\right)$. Similarly, $Y_{\ell i}\left(z_{i},d_{i}\right)$
with $\ell\in\\{1,2\\}$ denotes the corresponding potential outcome for
endpoint $\ell$, if the $i$-th subject were allocated to level $z_{i}$ of the
treatment and received level $d_{i}$. There are four potential outcomes. Since
each subject is randomised to one level of treatment, only one of the
potential outcomes per endpoint $\ell$, is observed, i.e. $Y_{\ell i}=Y_{\ell
i}(z_{i},D_{i}(z_{i}))=Y_{i}(z_{i})$.

The CACE for outcome $\ell$ can now be defined as

$\theta_{\ell}=E\left[\\{Y_{\ell i}(1)-Y_{\ell
i}(0)\\}\big{|}\\{D_{i}(1)-D_{i}(0)=1\\}\right].$ (2)

In addition to (i) SUTVA, the following assumptions are sufficient for
identification of the CACE, (Angrist et al., 1996): (ii) Ignorability of the
treatment assignment: $Z_{i}$ is independent of unmeasured confounders
(conditional on measured covariates) and the potential outcomes
$Z_{i}\perp\\!\\!\\!\perp U_{i},D_{i}(0),D_{i}(1),Y_{i}(0),Y_{i}(1).$
(iii) The random assignment predicts treatment received:
$Pr\\{D_{i}(1)=1\\}\neq Pr\\{D_{i}(0)=1\\}$.
(iv) Exclusion restriction: The effect of $Z$ on $Y_{\ell}$ must be via an
effect of $Z$ on $D$; $Z$ cannot affect $Y_{\ell}$ directly.
(v) Monotonicity: $D_{i}(1)\geq D_{i}(0)$.

The CACE can now be identified from equation (2) without any further
assumptions about the unobserved confounder; in fact, $U$ can be an effect
modifier of the relationship of $D$ and $Y$ (Didilez et al., 2010).

In the REFLUX study, the assumptions concerning the random assignment, (ii)
and (iii), are justified by design. The exclusion restriction assumption seems
plausible for the costs, since the costs of surgery are only incurred if the
patient actually has the procedure. We argue that it is also plausible that it
holds for QALYs, as the participants did not seem to have a preference for
either treatment, thus making the psychological effects of knowing to which
treatment one has been allocated minimal. The monotonicity assumption rules
out the presence of defiers. It seems fair to assume that there are no
participants who would refuse the REFLUX surgery when randomised to it, but
who would receive surgery when randomised to receive medical management.
Equation (2) implicitly assumes that receiving the intervention has the same
average effect in the linear scale, regardless of the level of $Z$ and $U$.
This average is however across different ‘versions’ of the intervention, as
the trial protocol did not prescribe a single surgical procedure, but allowed
for the surgeon to choose their preferred laparoscopy method, as would be the
case in routine clinical practice.

Since random allocation, $Z$, satisfies assumptions (ii)-(iv), we say it is an
instrument (or instrumental variable) for $D$. For binary instrument, the
simplest method of estimation of equation (2) in the IV framework is the Wald
estimator (Angrist et al., 1996):

$\hat{\theta}_{\ell,IV}=\frac{E(Y|Z=1)-E(Y|Z=0)}{E(D|Z=1)-E(D|Z=0)}$

Typically, estimation of these conditional expectations proceeds via an
approach known as two-stage least squares (2sls). The first stage fits a
linear regression to treatment received on treatment assigned. Then, in a
second stage, a regression model for the outcome on the predicted treatment
received is fitted:

$\displaystyle D_{i}$ $\displaystyle=$
$\displaystyle\alpha_{0}+\alpha_{1}Z_{i}+\omega_{1i}$ $\displaystyle Y_{\ell
i}$ $\displaystyle=$
$\displaystyle\beta_{0}+\beta_{IV}\hat{D_{i}}+\omega_{2i}$ (3)

where $\hat{\beta}_{IV}$ is an estimator for $\theta_{\ell}$. Covariates can
be used, by including them in both stages of the model. To obtain the correct
standard errors for the 2sls estimator, it is necessary to take into account
the uncertainty about the first stage estimates. The asymptotic standard error
for the 2sls CACE is given in Imbens and Angrist (1994), and implemented in
commonly used software packages.

OLS estimation produces first-stage residuals $\omega_{1i}$ that are
uncorrelated with the instrument, and this is sufficient to guarantee that the
2sls estimator is consistent for the CACE (Angrist et al., 2008). Therefore,
we restrict our attention here to models where the first-stage equation is
linear, even though the treatment received is binary. 333Non-linear versions
of the 2sls exist. See for example Clarke and Windmeijer (2012) for an
excellent review of methods for binary outcomes.

A key issue for settings such as CEA where there is interest in estimating the
CACE for bivariate outcomes, is that 2sls as implemented in most software
packages can only be readily applied to univariate outcomes. Ignoring the
correlation between the two endpoints is a concern for obtaining standard
errors of composite measures of the outcomes, e.g. INB, as this requires
accurate estimates of the covariance between the outcomes of interest (e.g.
costs and QALYs).

A simple way to address this problem would be to apply 2sls directly to the
composite measure, i.e. a net-benefit two-stage regression approach (Hoch et
al., 2006). However, it is known that net benefit regression is very sensitive
to outliers, and to distributional assumptions (Willan et al., 2004), and has
been recently shown to perform poorly when these assumptions are thought to be
violated (Mantopoulos et al., 2016). Moreover, such net benefit regression is
restrictive, in that it does not allow separate covariate adjustment for each
of the component outcomes, (e.g. baseline HRQoL, for the QALYs as opposed to
the costs). In addition, this simple approach would not be valid for
estimating the ICER, which is a non-linear function of the incremental costs
and QALYs. For these reasons, we do not consider this approach further.
Rather, we present here three flexible strategies for estimating a CACE of the
QALYs and the costs, jointly. The first approach combines SURs (equation 1)
and 2sls (equation 3) to obtain CACEs for both outcomes accounting for their
correlation. This simple approach is known in the econometrics literature as
three-stage least squares (3sls).

###  3.1 Three-stage least squares

Three-stage least squares (3sls) was developed for SUR systems with endogenous
regressors, i.e. any explanatory variables which are correlated with the error
term in equations (1) (Zellner and Theil, 1962). All the parameters appearing
in the system are estimated jointly, in three stages. The first two stages are
as for 2sls, but with the second stage applied to each of the outcomes.

$\displaystyle\mbox{(1st stage):}\qquad D_{i}$ $\displaystyle=$
$\displaystyle\alpha_{0}+\alpha_{1}Z_{i}+e_{0i}$ $\displaystyle\mbox{(2nd
stage):}\qquad Y_{1i}$ $\displaystyle=$
$\displaystyle\beta_{01}+\beta_{IV,1}\hat{D_{i}}+e_{1i}$ (4) $\displaystyle
Y_{2i}$ $\displaystyle=$
$\displaystyle\beta_{02}+\beta_{IV,2}\hat{D_{i}}+e_{2i}$ (5)

As with 2sls, the models can be extended to include baseline covariates. The
third stage is the same step used on a SUR with exogenous regressors (equation
(1)) for estimating the covariance matrix of the error terms from the two
equations (3.1) and (5). Thus, because we are assuming that $Z$ satisfies the
identification assumptions (i)-(v), Z is independent of the residuals at first
and second stage, i.e. $Z\perp\\!\\!\\!\perp e_{0i}$, $Z\perp\\!\\!\\!\perp
e_{1i}$, and $Z\perp\\!\\!\\!\perp e_{2i}$. Then, the 3sls procedure allows us
to obtain the covariance matrix between the residuals $e_{1i}$ and $e_{2i}$.
As with SURs, the 3sls approach does not require to make any distributional
assumptions, as estimation can be done by FGLS, and it is robust to
heteroscedasticity of the errors in the linear models for the outcomes
(Greene, 2002). We note that the 3sls estimator based on FGLS is consistent
only if the error terms in each equation of the system and the instrument are
independent, which is likely to hold here, as we are dealing with a randomised
instrument. In settings where this condition is not satisfied, other
estimation approaches such as generalised methods of moments (GMM) warrant
consideration (Schmidt, 1990). In the just-identified case, i.e. when there
are as many endogenous regressors as there are instruments, classical theory
about 3sls estimators shows that the GMM and the FGLS estimators coincide
(Greene, 2002). As the 3sls method uses an estimated variance-covariance
matrix, it is only asymptotically efficient (Greene, 2002).

###  3.2 Naive Bayesian estimators

Bayesian models have a natural appeal for cost-effectiveness analyses, as they
afford us the flexibility to estimate bivariate models on the expectations of
the two outcomes using different distributions, as proposed by (Nixon and
Thompson, 2005). These models are often specified by writing a marginal model
for one of the outcomes, e.g. the costs $Y_{1}$, and then, a model for
$Y_{2}$, conditional on $Y_{1}$.

For simplicity of exposition, we begin by assuming normality for both outcomes
and no adjustment for covariates. We have a marginal model for $Y_{1}$ and a
model for $Y_{2}$ conditional on $Y_{1}$ (Nixon and Thompson, 2005)

$\displaystyle Y_{1i}$ $\displaystyle\sim$ $\displaystyle
N(\mu_{1i},\sigma^{2}_{1})\qquad\qquad\quad\
\mu_{1i}=\beta_{0,1}+\beta_{1,1}\mbox{treat}_{i}$ (6) $\displaystyle
Y_{2i}\mid Y_{1i}$ $\displaystyle\sim$ $\displaystyle
N(\mu_{2i},\sigma^{2}_{2}(1-\rho^{2}))\qquad\mu_{2i}=\beta_{0,2}+\beta_{1,2}\mbox{treat}_{i}+\beta_{2,2}(y_{1i}-\mu_{1i}),$
(7)

where $\rho$ is the correlation between the outcomes. The linear relationship
between the two outcomes is represented by
$\beta_{2,2}=\rho\frac{\sigma_{2}}{\sigma_{1}}$.

Because of the non-compliance, in order to obtain a causal estimate of
treatment, we need to add a linear model for the treatment received, dependent
on randomisation $Z_{i}$, similar to the first equation of a 2sls. Formally,
this model (denoted uBN, for unadjusted Bayesian Normal) can be written with
three equations as follows:

$\begin{array}[]{l}D_{i}\sim N(\mu_{0i},\sigma^{2}_{0})\\\ Y_{1i}\sim
N(\mu_{1i},\sigma^{2}_{1})\\\ Y_{2i}\big{|}Y_{1i}\sim
N(\mu_{2i},\sigma^{2}_{2}(1-\rho^{2}))\\\ \end{array}\ \ \ \ \ \ \ \
\begin{array}[]{l}\mu_{0i}=\beta_{0,0}+\beta_{1,0}Z_{i}\\\
\mu_{1i}=\beta_{0,1}+\beta_{1,1}\mu_{0i}\\\
\mu_{2i}=\beta_{0,2}+\beta_{1,2}\mu_{0i}+\beta_{2,2}(y_{1i}-\mu_{1i})\end{array}$
(8)

This model is a bivariate version of the ‘unadjusted Bayesian’ method
previously proposed for Mendelian randomisation (Burgess and Thompson, 2012).
It is called unadjusted, because the variance structure of the outcomes is
assumed to be independent of the treatment received. The causal treatment
effect for outcome $Y_{\ell}$, with $\ell\in\\{1,2\\}$, is represented by
$\beta_{1,\ell}$ in equations (8). We use the Fisher’s z-transform of $\rho$,
i.e. $z={1\over 2}\log\left({1+\rho\over 1-\rho}\right)$, for which we assume
a vague normal prior, i.e. $z\sim N(0,10^{2})$. We also use vague multivariate
normal priors for the regression coefficient (with a precision of 0.01). For
standard deviations, we use $\sigma_{j}\sim\mbox{Unif}(0,10)$, for
$j\in\\{0,1,2\\}$. This is similar to the priors used in (Lancaster, 2004),
and are independent of the regression coefficient of treatment received on
treatment allocation $\beta_{1,0}$.

Cost data are notoriously right-skewed, and Gamma-distributions are often used
to model them. Thus, we can relax the normality assumption of equation (8),
and model $Y_{1}$ (i.e. cost) with a Gamma distribution, and treatment
received (binary), with a logistic regression. The health outcomes, $Y_{2}$,
are still modelled with a normal distribution, as is customary. Because we are
using a non-linear model for the treatment received, we use the predicted raw
residuals from this model as extra regressors in the outcome models, similar
to the 2-stage residual inclusion estimator (Terza et al., 2008). We model
$Y_{1}$ by its marginal distribution (Gamma) and $Y_{2}$ by a conditional
Normal distribution, given $Y_{1}$ (Nixon and Thompson, 2005). We call this
model ‘unadjusted Bayesian Gamma-Normal’ (uBGN) and write it as follows

$\displaystyle\begin{array}[]{l}D_{i}\sim\mbox{Bern}(\pi_{i})\\\
Y_{1i}\sim\mbox{Gamma}(\nu_{1i},\kappa_{1})\\\ Y_{2i}\mid Y_{1i}\sim
N(\mu_{2i},\sigma^{2}_{2}(1-\rho^{2}))\\\ \end{array}\ \ \ \ \ \ \ \
\begin{array}[]{l}\mbox{logit}(\pi_{i})=\alpha_{0}+\alpha_{1}Z_{i}\\\
r_{i}=D_{i}-\pi_{i}\\\
\mu_{1i}=\beta_{0,1}+\beta_{1,1}D_{i}+\beta_{1,r}r_{i}\\\
\mu_{2i}=\beta_{0,2}+\beta_{1,2}D_{i}+\beta_{2,r}r_{i}+\beta_{2,2}(y_{1i}-\mu_{1i})\end{array}$
(16)

where $\mu_{1}=\frac{\nu_{1}}{\kappa_{1}}$, is the mean of the Gamma
distributed costs, with shape $\nu_{1}$ and rate $\kappa_{1}$. Again, we
express $\beta_{2,2}=\rho\frac{\sigma_{2}}{\sigma_{1}}$, and assume a vague
Normal prior on the Fisher’s z-transform of $\rho$, $z\sim N(0,10^{2})$.The
prior distribution for $\nu_{1}$ is $\mbox{Gamma}(0.01,0.01)$. We also assume
a Gamma prior for the intercept term of the cost equation,
$\beta_{0,1}\sim\mbox{Gamma}(0.01,0.01)$. All the other regression parameters
have the same priors as those used in the uBN model.

The models introduced in this section, uBN and uBGN, are estimated in one
stage, allowing feedback between the regression equations and the propagation
of uncertainty. However, these ’unadjusted’ methods ignore the correlation
between the outcomes and the treatment received. This misspecification of the
covariance structure may result in biases in the causal effect, which are
likely to be more important at higher levels of non-compliance.

###  3.3 Bayesian simultaneous equations (BFL)

We now introduce an approach that models the covariance between treatment
received and outcomes appropriately, using a system of simultaneous equations.
This can be done via full or limited information maximum likelihood, or by
using MCMC to estimate the parameters in the model simultaneously allowing for
proper Bayesian feedback and propagation of uncertainty. Here, we propose a
Bayesian approach which is an extension of the methods presented in Burgess
and Thompson (2012); Kleibergen and Zivot (2003); Lancaster (2004).

This method treats the endogenous variable $D$ and the cost-effectiveness
outcomes as covariant and estimates the effect of treatment allocation, as
follows. Let $(D_{i},Y_{1i},Y_{2i})^{\top}$ be the transpose of vector of
outcomes, which now includes treatment received, as well as the bivariate
endpoints of interest. We treat all three variables as multivariate normally
distributed, so that

$\begin{array}[]{l}\left(\\!\\!\begin{array}[]{c}D_{i}\\\ Y_{1i}\\\
Y_{2i}\end{array}\\!\\!\right)\sim
N\left\\{\\!\\!\left(\begin{array}[]{c}\mu_{0i}\\\ \mu_{1i}\\\
\mu_{2i}\end{array}\right),\Sigma=\left(\\!\\!\\!\begin{array}[]{ccc}\sigma^{2}_{0}&s_{01}&s_{02}\\\
s_{01}&\sigma^{2}_{1}&s_{12}\\\
s_{02}&s_{12}&\sigma^{2}_{2}\end{array}\\!\\!\\!\right)\\!\\!\right\\};\end{array}\
\ \ \begin{array}[]{l}\mu_{0i}=\beta_{0,0}+\beta_{1,0}Z_{i}\\\
\mu_{1i}=\beta_{0,1}+\beta_{1,1}\beta_{1,0}Z_{i}\\\
\mu_{2i}=\beta_{0,2}+\beta_{1,2}\beta_{1,0}Z_{i}\end{array}$ (17)

where $s_{ij}=\mbox{cov}(Y_{i},Y_{j})$, and the causal treatment effect
estimates are $\beta_{1,1}$ and $\beta_{1,2}$ respectively. For the
implementation, we use vague normal priors for the regression coefficients,
i.e. $\beta_{m,j}\sim N(0,10^{2})$, for $j\in\\{0,1,2\\},m\in\\{0,1\\}$, and a
Wishart prior for the inverse of $\Sigma$ (Gelman and Hill, 2006).

##  4 Simulation study

We now use a factorial simulation study to assess the finite sample
performance of the alternative methods. The first factor is the proportion of
participants who do not comply with the experimental regime, when assigned to
it, expressed as a percentage of the total (one-sided non-compliance). Bias is
expected to increase with increasing levels of non-compliance. A systematic
review (Dodd et al., 2012) found that the percentage of non-compliance was
less than 30% in two-thirds of published RCTs, but greater than 50% in one-
tenth of studies. Here, two levels of non-compliance are chosen, 30% and 70%.
As costs are typically skewed, three different distributions (Normal, Gamma or
Inverse Gaussian – IG) are used to simulate cost data. As the 2sls approach
fails to accommodate the correlation between the endpoints, we examined the
impact of different levels of correlation on the methods’ performance; $\rho$
takes one of four values $\pm 0.4,\ \pm 0.8$. The final factor is the sample
size of the RCT, taking two settings $n=100,\text{and}\ 1000$. In total, there
are $2\times 3\times 4\times 2=48$ simulated scenarios.

To generate the data, we begin by simulating $U\sim N(0.50,0.25^{2})$,
independently from treatment allocation. $U$ represents a pre-randomisation
variable that is a common cause of both the outcomes and the probability of
non-compliance, i.e. it is a confounding variable, which we assume is
unobserved.

Now, let $S_{i}\sim\mbox{Bern}(\pi_{s})$ be the random variable denoting
whether the $i$th individual switches from allocated active treatment to
control. The probability $\pi_{s}$ of one-way non-compliance with allocated
treatment depends on $U$, in the following way,

$\pi_{s}=\begin{cases}p+0.1,\text{if\ }u>0.5,\\\
p-0.1,\text{otherwise}\end{cases}$ (18)

where $p$ denotes the corresponding average non-compliance percentage
expressed as a probability, i.e here $p\in\\{0.3,0.7\\}$. We now generate
$D_{i}$, the random variable of treatment received as

$D_{i}=\begin{cases}Z_{i},\text{if either\ }s_{i}=0\text{\ or\ }Z_{i}=0,\\\
1-Z_{i},\text{if\ }s_{i}=1\text{ and }Z_{i}=1\end{cases}$ (19)

where $Z_{i}$ denotes the random allocation for subject $i$.

Then, the means for both outcomes are assumed to depend linearly on treatment
received and the unobserved confounder $U$ as follows:

$\displaystyle\mu_{1}=E[Y_{1}]$ $\displaystyle=$ $\displaystyle
1.2+0.4D_{i}+0.16(u_{i}-0.5)$ (20) $\displaystyle\mu_{2}=E[Y_{2}]$
$\displaystyle=$ $\displaystyle 0.5+0.2D_{i}+0.04(u_{i}-0.5)$ (21)

Finally, the bivariate outcomes are generated using Gaussian copulas,
initially with normal marginals. In subsequent scenarios, we consider Gamma or
Inverse Gaussian marginals for $Y_{1}$ and normal for $Y_{2}$. The conditional
correlation between the outcomes, $\rho$, is set according to the
corresponding scenario.

For the scenarios where the endpoints are assumed to follow a bivariate normal
distribution, the variances of the outcomes are set to
$\sigma^{2}_{1}=0.2^{2},\sigma^{2}_{2}=0.1^{2}$ respectively, while for
scenarios with Gamma and IG distributed $Y_{1}$, the shape parameter is
$\eta=4$. For the Gamma case, this gives a variance for $Y_{1}$ equal to
$\sigma^{2}_{1}=0.36$ in control and $\sigma^{2}_{1}=0.64$ in the intervention
group. When $Y_{1}\sim\mbox{IG}(\mu_{1},\eta)$, the expected variance in the
control group is $\sigma^{2}_{1}=0.432$, and $\sigma^{2}_{1}=1.024$ in those
receiving the intervention.

The simulated endpoints represent cost-effectiveness variables that have been
rescaled for computational purposes, with costs divided by $1000$, and QALYs
by $0.1$, such that the true values are $\pounds 400$ (incremental costs), and
$0.02$ (incremental QALYs) and so with a threshold value of $\lambda=\pounds
30000$ per QALY, the true causal INB is $\pounds 200$.

For each simulated scenario, we obtained $M=2500$ sets. For the Bayesian
analyses, we use the median of the posterior distribution as the ‘estimate’ of
the parameter of interest, and the standard deviation of the posterior
distribution as the standard error. Equi-tailed 95% posterior credible
intervals are also obtained. We use the term confidence interval for the
Bayesian credible intervals henceforth, to have a unified terminology for both
Bayesian and frequentist intervals.

Once the corresponding causal estimate has been obtained in each of the 2500
replicated sets under each scenario in turn, we compute the median bias of the
estimates, coverage of 95% confidence intervals (CI), median CI width and root
mean square error (RMSE). We report median bias as opposed to mean bias,
because the BFL leads to a posterior distribution of the causal parameters
which is Cauchy-like (Kleibergen and Zivot, 2003). A method is ‘adequate’, if
it results in low levels of bias (median bias $\leq 5$%) with coverage rates
within 2.5% of the nominal value.

Implementation:
The 3sls was fitted using systemfit package in R using FGLS, and the Bayesian
methods were run using JAGS from R (r2jags). Two chains, each one with 5000
initial iterations and 1000 burn-in were used. The multiple chains allowed for
a check of convergence by the degree of their mixing and the initial
iterations enabled to estimate iteration autocorrelation. A variable number of
further 1000-iteration runs were performed until convergence was reached as
estimated by the absolute value of the Geweke statistics for the first 10% and
last 50% of iterations in a run being below 2.5. A final additional run of
5000 iterations was performed for each chain to achieve a total sample of
10000 iterations, and a MC error of about 1% of the parameter SE on which to
base the posterior estimates. For the uBGN, an offset of 0.01 is added to the
shape parameter $\nu_{1}$ for the Gamma distribution of the cost, to prevent
the sampled shape parameter to become too close to zero, which may result in
infinite densities. See the Appendix for the JAGS model code for BFL.

###  4.1 Simulation Study Results

Bias
Figure 1 shows the median bias corresponding to scenarios with 30% non-
compliance, by cost distributions (left to right) and levels of correlation
between the two outcomes, for sample sizes of $n=100$ (upper panel) and
$n=1000$ (lower panel). With the larger sample size, for all methods, bias is
negligible with normally distributed costs, and remain less than 5% when costs
are Gamma-distributed. However, when costs follow an Inverse Gaussian
distribution, and the absolute levels of correlation between the endpoints are
high ($\pm 0.8$), the uBGN approach results in biased estimates, around 10%
bias for the estimated incremental cost, and between 20 and 40% for the
estimated INB. With the small sample size and when costs follow a Gamma or
Inverse Gaussian distribution, both unadjusted Bayesian methods provide
estimates with moderate levels of bias. With 70% non-compliance (Figure A3 in
the Appendix), the unadjusted methods result in important biases which persist
even with large sample sizes, especially for scenarios with non-normal
outcomes. For small sample settings, uBN reports positive bias (10 to 20%) in
the estimation of incremental QALYs, and the resulting INB, irrespective of
the cost distribution. The uBGN method reports relatively unbiased estimates
of the QALYs, but large positive bias (up to 60%) in the estimation of costs,
and hence, there is substantial bias in the estimated INB (up to 200%).
Unadjusted Bayesian methods estimates are biased in the direction of the true
causal effect, because these methods ignore the positive correlation between
treatment received and each outcome introduced by the confounding variable. By
contrast, the BFL and the 3sls provide estimates with low levels of bias
across most settings.
CI coverage and width
Table 2 presents the results for CI coverage and width, for scenarios with a
sample size of $n=100$, absolute levels of correlation between the endpoints
of 0.4, and 30% non-compliance. All other results are shown in the Appendix.
The 2sls INB ignores the correlation between costs and QALYs, and thus,
depending on the direction of this correlation, 2sls reports CI coverage that
is above (positive correlation) or below (negative correlation) nominal
levels. This divergence from nominal levels increases with higher absolute
levels of correlation (see Appendix, Table A6).

The uBN approach results in over-coverage across many settings, with wide CIs.
For example, for both levels of non-compliance and either sample size, when
the costs are Normal, the CI coverage rates for both incremental costs and
QALYs exceed 0.98. The interpretation offered by Burgess and Thompson (2012)
is also relevant here: the uBN assumes that the treatment received and the
outcomes variance structures are uncorrelated, and so when the true
correlation is positive, the model overstates the variance and leads to wide
CIs. By contrast, the uBGN method results in low CI coverage rates for the
estimation of incremental costs, when costs follow an inverse Gaussian
distribution. This is because the model incorrectly assumes a Gamma
distribution, thereby underestimating the variance. The extent of the under-
coverage appears to increase with higher absolute values of the correlation
between the endpoints, with coverage as low as 0.68 (incremental costs) and
0.72 (INB) in scenarios where the absolute value of correlation between costs
and QALYs is 0.8. (see Appendix, Table A7).

The BFL approach reports estimates with CI coverage close to the nominal when
the sample size is large, but with excess coverage (greater than 0.975), and
relatively wide CI, when the sample size is $n=100$ (see Table 2 for 30%
noncompliance, and Table A4 in the Appendix for the result corresponding to
70% non-compliance). By contrast, the 3sls reports CI coverage within 2.5% of
nominal levels for each sample size, level of non-compliance, cost
distribution and level of correlation between costs and QALYs.
RMSE
Table 2 reports RMSE corresponding to 30% non-compliance, and $n=100$. The
least squares approaches result in lower RMSE than the other methods for the
summary statistic of interest, the INB. This pattern is repeated across other
settings, see the Appendix, Tables A10–A16.

##  5 Results for the motivating example

We now compare the methods in practice by applying them to the REFLUX dataset.
Only 48% of the individuals have completely observed cost-effectiveness
outcomes: there were 185 individuals with missing QALYs, 166 with missing
costs, and a further 13 with missing $\mbox{EQ5D}_{0}$ at baseline, with about
a third of those with missing outcomes having switched from their allocated
treatment. These missing data not only bring more uncertainty to our analysis,
but more importantly, unless the missing data are handled appropriately can
lead to biased causal estimates (Daniel et al., 2012). A complete case
analysis would be unbiased, albeit inefficient, if the missingness is
conditionally independent of the outcomes given the covariates in the model
(White et al., 2010), even when the covariates have missing data, as is the
case here.444This mechanism is a special case of missing not at random.
Alternatively, a more plausible assumption is to assume the missing data are
missing at random (MAR), i.e. the probability of missingness depends only on
the observed data, and use multiple imputation (MI) or a full Bayesian
analysis to obtain valid inferences.

Therefore, we perform MI prior to carrying out 2sls and 3sls analyses. We
begin by investigating all the possible associations between the covariates
available in the data set and the missingness, univariately for costs, QALYs
and baseline $\mbox{EQ5D}_{0}$. Covariates which are predictive of both, the
missing values and the probability of missing, are to be included in the
imputation model as auxiliary variables, as conditioning on more variables
helps make the MAR assumption more plausible. None of the fully observed
variables at our disposal (age, sex, and baseline BMI) satisfies these
criteria and therefore, we do not include any auxiliary variables in our
imputation models. Thus, we impute total cost, total QALYs and baseline
$\mbox{EQ5D}_{0}$, 50 times by chained equations, using predictive mean
matching (PMM), taking the 5 nearest neighbours as donors (White et al.,
2011), including treatment received in the imputation model and stratifying by
treatment allocation. We perform 2sls on costs and QALYs independently and
calculate (within MI) SE for the INB assuming independence between costs and
QALYs. For the 3sls approach, the model is fitted to both outcomes
simultaneously, and the post-estimation facilities are used to extract the
variance-covariance estimate, and compute the estimated INB and its
corresponding SE. We also use the CACE estimates of incremental cost and QALYs
to obtain the ICER. After applying each method to the 50 MI sets, we combine
the results using Rubin’s rules (Rubin, 1987).555Applying IV 2sls and 3sls
with multiply imputed datasets, and combining the results using Rubin’s rules
can be done automatically in Stata using mi estimate, cmdok: ivregress 2sls
and mi estimate, cmdok: reg3. In R, ivregress can be used with with.mids
command within mice, but systemfit cannot presently be combined with this
command so, Rubin’s rules have to be coded manually. Sample code is available
in the Appendix.

For the Bayesian approaches, the missing values become extra parameters to
model. Since baseline $\mbox{EQ5D}_{0}$ has missing observations, a model for
its distribution is added $\mbox{EQ5D}_{0}\sim N(\mu_{q0},\sigma_{q0}^{2})$,
with a vaguely informative prior for $\mu_{q0}\sim\mbox{Unif}(-0.5,1)$, and an
uninformative prior for $|\sigma_{q0}|\sim N(0,0.01)$. We add two extra lines
of code to the models to obtain posterior distributions for INB and ICERs. We
center the outcomes around the empirical mean (except for costs, when modelled
as Gamma) and re-scale the costs (dividing by 1000) to improve mixing and
convergence. We use two chains, initially running 15,000 iterations with 5,000
as burn-in. After checking visually for auto-correlation, an extra 10,000
iterations are needed to ensure that the density plots of the parameters
corresponding to the two chains are very similar, denoting convergence to the
stationary distribution. Enough iterations for each chain are kept to make the
total effective sample (after accounting for auto-correlation) equal to
10,000. 666Multivariate normal nodes cannot be partially observed in JAGS,
thus, we run BFL models on all available data within WinBUGs. An observation
with zero costs was set to missing when running the Bayesian Gamma model,
which requires strictly positive costs.

Table 1 shows the results for incremental costs, QALYs and INB for the
motivating example adjusted for baseline $\mbox{EQ5D}_{0}$. Bayesian posterior
distribution are summarised by their median value and 95% credible intervals.
The CACEs are similar across the methods, except for uBGN, where the
incremental QALYs CACE is nearly halved, resulting in a smaller INB with a CI
that includes 0. In line with the simulation results, this would suggest that,
where the uBGN is misspecified according to the assumed cost distribution, it
can provide a biased estimate of the incremental QALYs.

Comparing the CACEs to the ITT, we see that the incremental cost estimates
increases between an ITT and a CACE, as actual receipt of surgery carries with
it higher costs that the mere offering of surgery does not. Similarly, the
incremental QALYs are larger, meaning that amongst compliers, those receiving
surgery have a greater gain in quality of life, over the follow-up period. The
CACE for costs are relatively close to the per-protocol incremental costs
reported in the original study, $\pounds 2324\ (1780,2848)$. In contrast, the
incremental QALYs according to PP on complete-cases originally reported was
$0.3200\ (0.0837,0.5562)$, considerably smaller than our CACE estimates (Grant
et al., 2013). The ITT ICER obtained after MI was $\pounds 4135$, while using
causal incremental costs and QALYs, the corresponding estimates of the CACE
for ICER were $\pounds 4140$ (3sls), $\pounds 5189$ (uBN), $\pounds 5960$
(uBGN), and $\pounds 3948$ (BFL). The originally reported per-protocol ICER is
$\pounds 7932$ per extra QALY was obtained on complete cases only (Grant et
al., 2013).

These results may be sensitive to the modelling of the missing data. As a
sensitivity analysis to the MAR assumption, we present the complete case
analysis in Table A1 in the Appendix. The conclusions from complete-case
analysis are similar to those obtained under MAR.

We also explore the sensitivity to choices of priors, by re-running the BFL
analyses using different priors, first for the multivariate precision matrix,
keeping the priors for the coefficients normal, and then a second analysis,
with uniform priors for the regression coefficient, and an inverse Wishart
prior with 6 degrees of freedom and a identity scale matrix, for the
precision. The results are not materially changed (see Appendix, Table A2).

The results of the within-trial CEA suggest that amongst compliers,
laparoscopy is more cost-effective than medical management for patients
suffering from GORD. The results are robust to the choice of priors, and to
the assumptions about the missing data mechanism. The results for the uBGN
differ somewhat from the other models, and as our simulations show, the
concern is that such unadjusted Bayesian models are prone to bias from model
misspecification.

##  6 Discussion

This paper extends existing methods for CEA (Willan et al., 2003; Nixon and
Thompson, 2005), by providing IV approaches for obtaining causal cost-
effectiveness estimates for RCTs with non-compliance. The methods developed
here however are applicable to other settings with multivariate continuous
outcomes more generally, for example RCTs in education, with different
measures of attainment being combined into an overall score. To help
dissemination, we provide code in the Appendix.

We proposed exploiting existing 3sls methods and also considered IV Bayesian
models, which are extensions of previously proposed approaches for univariate
continuous outcomes. Burgess and Thompson (2012) found the BFL was median
unbiased and gave CI coverage close to nominal levels, albeit with wider CIs
than least-squares methods. Their ‘unadjusted Bayesian’ method, similar to our
uBN approach, assumes that the error term for the model of treatment received
on treatment allocated is uncorrelated with the error from the outcome models.
This results in bias in the direction of the true association, and affects the
CI coverage. Our simulation study shows that, in a setting with multivariate
outcomes, the bias can be substantial. A potential solution to this could be
using priors for the error terms that reflect this dependency explicitly. For
example, Rossi et al. (2012) propose using a prior for the errors that
explicitly depends on the coefficient $\beta_{1,0}$, the effect of treatment
allocation on treatment received, in equation (8). Kleibergen and Zivot (2003)
propose priors that also reflect this dependency explicitly, and replicate
better the properties of the 2sls. This is known as the “Bayesian two-stage
approach”.

The results of our simulations show that applying 2sls separately to the
univariate outcomes leads to inaccurate 95% CI around the INB, even with
moderate levels of correlation between costs and outcomes ($\pm 0.4$). Across
all the settings considered, the 3sls approach resulted in low levels of bias
for the INB and unlike 2sls, provided CI coverage close to nominal levels. BFL
performed well with large sample sizes, but produced standard deviations which
were too large when the sample size was small, as can be seen from the over-
coverage, with wide CIs.

The REFLUX study illustrated a common concern in CEA, in that the levels of
non-compliance in the RCT were different, in this case higher, to those in
routine practice. The CACEs presented provide the policy-maker with an
estimate of what the relative cost-effectiveness would be if all the RCT
participants had complied with their assigned treatment, which is
complementary to the ITT estimate. Since we judged the IV assumptions for
identification likely to hold in this case-study, we conclude that either 3sls
or BFL provide valid inferences for the CACE of INB. The re-analysis of the
REFLUX case study also provided the opportunity to investigate the sensitivity
to the choice of priors in practice. Here we found that our choice of weakly
informative priors, which were relatively flat in the region where the values
of the parameters were anticipated to be, together with samples of at least
size 100, had minimal influence on the posterior estimates. We repeated the
analysis using different vague priors for the parameters of interests and the
corresponding results were not materially changed.

We considered here relatively simple frequentist IV methods, namely 2sls and
3sls. One alternative approach to the estimation of CACEs for multivariate
responses, is to use linear structural equation modelling, estimated by
maximum-likelihood expectation-maximization (ML-EM) algorithm (Jo and Muthén,
2001). Further, we only considered those settings where a linear additive
treatment effect is of interest, and the assumptions for identification are
met. Where interest lies in systems of simultaneous non-linear equations with
endogeneous regressors, GMM or generalised structural equation models can be
used to estimate CACEs (Davidson and MacKinnon, 2004).

There are several options to study the sensitivity to departures from the
identification assumptions. For example, if the exclusion restriction does not
hold, a Bayesian parametric model can use priors on the non-zero direct effect
of randomisation on the outcome for identification (Conley et al., 2012;
Hirano et al., 2000). Since the models are only weakly identified, the results
would depend strongly on the parametric choices for the likelihood and the
prior distributions. In the frequentist IV framework, such modelling is also
possible, see Baiocchi et al. (2014) for an excellent tutorial on how to
conduct sensitivity analysis to violations of the ER and monotonicity
assumptions. Alternatively, violations of the ER can also be handled by using
baseline covariates to model the probability of compliance directly, within
structural equation modelling via ML-EM framework (Jo, 2002a, b).

Settings where the instrument is only weakly correlated with the endogenous
variable have not been considered here, because for binary non-compliance with
binary allocation, the percentage of one-way non-compliance would need to be
in excess of 85%, for the F-statistic of the randomization instrument to be
less than 10, the traditional cutoff beneath which an instrument is regarded
as ‘weak’. Such levels of non-compliance are not realistic in practice, with
the reported median non-compliance equal to 12% (Zhang et al., 2014).
Nevertheless, Bayesian IV methods have been shown to perform better than 2sls
methods when the instrument is weak (Burgess and Thompson, 2012).

Also, for simplicity, we restricted our analysis of the case study to MAR and
complete cases assumptions. Sensitivity to departures from these assumptions
is beyond the scope of this paper, but researchers should be aware of the need
to think carefully about the possible causes of missingness, and conduct
sensitivity analysis under MNAR, assuming plausible differences in the
distributions of the observed and the missing data. When addressing the
missing data through Bayesian methods, the posterior distribution can be
sensitive to the choice of prior distribution, especially with a large amount
of missing data (Hirano et al., 2000).

Future research directions could include exploiting the additional flexibility
of the Bayesian framework to incorporate informative priors, perhaps as part
of a comprehensive decision modelling approach. The methods developed here
could also be extended to handle time-varying non-compliance.

### Acknowledgements

We thank M. Sculpher, R. Faria, D. Epstein, C. Ramsey and the REFLUX study
team for access to the data.
K. DiazOrdaz was supported by UK Medical Research Council Career development
award in Biostatistics MR/L011964/1. This report is independent research
supported by the National Institute for Health Research (Senior Research
Fellowship, R. Grieve, SRF-2013-06-016). The views expressed in this
publication are those of the author(s) and not necessarily those of the NHS,
the National Institute for Health Research or the Department of Health.

Figure 1: Median Bias for scenarios with 30% non-compliance and sample sizes
of (a) $n=100$ (top) and (b) $n=1000$. Results are stratified by cost
distribution, and correlation between cost and QALYs. The dotted line
represents zero bias. Results for 2sls (not plotted) are identical to those
for 3sls; uBGN was not applied to Normal cost data.

(a) n=100 (b) n=1000

Table 1: The REFLUX study: descriptive statistics and cost-effectiveness according to ITT and alternative methods for estimating the CACE. Follow-up period is five years, and treatment switches are defined within the first year post randomisation. Costs and INB numbers rounded to the nearest integer. | Medical management | Laparoscopic surgery
---|---|---
N Assigned | $179$ | $178$
N (%) Switched |  $10$ $(8.3)$ |  $67$ $(28.4)$
N (%) missing costs | 83 (46) | 83 (47)
Mean (SD) observed cost in £ | $1258\ (1687)$ | $2971\ (1828)$
N (%) missing QALYs | 91 (51) | 94 (53)
Mean (SD) observed QALYs | $3.52\ (0.99)$ | $3.74\ (0.90)$
Baseline variables |  |
N (%) missing $\mbox{EQ5D}_{0}$ | 6 (3) | 7 (4)
Mean (SD) observed $\mbox{EQ5D}_{0}$ | $0.72\ (0.25)$ | $0.71\ (0.26)$
Correlation between costs and QALYs | $-0.42$ | $-0.07$
Correlation of costs and QALYs | $-0.36$ | $-0.18$
by treatment received |  |
Incremental costs, QALYs and INB of surgery vs medicine
Outcome | Method | estimate ( 95% CI )
Incremental cost |  |
| ITT | $1103\ \ (593,1613)$
| 2sls | $1899\ (1073,2724)$
| 3sls | $1899\ (1073,2724)$
| uBN | $2960\ (2026,3998)$
| uBGN | $2176\ (1356,3031)$
| BFL | $2030\ (1170,2878)$
Incremental QALYs |  |
| ITT | $0.295\ (0.002,0.589)$
| 2sls | $0.516\ (0.103,0.929)$
| 3sls | $0.516\ (0.103,0.929)$
| uBN | $0.568\ (0.181,0.971)$
| uBGN | $0.268\ (-0.229,0.759)$
| BFL | $0.511\ (0.121,0.947)$
INB |  |
| ITT | $7763\ (-1059,16585)$
| 2sls | $13587\ (1101,26073)$
| 3sls | $13587\ (1002,26173)$
| uBN | $14091\ (2485,26086)$
| uBGN | $5869\ (-9204,20740)$
| BFL | $13340\ (1406,26315)$
Table 2: CI Coverage rates (CR) and median width for incremental cost, QALYs, and INB, across scenarios with 30% non-compliance, sample size $n=100$ and moderate correlation between outcomes (with odd rows corresponding to positive correlation $\rho$ and even rows to negative). uBGN was not applied in settings with normal cost data. | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---
$Y_{1}\sim N$ | CR | CIW | CR | CIW | CR | CIW | CR | CIW | CR | CIW
Cost | .952 | .228 | .952 | .228 | .992 | .312 |  |  | .988 | .299
$\rho<0$ | .952 | .229 | .952 | .229 | .993 | .325 |  |  | .986 | .297
QALYs | .946 | .112 | .946 | .112 | .988 | .155 |  |  | .950 | .121
$\rho<0$ | .950 | .113 | .950 | .113 | .992 | .163 |  |  | .950 | .121
INB | .988 | 405 | .953 | 319 | .982 | 398 |  |  | .966 | 376
$\rho<0$ | .900 | 409 | .948 | 475 | .951 | 509 |  |  | .962 | 525
$Y_{1}\sim G$ |  |  |  |  |  |  |  |  |  |
Cost | .952 | .756 | .952 | .756 | .955 | .815 | .941 | .818 | .954 | .823
$\rho<0$ | .942 | .759 | .942 | .759 | .949 | .828 | .936 | .822 | .945 | .811
QALYs | .959 | .113 | .959 | .113 | .993 | .160 | .960 | .122 | .960 | .122
$\rho<0$ | .959 | .113 | .949 | .113 | .995 | .163 | .954 | .122 | .954 | .122
INB | .982 | 829 | .948 | 696 | .958 | 764 | .942 | 748 | .956 | 760
$\rho<0$ | .914 | 833 | .948 | 943 | .930 | 921 | .941 | 1019 | .951 | 1014
$Y_{1}\sim IG$ |  |  |  |  |  |  |  |  |  |
Cost | .951 | .880 | .951 | .880 | .958 | .949 | .904 | .866 | .956 | .945
$\rho<0$ | .950 | .878 | .950 | .878 | .958 | .951 | .905 | .864 | .954 | .932
QALYs | .945 | .112 | .945 | .112 | .991 | .161 | .944 | .120 | .999 | .206
$\rho<0$ | .954 | .112 | .954 | .112 | .993 | .161 | .952 | .120 | .999 | .204
INB | .980 | 944 | .954 | 818 | .959 | 889 | .917 | 814 | .984 | 1001
$\rho<0$ | .917 | 942 | .947 | 1049 | .934 | 1034 | .911 | 1041 | .971 | 1203
Table 3: RMSE for incremental Cost, QALYs and INB across scenarios with 30% non-compliance, moderate correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.4 | 0.058 | 0.060 |  | 0.059
|  | -0.4 | 0.060 | 0.062 |  | 0.061
| QALYs |  |  |  |  |
|  | 0.4 | 0.029 | 0.030 |  | 0.030
|  | -0.4 | 0.029 | 0.030 |  | 0.030
| INB |  |  |  |  |
|  | 0.4 | 83 | 84 |  | 87
|  | -0.4 | 125 | 127 |  | 125
Gamma | Cost |  |  |  |  |
|  | 0.4 | 0.198 | 0.202 | 0.212 | 0.202
|  | -0.4 | 0.200 | 0.204 | 0.212 | 0.203
| QALYs |  |  |  |  |
|  | 0.4 | 0.030 | 0.030 | 0.030 | 0.029
|  | -0.4 | 0.029 | 0.030 | 0.030 | 0.030
| INB |  |  |  |  |
|  | 0.4 | 181 | 184 | 193 | 184
|  | -0.4 | 246 | 251 | 261 | 252
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.4 | 0.230 | 0.232 | 0.252 | 0.232
|  | -0.4 | 0.230 | 0.232 | 0.250 | 0.232
| QALYs |  |  |  |  |
|  | 0.4 | 0.029 | 0.030 | 0.030 | 0.030
|  | -0.4 | 0.029 | 0.030 | 0.030 | 0.030
| INB |  |  |  |  |
|  | 0.4 | 211 | 214 | 231 | 214
|  | -0.4 | 273 | 278 | 296 | 278

##  Appendix A  Sensitivity analyses for the REFLUX RCT

Our motivating example trial had over 50% missing cost-effectiveness outcomes,
with a few missing data in the baseline covariate. As primary analysis, we
performed MI and full Bayesian analyses, that intrinsically obtain a
predictive posterior distribution for these missing values, and reported the
results of these analyses in the main text.

As a sensitivity to the MAR assumption, we present here the complete case
analysis (N=172) of the REFLUX RCT, adjusting for baseline EQ5D. The results
are valid if conditional on baseline EQ5D, the missing data mechanism is
independent of the outcomes.

Table A1: Estimated cost-effectiveness for the REFLUX RCT, over a time horizon
of five years according to ITT and alternative methods for estimating the CACE
based on complete cases. Treatment switches are defined within the first year
post randomisation. Costs and INB numbers rounded to the nearest integer.
Incremental costs, QALYs and INB of surgery vs medicine
---
Outcome | Method | estimate ( 95% CI )
Incremental cost |  |
| ITT | $1626\ (1098,2153)$
| 2sls | $2102\ (1484,2720)$
| 3sls | $2102\ (1484,2720)$
| uBN | $2112\ (1401,2879)$
| uBGN | $2195\ (1481,3065)$
| BFL | $2081\ (1457,2705)$
Incremental QALYs |  |
| ITT | $0.319\ (-0.055,0.694)$
| 2sls | $0.412\ (0.137,0.688)$
| 3sls | $0.412\ (0.137,0.688)$
| uBN | $0.413\ (0.136,0.704)$
| uBGN | $0.285\ (-0.066,0.629)$
| BFL | $0.414\ (0.135,0.703)$
INB |  |
| ITT | $7948.\ (-3259,19155)$
| 2sls | $10275\ (1997,18554)$
| 3sls | $10275\ (1828,18724)$
| uBN | $10280,\ (1947,18980)$
| uBGN | $6362\ -4497,16941)$
| BFL | $10353\ (1789,19165)$

As can be seen in Table A1, the conclusions each model reaches are not
substantially changed from those under MAR (either using MI or full Bayesian),
presented in the main text.

###  A.1 Sensitivity analysis of the BFL to choice of priors

We now study the sensitivity of the BFL model to the choice of priors for the
parameters in the model. Recall that this model can be written as

$\begin{array}[]{l}\left(\begin{array}[]{c}D_{i}\\\ Y_{1i}\\\
Y_{2i}\end{array}\right)\sim N\left\\{\left(\begin{array}[]{c}\mu_{0i}\\\
\mu_{1i}\\\ \mu_{2i}\end{array}\right),\Sigma\right\\}\end{array}\ \ \ \ \ \ \
\
\begin{array}[]{l}\Sigma=\left(\begin{array}[]{ccc}\sigma^{2}_{0}&s_{01}&s_{02}\\\
s_{01}&\sigma^{2}_{1}&s_{12}\\\
s_{02}&s_{12}&\sigma^{2}_{2}\end{array}\right)\end{array}$ (22)

where $s_{ij}=\mbox{cov}(Y_{i},Y_{j})$, and

$\displaystyle\mu_{0i}$ $\displaystyle=$
$\displaystyle\beta_{0,0}+\beta_{1,0}Z_{i}$ $\displaystyle\mu_{1i}$
$\displaystyle=$ $\displaystyle\beta_{0,1}+\beta_{1,1}\beta_{1,0}Z_{i}$
$\displaystyle\mu_{2i}$ $\displaystyle=$
$\displaystyle\beta_{0,2}+\beta_{1,2}\beta_{1,0}Z_{i}$ (23)

Originally, we chose normal priors for the regression coefficients,
$\beta_{m,j}\sim N(0,10^{2})$, for $j\in\\{0,1,2\\},m\in\\{0,1\\}$, and a
Wishart prior for the inverse of $\Sigma$ (Gelman and Hill, 2006).

Here, we first vary the prior distribution for $\Sigma$, and assume a
structured covariance matrix (see Congdon (2007) example 5.8, and Lunn et al.
(2012), example 9.1.4). Thus, we write:

$\begin{array}[]{l}\Sigma=\left(\begin{array}[]{ccc}s_{00}s_{00}&\rho_{1}s_{00}s_{11}&\rho_{2}s_{00}s_{22}\\\
\rho_{1}s_{00}s_{11}&s_{11}s_{11}&\rho_{3}s_{11}s_{22}\\\
\rho_{2}s_{00}s_{22}&\rho_{3}s_{11}s_{22}&s_{22}s_{22}\end{array}\right)\end{array}$
(24)

and assumed the following priors

$s_{ij}\sim N(0,10^{2}),\qquad\rho_{j}\sim\mbox{Unif}[-1,1]$ (25)

In a secondary sensitivity analysis, we assume a Wishart prior for $\Sigma$ as
before, but used uniform priors for the regression coefficients,
$\beta_{m,j}\sim\mbox{Unif}[-10,10]$, for $j\in\\{0,1,2\\},m\in\\{0,1\\}$. We
did these changes on both available cases and complete cases. The results
corresponding to INB are reported in Table A2.

Table A2: Results from sensitivity analysis of the BFL to prior specifications
for the INB of the motivating example: CEA within the REFLUX trial. INB
reported to the nearest integer, in Great British Pounds (GBP).

Data | Priors | INB estimate | 95% CI
---|---|---|---
Available cases | Structured precision matrix1 | $10421$ | $(1878,19187)$
| normal priors for $\beta_{m,j}$ |  |
Available cases | Wishart prior for the precision | $10405$ | $(1733,19330)$
| Unif prior for $\beta_{m,j}$ |  |
Complete cases | Structured precision matrix | 10425 | $(1792,19156)$
| normal priors for $\beta_{m,j}$ |  |
Complete cases | Wishart prior for the precision | $10428$ | $(1894,19279)$
| Unif prior for $\beta_{m,j}$ |  |

  * 1

Following the specifications of equation (25) in this Appendix.

##  Appendix B Supplementary Results from the simulations

Figure 2: Median Bias for scenarios with 70% non-compliance and sample sizes
of (a) $n=100$ (top) and (b) $n=1000$. Results are stratified by cost
distribution, and correlation between cost and QALYs. The dotted line
represents zero bias. Results for 2sls (not plotted) are identical to those
for 3sls; uBGN was not applied to Normal cost data.

(a) n=100 (b) n=1000

Table A3: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 30% non-compliance, moderate correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.949 | 0.072 | 0.949 | 0.072 | 0.989 | 0.094 |  |  | 0.956 | 0.075
|  | -0.4 | 0.957 | 0.072 | 0.957 | 0.072 | 0.992 | 0.097 |  |  | 0.963 | 0.075
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.942 | 0.036 | 0.942 | 0.036 | 0.982 | 0.046 |  |  | 0.937 | 0.035
|  | -0.4 | 0.947 | 0.036 | 0.947 | 0.036 | 0.991 | 0.048 |  |  | 0.944 | 0.035
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.986 | 129 | 0.951 | 102 | 0.977 | 120 |  |  | 0.950 | 101
|  | -0.4 | 0.920 | 129 | 0.956 | 150 | 0.962 | 155 |  |  | 0.949 | 149
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.949 | 0.240 | 0.949 | 0.240 | 0.955 | 0.246 | 0.943 | 0.235 | 0.951 | 0.242
|  | -0.4 | 0.953 | 0.240 | 0.953 | 0.240 | 0.958 | 0.248 | 0.945 | 0.235 | 0.954 | 0.243
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.947 | 0.036 | 0.947 | 0.036 | 0.990 | 0.048 | 0.946 | 0.036 | 0.943 | 0.035
|  | -0.4 | 0.950 | 0.036 | 0.950 | 0.036 | 0.993 | 0.048 | 0.949 | 0.036 | 0.950 | 0.036
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.985 | 262 | 0.950 | 221 | 0.957 | 229 | 0.944 | 218 | 0.949 | 221
|  | -0.4 | 0.920 | 262 | 0.952 | 297 | 0.932 | 276 | 0.946 | 292 | 0.954 | 301
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.4 | 0.949 | 0.283 | 0.949 | 0.283 | 0.958 | 0.290 | 0.902 | 0.248 | 0.947 | 0.281
|  | -0.4 | 0.952 | 0.283 | 0.952 | 0.283 | 0.957 | 0.291 | 0.901 | 0.249 | 0.946 | 0.277
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.954 | 0.036 | 0.954 | 0.036 | 0.990 | 0.048 | 0.942 | 0.035 | 0.967 | 0.038
|  | -0.4 | 0.949 | 0.036 | 0.949 | 0.036 | 0.989 | 0.048 | 0.942 | 0.035 | 0.961 | 0.038
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.976 | 303 | 0.956 | 263 | 0.958 | 268 | 0.914 | 235 | 0.958 | 268
|  | -0.4 | 0.922 | 302 | 0.950 | 336 | 0.932 | 315 | 0.906 | 300 | 0.947 | 333

Table A4: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 70% non-compliance, moderate correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.968 | 0.550 | 0.968 | 0.550 | 0.988 | 1.33 |  |  | 0.991 | 1.13
|  | -0.4 | 0.966 | 0.545 | 0.966 | 0.545 | 0.990 | 1.30 |  |  | 0.985 | 1.08
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.968 | 0.271 | 0.968 | 0.271 | 0.988 | 0.665 |  |  | 0.964 | 0.445
|  | -0.4 | 0.964 | 0.269 | 0.964 | 0.269 | 0.986 | 0.665 |  |  | 0.954 | 0.435
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.992 | 981 | 0.966 | 771 | 0.985 | 1453 |  |  | 0.981 | 1411
|  | -0.4 | 0.917 | 973 | 0.966 | 1131 | 0.957 | 1788 |  |  | 0.958 | 1912
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.966 | 1.68 | 0.966 | 1.68 | 0.964 | 2.640 | 0.944 | 2.21 | 0.960 | 2.73
|  | -0.4 | 0.969 | 1.69 | 0.969 | 1.69 | 0.965 | 2.690 | 0.950 | 2.21 | 0.964 | 2.71
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.972 | 0.267 | 0.972 | 0.267 | 0.994 | 0.624 | 0.959 | 0.348 | 0.965 | 0.42
|  | -0.4 | 0.965 | 0.269 | 0.965 | 0.269 | 0.987 | 0.625 | 0.954 | 0.351 | 0.954 | 0.42
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.988 | 1870 | 0.966 | 1557 | 0.964 | 2456 | 0.947 | 2041 | 0.970 | 2503
|  | -0.4 | 0.939 | 1876 | 0.971 | 2138 | 0.947 | 2984 | 0.958 | 2774 | 0.962 | 3344
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.4 | 0.965 | 1.92 | 0.965 | 1.92 | 0.967 | 3.04 | 0.912 | 2.29 | 0.972 | 3.31
|  | -0.4 | 0.958 | 1.90 | 0.958 | 1.90 | 0.966 | 3.05 | 0.910 | 2.30 | 0.968 | 3.37
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.967 | 0.268 | 0.967 | 0.268 | 0.988 | 0.637 | 0.948 | 0.355 | 1 | 0.766
|  | -0.4 | 0.966 | 0.271 | 0.966 | 0.271 | 0.990 | 0.631 | 0.945 | 0.359 | 0.999 | 0.774
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.990 | 2082 | 0.961 | 1777 | 0.963 | 2827 | 0.912 | 2157 | 0.991 | 3554
|  | -0.4 | 0.936 | 2071 | 0.970 | 2335 | 0.912 | 3389 | 0.924 | 2855 | 0.979 | 4819

Table A5: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 70% non-compliance, moderate correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.954 | 0.169 | 0.954 | 0.169 | 0.986 | 0.223 |  |  | 0.958 | 0.178
|  | -0.4 | 0.953 | 0.169 | 0.953 | 0.169 | 0.991 | 0.230 |  |  | 0.960 | 0.179
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.946 | 0.083 | 0.946 | 0.083 | 0.986 | 0.111 |  |  | 0.942 | 0.084
|  | -0.4 | 0.953 | 0.083 | 0.953 | 0.083 | 0.992 | 0.115 |  |  | 0.947 | 0.083
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.986 | 301 | 0.948 | 238 | 0.976 | 287 |  |  | 0.944 | 240
|  | -0.4 | 0.909 | 301 | 0.944 | 349 | 0.953 | 368 |  |  | 0.942 | 356
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.952 | 0.526 | 0.952 | 0.526 | 0.958 | 0.556 | 0.939 | 0.523 | 0.952 | 0.541
|  | -0.4 | 0.952 | 0.525 | 0.952 | 0.525 | 0.958 | 0.556 | 0.943 | 0.524 | 0.952 | 0.541
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.951 | 0.083 | 0.951 | 0.0832 | 0.995 | 0.115 | 0.947 | 0.085 | 0.949 | 0.084
|  | -0.4 | 0.955 | 0.083 | 0.955 | 0.0832 | 0.995 | 0.115 | 0.954 | 0.084 | 0.955 | 0.085
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.980 | 582 | 0.954 | 486 | 0.932 | 626 | 0.946 | 487 | 0.951 | 496
|  | -0.4 | 0.916 | 581 | 0.953 | 661 | 0.932 | 626 | 0.944 | 662 | 0.953 | 678
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.4 | 0.952 | 0.595 | 0.952 | 0.595 | 0.959 | 0.621 | 0.909 | 0.538 | 0.950 | 0.603
|  | -0.4 | 0.946 | 0.955 | 0.599 | 0.955 | 0.599 | 0.629 | 0.910 | 0.542 | 0.948 | 0.599
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.943 | 0.083 | 0.943 | 0.083 | 0.991 | 0.114 | 0.938 | 0.084 | 0.959 | 0.0918
|  | -0.4 | 0.083 | 0.946 | 0.083 | 0.96 | 0.988 | 0.115 | 0.942 | 0.084 | 0.957 | 0.0905
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.4 | 0.975 | 645 | 0.952 | 552 | 0.954 | 577 | 0.918 | 510 | 0.953 | 571
|  | -0.4 | 0.920 | 649 | 0.952 | 728 | 0.93 | 692 | 0.914 | 668 | 0.950 | 736

Table A6: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 30% non-compliance, high correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.951 | 0.229 | 0.951 | 0.229 | 0.984 | 0.299 |  |  | 0.987 | 0.299
|  | -0.8 | 0.946 | 0.228 | 0.946 | 0.228 | 0.992 | 0.324 |  |  | 0.986 | 0.295
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.954 | 0.113 | 0.954 | 0.113 | 0.985 | 0.148 |  |  | 0.958 | 0.122
|  | -0.8 | 0.949 | 0.112 | 0.949 | 0.112 | 0.992 | 0.162 |  |  | 0.952 | 0.121
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.999 | 408 | 0.953 | 207 | 0.998 | 373 |  |  | 0.990 | 279
|  | -0.8 | 0.867 | 407 | 0.947 | 533 | 0.93 | 514 |  |  | 0.957 | 585
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.955 | 0.757 | 0.955 | 0.757 | 0.947 | 0.764 | 0.922 | 0.681 | 0.961 | 0.822
|  | -0.8 | 0.946 | 0.757 | 0.946 | 0.757 | 0.954 | 0.825 | 0.947 | 0.820 | 0.951 | 0.817
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.952 | 0.113 | 0.952 | 0.113 | 0.986 | 0.152 | 0.928 | 0.109 | 0.957 | 0.123
|  | -0.8 | 0.947 | 0.113 | 0.947 | 0.113 | 0.99 | 0.162 | 0.950 | 0.121 | 0.954 | 0.123
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.999 | 828 | 0.955 | 539 | 0.99 | 707 | 0.942 | 528 | 0.966 | 596
|  | -0.8 | 0.884 | 831 | 0.945 | 1041 | 0.901 | 920 | 0.882 | 902 | 0.95 | 1118
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.8 | 0.954 | 0.883 | 0.954 | 0.883 | 0.948 | 0.894 | 0.836 | 0.745 | 0.961 | 0.956
|  | -0.8 | 0.953 | 0.884 | 0.953 | 0.884 | 0.961 | 0.956 | 0.826 | 0.750 | 0.955 | 0.938
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.954 | 0.113 | 0.954 | 0.113 | 0.989 | 0.154 | 0.900 | 0.108 | 1 | 0.204
|  | -0.8 | 0.951 | 0.113 | 0.951 | 0.113 | 0.993 | 0.163 | 0.916 | 0.108 | 1 | 0.204
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.994 | 947 | 0.952 | 674 | 0.978 | 825 | 0.837 | 591 | 0.990 | 877
|  | -0.8 | 0.893 | 948 | 0.953 | 1155 | 0.910 | 1042 | 0.842 | 993 | 0.969 | 1307

Table A7: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 30% non-compliance, high correlation between outcomes and sample size $n=1000$.uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | .949 | 0.072 | 0.949 | 0.072 | 0.985 | 0.089 |  |  | 0.951 | 0.074
|  | -0.8 | 0.954 | 0.072 | 0.954 | 0.072 | 0.996 | 0.097 |  |  | 0.955 | 0.074
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.943 | 0.036 | 0.943 | 0.036 0 | 0.988 | 0.044 |  |  | 0.938 | 0.035
|  | -0.8 | 0.952 | 0.036 | 0.952 | 0.036 | 0.992 | 0.048 |  |  | 0.951 | 0.036
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 1 | 129 | 0.95 | 65.4 | 1 | 112 |  |  | 0.946 | 66.3
|  | -0.8 | 0.864 | 129 | 0.952 | 168.4 | 0.928 | 154 |  |  | 0.946 | 170.2
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.946 | 0.240 | 0.946 | 0.240 | 0.938 | 0.230 | 0.897 | 0.199 | 0.944 | 0.239
|  | -0.8 | 0.949 | 0.240 | 0.949 | 0.240 | 0.955 | 0.249 | 0.896 | 0.202 | 0.950 | 0.241
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.948 | 0.036 | 0.948 | 0.036 | 0.983 | 0.045 | 0.921 | 0.032 | 0.946 | 0.035
|  | -0.8 | 0.952 | 0.036 | 0.952 | 0.036 | 0.990 | 0.048 | 0.925 | 0.032 | 0.950 | 0.036
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.996 | 263 | 0.950 | 172 | 0.982 | 213 | 0.928 | 153 | 0.948 | 172
|  | -0.8 | 0.883 | 263 | 0.946 | 329 | 0.899 | 277 | 0.898 | 275 | 0.946 | 330
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.8 | 0.946 | 0.283 | 0.946 | 0.283 | 0.936 | 0.271 | 0.684 | 0.217 | 0.942 | 0.279
|  | -0.8 | 0.952 | 0.283 | 0.952 | 0.283 | 0.955 | 0.290 | 0.701 | 0.221 | 0.950 | 0.280
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.950 | 0.0356 | 0.95 | 0.0356 | 0.986 | 0.0457 | 0.833 | 0.032 | 0.964 | 0.038
|  | -0.8 | 0.944 | 0.0356 | 0.944 | 0.0356 | 0.990 | 0.048 | 0.849 | 0.032 | 0.965 | 0.038
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.994 | 303 | 0.952 | 218 | 0.972 | 250 | 0.722 | 173 | 0.954 | 221
|  | -0.8 | 0.89 | 302 | 0.947 | 367 | 0.904 | 315 | 0.722 | 292 | 0.946 | 364

Table A8: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 70% non-compliance, high correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.972 | 0.546 | 0.972 | 0.546 | 0.982 | 1.29 |  |  | 0.988 | 1.09
|  | -0.8 | 0.970 | 0.550 | 0.970 | 0.550 | 0.99 | 1.36 |  |  | 0.988 | 1.07
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.97 | 0.269 | 0.97 | 0.269 | 0.982 | 0.64 |  |  | 0.966 | 0.443
|  | -0.8 | 0.966 | 0.272 | 0.966 | 0.272 | 0.989 | 0.662 |  |  | 0.960 | 0.446
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 1 | 975 | 0.973 | 495 | 1 | 1322 |  |  | 0.993 | 1032
|  | -0.8 | 0.875 | 981 | 0.966 | 1284 | 0.936 | 1858 |  |  | 0.957 | 2096
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.967 | 1.68 | 0.967 | 1.68 | 0.95 | 2.51 | 0.921 | 1.82 | 0.967 | 2.77
|  | -0.8 | 0.965 | 1.70 | 0.965 | 1.70 | 0.965 | 2.75 | 0.904 | 1.88 | 0.964 | 2.80
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.969 | 0.268 | 0.969 | 0.268 | 0.984 | 0.609 | 0.935 | 0.315 | 0.967 | 0.430
|  | -0.8 | 0.966 | 0.272 | 0.966 | 0.272 | 0.989 | 0.626 | 0.930 | 0.324 | 0.958 | 0.440
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 1 | 1871 | 0.963 | 1175 | 0.992 | 2233 | 0.947 | 1397 | 0.975 | 1964
|  | -0.8 | 0.895 | 1885 | 0.965 | 2383 | 0.908 | 3096 | 0.903 | 2631 | 0.954 | 3867
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.8 | 0.96 | 1.92 | 0.96 | 1.92 | 0.947 | 2.92 | 0.839 | 2.03 | 0.966 | 3.42
|  | -0.8 | 0.964 | 1.91 | 0.964 | 1.91 | 0.966 | 3.05 | 0.855 | 2.03 | 0.968 | 3.37
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.963 | 0.271 | 0.963 | 0.271 | 0.983 | 0.645 | 0.906 | 0.322 | 0.999 | 0.789
|  | -0.8 | 0.973 | 0.27 | 0.973 | 0.270 | 0.99 | 0.631 | 0.918 | 0.325 | 0.999 | 0.774
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.998 | 2091 | 0.958 | 1417 | 0.981 | 2654 | 0.863 | 1586 | 0.990 | 3145
|  | -0.8 | 0.902 | 2075 | 0.966 | 2568 | 0.912 | 3389 | 0.866 | 2755 | 0.979 | 4819

Table A9: CI Coverage rates and median width for incremental Cost, QALYs, and INB, across scenarios with 70% non-compliance, high correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. |  |  | 2sls | 3sls | uBN | uBGN | BFL
---|---|---|---|---|---|---|---
Cost distribution |  | $\rho$ | coverage | width | coverage | width | coverage | width | coverage | width | coverage | width
Normal | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.959 | 0.169 | 0.959 | 0.169 | 0.983 | 0.211 |  |  | 0.959 | 0.175
|  | -0.8 | 0.952 | 0.169 | 0.952 | 0.169 | 0.991 | 0.230 |  |  | 0.954 | 0.177
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.957 | 0.083 | 0.957 | 0.083 | 0.983 | 0.106 |  |  | 0.951 | 0.084
|  | -0.8 | 0.954 | 0.083 | 0.954 | 0.083 | 0.991 | 0.115 |  |  | 0.951 | 0.085
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 1 | 302 | 0.952 | 153 | 0.999 | 269 |  |  | 0.952 | 158
|  | -0.8 | 0.873 | 302 | 0.954 | 395 | 0.931 | 369 |  |  | 0.953 | 406
Gamma | Cost |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.952 | 0.527 | 0.952 | 0.527 | 0.938 | 0.514 | 0.890 | 0.442 | 0.948 | 0.535
|  | -0.8 | 0.950 | 0.525 | 0.95 | 0.525 | 0.954 | 0.556 | 0.904 | 0.449 | 0.949 | 0.542
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.958 | 0.0833 | 0.958 | 0.0833 | 0.986 | 0.108 | 0.923 | 0.076 | 0.952 | 0.0841
|  | -0.8 | 0.951 | 0.0832 | 0.951 | 0.0832 | 0.992 | 0.115 | 0.930 | 0.077 | 0.947 | 0.0843
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.998 | 583 | 0.953 | 368 | 0.986 | 479 | 0.923 | 336 | 0.950 | 374
|  | -0.8 | 0.888 | 582 | 0.953 | 734 | 0.902 | 627 | 0.904 | 627 | 0.951 | 753
Inverse | Cost |  |  |  |  |  |  |  |  |  |  |
Gaussian |  | 0.8 | 0.952 | 0.599 | 0.952 | 0.599 | 0.938 | 0.583 | 0.785 | 0.474 | 0.946 | 0.602
|  | -0.8 | 0.959 | 0.600 | 0.959 | 0.600 | 0.962 | 0.630 | 0.786 | 0.482 | 0.955 | 0.604
| QALYs |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.948 | 0.083 | 0.948 | 0.083 | 0.992 | 0.109 | 0.891 | 0.076 | 0.96 | 0.091
|  | -0.8 | 0.956 | 0.083 | 0.956 | 0.083 | 0.992 | 0.115 | 0.896 | 0.077 | 0.968 | 0.091
| INB |  |  |  |  |  |  |  |  |  |  |
|  | 0.8 | 0.996 | 650 | 0.946 | 447 | 0.978 | 542 | 0.800 | 370 | 0.955 | 468
|  | -0.8 | 0.899 | 650 | 0.96 | 801 | 0.915 | 694 | 0.806 | 654 | 0.958 | 807

Table A10: RMSE for incremental Cost, QALYs and INB across scenarios with 30% non-compliance, high correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.8 | 0.059 | 0.060 |  | 0.060
|  | -0.8 | 0.059 | 0.061 |  | 0.060
| QALYs |  |  |  |  |
|  | 0.8 | 0.029 | 0.030 |  | 0.030
|  | -0.8 | 0.029 | 0.030 |  | 0.030
| INB |  |  |  |  |
|  | 0.8 | 54 | 55 |  | 56
|  | -0.8 | 137 | 140 |  | 139
Gamma | Cost |  |  |  |  |
|  | 0.8 | 0.192 | 0.196 | 0.194 | 0.196
|  | -0.8 | 0.196 | 0.200 | 0.205 | 0.200
| QALYs |  |  |  |  |
|  | 0.8 | 0.029 | 0.030 | 0.030 | 0.030
|  | -0.8 | 0.029 | 0.030 | 0.030 | 0.030
| INB |  |  |  |  |
|  | 0.8 | 138 | 141 | 137 | 141
|  | -0.8 | 271 | 276 | 281 | 276
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.8 | 0.226 | 0.229 | 0.272 | 0.229
|  | -0.8 | 0.225 | 0.228 | 0.272 | 0.228
| QALYs |  |  |  |  |
|  | 0.8 | 0.029 | 0.030 | 0.032 | 0.030
|  | -0.8 | 0.028 | 0.029 | 0.031 | 0.029
| INB |  |  |  |  |
|  | 0.8 | 175 | 178 | 210 | 178
|  | -0.8 | 291 | 296 | 345 | 296
Table A11: RMSE for incremental Cost, QALYs and INB across scenarios with 30% non-compliance, moderate correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.4 | 0.018 | 0.018 |  | 0.018
|  | -0.4 | 0.018 | 0.018 |  | 0.018
| QALYs |  |  |  |  |
|  | 0.4 | 0.009 | 0.009 |  | 0.009
|  | -0.4 | 0.009 | 0.009 |  | 0.009
| INB |  |  |  |  |
|  | 0.4 | 26 | 26 |  | 26
|  | -0.4 | 38 | 38 |  | 39
Gamma | Cost |  |  |  |  |
|  | 0.4 | 0.061 | 0.061 | 0.061 | 0.061
|  | -0.4 | 0.061 | 0.061 | 0.061 | 0.061
| QALYs |  |  |  |  |
|  | 0.4 | 0.009 | 0.009 | 0.009 | 0.009
|  | -0.4 | 0.009 | 0.009 | 0.009 | 0.009
| INB |  |  |  |  |
|  | 0.4 | 55 | 55 | 56 | 55
|  | -0.4 | 75 | 75 | 76 | 75
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.4 | 0.072 | 0.072 | 0.075 | 0.072
|  | -0.4 | 0.072 | 0.072 | 0.076 | 0.073
| QALYs |  |  |  |  |
|  | 0.4 | 0.009 | 0.009 | 0.035 | 0.009
|  | -0.4 | 0.009 | 0.009 | 0.035 | 0.009
| INB |  |  |  |  |
|  | 0.4 | 66 | 66 | 69 | 66
|  | -0.4 | 85 | 86 | 89 | 86
Table A12: RMSE for incremental Cost, QALYs and INB across scenarios with 30% non-compliance, high correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.8 | 0.019 | 0.019 |  | 0.019
|  | -0.8 | 0.018 | 0.018 |  | 0.018
| QALYs |  |  |  |  |
|  | 0.8 | 0.009 | 0.009 |  | 0.009
|  | -0.8 | 0.009 | 0.009 |  | 0.009
| INB |  |  |  |  |
|  | 0.8 | 17 | 17 |  | 17
|  | -0.8 | 43 | 43 |  | 44
Gamma | Cost |  |  |  |  |
|  | 0.8 | 0.062 | 0.062 | 0.061 | 0.062
|  | -0.8 | 0.062 | 0.062 | 0.062 | 0.062
| QALYs |  |  |  |  |
|  | 0.8 | 0.009 | 0.009 | 0.009 | 0.009
|  | -0.8 | 0.009 | 0.009 | 0.009 | 0.009
| INB |  |  |  |  |
|  | 0.8 | 44 | 44 | 43 | 44
|  | -0.8 | 85 | 85 | 84 | 85
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.8 | 0.073 | 0.073 | 0.105 | 0.073
|  | -0.8 | 0.072 | 0.072 | 0.103 | 0.072
| QALYs |  |  |  |  |
|  | 0.8 | 0.009 | 0.009 | 0.011 | 0.009
|  | -0.8 | 0.009 | 0.009 | 0.011 | 0.009
| INB |  |  |  |  |
|  | 0.8 | 56 | 56 | 80 | 56
|  | -0.8 | 94 | 94 | 132 | 94
Table A13: RMSE for incremental Cost, QALYs and INB across scenarios with 70% non-compliance, moderate correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.4 | 0.145 | 0.207 |  | 0.181
|  | -0.4 | 0.146 | 0.203 |  | 0.182
| QALYs |  |  |  |  |
|  | 0.4 | 0.074 | 0.104 |  | 0.089
|  | -0.4 | 0.074 | 0.109 |  | 0.088
| INB |  |  |  |  |
|  | 0.4 | 212 | 272 |  | 264
|  | -0.4 | 309 | 396 |  | 383
Gamma | Cost |  |  |  |  |
|  | 0.4 | 0.467 | 0.565 | 0.682 | 0.559
|  | -0.4 | 0.456 | 0.553 | 0.636 | 0.549
| QALYs |  |  |  |  |
|  | 0.4 | 0.072 | 0.097 | 0.082 | 0.085
|  | -0.4 | 0.073 | 0.099 | 0.083 | 0.086
| INB |  |  |  |  |
|  | 0.4 | 426 | 520 | 429 | 500
|  | -0.4 | 570 | 697 | 758 | 687
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.4 | 0.521 | 0.622 | 0.946 | 0.622
|  | -0.4 | 0.516 | 0.635 | 1.096 | 0.635
| QALYs |  |  |  |  |
|  | 0.4 | 0.073 | 0.094 | 0.092 | 0.094
|  | -0.4 | 0.074 | 0.092 | 0.110 | 0.092
| INB |  |  |  |  |
|  | 0.4 | 484 | 583 | 632 | 583
|  | -0.4 | 630 | 865 | 967 | 865
Table A14: RMSE for incremental Cost, QALYs and INB across scenarios with 70% non-compliance, high correlation between outcomes and sample size $n=100$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.8 | 0.147 | 0.201 |  | 0.180
|  | -0.8 | 0.147 | 0.209 |  | 0.182
| QALYs |  |  |  |  |
|  | 0.8 | 0.072 | 0.100 |  | 0.086
|  | -0.8 | 0.075 | 0.109 |  | 0.090
| INB |  |  |  |  |
|  | 0.8 | 131 | 169 |  | 169
|  | -0.8 | 351 | 447 |  | 425
Gamma | Cost |  |  |  |  |
|  | 0.8 | 0.450 | 0.543 | 0.604 | 0.531
|  | -0.8 | 0.454 | 0.546 | 0.563 | 0.538
| QALYs |  |  |  |  |
|  | 0.8 | 0.072 | 0.097 | 0.081 | 0.086
|  | -0.8 | 0.073 | 0.096 | 0.084 | 0.086
| INB |  |  |  |  |
|  | 0.8 | 313 | 377 | 429 | 368
|  | -0.8 | 638 | 768 | 758 | 773
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.8 | 0.532 | 0.637 | 0.902 | 0.637
|  | -0.8 | 0.526 | 0.635 | 1.010 | 0.635
| QALYs |  |  |  |  |
|  | 0.8 | 0.074 | 0.092 | 0.100 | 0.092
|  | -0.8 | 0.074 | 0.092 | 0.096 | 0.092
| INB |  |  |  |  |
|  | 0.8 | 395 | 468 | 615 | 468
|  | -0.8 | 704 | 865 | 1022 | 865
Table A15: RMSE for incremental Cost, QALYs and INB across scenarios with 70% non-compliance, moderate correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.4 | 0.043 | 0.043 |  | 0.043
|  | -0.4 | 0.043 | 0.043 |  | 0.043
| QALYs |  |  |  |  |
|  | 0.4 | 0.022 | 0.022 |  | 0.022
|  | -0.4 | 0.022 | 0.022 |  | 0.022
| INB |  |  |  |  |
|  | 0.4 | 61 | 62 |  | 62
|  | -0.4 | 90 | 91 |  | 91
Gamma | Cost |  |  |  |  |
|  | 0.4 | 0.135 | 0.135 | 0.138 | 0.136
|  | -0.4 | 0.134 | 0.135 | 0.138 | 0.135
| QALYs |  |  |  |  |
|  | 0.4 | 0.021 | 0.021 | 0.021 | 0.021
|  | -0.4 | 0.021 | 0.021 | 0.021 | 0.021
| INB |  |  |  |  |
|  | 0.4 | 124 | 169 | 126 | 125
|  | -0.4 | 168 | 169 | 172 | 170
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.4 | 0.151 | 0.153 | 0.158 | 0.153
|  | -0.4 | 0.152 | 0.153 | 0.158 | 0.153
| QALYs |  |  |  |  |
|  | 0.4 | 0.021 | 0.022 | 0.022 | 0.022
|  | -0.4 | 0.021 | 0.022 | 0.022 | 0.022
| INB |  |  |  |  |
|  | 0.4 | 141 | 142 | 147 | 142
|  | -0.4 | 185 | 187 | 193 | 187
Table A16: RMSE for incremental Cost, QALYs and INB across scenarios with 70% non-compliance, high correlation between outcomes and sample size $n=1000$. uBGN was not applied in settings with normal cost data. Numbers for INB have been rounded to the nearest integer. Cost distribution |  | $\rho$ | 3sls (2sls) | uBN | uBGN | BFL
---|---|---|---|---|---|---
Normal | Cost |  |  |  |  |
|  | 0.8 |  | 0.043 |  | 0.043
|  | -0.8 | 0.043 | 0.043 |  | 0.043
| QALYs |  |  |  |  |
|  | 0.8 | 0.021 | 0.021 |  | 0.021
|  | -0.8 | 0.021 | 0.022 |  | 0.021
| INB |  |  |  |  |
|  | 0.8 | 39 | 39 |  | 39
|  | -0.8 | 100 | 101 |  | 101
Gamma | Cost |  |  |  |  |
|  | 0.8 | 0.135 | 0.137 | 0.135 | 0.137
|  | -0.8 | 0.133 | 0.135 | 0.133 | 0.135
| QALYs |  |  |  |  |
|  | 0.8 | 0.021 | 0.022 | 0.022 | 0.021
|  | -0.8 | 0.021 | 0.021 | 0.021 | 0.021
| INB |  |  |  |  |
|  | 0.8 | 95 | 96 | 94 | 96
|  | -0.8 | 187 | 188 | 186 | 188
Inverse | Cost |  |  |  |  |
Gaussian |  | 0.8 | 0.155 | 0.157 | 0.190 | 0.157
|  | -0.8 | 0.151 | 0.153 | 0.192 | 0.153
| QALYs |  |  |  |  |
|  | 0.8 | 0.021 | 0.021 | 0.023 | 0.021
|  | -0.8 | 0.021 | 0.021 | 0.023 | 0.021
| INB |  |  |  |  |
|  | 0.8 | 118 | 119 | 144 | 119
|  | -0.8 | 201 | 203 | 250 | 203

###  B.1 Software code for performing MI and 3sls in Stata

The code belows performs MI of total costs, total QALYs and the baseline EQ5D
summary score eq5d_b, by chained equations using predictive mean matching with
the 5 nearest neighbours (specified by the option knn(5). We include treatment
received  txreceived as a variable in the imputation models, and perform the
imputations separately by randomised treatment arm, rnd (done by specifying
the option by (rnd) .

After obtaining 50 imputations, 3sls estimates are obtained, and pooled
following Rubin’s rules.

    *MI
    sort rnd
    mi set mlong
    mi register imputed total_qalys total_costs eq5d_b
    mi impute chained (pmm, knn(5) ) total_qalys  eq5d_b total_costs=   txreceived,///
    ¯ by (rnd) add(50)   rseed (110)

    *3sls
    mi estimate, cmdok: reg3 (total_costs total_qalys== txreceived), exog(rnd) ireg
    mi estimate(inb: _b[total_qalys:txreceived] *30000 - _b[total_costs:txreceived]), ///
    ¯cmdok: reg3 (total_costs total_qalys==txreceived), exog(rnd)

    mi estimate (icer: _b[total_costs:txreceived]/_b[total_qalys:txreceived] ), ///
    ¯cmdok: reg3 (total_costs total_qalys== txreceived), exog(rnd)

###  B.2 Software code for performing MI and 3sls in R

The following code performs MI by chain equations using PMM. The default
option for the number of donors is 5, so we do not need to specify this. The
code also carries out 3sls using systemfit in R in each imputed data set.
Rubin’s rules for pooling the estimates across imputed sets are coded
manually.

    rm(list=ls())
    library(MASS,"AER",systemfit,mice,foreign)
    data <- read.dta("../data.dta")

    #the data must be ordered as follows total_qalys  total_costs   txreceived eq5d_b

    #MI using MICE in R
    M<-50  #50 imputations
    data<-data[order(data$rnd),]       # order data by arm

    #we impute separate in each randomised arm
    data0<-subset(data,rnd==0)
    data1<-subset(data,rnd==1)
    #the variable rnd must be dropped from data1 and data0

    #rnd=1
    imp1 <-mice(data1, m = M, method =c("pmm","pmm", "","pmm"),print=FALSE)
    #rnd=0
    imp0 <-mice(data0, m = M, method =c("pmm","pmm", "","pmm"),print=FALSE)

    ##merging data
    mj<-imput_no<- rep(1:M, each=length(data$rnd))
    data_mi <-data[rep(1:nrow(data), M), ]
    data_mi$mj<-imput_no

    data_mi$baseq_imput<-data_mi$qaly_imput<-<-array(NA, dim=M*length(data$rnd))
    data_mi$cost_imput<-array(NA, dim=M*length(data$rnd))

    for (j in 1:M){
      data_mi$baseq_imput[data_mi$mj==j]<-c(complete(imp0,j)$baseq , complete(imp1,j)$baseq)
      data_mi$qaly_imput[data_mi$mj==j]<-c(complete(imp0,j)$qaly  , complete(imp1,j)$qaly )
      data_mi$cost_imput[data_mi$mj==j]<-c(complete(imp0,j)$cost , complete(imp1,j)$cost)
    }
    ######################################################
    ####  3sls
    ######################################################
    coef.c<-var.coef.c<-matrix(NA, nrow=M, ncol=1)
    coef.q<-var.coef.q<-matrix(NA, nrow=M, ncol=1)
    coef.inb<-var.coef.inb<-matrix(NA, nrow=M, ncol=1)
    coef.icer.3sls<-matrix(NA, nrow=M, ncol=1)

    for (j in 1:M){
      eq.cost <-  data_mi$cost_imput[data_mi$mj==j] ~
    ¯¯¯¯¯data_mi$exp[data_mi$mj==j]+ data_mi$baseq_imput[data_mi$mj==j]
      eq.qaly <-  data_mi$qaly_imput[data_mi$mj==j] ~
    ¯¯¯¯¯ data_mi$exp[data_mi$mj==j]+ data_mi$baseq_imput[data_mi$mj==j]
      eq.Exp.Rand <- data_mi$exp[data_mi$mj==j] ~
    ¯¯¯¯¯data_mi$rnd[data_mi$mj==j]+ data_mi$baseq_imput[data_mi$mj==j]

      System3eq <- list( cost=eq.cost, qaly=eq.qaly, trt.received = eq.Exp.Rand )
      reg3 <- systemfit(System3eq, method = "3SLS",
                        inst = ~ data_mi$rnd[data_mi$mj==j] + data_mi$baseq_imput[data_mi$mj==j])

      coef.q[j]<-reg3$coefficients["qaly_data_mi$exp[data_mi$mj == j]"]
      coef.c[j]<-reg3$coefficients["cost_data_mi$exp[data_mi$mj == j]"]
      var.coef.q[j]<-diag(reg3$coefCov)["qaly_data_mi$exp[data_mi$mj == j]"]
      var.coef.c[j]<-diag(reg3$coefCov)["cost_data_mi$exp[data_mi$mj == j]"]
      cov<-reg3$coefCov["qaly_data_mi$exp[data_mi$mj == j]","cost_data_mi$exp[data_mi$mj == j]"]
      coef.inb[j]<- coef.q[j]*30000-coef.c[j]
      var.coef.inb[j]<-30000^2*(var.coef.q[j])+  var.coef.c[j] -(2*30000*cov)
      coef.icer.3sls[j]<-(coef.c[j])/(coef.q[j])
     ¯}
    ###### Apply Rubin’s rules ####
    mean.c <-mean(coef.c)
    var.c <-mean(var.coef.c)
    with.var.c<- var.c   # within-imputation variance
    bet.var.c<-(1/(M-1))*sum((coef.c-mean.c)^2)
    se.c<- sqrt(with.var.c+(1+1/M)*bet.var.c)

    ##Q
    mean.q <-mean(coef.q)
    var.q <-mean(var.coef.q)
    with.var.q<- var.q   # within-imputation variance
    bet.var.q<-(1/(M-1))*sum((coef.q-mean.q)^2)
    se.q<- sqrt(with.var.q+(1+1/M)*bet.var.q)

    ##INB
    mean.inb <-mean(coef.inb)
    var.inb <-mean(var.coef.inb)
    with.var.inb<- var.inb   # within-imputation variance
    bet.var.inb<-(1/(M-1))*sum((coef.inb-mean.inb)^2)
    se.inb<- sqrt(with.var.inb+(1+1/M)*bet.var.inb)

    ##ONLY Mean ICER
    mean.icer <-mean(coef.icer.3sls)

    res.3sls<-c(mean.c, se.c, mean.q, se.q, mean.inb, se.inb, mean.icer)

###  B.3 Model code for BFL in JAGS

The following model can also be run in WinBUGs. This may be necessary if there
is missing outcome data in the analysis at hand, since multivariate normal
nodes cannot be partially observed in the current version of JAGS.

The vector y=(txreceived, total_qalys,total_costs). The other variables are as
before.

    model{
    for (i in 1:n) {
       ¯  eq5d_b[i] ~  dnorm(qb, pr.qb)

        y[i,1:3] ~ dmnorm(mu[i,], pr[,])
        mu[i,1] <- b[1,1] + b[1,2] * rnd[i]
         mu[i,2] <- b[2,1] + b[2,2]*b[1,2] * rnd[i] +  c*eq5d_b[i]
         mu[i,3] <- b[3,1] + b[3,2]* b[1,2] *rnd[i] +  g*eq5d_b[i]
    ¯¯¯}
    #¯¯Priors:
    pr[1:3,1:3]~dwish(R,3)
          R<-diag(1,3)

        for (k in 1:3) {
              b[k,1:2] ~ dmnorm(vm[],vp[,])
        }
        for (k in 1:2) {
              vm[k] <- 0
              vp[k,k] <- 0.1
        }
        vp[1,2] <- 0
        vp[2,1] <- 0
                      pr.qb<- 1 / (qb.sd * qb.sd)
    ¯¯qb.sd <- exp(log.qb.sd)
    ¯¯log.qb.sd ~ dnorm(0, 0.01)
    ¯¯qb ~ dunif(-0.5,1)
    ¯¯c ~ dnorm(0, 0.01)
    ¯¯g ~ dnorm(0, 0.01)
    ¯¯
      #functionals of the parameters
      inb<- 30000*b[2,2]-b[3,2]
      icer<- b[3,2]/(b[2,2]+0.000001) #to avoid dividing by a zero

    }

## References

  * Angrist et al. (1996) Angrist, J. D., Imbens, G. W. and Rubin, D. B. (1996), ‘Identification of causal effects using instrumental variables’, Journal of the American Statistical Association 91(434), pp. 444–455.
  * Angrist et al. (2008) Angrist, J. D., Pischke, J. S. (2008), Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton University Press.
  * Baiocchi et al. (2014) Baiocchi, M., Cheng, J. and Small, D. S. (2014), ‘Instrumental variable methods for causal inference’, Statistics in Medicine 33(13), 2297–2340.
  * Brilleman et al. (2015) Brilleman, S., Metcalfe, C., Peters, T. and Hollingsworth, W. (2015), ‘The reporting of treatment non-adherence and its associated impact on economic evaluations conducted alongside randomised trials: a systematic review.’, Value in Health 19 (1), pp., 99 – 108.
  * Burgess and Thompson (2012) Burgess, S. and Thompson, S. G. (2012), ‘Improving bias and coverage in instrumental variable analysis with weak instruments for continuous and binary outcomes’, Statistics in Medicine 31(15), 1582–1600.
  * Clarke and Windmeijer (2012) Clarke, P. S. and Windmeijer, F. (2012), ‘Instrumental Variable Estimators for Binary Outcomes.’ Journal of the American Statistical Association 107(500):1638–1652.
  * Congdon (2007) Congdon, P. (2007), Bayesian Statistical Modelling, 2nd ed., Wiley Series in Probability and Statistics.
  * Conley et al. (2012) Conley, T. G., Hansen, C. B. and Rossi, P. E. (2012), ‘Plausibly exogenous’, Review of Economics and Statistics 94(1), 260–272.
  * Daniel et al. (2012) Daniel, R. M., Kenward, M. G., Cousens, S. N. and De Stavola, B. L. (2012), ‘Using causal diagrams to guide analysis in missing data problems’, Statistical Methods in Medical Research: 21(3), 243–256.
  * Davidson and MacKinnon (2004) Davidson, R. and MacKinnon, J. G. (2004). Economic theory and methods, New York: Oxford UniversityPress.
  * Didilez et al. (2010) Didelez, V., Meng, S. and Sheehan, N. (2010) ‘ Assumptions of IV Methods for Observational Epidemiology’. Statistical Science,25(1), 22–40.
  * Dodd et al. (2012) Dodd, S., White, I. and Williamson, P. (2012), ‘Nonadherence to treatment protocol in published randomised controlled trials: a review’, Trials 13(1), 84.
  * Gelman and Hill (2006) Gelman, A. and Hill, J. (2006), Data Analysis Using Regression and Multilevel/Hierarchical Models, Analytical Methods for Social Research, Cambridge University Press.
  * Grant et al. (2008) Grant, A. M., Boachie, C., Cotton, S. C., Faria, R., Bojke, L. and Epstein, D. ( 2013)  ‘Clinical and economic evaluation of laparoscopic surgery compared with medical management for gastro-oesophageal reflux disease: a 5-year follow-up of multicentre randomised trial (the REFLUX trial).’, Health Technology Assessment 17(22).
  * Grant et al. (2013) Grant, A., Wileman, S., Ramsay, C., Boyke, L., Epstein, D. and Sculpher, M. (2008), ‘The effectiveness and cost-effectiveness of minimal access surgery amongst people with gastro-oesophageal reflux disease- a uk collaborative study. the REFLUX trial.’ Health Technology Assessment 12 (31).
  * Greene (2002) Greene, W. (2002). Econometric Analysis, Prentice-Hall international editions, Prentice Hall.
  * Hernán and Robins (2006) Hernán, M. A. and Robins, J. M. ‘Instruments for causal inference: an epidemiologist’s dream?’, Epidemiology 17 (4), 360–372.
  * Hirano et al. (2000) Hirano, K. , Imbens, G. W. , Rubin, D. B. and Zhou, X. H. ‘Assessing the effecr of an influenza vaccine in an encouragement design’, Biostatistics 1, 69 –88.
  * Hoch et al. (2006) Hoch, J. S., Briggs, A. H. and Willan, A. R., (2002). ‘Something old, something new, something borrowed, something blue: A framework for the marriage of health econometrics and cost-effectiveness analysis’. Health economics, 11 (5) 415–430.
  * Hughes et al. (2016) Hughes, D., Charles, J., Dawoud, D., Edwards, R. T., Holmes,E. , Jones, C. , Parham, P. , Plumpton, C. , Ridyard, C., Lloyd-Williams, H., Wood, E., and Yeo,S. T. (2016). ‘Conducting economic evaluations alongside randomised trials: Current methodological issues and novel approaches.’ PharmacoEconomics, 1–15.
  * Imbens and Angrist (1994) Imbens, G. W. and Angrist, J. D. (1994), ‘Identification and estimation of local average treatment effects’, Econometrica 62(2), 467–475.
  * Imbens and Rubin (1997) Imbens, G. W. and Rubin, D. B. (1997), ‘Bayesian inference for causal effects in randomized experiments with noncompliance’, The Annals of Statistics 25(1), 305–327.
  * Jo (2002a) Jo B. (2002a) ‘Estimating intervention effects with noncompliance: Alternative model specifications’. Journal of Educational and Behavioral Statistics. 27:385 –420.
  * Jo (2002b) Jo B. (2002b) ‘Model misspecification sensitivity analysis in estimating causal effects of interventions with noncompliance.’ Statistics in Medicine. 21: 3161– 3181.
  * Jo and Muthén (2001) Jo B. and Muthén B. O. (2001) ‘Modeling of intervention effects with noncompliance: a latent variable approach for randomised trials,’ In Marcoulides GA, Schumacker RE, eds. New developments and techniques in structural equation modeling. Lawrence Erlbaum Associates, Mahwah, New Jersey: 57–87.
  * Kleibergen and Zivot (2003) Kleibergen, F. and Zivot, E. (2003), ‘Bayesian and classical approaches to instrumental variable regression’, Journal of Econometrics 114(1), 29 – 72.
  * Lancaster (2004) Lancaster, T. (2004), Introduction to Modern Bayesian Econometrics, Wiley.
  * Latimer et al. (2014) Latimer, N. R., Abrams, K., Lambert, P., Crowther, M., Wailoo, A., Morden, J., Akehurst, R. and Campbell, M. (2014), ‘Adjusting for treatment switching in randomised controlled trials - a simulation study and a simplified two-stage method’, Statistical Methods in Medical Research: 0962280214557578.
  * Lunn et al. (2012) Lunn, D. and Jackson, C. and Best, N. and Thomas, A. and Spiegelhalter, D., (2012) The BUGS Book: A Practical Introduction to Bayesian Analysis, Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis.
  * Mantopoulos et al. (2016) Mantopoulos, T., Mitchell, P.M., Welton, N.J. McManus, R. and Andronis, L. (2016) ‘Choice of statistical model for cost-effectiveness analysis and covariate adjustment: empirical application of prominent models and assessment of their results’, The European Journal of Health Economics 17:(8) 927–938.
  * NICE (2013) NICE (2013) , Guide to the Methods of Technology Appraisal, National Institute for Health and Care Excellence, London, UK.
  * Nixon and Thompson (2005) Nixon, R. M. and Thompson, S. G. (2005), ‘Methods for incorporating covariate adjustment, subgroup analysis and between-centre differences into cost-effectiveness evaluations.’ Health economics 14(12), 1217–29.
  * Rossi et al. (2012) Rossi, P., Allenby, G. and McCulloch, R. (2012), Bayesian Statistics and Marketing, Wiley Series in Probability and Statistics, Wiley.
  * Rubin (1987) Rubin, D. (1987). Multiple Imputation for Nonresponse in Surveys. Chichester: Wiley.
  * Schmidt (1990) Schmidt, P. (1990). ‘Three-stage least squares with different instruments for different equations’, Journal of Econometrics 43(3), 389 – 394.
  * Terza et al. (2008) Terza, J. V., Basu, A. and Rathouz, P. J. (2008). ‘Two-stage residual inclusion estimation: Addressing endogeneity in health econometric modeling’, Journal of Health Economics 27(3), 531 – 543.
  * White et al. (2010) White, I. R. and Carlin, J. B. (2010). ‘Bias and efficiency of multiple imputation compared with complete-case analysis for missing covariate values’. Statistics in Medicine, 28: 2920–2931.
  * White et al. (2011) White, I. R., Royston, P. and Wood, A. M. (2011). ‘Multiple imputation using chained equations: issues and guidance for practice’, Statistics in Medicine, 30 (4), 377–399.
  * Willan (2006) Willan,A. R. (2006). ‘Statistical Analysis of cost-effectiveness data from randomised clinical trials’. Expert Revision Pharmacoeconomics Outcomes Research 6, 337–346.
  * Willan et al. (2004) Willan, A. R., Briggs, A. and Hoch,J. (2004). ‘Regression methods for covariate adjustment and subgroup analysis for non-censored cost-effectiveness data’. Health Econ. 13(5), 461–475.
  * Willan et al. (2003) Willan, A. R., Chen, E., Cook, R. and Lin, D. (2003). ‘Incremental net benefit in randomized clinical trials with qualify-adjusted survival’. Statistics in Medicine 22, 353–362.
  * Zellner (1962) Zellner, A. (1962). ‘An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias’. Journal of the American Statistical Association 57(298), pp. 348–368.
  * Zellner and Huang (1962) Zellner, A. and Huang D. S. (1962). ‘Further Properties of Efficient Estimators for Seemingly Unrelated Regression Equations’. International Economic Review 3(3), pp. 300–313.
  * Zellner and Theil (1962) Zellner, A. and Theil, H. (1962), ‘Three-stage least squares: Simultaneous estimation of simultaneous equations’. Econometrica 30(1), pp. 54–78.
  * Zhang et al. (2014) Zhang, Z., Peluso, M. J., Gross, C. P., Viscoli, C. M. and Kernan, W. N. (2014). ‘Adherence reporting in randomized controlled trials’. Clinical Trials 11(2), 195–204.