Longitudinal Models in missingHE
Currently, only longitudinal selection models can be
fitted using missingHE, for which the package provides a
series of customisation options to allow a flexible specification of the
models in terms of modelling assumptions and prior choices. These can be
extremely useful for handling the typical features of trial-based CEA
data, such as non-normality, clustering, and type of missingness
mechanism. This tutorial shows how it is possible to customise different
aspects of longitudinal models using the arguments of each type of
function in the package. Throughout, we will use the built-in dataset
called PBS as a toy example, which is directly available
when installing and loading missingHE in your
R workspace. See the vignette called Introduction to
missingHE for an introductory tutorial of each function in
missingHE and a general presentation of the data used in
trial-based economics evaluations.
If you would like to have more information on the package, or would
like to point out potential issues in the current version, feel free to
contact the maintainer at a.gabrio@maastrichtuniversity.nl. Suggestions on how
to improve the package are also very welcome.
Changing the distributional assumptions
A general concern when analysing trial-based CEA data is that, in
many cases, both effectiveness and costs are characterised by highly
skewed distributions, which may cause standard normality modelling
assumptions to be difficult to justify, especially for small samples.
missingHE allows to chose among a range of parametric
distributions for modelling both outcome variables, which were selected
based on the most common choices in standard practice and the
literature. In the context of a trial-based analysis, health economic
outcomes are typically collected at baseline and a series of follow-up
time points. Traditionally, the analysis is conducted at the aggregate
level using cross-sectional effectiveness and cost outcomes obtained by
combining the outcomes collected at different times as this considerably
simplify the analysis task which does not require a longitudinal model
specification. However, such an approach can be justified only when very
limited missing outcome values occur throughout the time period of the
study. When this is not true, then focussing at the aggregate level may
considerably hinder the validity and reliability of the results in terms
of either bias or efficiency. To deal with this problem, then a
longitudinal model specification is reuiquired which allows to make full
use of all observed outcome values collected in the trial while also
being able to specify different missingness assumptions based on the
different missingness patterns observed across the trial period.
Different approaches are available for modelling longitudinal or
repeated measurement binary outcomes under different missingness
assumptions, such as longitudinal selection models, which are
implemented in missingHE together with a series of
customisation options in terms of model specification. Within the model,
the specific type of distributions for the effectiveness or utility
(\(u\)) and cost (\(c\)) outcome at each time point of the
analysis can be selected by setting the arguments dist_u
and dist_c to specific character names. Available choices
include: Normal ("norm") and Beta ("beta")
distributions for \(u\) and Normal
("norm") and Gamma ("gamma") for \(c\). Distributions for modelling both
discrete and binary effectiveness variables are also available, such as
Poisson ("pois") and Bernoulli (bern)
distributions. The full list of available distributions for each type of
outcome can be seen by using the help function on each
function of the package.
In the PBS dataset the longitudinal health economic
variables are the utilities (from EQ5D-3L questionnaires) and costs
(from clinic records) collected at baseline (time = 1) and two follow-up
times (6 months, time = 2; and 12 months, time = 3). To have an idea of
what the data look like, we can inspect the first entries for each
variable at a specific time (e.g. time = 1) by typing
> #first 10 entries of u and c at time 1
> head(PBS$u[PBS$time == 1], n = 10)
[1] 0.17300001 0.85000002 -0.16599999 0.25800008 0.03800001 1.00000000
[7] 0.41400003 0.71000004 0.88300002 0.81500006
> head(PBS$c[PBS$time == 1], n = 10)
[1] 9214.0 2492.5 15283.0 1858.0 797.0 358.0 803.0 336.0 472.0
[10] 959.0
We can check the empirical histograms of the data at different times,
for example \(u\) at time \(1\) and \(2\) by treatment group by typing
> par(mfrow=c(2, 2))
> hist(PBS$u[PBS$t==1 & PBS$time == 1], main = "Utility at time 1 - Control")
> hist(PBS$u[PBS$t==2 & PBS$time == 1], main = "Utility at time 1 - Intervention")
> hist(PBS$u[PBS$t==1 & PBS$time == 2], main = "Utility at time 2 - Control")
> hist(PBS$u[PBS$t==2 & PBS$time == 2], main = "Utility at time 2 - Intervention")
We can also see that the proportion of missing values in \(u\) at each time is moderate in both
treatment groups.
> #proportions of missing values in the control group
> sum(is.na(PBS$u[PBS$t==1 & PBS$time == 1])) / length(PBS$u[PBS$t==1 & PBS$time == 1])
[1] 0.06617647
> sum(is.na(PBS$u[PBS$t==1 & PBS$time == 2])) / length(PBS$u[PBS$t==1 & PBS$time == 2])
[1] 0.125
> sum(is.na(PBS$u[PBS$t==1 & PBS$time == 3])) / length(PBS$u[PBS$t==1 & PBS$time == 3])
[1] 0.08088235
>
> #proportions of missing values in the intervention group
> sum(is.na(PBS$u[PBS$t==2 & PBS$time == 1])) / length(PBS$u[PBS$t==2 & PBS$time == 1])
[1] 0.0462963
> sum(is.na(PBS$u[PBS$t==2 & PBS$time == 2])) / length(PBS$u[PBS$t==2 & PBS$time == 2])
[1] 0.05555556
> sum(is.na(PBS$u[PBS$t==2 & PBS$time == 3])) / length(PBS$u[PBS$t==2 & PBS$time == 3])
[1] 0.0462963
As an example, we fit a longitudinal selection model assuming Normal
distributions to handle \(u\), and we
choose Gamma distributions to capture the skewness in the costs. We note
that, in case some of individuals have costs that are equal to zero (as
in the PBS dataset), standard parametric distributions with
a positive support are not typically defined at \(0\) (e.g. the Gamma distributions), making
their implementation impossible. Thus, in these cases, it is necessary
to use a trick to modify the boundary values before fitting the model. A
common approach is to add a small constant to the cost variables. These
can be done by typing
> PBS$c <- PBS$c + 0.01
We note that, although simple, this strategy has the potential
drawback that results may be affected by the choice of the constant
added and therefore sensitivity analysis to the value used is typically
recommended.
We are now ready to fit our longitudinal selection model to the
PBS dataset using the following command
> NN.sel=selection_long(data = PBS, model.eff = u ~ 1, model.cost = c ~ u,
+ model.mu = mu ~ gender,
+ model.mc = mc ~ gender, type = "MAR",
+ n.iter = 500, dist_u = "norm", dist_c = "norm", time_dep = "none")
The arguments dist_u = "norm" and
dist_c = "norm" specify the distributions assumed for the
outcome variables and, in the model of \(c\), we also take into account the possible
association between the two outcomes at each time by uncluding the
utility as a covariate (u). According to the type of
distributions chosen, missingHE automatically models the
dependence between covariates and the mean outcome on a specific scale
to reduce the chance of incurring into numerical problems due to the
constraints of the distributions. For example, for both Poisson and
Gamma distributions means are modelled on the log scale, while for Beta
and Bernoulli distirbutions they are modelled on the logit scale. To see
the modelling scale used by missingHE according to the type
of distribution selected, use the help command on each
function of the package. The optional argument time_dep
allows to specify the type of time dependence structure assumed by the
model, here corresponding to no temporal dependence
("none"). An alternative choice available in
missingHE is to set time_dep = "AR1", which
assumes an autoregressive of order one time structure.
The model assumes MAR conditional on gender as auxiliary
variable for predicting missingness in both outcomes. We can look at how
the model generate imputations for the outcomes by treatment group using
the generic function plot. For example, we can look at how
the missing \(u\) at time \(3\) are imputed by typing
> plot(NN.sel, outcome = "effects", time_plot = 3)
Summary results of our model from a statistical perspective can be
inspected using the command coef, which extracts the
estimates of the mean regression coefficients for \(u\) and \(c\) by treatment group and time point. By
default, the lower and upper bounds provide the \(95\%\) credibile intervals for each
estimate (based on the \(0.025\) and
\(0.975\) quantiles of the posterior
distribution). For example, we can inspect the coefficients in the
utility and cost models at time \(2\)
by typing
> coef(NN.sel, time = 2)
$Comparator
$Comparator$Effects
mean sd lower upper
(Intercept) 0.502 0.032 0.441 0.562
$Comparator$Costs
mean sd lower upper
(Intercept) 1443.615 165.982 1123.037 1754.898
u -2062.772 437.469 -2898.600 -1219.158
$Reference
$Reference$Effects
mean sd lower upper
(Intercept) 0.641 0.033 0.575 0.7
$Reference$Costs
mean sd lower upper
(Intercept) 2822.054 234.262 2372.573 3285.865
u -1903.106 684.131 -3346.278 -598.933
The entire posterior distribution for each parameter of the model can
also be extracted from the output of the model by typing
NN.sel$model_output, which returns a list object containing
the posterior estimates for each model parameter. An overall summary of
the economic analysis based on the model estimates can be obtained using
the summary command
> summary(NN.sel)
Cost-effectiveness analysis summary
Comparator intervention: intervention 1
Reference intervention: intervention 2
Parameter estimates under MAR assumption
Comparator intervention
mean sd LB UB
mean effects (t = 1) 0.491 0.019 0.46 0.522
mean costs (t = 1) 2888.862 388.703 2197.963 3488.828
Reference intervention
mean sd LB UB
mean effects (t = 2) 0.615 0.021 0.58 0.645
mean costs (t = 2) 5677.399 292.016 5204.342 6142.506
Incremental results
mean sd LB UB
delta effects 0.123 0.029 0.075 0.17
delta costs 2788.537 506.129 1942.413 3661.652
ICER 22631.358
which shows summary statistics for the mean effectiveness and costs
in each treatment group, for the mean differentials and the estimate of
the ICER. It is important to clarify how these results are obtained in
that they pertain summary statistics for aggregated health economic
outcomes (e.g. QALYs and Total costs) which are retrieved by combining
the posterior results from the longitudinal model for each type of
outcome and time point of the analysis. More specifically, the function
selection_long, after deriving the posterior results from
the model, applies common CEA techniques for computing aggregate
variables based on the posterior mean results of the effectiveness and
cost variables at each time point of the analysis. These include: Area
Under the Curve approach for \(u\) and
total sum over the follow-up period for \(c\). By default, missingHE
assumes a value of \(0.5\) and \(1\) for the weights used to respectively
calculate the aggregate mean QALYs and Total costs over the time period
of the study. When desired, users can provide their own weights as
additional argument named qaly_calc and
tcost_calc that can be passed into the function
selection_long when fitting the model.
Including random effects terms
For each type of model, missingHE allows the inclusion
of random effects terms to handle clustered data. To be precise, the
term random effects does not have much meaning within a
Bayesian approach since all parameters are in fact random variables.
However, this terminology is quite useful to explain the structure of
the model.
We show how random effects can be added to the model of \(u\) and \(c\) within a longitudinal selection model
fitted to the PBS dataset using the function
selection_long. The clustering variable over which the
random effects are specified is the factor variable site,
representing the centres at which data were collected in the trial.
Using the same distributional assumptions of the selection model, we fit
the pattern mixture model by typing
> PBS$site <- factor(PBS$site)
> NN.re=selection_long(data = PBS, model.eff = u ~ 1 + (1 | site),
+ model.cost = c ~ u + (1 | site), model.mu = mu ~ gender,
+ model.mc = mc ~ gender, type = "MAR",
+ n.iter = 500, dist_u = "norm", dist_c = "norm", time_dep = "none")
The function fits a random intercept only model for \(u\) and \(c\), as indicated by the notation
(1 | site)and (1 | site). In both models,
site is the clustering variable over which the random
coefficients are estimated. missingHE allows the user to
choose among different clustering variables for the model of \(u\) and \(c\) if these are available in the dataset.
Random intercept and slope models can be specified using the notation
(1 + gender | site), where gender is the name
of a covariate which should be inlcuded as fixed effects in the
corresponding outcome model. Finally, it is also possible to specify
random slope only models using the notation
(0 + gender | site), where 0 indicates the
removal of the random intercept.
Coefficient estimates for the random effects at each time can be
extracted using the coef function and setting the argument
random = TRUE (if set to FALSE then the fixed
effects estimates are displayed).
> coef(NN.re, time = 3, random = TRUE)
$Comparator
$Comparator$Effects
mean sd lower upper
(Intercept) 1 4.777 0.422 4.137 5.344
(Intercept) 2 4.797 0.430 4.159 5.386
(Intercept) 3 4.715 0.429 4.065 5.301
(Intercept) 4 4.839 0.421 4.209 5.440
(Intercept) 5 4.747 0.427 4.102 5.316
(Intercept) 6 4.752 0.427 4.111 5.349
(Intercept) 7 4.699 0.426 4.038 5.273
(Intercept) 8 4.749 0.424 4.091 5.331
(Intercept) 9 4.782 0.430 4.111 5.388
(Intercept) 10 4.916 0.439 4.286 5.601
(Intercept) 11 4.754 0.428 4.103 5.336
(Intercept) 12 4.754 0.422 4.133 5.349
$Comparator$Costs
mean sd lower upper
(Intercept) 1 -8.493 74.892 -177.766 131.328
(Intercept) 2 -5.464 75.301 -173.826 144.953
(Intercept) 3 -2.938 74.825 -162.873 157.590
(Intercept) 4 -4.574 76.058 -163.989 159.071
(Intercept) 5 6.555 71.536 -130.067 168.417
(Intercept) 6 -3.768 72.618 -153.317 160.399
(Intercept) 7 -4.558 73.906 -186.342 145.056
(Intercept) 8 -1.918 71.450 -147.770 157.473
(Intercept) 9 -4.665 71.335 -150.948 143.872
(Intercept) 10 1.420 74.571 -144.955 157.894
(Intercept) 11 -11.235 75.023 -186.677 128.839
(Intercept) 12 3.091 70.912 -144.800 142.586
$Reference
$Reference$Effects
mean sd lower upper
(Intercept) 1 0.712 0.422 -0.014 1.294
(Intercept) 2 0.767 0.436 -0.008 1.393
(Intercept) 3 0.702 0.421 -0.017 1.272
(Intercept) 4 0.720 0.423 -0.005 1.318
(Intercept) 5 0.680 0.418 -0.033 1.252
(Intercept) 6 0.678 0.415 -0.026 1.247
(Intercept) 7 0.745 0.431 0.004 1.322
(Intercept) 8 0.759 0.432 0.001 1.344
(Intercept) 9 0.651 0.407 -0.032 1.198
(Intercept) 10 0.668 0.412 -0.035 1.234
(Intercept) 11 0.755 0.432 0.002 1.381
$Reference$Costs
mean sd lower upper
(Intercept) 1 2.368 59.647 -122.734 137.984
(Intercept) 2 -6.763 58.955 -130.632 112.049
(Intercept) 3 -3.104 58.607 -124.665 122.053
(Intercept) 4 -9.590 59.754 -146.144 121.782
(Intercept) 5 -14.998 60.737 -158.812 110.097
(Intercept) 6 5.453 62.314 -123.702 138.898
(Intercept) 7 -11.447 62.219 -142.054 131.603
(Intercept) 8 -1.988 63.197 -132.991 134.336
(Intercept) 9 -5.102 63.259 -144.608 129.732
(Intercept) 10 5.547 58.731 -111.926 132.147
(Intercept) 11 -10.979 59.127 -140.742 98.179
For both \(u\) and \(c\) models, summary statistics for the
random coefficient estimates are displayed for each treatment group and
each of the clusters in site.
Changing the priors
By default, all models in missingHE are fitted using
vague prior distributions so that posterior results are essentially
derived based on the information from the observed data alone. This
ensures a rough approximation to results obtained under a frequentist
approach based on the same type of models.
However, in some cases, it may be reasonable to use more informative
priors to ensure a better stability of the posterior estimates by
restricting the range over which estimates can be obtained. For example
if, based on previous evidence or knowledge, the range over which a
specific parameter is likely to vary is known, then priors can be
specified so to give less weight to values outside that range when
deriving the posterior estimates. However, unless the user is familiar
with the choice of informative priors, it is generally recommended not
to change the default priors of missingHE as the unintended
use of informative priors may substantially drive posterior estimates
and lead to incorrect results.
For each type of model in missingHE, priors can be
modified using the argument prior, which allows to change
the hyperprior values for each model parameter. The interpretation of
the prior values change according to the type of parameter and model
considered. For example, we can fit a longitudinal selection model to
the PBS dataset using more informative priors on some
parameters.
Prior values can be modified by first creating a list object which,
for example, we call my.prior. Within this list, we create
a number of elements (vectors of length two) which should be assigned
specific names based on the type of parameters which priors we want to
change.
> my.prior <- list(
+ "alpha0.prior" = c(0 , 0.0000001),
+ "beta0.prior" = c(0, 0.0000001),
+ "gamma0.prior.c"= c(0, 1),
+ "gamma.prior.c" = c(0, 0.01),
+ "sigma.prior.c" = c(0, 100)
+ )
The names above have the following interpretations in terms of the
model parameters:
"alpha0.prior" is the intercept of the model of
\(u\). The first and second elements
inside the vector for this parameter are the mean and precision (inverse
of variance) that should be used for the normal prior given to this
parameter by missingHE.
"beta0.prior" is the intercept of the model of \(c\). The first and second elements inside
the vector for this parameter are the mean and precision (inverse of
variance) that should be used for the normal prior given to this
parameter by missingHE.
"gamma0.prior.c" is the intercept of the model of
\(mc\). The first and second elements
inside the vector for this parameter are the mean and precision (inverse
of variance) that should be used for the logistic prior given to this
parameter by missingHE.
"gamma.prior.c" are the regression coefficients
(exclusing the intercept) of the model of \(mc\). The first and second elements inside
the vector for this parameter are the mean and precision (inverse of
variance) that should be used for the normal priors given to each
coefficient by missingHE.
"sigma.prior.c" is the standard deviation of the
model of \(c\). The first and second
elements inside the vector for this parameter are the lower and upper
bounds that should be used for the uniform prior given to this parameter
by missingHE.
The values shown above are the default values set in
missingHE for each of these parameters. It is possible to
change the priors by providing different values, for example by
increasing the precision for some of the coefficient estimates or
decreasing the upper bound for standard deviation parameters. Different
names should be used to indicate for which parameter the prior should be
modified, keeping in mind that the full list of names that can be used
varies depending on the type of models and modelling assumptions
specified. The full list of parameter names for each type of model can
be assessed using the helpcommand on each function of
missingHE.
We can now fit the hurdle model using our priors by typing
> NN.prior=selection_long(data = PBS, model.eff = u ~ 1 + (1 | site),
+ model.cost = c ~ u + (1 | site), model.mu = mu ~ gender,
+ model.mc = mc ~ gender, type = "MAR",
+ n.iter = 500, dist_u = "norm", dist_c = "norm", time_dep = "none",
+ prior = my.prior)
Finally, we can inspect the statistical results from the model at
time \(3\) by typing
> coef(NN.prior, random = FALSE, time = 3)
$Comparator
$Comparator$Effects
mean sd lower upper
(Intercept) 2.913 0.382 1.828 3.469
$Comparator$Costs
mean sd lower upper
(Intercept) 1403.267 330.903 734.967 2033.578
u -1412.674 1009.988 -3431.254 679.444
$Reference
$Reference$Effects
mean sd lower upper
(Intercept) -0.305 1.183 -1.673 1.078
$Reference$Costs
mean sd lower upper
(Intercept) 2853.720 192.780 2468.538 3242.845
u -1061.623 539.509 -2146.405 -28.792
and
> coef(NN.prior, random = TRUE, time = 3)
$Comparator
$Comparator$Effects
mean sd lower upper
(Intercept) 1 -2.423 0.388 -2.975 -1.343
(Intercept) 2 -2.410 0.392 -2.973 -1.324
(Intercept) 3 -2.488 0.385 -3.039 -1.380
(Intercept) 4 -2.363 0.392 -2.948 -1.328
(Intercept) 5 -2.456 0.385 -2.995 -1.374
(Intercept) 6 -2.446 0.388 -2.999 -1.362
(Intercept) 7 -2.496 0.376 -3.035 -1.407
(Intercept) 8 -2.451 0.384 -3.009 -1.375
(Intercept) 9 -2.425 0.393 -3.002 -1.331
(Intercept) 10 -2.292 0.421 -2.935 -1.252
(Intercept) 11 -2.446 0.384 -2.998 -1.377
(Intercept) 12 -2.445 0.390 -3.009 -1.365
$Comparator$Costs
mean sd lower upper
(Intercept) 1 2.214 60.465 -133.780 94.893
(Intercept) 2 2.264 64.379 -125.602 129.749
(Intercept) 3 3.289 58.603 -120.273 99.872
(Intercept) 4 6.973 60.166 -117.305 112.142
(Intercept) 5 10.470 61.587 -107.603 131.094
(Intercept) 6 0.032 65.208 -139.383 113.221
(Intercept) 7 4.713 60.664 -114.812 111.643
(Intercept) 8 3.316 60.012 -123.457 121.502
(Intercept) 9 3.300 63.544 -126.000 106.322
(Intercept) 10 3.555 59.645 -122.852 105.364
(Intercept) 11 1.197 61.502 -125.707 117.566
(Intercept) 12 2.789 60.398 -128.851 111.930
$Reference
$Reference$Effects
mean sd lower upper
(Intercept) 1 0.921 1.187 -0.469 2.309
(Intercept) 2 0.982 1.182 -0.434 2.365
(Intercept) 3 0.917 1.188 -0.514 2.302
(Intercept) 4 0.928 1.182 -0.472 2.319
(Intercept) 5 0.893 1.189 -0.504 2.294
(Intercept) 6 0.884 1.186 -0.516 2.267
(Intercept) 7 0.954 1.189 -0.448 2.322
(Intercept) 8 0.968 1.181 -0.430 2.347
(Intercept) 9 0.861 1.184 -0.576 2.254
(Intercept) 10 0.872 1.182 -0.550 2.257
(Intercept) 11 0.966 1.185 -0.437 2.318
$Reference$Costs
mean sd lower upper
(Intercept) 1 16.437 72.968 -117.404 173.488
(Intercept) 2 -6.281 69.757 -153.603 126.842
(Intercept) 3 2.854 70.069 -128.847 141.100
(Intercept) 4 7.369 63.861 -131.735 145.198
(Intercept) 5 -8.363 65.953 -159.498 122.817
(Intercept) 6 5.941 66.065 -123.337 139.964
(Intercept) 7 0.202 71.886 -147.109 138.875
(Intercept) 8 2.823 71.076 -154.977 144.029
(Intercept) 9 -1.760 69.268 -162.141 129.411
(Intercept) 10 16.798 73.703 -132.601 175.583
(Intercept) 11 -0.431 70.695 -149.705 146.676