# Errors-in-variables regression
library(INLA); library(brinla)
# simulation a dataset with heteroscedatic measurment error
set.seed(5)
n = 100
beta0 = 1
beta1 = 5
prec.y = 1
prec.u = 1
prec.x = 1
## generate the true unobserved predictor x
x <- rnorm(n, sd = 1/sqrt(prec.x))
## generate the observed predictor w with heteroscedastic error
d <- runif(n, min = 0.5, max = 1.5)
w <- x + rnorm(n, sd = 1/sqrt(prec.u * d))
## generate the response from the true model with the unobserved x
y <- beta0 + beta1*x + rnorm(n, sd = 1/sqrt(prec.y))
sim.data <- data.frame(y, w, d)
## naive method using INLA
sim.inla <- inla(y ~ w, data = sim.data, family = "gaussian")
# summary(sim.inla)
round(sim.inla$summary.fixed, 4)
## set initial values of hyperparameters
init.prec.u <- prec.u
init.prec.x <- var(w) - 1/prec.u
init.prec.y <- sim.inla$summary.hyperpar$mean
## set prior parameters
prior.beta = c(0, 0.0001)
prior.prec.u = c(10, 10)
prior.prec.x = c(10, 10)
prior.prec.y = c(10, 10)
# define the mec model formula
formula = y ~ f(w, model = "mec", scale = d, values = w,
hyper = list(
beta = list(prior = "gaussian", param = prior.beta, fixed = FALSE),
prec.u = list(prior = "loggamma", param = prior.prec.u, initial = log(init.prec.u), fixed = FALSE),
mean.x = list(prior = "gaussian", initial = 0, fixed = TRUE),
prec.x = list(prior = "loggamma", param = prior.prec.x, initial = log(init.prec.x), fixed = FALSE)))
# fit the mec model
sim.mec.inla <- inla(formula, data = sim.data, family = "gaussian",
control.family = list(hyper = list(prec = list(param = prior.prec.y, initial = log(init.prec.y), fixed=FALSE))),
control.fixed = list(mean.intercept = prior.beta[1], prec.intercept = prior.beta[2]))
# summary(sim.mec.inla)
round(sim.mec.inla$summary.fixed, 4)
round(sim.mec.inla$summary.hyperpar, 4)
# Plot the results
plot(w, y, xlim=c(-4, 4), col= "grey")
curve(1 + 5*x, -4,4, add=TRUE, lwd=2)
curve(sim.mec.inla$summary.fixed$mean + sim.mec.inla$summary.hyperpar$mean[2]*x, -4,4, add=TRUE, lwd=2, lty = 5)
curve(sim.inla$summary.fixed$mean[1] + sim.inla$summary.fixed$mean[2]*x, -4,4, add=TRUE, lwd=2, lty = 4)
# ---------------------------
# framingham study from Carroll
# Read the data
data(framingham, package = "brinla")
fram <- subset(framingham, select = c("AGE", "SMOKE", "FIRSTCHD"))
fram$W1 <- log((framingham$SBP21 + framingham$SBP22)/2 - 50)
fram$W2 <- log((framingham$SBP31 + framingham$SBP32)/2 - 50)
fram$W <- (fram$W1 + fram$W2)/2
# Set priors
prior.beta <- c(0, 0.01)
prior.lambda <- c(0, 0.01)
# estimate initial values of precision parameters
prec.u <- 1/(var(fram$W1 - fram$W2)/2)
prec.x <- 1/(var((fram$W1 + fram$W2)/2) - (1/4)*(var((fram$W1 - fram$W2)/2)))
prior.prec.x <- c(prec.x, 1)
prior.prec.u <- c(prec.u, 1)
# Create the data frame for the joint modeling
n <- nrow(fram)
Y <- matrix(NA, 4*n, 3)
Y[1:n, 1] <- fram$FIRSTCHD
Y[n+(1:n), 2] <- rep(0, n)
Y[2*n+(1:n), 3] <- fram$W1
Y[3*n+(1:n), 3] <- fram$W2
beta.0 <- c(rep(1, n), rep(NA, n), rep(NA, n), rep(NA, n))
beta.x <- c(1:n, rep(NA, n), rep(NA, n), rep(NA, n))
idx.x <- c(rep(NA, n), 1:n, 1:n, 1:n)
weight.x <- c(rep(1, n), rep(-1, n), rep(1, n), rep(1,n))
beta.smoke <- c(fram$SMOKE, rep(NA, n), rep(NA, n), rep(NA,n))
beta.age <- c(fram$AGE, rep(NA, n), rep(NA, n), rep(NA,n))
lambda.0 <- c(rep(NA, n), rep(1, n), rep(NA, n), rep(NA, n))
lambda.smoke <- c(rep(NA, n), fram$SMOKE, rep(NA, n), rep(NA, n))
lambda.age <- c(rep(NA, n), fram$AGE, rep(NA, n), rep(NA, n))
Ntrials <- c(rep(1, n), rep(NA, n), rep(NA, n), rep(NA, n))
fram.jointdata <- data.frame(Y=Y,
beta.0=beta.0, beta.x=beta.x, beta.smoke=beta.smoke, beta.age=beta.age,
idx.x=idx.x, weight.x=weight.x,
lambda.0=lambda.0, lambda.smoke=lambda.smoke, lambda.age=lambda.age,
Ntrials=Ntrials)
fram.formula <- Y ~ f(beta.x, copy = "idx.x", hyper = list(beta = list(param = prior.beta, fixed = FALSE))) + f(idx.x, weight.x, model = "iid", values = 1:n, hyper = list(prec = list(initial = -15, fixed = TRUE))) + beta.0 - 1 + beta.smoke + beta.age + lambda.0 + lambda.smoke + lambda.age
fram.mec.inla <- inla(fram.formula, Ntrials = Ntrials, data = fram.jointdata, family = c("binomial", "gaussian", "gaussian"),
control.family = list(
list(hyper = list()),
list(hyper = list(
prec = list(initial = log(prec.x),
param = prior.prec.x,
fixed = FALSE))),
list(hyper = list(
prec = list(initial=log(prec.u),
param = prior.prec.u,
fixed = FALSE)))),
control.fixed = list(
mean = list(beta.0=prior.beta[1], beta.smoke=prior.beta[1], beta.age=prior.beta[1],
lambda.0=prior.lambda[1], lambda.smoke=prior.lambda[1], lambda.age=prior.lambda[1]),
prec = list(beta.0=prior.beta[2], beta.smoke=prior.beta[2], beta.age=prior.beta[2],
lambda.0=prior.lambda[2], lambda.smoke=prior.lambda[2], lambda.age=prior.lambda[2]))
)
fram.mec.inla <-inla.hyperpar(fram.mec.inla)
# summary(fram.mec.inla)
round(fram.mec.inla$summary.fixed, 4)
round(fram.mec.inla$summary.hyperpar, 4)
plot(fram.mec.inla)
# compare the result with the model ignoring measurement error
fram.naive.inla <- inla(FIRSTCHD~ SMOKE + AGE + W, family = "binomial", data = fram)
# summary(fram.naive.inla)
round(fram.naive.inla$summary.fixed, 4)
# ---------------------------
# Bronchitis study from Carroll: Berkson Error
# Read the data
data(bronch, package = "brinla")
bronch1 <- subset(bronch, dust >=1.28)
round(prop.table(table(bronch1$cbr)),4)
round(prop.table(table(bronch1$smoking)),4)
round(c(mean(bronch1$dust), sd(bronch1$dust)), 4)
round(cor(bronch1$dust, bronch1$expo), 4)
round(by(bronch1$dust, bronch1$smoking, mean), 4)
# hist(bronch1$dust)
prior.beta <- c(0, 0.01)
prec.u <- 1/1.3
prior.prec.u <- c(1/1.3, 0.01)
bronch.formula <- cbr ~ smoking + expo + f(dust, model="meb", hyper = list(beta = list(param = prior.beta, fixed = FALSE), prec.u = list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))
bronch.meb.inla <- inla(bronch.formula, data = bronch1, family = "binomial",
control.fixed = list(mean.intercept = prior.beta[1], prec.intercept = prior.beta[2], mean = prior.beta[1], prec = prior.beta[2]))
bronch.meb.inla <- inla.hyperpar(bronch.meb.inla)
# summary(bronch.meb.inla)
round(bronch.meb.inla$summary.fixed, 4)
round(bronch.meb.inla$summary.hyperpar, 4)
plot(bronch.meb.inla)
# compare with the model ignoring measurement error
bronch.naive.inla <- inla(cbr ~ dust + smoking + expo, data = bronch1, family = "binomial",
control.fixed = list(mean.intercept = prior.beta[1], prec.intercept = prior.beta[2], mean = prior.beta[1], prec = prior.beta[2]))
# summary(bronch.naive.inla)
round(bronch.naive.inla$summary.fixed, 4)
# Using copy feature to fit the model
# Create the data frame for the joint modeling
n <- nrow(bronch1)
Y <- matrix(NA, 2*n, 2)
Y[1:n, 1] <- bronch1$cbr
Y[n+(1:n), 2] <- -1* bronch1$dust
beta.0 <- c(rep(1, n), rep(NA, n))
beta.x <- c(1:n, rep(NA, n))
idx.x <- c(rep(NA, n), 1:n)
weight.x <- c(rep(NA, n), -rep(1, n))
beta.smoking <- c(bronch1$smoking, rep(NA, n))
beta.expo <- c(bronch1$expo, rep(NA, n))
Ntrials <- c(rep(1, n), rep(NA, n))
bronch.jointdata <- data.frame(Y=Y,
beta.0=beta.0, beta.x=beta.x, beta.smoking=beta.smoking, beta.expo=beta.expo,
idx.x=idx.x, weight.x=weight.x, Ntrials = Ntrials)
bronch.formula2 <- Y ~ beta.0 - 1 + beta.smoking + beta.expo + f(beta.x, copy = "idx.x", hyper = list(beta = list(param = prior.beta, fixed = FALSE))) + f(idx.x, weight.x, model = "iid", values = 1:n, hyper = list(prec = list(initial = -15, fixed = TRUE)))
bronch.mec.inla2 <- inla(bronch.formula2, Ntrials = Ntrials, data = bronch.jointdata, family = c("binomial", "gaussian"),
control.family = list(
list(hyper = list()),
list(hyper = list(
prec = list(initial=log(prec.u),
param = prior.prec.u,
fixed = FALSE)))),
control.fixed = list(mean.intercept = prior.beta[1], prec.intercept = prior.beta[2], mean = prior.beta[1], prec = prior.beta[2]))
bronch.mec.inla2 <- inla.hyperpar(bronch.mec.inla2)
summary(bronch.mec.inla2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.