Fitting MNAR models in missingHE

When analysing partially-observed data, any statistical method makes either explicit or implicit assumptions about the missing values which can never be verified from the data at hand. Typically, analyses rely on a missing at random (MAR) assumption, that is they assume that the observed data are able to fully explain missingness. However, there is always the chance that this assumption is not correct and that missingness also depends on some values which are not observed, leading to a missing not at random (MNAR) assumption. Thus, it is extremely important that the robustness of the results to a range of alternative missingness assumptions is assessed in sensitivity analysis, including MNAR.

Each of the three types of missingness models implemented in missingHE through the functions selection, pattern, and hurdle can be fitted under MNAR for either or both the effectiveness and cost outcomes. In each model, MNAR assumptions are expressed in terms of some suitably-defined departures from a MAR assumption about the mean effectiveness and cost parameters. This tutorial illustrates how to specify a MNAR model in missingHE using the built-in data set called MenSS, which is directly available when installing and loading missingHE in your R workspace, as a toy example. See Introduction to missingHE for an introductory tutorial of each type of missingness model in missingHE and a general presentation of the data from the MenSS data set.

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 [email protected]. Suggestions on how to improve the package are also very welcome. A development version of the package is also available from the author’s main GitHub repository at https://github.com/AnGabrio/missingHE.

Handling MNAR using selection models

Selection models encode MNAR assumptions by directly modelling the missingness mechanism, that is the model for the missing data indicator, as a function of some partially-observed variables. In missingHE, MNAR is specified by including the outcome variables (\(e\) and/or \(c\)) inside the model formulae of the corresponding missing indicators (model.me and model.mc). For example, a model encoding a MNAR assumption for the missing effectiveness mechanism and a MAR assumption for the missing cost mechanism can be specified by setting

model.me <- me ~  e + ...
model.mc <- mc ~  1 + ...

where me and mc indicate the effectiveness and cost missingness indicator variables, while ... is a suitable form for the combination of covariates to be used in each model. Inclusion of e into model.me instructs missingHE to include the partially-observed effectiveness variable into the logistic regression model for its missingness indicator. Since the regression coefficient capturing the dependence between missingness and the outcome (\(\delta_e\)) cannot be fully estimated from the observed data, its estimation implies some form of MNAR assumption for the missing effectiveness mechanism. In selection models, two sources of external information are used to identify these parameters: informative prior distributions and distributional outcome assumptions.

Informative priors on \(\delta_e\) (or \(\delta_c\) for costs) can be specified in the form of a R list object containing information about the desired prior specification, including name of the parametric distribution assumed (among a set of pre-defined choices) and its corresponding hyperparameter values. This list can then be passed to the selection function using the optional argument prior to overwrite the default priors (standard normals) assumed by missingHE for these parameters.
The distribution assumed for each outcome can be selected among a pre-defined set of parametric choices (selected among popular choices in the literature) using the arguments dist_e and dist_c inside the selection function. Type help(selection) to access the full list of available distributions for \(e\) and \(c\).

As an example, let us fit a MNAR selection model for the effectiveness variables (QALYs) available in the MenSS data set using the following command

> NN.sel1 = selection(data = MenSS, model.eff = e ~ trt + u.0, 
+                     model.cost = c ~ trt, model.me = me ~ e, 
+                     model.mc = mc ~ 1, type = "MNAR", 
+                     n.iter = 1000, ref = 2, 
+                     dist_e = "norm", dist_c = "norm")

The MNAR assumption is encoded into the model by: setting the argument type = "MNAR"; adding the term e inside model.me. As for the outcome model, both effectiveness and costs are assumed to follow normal distributions, with the inclusion of the baseline utilities (u.0) as covariates into the model of \(e\). Since we did not provide custom priors, default standard normals are specified for \(\delta_e\).

In general, missingHE allows a separate specification of assumptions for the effectiveness and cost missingness mechanisms through the arguments model.me and model.mc. When desired, other fully-observed variables can also be included into the missingness models (either under MAR or MNAR assumptions), to improve the estimation of the model paramters. Estimates of quantities of interest under the selected missingness assumptions can then be retrieved through the print function.

> print(NN.sel1)
#>                 mean          sd        2.5%         50%       97.5% Rhat n.eff
#> alpha[1]   0.5599533  0.06907426   0.4575652   0.5495287   0.7207763 1.00  1000
#> alpha[2]   0.0306913  0.03126271  -0.0339923   0.0314613   0.0901279 1.23    11
#> alpha[3]   0.3718716  0.07022373   0.2005000   0.3823835   0.4697024 1.01  1000
#> beta[1]   81.3972696 26.89877688  27.6844562  80.9085134 132.0675030 1.00   880
#> beta[2]   27.2147242 28.18871012 -26.4252899  27.5554026  84.2765172 1.00  1000
#> delta_e   -0.2007956  1.00513938  -1.9452491  -0.3234225   1.7632070 1.01   130
#> gamma_c    0.9016351  0.17221357   0.5780955   0.9007196   1.2308978 1.01   350
#> gamma_e    1.0924554  0.92128733  -0.7363396   1.2038505   2.6780230 1.01   120
#> p_c        0.7100066  0.03508521   0.6406291   0.7110974   0.7739756 1.00   380
#> p_e        0.7131744  0.03471405   0.6442108   0.7134823   0.7820707 1.00  1000
#> s_c      246.5582447 28.36901150 198.7271974 244.9672838 305.0147032 1.00   450
#> s_e        0.0915624  0.00985752   0.0739536   0.0903899   0.1150824 1.00   510
#> tau_c      0.0000171  0.00000387   0.0000107   0.0000167   0.0000253 1.00   450
#> tau_e    123.3572028 26.08544154  75.5060777 122.3942175 182.8445814 1.00   510
#> tmu_c     95.7748598 28.40123049  39.3187191  95.8555556 149.2801288 1.01   390
#> tmu_e      0.9041271  0.01583446   0.8755038   0.9030595   0.9374956 1.02   200

We now consider an alternative MNAR specification where we encode informative priors on \(\delta_e\). In general, it is difficult to attach a specific interpretation to the values for this parameter (and therefore also to its priors) because its effect may vary depending on the type of distributional assumptions made on the outcome. In this applied example we will therefore only illustrate how some already-defined suitable prior information can be passed into the model, rather than focussing on the plausibility of the chosen values within the current context.

We first define our priors by creating a list object called my.prior. Within this list, we create a vector of length three called "delta.prior.e" containing the following elements: name of the prior distribution ("norm"); the hyperprior values for the selected distribution (mean of \(10\) and precision of \(1\)).

> my.prior <- list(
+   "delta.prior.e" = c("norm", 10, 1)
+ )

As a simple exercise, we increase the prior mean of \(\delta_e\) to \(10\) to assess the impact on posterior estimates of a quite informative prior about this parameter (relatively high positive value on the logit scale). We then fit the model using the new prior specification by passing the object my.prior to the argument prior inside selection.

> NN.sel2 = selection(data = MenSS, model.eff = e ~ trt + u.0, 
+                     model.cost = c ~ trt, model.me = me ~ e, 
+                     model.mc = mc ~ 1, type = "MNAR", 
+                     n.iter = 1000, ref = 2, prior = my.prior,
+                     dist_e = "norm", dist_c = "norm")

We can now check the results and compare them to those obtained under the initial MNAR model.

