Simple example with missingness and two types of measurement error
| Error types | Likelihood | Response | Covariate with error | Other covariate(s) |
|---|---|---|---|---|
| Berkson, classical, missing values | Gaussian | \(y\) | \(x\) | \(z\) |
This is a simple simulation with Berkson and classical error as well as missing data, to check that the package works as expected in that scenario.
Generating the data
set.seed(2024)
n <- 1000
# Covariate without error:
z <- rnorm(n, mean = 0, sd = 1)
# Berkson error:
u_b <- rnorm(n, sd = 1)
alpha.0 <- 1; alpha.z <- 2
r <- rnorm(n, mean = alpha.0 + alpha.z*z, sd = 1)
x <- r + u_b # Turn off Berkson by commenting out "+ u_b"
# Response:
beta.0 <- 1; beta.x <- 2; beta.z <- 2
y <- beta.0 + beta.x*x + beta.z*z + rnorm(n)
# Classical error:
u_c <- rnorm(n, sd = 1)
x_obs <- r + u_c
# Missingness:
m_pred <- -1.5 - 0.5*z # This gives a mean probability of missing of ca 0.2.
m_prob <- exp(m_pred)/(1 + exp(m_pred))
m_index <- as.logical(rbinom(n, 1, prob = m_prob)) # MAR
# m_index <- sample(1:n, 0.2*n, replace = FALSE) # MCAR
x_obs[m_index] <- NA
simple_data <- data.frame(y = y, x = x_obs, z = z)Fitting the model
# Fit the model
simple_model <- fit_inlamemi(data = simple_data,
formula_moi = y ~ x + z,
formula_imp = x ~ z,
family_moi = "gaussian",
error_type = c("berkson", "classical"),
prior.prec.moi = c(10, 9), # Gamma(10, 9)
prior.prec.berkson = c(10, 9), # Gamma(10, 9)
prior.prec.classical = c(10, 9), # Gamma(10, 9)
prior.prec.imp = c(10, 9), # Gamma(10, 9)
prior.beta.error = c(0, 1/1000), # N(0, 10^3)
initial.prec.moi = 1,
initial.prec.berkson = 1,
initial.prec.classical = 1,
initial.prec.imp = 1)
summary(simple_model)
#> Formula for model of interest:
#> y ~ x + z
#>
#> Formula for imputation model:
#> x ~ z
#>
#> Error types:
#> [1] "berkson" "classical"
#>
#> Fixed effects for model of interest:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> beta.0 1.033298 0.2179847 0.6149367 1.033448 1.439182 1.028936
#> beta.z 1.918095 0.3864149 1.2320486 1.911048 2.564599 1.916255
#>
#> Coefficient for variable with measurement error and/or missingness:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> beta.x 1.973451 0.2021624 1.572865 1.974339 2.36885 1.978043
#>
#> Fixed effects for imputation model:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> alpha.x.0 1.033078 0.05059793 0.9338142 1.033086 1.132297 1.033086
#> alpha.x.z 2.024712 0.05226033 1.9222594 2.024694 2.127263 2.024694
#>
#> Model hyperparameters (apart from beta.x):
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> Precision for model of interest 1.1282039 0.3610639 0.5673544 1.0788282 1.973095 0.9880770
#> Precision for x berkson model 1.1283327 0.3464080 0.5947152 1.0795115 1.944788 0.9885253
#> Precision for x classical model 0.9256935 0.1092025 0.7301264 0.9190625 1.159482 0.9054505
#> Precision for x imp model 0.9777885 0.1259149 0.7530336 0.9699057 1.247940 0.9545509simple.truth <- tibble::tribble(
~"variable", ~"value",
"beta.x", beta.x,
"beta.z", beta.z,
# "beta.0", beta.0,
"alpha.x.z", alpha.z,
# "alpha.0", alpha.0
)
plot(simple_model, plot_intercepts = FALSE) +
geom_point(data = simple.truth, aes(x = value))
plot of chunk unnamed-chunk-4
