# NOT RUN {
# example of Ybarra and Lohr (2008)
m <- 50
X <- rnorm(m, mean=5, sd=3) # true covariate values
v <- rnorm(m, sd=2)
theta <- 1 + 3*X + v # true values
psi <- rgamma(m, shape=4.5, scale=2)
e <- rnorm(m, sd=sqrt(psi)) # sampling error
y <- theta + e # direct estimates
C <- c(rep(3, 10), rep(0, 40)) # measurement error for first 10 values
W <- X + rnorm(m, sd=sqrt(C)) # covariate subject to measurement error
# fit Ybarra-Lohr model
sampler <- create_sampler(
y ~ 1 + mec(~ 0 + W, V=C) + gen(factor=~local_),
Q0=1/psi, sigma.fixed=TRUE, linpred="fitted"
)
sim <- MCMCsim(sampler, n.iter=800, n.chain=2, store.all=TRUE, verbose=FALSE)
(summ <- summary(sim))
plot(X, W, xlab="true X", ylab="inferred X")
points(X, summ$mec2_X[, "Mean"], col="green")
abline(0, 1, col="red")
legend("topleft", legend=c("prior mean", "posterior mean"), col=c("black", "green"), pch=c(1,1))
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab