inst/doc/Model_customisation_in_missingHE.R

In missingHE: Missing Outcome Data in Health Economic Evaluation

## ----echo = FALSE, message = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(prompt = TRUE, highlight = F, background = '#FFFFFF',
                      collapse = T, comment = "#>")
library(missingHE)
library(bookdown)
library(ggplot2)

options(width = 300)
set.seed(1014)

## ----seldist, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE, results='hide'-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#rename trt levels
MenSS$trt <- factor(MenSS$trt)
levels(MenSS$trt) <- c("SoC", "MenSS")
MenSS$e <- MenSS$e - 0.01 #ensure no ones QALYs occur
MenSS$c <- MenSS$c + 0.01 #ensure no zero costs

#fit models with different distributions for outcomes
#1=Normal-Normal
sm1_nn <- selection(data = MenSS, dist_e = "norm", dist_c = "norm",
                    model.eff = e ~ trt + u.0, model.cost = c ~ trt + e,
                    model.me = me ~ 1, model.mc = mc ~ 1,
                    type = "MAR", n.iter = 1000, ref = 2)
#2=Normal-Gamma
sm1_ng <- selection(data = MenSS, dist_e = "norm", dist_c = "gamma",
                    model.eff = e ~ trt + u.0, model.cost = c ~ trt + e,
                    model.me = me ~ 1, model.mc = mc ~ 1,
                    type = "MAR", n.iter = 1000, ref = 2)
#3=Beta-Gamma
sm1_bg <- selection(data = MenSS, dist_e = "beta", dist_c = "gamma",
                    model.eff = e ~ trt + u.0, model.cost = c ~ trt + e,
                    model.me = me ~ 1, model.mc = mc ~ 1,
                    type = "MAR", n.iter = 1000, ref = 2)

## ----selpic, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#estimate looic for all models based on complete cases
looic_m123 <- pic(x = list(sm1_nn, sm1_ng, sm1_bg),
                  criterion = "looic", cases = "cc")
#print criteria for each model
looic_m123$pic

## ----selpeff, echo=TRUE, eval=FALSE, message=FALSE, warning=FALSE, error=FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# #print peff values
# looic_m123$peff

## ----selppc, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#compare densities of observed vs replicated outcome data under each model
ppc_dens_m1 <- ppc(x = sm1_nn, type = "dens_overlay", outcome = "both",
                   ndisplay = 20, trt = "none")
ppc_dens_m2 <- ppc(x = sm1_ng, type = "dens_overlay", outcome = "costs",
                   ndisplay = 20, trt = "none")
ppc_dens_m3 <- ppc(x = sm1_bg, type = "dens_overlay", outcome = "effects",
                   ndisplay = 20, trt = "none")

## ----figplotppc1, echo=FALSE, eval=TRUE, tidy=TRUE, error=FALSE, warning=FALSE, dpi=300, fig.show='hold',fig.cap="Posterior densities of 20 replicated data under model 1 compared to the empirical density of the original data.", out.width='100%', fig.pos='h', out.extra=''---------------------------
ppc_dens_m1

## ----figplotppc2, echo=FALSE, eval=TRUE, tidy=TRUE, error=FALSE, warning=FALSE, dpi=300, fig.show='hold',fig.cap="Posterior densities of 20 replicated data under model 2 compared to the empirical density of the original data.", out.width='100%', fig.pos='h', out.extra=''---------------------------
ppc_dens_m2

## ----figplotppc3, echo=FALSE, eval=TRUE, tidy=TRUE, error=FALSE, warning=FALSE, dpi=300, fig.show='hold',fig.cap="Posterior densities of 20 replicated data under model 3 compared to the empirical density of the original data.", out.width='100%', fig.pos='h', out.extra=''---------------------------
ppc_dens_m3

## ----selmar, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE, results='hide'--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#fit selection model 1 (Normal-Normal) under MAR conditional on u.0 and age
sm1_nn_mar <- selection(data = MenSS, dist_e = "norm", dist_c = "norm",
                    model.eff = e ~ trt + u.0, model.cost = c ~ trt + e,
                    model.me = me ~ u.0 + age, model.mc = mc ~ u.0 + age,
                    type = "MAR", n.iter = 1000, ref = 2)

## ----selprint, echo=TRUE, eval=TRUE, tidy=TRUE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#print posterior summaries for fixed effects
print(x = sm1_nn_mar, display = "fixed")

## ----selpe, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#store posterior samples of p_e and compute some custom summaries
p_e <- sm1_nn_mar$model_output$model$BUGSoutput$sims.list$p_e
summary(c(p_e))

## ----selcoef, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#print coefficient estimates from model.eff and model.cost
coef(sm1_nn_mar, random = FALSE, digits = 2)

## ----selcov, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE, results='hide'--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#fit selection model 1 (Normal-Normal) with covariates into both outcome and missingness models
sm1_nn_cov <- selection(data = MenSS, dist_e = "norm", dist_c = "norm",
                    model.eff = e ~ trt + u.0 + age, 
                    model.cost = c ~ trt + age + employment + e,
                    model.me = me ~ u.0 + age, model.mc = mc ~ age + employment,
                    type = "MAR", n.iter = 1000, ref = 2)

## ----selcovre, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE, results='hide'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#fit selection model 1 (Normal-Normal) with covariates and random intercepts
#into the model
sm1_nn_cov_re <- selection(data = MenSS, dist_e = "norm", dist_c = "norm",
                    model.eff = e ~ trt + u.0 + age + (1 | site), 
                    model.cost = c ~ trt + age + employment + e + (1 | site),
                    model.me = me ~ u.0 + age, model.mc = mc ~ age + employment,
                    type = "MAR", n.iter = 1000, ref = 2)

## ----figplotdiag, echo=FALSE, eval=TRUE, tidy=TRUE, error=FALSE, warning=FALSE, dpi=300, fig.show='hold',fig.cap="Checking convergence using the diagnostic function for a family of model parameters estimated in missingHE, for example through inspection of the autocorrelation plots.", out.width='100%', fig.pos='h', out.extra=''----
diagnostic(x = sm1_nn_cov, type = "acf", param = "alpha")

## ----selprior, echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE, error=FALSE, results='hide'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#update priors on outcome regression coefficients and standard deviations
myprior <- list("beta.prior" = c("norm", 0, 0.01), 
                "beta_f.prior" = c("norm", 0, 0.001),
                "alpha.prior" = c("norm", 0, 0.01),
                "sigma.prior.e" = c("unif", 0, 5),
                "sigma.prior.c" = c("unif", 0, 1000))
#update model with new priors
sm2_nn_cov <- selection(data = MenSS, dist_e = "norm", dist_c = "norm",
                    model.eff = e ~ trt + u.0 + age, 
                    model.cost = c ~ trt + age + employment + e,
                    model.me = me ~ u.0 + age, model.mc = mc ~ age + employment,
                    type = "MAR", n.iter = 1000, ref = 2,
                    prior = myprior)