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 a.gabrio@maastrichtuniversity.nl. 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.4645247 0.05805454 0.3771471 0.45264 0.5915139 1.20 12
#> alpha[2] 0.0396586 0.02866330 -0.0150031 0.03978 0.0926211 1.00 1000
#> alpha[3] 0.4707027 0.06003276 0.3393986 0.48226 0.5573449 1.12 18
#> beta[1] 209.1735730 45.60870612 121.6519722 210.23479 294.2662604 1.02 110
#> beta[2] -19.4350843 69.88344424 -154.1825108 -22.14496 114.8012948 1.01 280
#> delta_e -0.5383484 0.93689777 -2.0608956 -0.64034 1.2667487 1.32 9
#> gamma_c 0.9086904 0.17387289 0.5765299 0.90861 1.2421559 1.00 1000
#> gamma_e 1.4008613 0.86014021 -0.2873482 1.49168 2.7933193 1.33 8
#> p_c 0.7114310 0.03553662 0.6402685 0.71272 0.7759390 1.00 1000
#> p_e 0.7142690 0.03558383 0.6452954 0.71594 0.7822661 1.00 320
#> s_c 224.4636258 25.34699051 177.7613398 223.38018 280.0206791 1.00 1000
#> s_e 0.0938185 0.01097212 0.0751840 0.09305 0.1179727 1.03 72
#> tau_c 0.0000206 0.00000466 0.0000128 0.00002 0.0000316 1.00 1000
#> tau_e 118.2719897 27.35250756 71.8517200 115.49780 176.9087561 1.03 72
#> tmu_c 198.9059813 35.04604472 130.2904546 198.21469 268.7510775 1.01 280
#> tmu_e 0.9005966 0.01617143 0.8687151 0.90136 0.9300405 1.04 44
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.6130448 0.06828868 0.4988758 0.6043103 0.7568512 1.29 9
#> alpha[2] 0.0532045 0.02792078 -0.0044570 0.0537036 0.1058863 1.04 48
#> alpha[3] 0.3751754 0.06534372 0.2452718 0.3787663 0.4831373 1.28 9
#> beta[1] 207.7461065 42.14475381 126.5171425 207.7454173 287.0265731 1.00 880
#> beta[2] -18.1003446 67.38455655 -145.5341559 -16.7681368 109.1493530 1.00 1000
#> delta_e 9.2824744 0.79068947 7.8656903 9.2175135 10.8477959 1.23 11
#> gamma_c 0.8995441 0.17088771 0.5863341 0.8960545 1.2584333 1.01 270
#> gamma_e -7.9196580 0.81314691 -9.5119306 -7.8520340 -6.4915169 1.21 12
#> p_c 0.7095999 0.03491746 0.6425236 0.7101380 0.7787563 1.01 260
#> p_e 0.7501926 0.03901182 0.6718877 0.7505270 0.8229357 1.00 1000
#> s_c 224.7536158 25.28077341 184.3112195 222.2413240 285.4473227 1.00 1000
#> s_e 0.1054995 0.01982358 0.0798827 0.1009040 0.1623999 1.02 1000
#> tau_c 0.0000205 0.00000436 0.0000123 0.0000202 0.0000294 1.00 1000
#> tau_e 97.8572500 29.90289459 37.9165946 98.2162108 156.7090571 1.02 1000
#> tmu_c 198.1836603 34.63415667 129.4307678 197.9788028 259.8079372 1.00 360
#> tmu_e 0.9720261 0.02864903 0.9312184 0.9659801 1.0496229 1.01 1000
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.
The comparison between the graphs in Figure 1 and Figure 2 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,
\tag{1}
\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 1 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.5479444 0.07259214 0.3954307 0.5483905 0.6822834 1.00 1000
#> alpha_p[2,1] 0.0315237 0.02802849 -0.0216950 0.0329868 0.0862579 1.00 640
#> alpha_p[3,1] 0.3871789 0.07653578 0.2443980 0.3856356 0.5519797 1.00 1000
#> alpha_p[1,2] 0.5479444 0.07259214 0.3954307 0.5483905 0.6822834 1.00 1000
#> alpha_p[2,2] 0.0315237 0.02802849 -0.0216950 0.0329868 0.0862579 1.00 640
#> alpha_p[3,2] 0.3871789 0.07653578 0.2443980 0.3856356 0.5519797 1.00 1000
#> delta_e -0.2509689 0.02881218 -0.2968230 -0.2509113 -0.2031622 1.00 1000
#> p_prob[1] 0.2905850 0.03567906 0.2238209 0.2902693 0.3633415 1.00 1000
#> p_prob[2] 0.7094150 0.03567906 0.6366585 0.7097307 0.7761791 1.00 1000
#> s_c_p[1] 224.9232498 23.73498334 185.3220233 223.1927938 273.9136114 1.00 1000
#> s_c_p[2] 224.9232498 23.73498334 185.3220233 223.1927938 273.9136114 1.00 1000
#> s_e_p[1] 0.0917887 0.01058907 0.0736300 0.0908433 0.1150196 1.01 1000
#> s_e_p[2] 0.0917887 0.01058907 0.0736300 0.0908433 0.1150196 1.01 1000
#> tau_c_p[1] 0.0000204 0.00000426 0.0000133 0.0000201 0.0000291 1.00 1000
#> tau_c_p[2] 0.0000204 0.00000426 0.0000133 0.0000201 0.0000291 1.00 1000
#> tau_e_p[1] 123.3699895 27.90520431 75.5886091 121.1752491 184.4548702 1.01 1000
#> tau_e_p[2] 123.3699895 27.90520431 75.5886091 121.1752491 184.4548702 1.01 1000
#> tmu_c 195.4763923 33.37473195 127.3752258 194.7342015 261.9939599 1.00 1000
#> tmu_e 0.7279893 0.02654854 0.6754037 0.7278835 0.7779161 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.
The graphs in Figure 3 and Figure 4 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.4932329 0.4615 -0.2568193 0.4784702 1.4479435 1.09 22
#> alpha[2,1] 0.1823860 0.2674 -0.3256950 0.1664204 0.7290027 1.02 88
#> alpha[3,1] 1.4309554 0.5103 0.3484756 1.4601082 2.2846041 1.10 22
#> alpha[1,2] 16.1180956 0.0000 16.1180956 16.1180956 16.1180956 1.00 1
#> alpha[2,2] 0.0000000 0.0000 0.0000000 0.0000000 0.0000000 1.00 1
#> alpha[3,2] 0.0000000 0.0000 0.0000000 0.0000000 0.0000000 1.00 1
#> beta[1,1] 5.4560833 0.2008 5.0577985 5.4615020 5.8493035 1.01 940
#> beta[2,1] 0.3222346 0.1813 -0.0564962 0.3242877 0.6683112 1.00 1000
#> beta[1,2] -16.1180957 0.0000 -16.1180957 -16.1180957 -16.1180957 1.00 1
#> beta[2,2] 0.0000000 0.0000 0.0000000 0.0000000 0.0000000 1.00 1
#> gamma_c -1.0598672 0.3353 -1.7630070 -1.0491395 -0.4375913 1.00 1000
#> gamma_e 0.3551871 0.2427 -0.0918544 0.3541665 0.8220481 1.01 230
#> p_c 0.2624047 0.0629 0.1464141 0.2593904 0.3923151 1.00 1000
#> p_e 0.5866641 0.0581 0.4770525 0.5876275 0.6946709 1.01 240
#> s_c[1] 244.1658669 53.2949 163.6391302 235.8339161 378.9198524 1.01 170
#> s_c[2] 0.0000001 0.0000 0.0000001 0.0000001 0.0000001 1.00 1
#> s_e[1] 0.0932926 0.0131 0.0716440 0.0922705 0.1243735 1.00 1000
#> s_e[2] 0.0000001 0.0000 0.0000001 0.0000001 0.0000001 1.00 1
#> tmu_c 207.4819285 38.4730 141.8072848 204.2785075 293.6964992 1.00 350
#> tmu_e 0.9434835 0.0112 0.9189496 0.9445690 0.9622638 1.02 120
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
As Figure 5 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.
According to Figure 6 and Figure 7, 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\)).
