Example in inlabru

library(INLA)
library(inlabru)

inla.setOption(num.threads = "1:1")

We here show how to fit the measurement error and missing data model in inlabru. The inlabru code is written by Sara Martino, and commented by Emma Skarstein.

Loading the data

For this demonstration we will use the simulated data from the Simulation example. Just to refresh: this is a situation where we have one covariate (\(\boldsymbol{x}\)) that has both classical error, Berkson error and missing data. We also observe another covariate (\(\boldsymbol{z}\)) that has no error.

data <- read.csv("data/simulated_data.csv")
n <- nrow(data)

Priors

We use the exact same priors as in the original simulation study.

# Priors for model of interest coefficients
prior.beta = c(0, 1/1000) # N(0, 10^3)

# Priors for exposure model coefficients
prior.alpha <- c(0, 1/10000) # N(0, 10^4) 

# Priors for y, measurement error and true x-value precision
prior.prec.y <- c(10, 9) # Gamma(10, 9)
prior.prec.u_b <- c(10, 9) # Gamma(10, 9)
prior.prec.u_c <- c(10, 9) # Gamma(10, 9)
prior.prec.r <- c(10, 9) # Gamma(10, 9) 


# Initial values
prec.y <- 1
prec.u_b <- 1
prec.u_c <- 1
prec.r <- 1

Specifying the data

# Regression model of interest 
data_moi <- data.frame(y = data$y, z = data$z, weight.x = 1, weight.r = 1, r = 1:n) 

# Berkson ME model
data_berkson <- data.frame(zero = 0, weight.x = -1, weight.r = 1, r = 1:n) 

# Classical ME model
data_classical <- data.frame(w = data$w, weight.x = 1, weight.r = 1, r = 1:n)

# Imputation model
data_imputation <- data.frame(zero = 0, z = data$z, weight.x = 1, weight.r = -1, r = 1:n)

Formulas

cmp_new <- ~ Intercept(1, model = "linear", prec.linear = prior.beta[2]) +
  beta_z(main = z, model = "linear", prec.linear = prior.beta[2]) +
  x_eff(r, weight.x, model = "iid",  hyper = list(prec = list(initial = -15, fixed=TRUE))) +
  x_eff_copy(r, copy="x_eff", hyper = list(beta = list(param = prior.beta, fixed=FALSE))) +
  r_eff(r, weight.r, model = "iid",  hyper = list(prec = list(initial = -15, fixed=TRUE))) +
  alpha_0(main = 1, model = "linear", prec.linear = prior.alpha[2]) +
  alpha_z(main = z, model = "linear", prec.linear = prior.alpha[2])


lik_moi <- like(formula = y ~ .,
                family = "gaussian",
                include = c("Intercept", "beta_z", "x_eff_copy"),
                control.family = list(
                  hyper = list(prec = list(initial = log(prec.y), 
                                           param = prior.prec.y, 
                                           fixed = FALSE))),
            data = data_moi)

lik_berkson <- like(formula = zero ~ .,
                    family = "gaussian",
                    include = c("x_eff", "r_eff"),
                    control.family = list(
                      hyper = list(prec = list(initial = log(prec.u_b),
                                               param = prior.prec.u_b,
                                               fixed = FALSE))),
                    data = data_berkson)

lik_classical <- like(formula = w ~ .,
                      family = "gaussian",
                      include = c("r_eff"),
                      control.family =  list(
                        hyper = list(prec = list(initial = log(prec.u_c), 
                                                 param = prior.prec.u_c, 
                                                 fixed = FALSE))),
            data  = data_classical)

lik_imputation <- like(formula = zero ~ .,
                       family = "gaussian",
                       include = c("alpha_0", "alpha_z","r_eff"),
                       control.family =   list(
                         hyper = list(prec = list(initial = log(prec.r), 
                                                  param = prior.prec.r, 
                                                  fixed = FALSE))),
             data = data_imputation)

Fitting the model

set.seed(1)
bru_options_set(bru_verbose = 1)
fit <- bru(components = cmp_new,
           lik_moi,
           lik_berkson,
           lik_classical,
           lik_imputation,
           options = list(verbose = F,
                          bru_max_iter = 20,
                          inla.mode  = "experimental"))
fit$summary.fixed

               mean         sd 0.025quant  0.5quant 0.975quant mode
Intercept 0.8705529 0.21579985  0.4541298 0.8682062   1.289494   NA
beta_z    1.8088724 0.38738755  1.0837732 1.8034426   2.527867   NA
alpha_0   1.0053028 0.04731152  0.9125007 1.0053013   1.098114   NA
alpha_z   1.9876737 0.04890610  1.8917558 1.9876680   2.083624   NA
                   kld
Intercept 9.468548e-08
beta_z    1.200887e-04
alpha_0   2.479969e-13
alpha_z   1.624621e-12

fit$summary.hyperpar

                                               mean        sd 0.025quant
Precision for the Gaussian observations    1.129095 0.3352055  0.5918354
Precision for the Gaussian observations[2] 1.149144 0.3118757  0.6520526
Precision for the Gaussian observations[3] 1.097736 0.1285802  0.8670084
Precision for the Gaussian observations[4] 1.025116 0.1185145  0.8099070
Beta for x_eff_copy                        2.078488 0.1895427  1.7012042
                                           0.5quant 0.975quant mode
Precision for the Gaussian observations    1.089183   1.901142   NA
Precision for the Gaussian observations[2] 1.110314   1.873069   NA
Precision for the Gaussian observations[3] 1.089934   1.373323   NA
Precision for the Gaussian observations[4] 1.018916   1.276604   NA
Beta for x_eff_copy                        2.079849   2.448850   NA