> print(NN.sel2)
#>                 mean         sd        2.5%         50%       97.5% Rhat n.eff
#> alpha[1]   0.6600167  0.0781434   0.5238052   0.6597518   0.8014864 1.02   210
#> alpha[2]   0.0502063  0.0280245  -0.0012221   0.0496505   0.1100923 1.00   740
#> alpha[3]   0.3259157  0.0772075   0.1823807   0.3231588   0.4701350 1.03   340
#> beta[1]   82.3701863 26.5665264  30.9678959  82.9121768 135.4433126 1.01   270
#> beta[2]   25.4637825 28.2280319 -29.3187886  24.6450643  79.7246617 1.00   750
#> delta_e    9.3508834  1.0795372   7.1374088   9.4664208  11.3645440 1.15    15
#> gamma_c    0.9028063  0.1653424   0.5809345   0.9033879   1.2351742 1.00   430
#> gamma_e   -7.9947504  1.1115009  -9.9723971  -8.1114270  -5.7377937 1.15    14
#> p_c        0.7103482  0.0337771   0.6412824   0.7116452   0.7747229 1.00   400
#> p_e        0.7523614  0.0381375   0.6734435   0.7548277   0.8199163 1.00   800
#> s_c      246.3355475 27.3402950 197.8590928 244.9103220 303.8526205 1.00  1000
#> s_e        0.1082358  0.0160752   0.0835517   0.1060888   0.1449185 1.01   120
#> tau_c      0.0000171  0.0000038   0.0000108   0.0000167   0.0000255 1.00  1000
#> tau_e     90.8517357 25.6821358  47.6159450  88.8507670 143.2485224 1.01   120
#> tmu_c     95.8227507 27.5312560  41.4218241  96.4580300 149.6458972 1.01   540
#> tmu_e      0.9739711  0.0243093   0.9333921   0.9711420   1.0265348 1.02   120

We see that, with respect to the results from NN.sel1, the mean effectiveness estimates across both groups (tmu_e) are on average higher. Average mean estimates for each treatment group and associated incremental results for each model can also be inspected by using the summary function. To have a better idea of the impact in terms of cost-effectiveness conclusions, we can use the function ceac.plot, from the package BCEA, to display the cost-effectiveness acceptability curves based on the results from each model.

Cost-effectiveness acceptability curves (CEAC) derived from model 1 fitted under default MNAR prior specification for the missing effectiveness mechanism.

Cost-effectiveness acceptability curves (CEAC) derived from model 2 fitted under custom MNAR prior specification for the missing effectiveness mechanism.

The comparison between the graphs in Figure @ref(fig:figplotceamnarsel1) and Figure @ref(fig:figplotceamnarsel2) shows that CEA conclusions are somewhat affected by the specific assumptions made about the missingness mechanism, therefore suggesting that the results of the model are very robust to the MNAR assumptions considered. It is very important that the scenarios explored are informed based on some external information (e.g. expert opinion) to provide a range of plausible assumptions to assess in sensitivity analysis.

Handling MNAR using pattern mixture models

Pattern mixture models encode MNAR assumptions through the specification of two elements: identifying restrictions and sensitivity parameters. Identifying restrictions consist in modelling restrictions imposed to identify those parameters in each missingness pattern that cannot otherwise be identified from the data. In the context of a bivariate cross-sectional outcome (\(e,c\)), a total of four missingness patterns can be observed based on the combination of possible values taken by the outcome-specific missingness indicators \(d^m=(m_i^c,m^e_i)\). These are: complete cases \((0,0)\); missing costs only \((1,0)\); missing effects only \((0,1)\); fully missing cases \((1,1)\). A typical example of these restrictions is when unidentifiable parameters in a given pattern (e.g. fully missing cases) are identified by matching those from a corresponding identifiable pattern (e.g. complete cases). Under MAR, these restrictions are the only element used to achieve the identification of the model. However, when MNAR assumptions are specified, identifying restrictions are often combined with additional parameters which are used to identify the model under a range of departures from MAR. These so-called “sensitivity parameters” are entirely identified using informative priors encoding evidence that is external to the observed data. Unlike the priors for the MNAR parameters in selection models, sensitivity parameters in pattern mixture models can have more intuitive and easier-to-assess interpretations in terms of their impact on posterior results.

missingHE allows the specification of MNAR assumptions through the arguments restriction and type inside the function pattern. The first argument captures the type of identifying restrictions imposed: either complete case ("CC") or available case ("AC") restrictions. The second argument specifies the type of mechanism assumed for the missing values: either MAR or MNAR. Under MAR, the model is identified only based on the given identifying restrictions imposed at the location level in each unidentified pattern for both outcomes. For example, when "CC" restrictions are selected, the mean effectiveness and cost parameters in the fully missing cases \(\mu_e^{(1,1)}\) and \(\mu_c^{(1,1)}\) are set equal to those from the complete cases \(\mu_e^{(0,0)}\) and \(\mu_c^{(0,0)}\), respectively. Under MNAR, outcome-specific sensitivity parameters (\(\delta_e\), \(\delta_c\)) are additively included within the chosen identification strategy. Using the same example, under MNAR, the mean outcomes in the fully missing cases are then identified as:

\[\begin{equation} \mu_e^{(1,1)}=\mu_e^{(0,0)}+\delta_e \;\;\; \text{and} \;\;\; \mu_c^{(1,1)}=\mu_c^{(0,0)}+\delta_c, (\#eq:pmsp) \end{equation}\]

where \(\delta_e\) and \(\delta_c\) are by default assigned with minimally informative uniform priors between \(0\) ans \(1\) by missingHE (both set to \(0\) under MAR). Sensitivity analysis to alternative MNAR assumption scenarios about the mean outcomes can then be conducted by:

Choosing alternative identifying restriction strategies to be passed to the argument restriction, i.e. either "CC" or "AC".
Specify informative priors on \(\delta_e\) and/or \(\delta_c\) in the form of a R list object containing information about the desired prior specification, including name of the parametric distribution assumed (among a set of pre-defined choices) and its corresponding hyperparameter values. This list can then be passed to the pattern function using the optional argument prior to overwrite the default priors (minimally informative uniforms) assumed by missingHE for these parameters under MNAR.

Going back to our example, based on external information, assume that we expect the mean effectiveness among the missing cases to be on average between \(0.2\) and \(0.3\) lower compared to that among the observed cases across both treatment groups. This information can be encoded into a prior distribution for \(\delta_e\) within a pattern mixture model in a similar way to what shown in Section @ref(sec-sm) for selection models. The following code shows how to instruct missingHE to specify a uniform distribution ("unif"), with hyperprior boundary values of \(0.2\) and \(0.3\), for \(\delta_e\).

> my.prior <- list(
+   "delta.prior.e" = c("unif", -0.3, -0.2)
+ )

Next, we proceed to fit the desired MNAR pattern mixture model using the function pattern, setting the arguments restriction = "CC", type = "MNAR" and overwriting default prior specifications for \(\delta_e\) with those contained in the object my.prior.

> NN.pat2 = pattern(data = MenSS, model.eff = e ~ trt + u.0, 
+                 model.cost = c ~ trt, type = "MNAR", 
+                 restriction = "CC", n.iter = 1000, 
+                 ref = 2, prior = my.prior,
+                 dist_e = "norm", dist_c = "norm")

The function includes the baseline utilities in the model for \(e\) and achieves identification under MNAR using complete case restrictions (restriction = "CC") and informative priors on the sensitivity parameters for the mean \(e\). Economic results in terms of posterior summaries about the mean parameters from the model can be seen using the print function

> print(NN.pat2)
#>                     mean          sd        2.5%         50%      97.5% Rhat n.eff
#> alpha_p[1,1]   0.5503279  0.07384925   0.4029207   0.5479382   0.700684 1.01   220
#> alpha_p[2,1]   0.0327223  0.02744614  -0.0175976   0.0313130   0.085599 1.00   660
#> alpha_p[3,1]   0.3841949  0.07837688   0.2264881   0.3862823   0.537581 1.00   330
#> alpha_p[1,2]   0.5503279  0.07384925   0.4029207   0.5479382   0.700684 1.01   220
#> alpha_p[2,2]   0.0327223  0.02744614  -0.0175976   0.0313130   0.085599 1.00   660
#> alpha_p[3,2]   0.3841949  0.07837688   0.2264881   0.3862823   0.537581 1.00   330
#> delta_e       -0.2508527  0.02890294  -0.2977674  -0.2502550  -0.202738 1.00  1000
#> p_prob[1]      0.2905370  0.03578220   0.2241147   0.2901687   0.361951 1.00  1000
#> p_prob[2]      0.7094630  0.03578220   0.6380491   0.7098313   0.775885 1.00  1000
#> s_c_p[1]     246.3257315 27.12331900 200.1628716 245.0536988 304.749204 1.01   190
#> s_c_p[2]     246.3257315 27.12331900 200.1628716 245.0536988 304.749204 1.01   190
#> s_e_p[1]       0.0914307  0.01013916   0.0729561   0.0906417   0.112686 1.00  1000
#> s_e_p[2]       0.0914307  0.01013916   0.0729561   0.0906417   0.112686 1.00  1000
#> tau_c_p[1]     0.0000171  0.00000372   0.0000108   0.0000167   0.000025 1.01   190
#> tau_c_p[2]     0.0000171  0.00000372   0.0000108   0.0000167   0.000025 1.01   190
#> tau_e_p[1]   123.9650296 26.95010786  78.7513465 121.7150425 187.878502 1.00  1000
#> tau_e_p[2]   123.9650296 26.95010786  78.7513465 121.7150425 187.878502 1.00  1000
#> tmu_c         96.7446667 28.15088041  44.5727956  97.1805039 152.863887 1.00   750
#> tmu_e          0.7284244  0.02659706   0.6764983   0.7289712   0.778808 1.00  1000

For comparison, we also fit the same model under an alternative MNAR specification where we assume a normal distribution, roughly defined on the same value range, for the sensitivity parameters

> my.prior2 <- list(
+   "delta.prior.e" = c("norm", -0.25, 1/(0.05^2))
+ )

which is then passed to the argument prior in the pattern function

> NN.pat3 = pattern(data = MenSS, model.eff = e ~ trt + u.0, 
+                 model.cost = c ~ trt, type = "MNAR", 
+                 restriction = "CC", n.iter = 1000, 
+                 ref = 2, prior = my.prior2,
+                 dist_e = "norm", dist_c = "norm")

We assess the impact on the cost-effectiveness results between the two MNAR specifications by looking at the acceptability curves derived from each model using again the function ceac.plot from the package BCEA.

Cost-effectiveness acceptability curves (CEAC) derived from model 1 fitted under uniform MNAR prior specification for the effectiveness sensitivity parameter.

Cost-effectiveness acceptability curves (CEAC) derived from model 2 fitted under normal MNAR prior specification for the effectiveness sensitivity parameter.

The graphs in Figure @ref(fig:figplotceamnarpat1) and Figure @ref(fig:figplotceamnarpat2) indicate a slightly lower chance for the new intervention to be cost-effective when comparing the results under an informative normal distribution for \(\delta_e\) (NN.pat3) compared to when a uniform distribution is assumed (NN.pat3). Note that, in real applications, the range of values and choice of distribution for the sensitivity parameters should be informed based on some external source of information (e.g. expert opinion) to guide the elicitation process.

Handling MNAR using hurdle models

Hurdle models consist of two-part mixture models, where the outcome variable (i.e. \(e\), \(c\) or both) are split into two components: one associated with identical values (called “structural”) which are assumed to be fixed; and one associated with a set of values (called non-structural) which are assumed to be generated from a given distribution. missingHE models the non-structural outcome values using some parametric distributions and uses logistic regressions to estimate the probability of incurring in a structural value. In particular, the function hurdle requires the user to specify separate formulae for: the models estimating the mean outcome for the non-structural component through the usual arguments model.eff and model.cost; the logistic regressions estimating the structural value probabilities using the arguments model.se and model.sc:

model.se <- se ~  1 + ...
model.sc <- sc ~  1 + ...

where se and sc indicate the effectiveness and cost structural value indicator variables, while ... is a suitable form for the combination of covariates to be used in each logistic regression model. Adopting a nomenclature similar to the one used for describing the missingness mechanism, it is possible to argue that exclusion of any covariates from these models corresponds to a “Structural Completely At Random” (SCAR) assumption, since the probability of a structural value does not depend on any observed variable. Similarly, inclusion of some covariates into these models suggests a “Structural At Random” (SAR) assumption, where the structural probability can only be correctly estimated conditional on some other observed variables. In contrast to selection or pattern mixture models, hurdle models do not provide an explicit framework for encoding informative missingness assumptions given that the missingness process is not directly modelled. However, they can be used to assess the impact of specific MNAR scenarios in terms of arbitrary assumptions about the number of missing individuals assigned to a structural value in the model.

To illustrate how this is implemented in missingHE, let us first consider a standard hurdle model specification under MAR. We specify the model using hurdle to handle both structural ones and zeros in \(e\) and \(c\) from the MenSS data set, where values of one and zero are specified by setting the arguments se = 1 and sc = 0 inside hurdle, respectively. Since these values are assumed to be “fixed”, their occurrence does not affect the choice of a distributional form for the non-structural component of the outcomes. As a result, we can assign distributions that are not defined at \(1\) and \(0\) to the non-structural outcome values, such as Beta and Gamma distributions for \(e\) and \(c\), respectively.

> NN.hur1 = hurdle(data = MenSS, model.eff = e ~ trt + u.0, 
+                  model.cost = c ~ trt, model.se = se ~ 1, 
+                  model.sc = sc ~ 1, type = "SCAR", 
+                  se = 1, sc = 0, n.iter = 1000, ref = 2, 
+                  dist_e = "beta", dist_c = "gamma")

The model assumes that the mechanisms of the structural values in both outomes do not depend on any observed covariate, i.e. a structural completely at random (SCAR) mechanism. The function hurdle automatically assigns all individuals with an observed one effect and zero cost to the structural components of the effectiveness and cost mixture distributions, respectively. In general, we do not know to which component of the mixture the individuals with a missing value should be assigned, as this information cannot be obtained from the data. However, based on some external information that we may have, we arbitrarily impose this assignment (which effectively corresponds to a MNAR mechanism).

In missingHE, this can be achieved by first creating an indicator variable (named s_e), telling for each individual whether a structural value is observed (s_e = 1), not observed (s_e = 0) or unknown (s_e = NA). For example, focussing on the effectiveness variables, we can obtain this indicator by typing

> s_e <- ifelse(MenSS$e == 1, 1, 0)
> 
> #number of ones
> sum(s_e == 1, na.rm = T)
#> [1] 17

Next, for all or some of the individuals with a missing effect value, we set the value of s_e = 1 to assign them to the structural component of the hurdle model. For example, we may believe that it is likely for all individuals aged \(< 22\) to be associated with a perfect health status (i.e. \(e = 1\)). We can encode this information into a new object called mys_e by typing

> mys_e <- ifelse(is.na(s_e) & MenSS$age < 22, 1, s_e)
> 
> #number of ones
> sum(mys_e == 1, na.rm = T)
#> [1] 41

The number of individuals associated with \(e = 1\) has considerably increased with respect to the observed data alone. We can now proceed to fit the hurdle model but replacing the original structural indicator variable of \(e\) with the one newly created by setting the optional argument s_e = mys_e.

> NN.hur2 = hurdle(data = MenSS, model.eff = e ~ trt + u.0, 
+                  model.cost = c ~ trt, model.se = se ~ 1, 
+                  model.sc = sc ~ 1, type = "SCAR", s_e = mys_e,
+                  se = 1, sc = 0, n.iter = 1000, ref = 2, 
+                  dist_e = "beta", dist_c = "gamma")

We can then inspect the posterior results by typing

> print(NN.hur2)
#>                 mean      sd      2.5%       50%     97.5% Rhat n.eff
#> alpha[1,1]   0.66381  0.4701  -0.16582   0.60570   1.68012 1.43     7
#> alpha[2,1]   0.11122  0.2592  -0.36718   0.09521   0.64794 1.16    14
#> alpha[3,1]   1.25280  0.5095   0.14532   1.30753   2.14493 1.43     7
#> alpha[1,2]  11.51292  0.0000  11.51292  11.51292  11.51292 1.00     1
#> alpha[2,2]   0.00000  0.0000   0.00000   0.00000   0.00000 1.00     1
#> alpha[3,2]   0.00000  0.0000   0.00000   0.00000   0.00000 1.00     1
#> beta[1,1]    5.44192  0.2023   5.05263   5.44759   5.84796 1.04    48
#> beta[2,1]    0.33657  0.1956  -0.01148   0.33042   0.73549 1.02    87
#> beta[1,2]  -11.51293  0.0000 -11.51293 -11.51293 -11.51293 1.00     1
#> beta[2,2]    0.00000  0.0000   0.00000   0.00000   0.00000 1.00     1
#> gamma_c     -1.06444  0.3413  -1.78460  -1.05053  -0.44883 1.00  1000
#> gamma_e      0.34367  0.2486  -0.13107   0.33933   0.82242 1.00  1000
#> p_c          0.26171  0.0639   0.14374   0.25912   0.38964 1.00  1000
#> p_e          0.58384  0.0595   0.46728   0.58403   0.69475 1.00  1000
#> s_c[1]     241.18534 51.3715 166.79319 233.20638 372.26140 1.02    71
#> s_c[2]       0.00001  0.0000   0.00001   0.00001   0.00001 1.00     1
#> s_e[1]       0.09421  0.0144   0.07060   0.09296   0.12623 1.03    70
#> s_e[2]       0.00001  0.0000   0.00001   0.00001   0.00001 1.00     1
#> tmu_c      206.13097 37.2175 144.53745 202.49336 291.00370 1.02    85
#> tmu_e        0.94191  0.0114   0.91744   0.94289   0.96247 1.02    90

and we can look at how imputations in each treatment group are carried out based on the updated model using the generic plot function.

#> Loading required namespace: ggthemes

Scatter plots of observed and imputed effectiveness data in the SoC (trt=1) arm based on an hurdle model assuming all individuals younger than 22 years old are associated with a perfect healf status.

As Figure @ref(fig:figplothm2) clearly shows, the model imputes missing effect values in "SoC" as one with very small credible intervals. Finally, we compare the economic results from the two alternative hurdle models using the ceac.plot function from the BCEA package.

Cost-effectiveness acceptability curves (CEAC) derived from the hurdle model fitted under a MAR assumption about the number of structural ones in SoC.

Cost-effectiveness acceptability curves (CEAC) derived from the hurdle models fitted under a MNAR assumption about the number of structural ones in SoC.

According to Figure @ref(fig:figplotceamnarhm1) and Figure @ref(fig:figplotceamnarhm2), the probability of cost-effectiveness for the hurdle model under MAR (NN.hur1) shows a moderate increase alongside the willingness to pay threshold up to \(0.7\), where a plateau is observed. Conversely, results from the MNAR model (NN.hur2) indicate an higher chance of cost-effectiveness up to about \(0.8\) for threshold values above \(30000\). This is because the difference in the number of individuals assigned to a structural one between the treatments is more in favour of "MenSS" (\(22\) vs \(19\)) under the specified MNAR assumptions compared with the number observed under MAR (\(8\) vs \(9\)).