micar: Multiple imputitaon with INLA using a CAR Spatial Model

In becarioprecario/MIINLA: Multiple Imputation with INLA

Multiple imputitaon with INLA using a CAR Spatial Model

Description

Multiple imputation using a spatial CAR regression model. The response may have missing values to be imputed, but the covariates must be fully observed.

Usage

inla.rgeneric.micar.model(
  cmd = c("graph", "Q", "mu", "initial", "log.norm.const", "log.prior", "quit"),
  theta = NULL
)

Arguments

cmd	Arguments used by latent effects defined using the 'rgeneric' latent effect.
theta	Vector of hyperparameters.

Details

This function used is used to define a latent effect that is a linear term on some covariates with missing observations. However, multiple imputation is performed on the missing values of the covariates internally using another linear model. For this reason, this package requires the following arguments when defining the latent effect:

• x Vector of covariates (with missing observations).

• W SCALED (i.e., divided by its maximum eigenvalue) adjacency SPARSE matrix for spatial effect.

• n Number of observations.

• idx.na Index with the positions of the missing observations.

This model is defined using the 'f()' function and an index of NA's is set so that imputation is done but the covariate not included in the actual model. Then, this latent effect is 'copied', which makes the covariates (observed and imputed values) into a linear term in the liner predictor. See the example.

Value

This is used internally by the 'INLA::inla()'.

Examples

library(spdep)
library(rgdal)
library(sp)
library(INLA)

# Load data
nc.sids <- readOGR(system.file("shapes/sids.shp", package="spData")[1])
proj4string(nc.sids) <- CRS("+proj=longlat +ellps=clrk66")
# Compute covariate and expected counts
nc.sids$NWPROP74 <- (nc.sids$NWBIR74 / nc.sids$BIR74)
nc.sids$EXP74 <- nc.sids$BIR74 * sum(nc.sids$SID74) / sum(nc.sids$BIR74)

# Create covariate with missing observations
nc.sids$NWPROP74M <- nc.sids$NWPROP74
idx.na <- sample(1:100, 20) #1:10 #seq(1, 100, by = 10)
nc.sids$NWPROP74M [idx.na] <- NA

# Standard model
m0 <- inla(SID74 ~ NWPROP74, family = "poisson", data = as.data.frame(nc.sids),
 E = EXP74)
summary(m0)

# Imputation model
adj <- poly2nb(nc.sids)
W <- as(nb2mat(adj, style = "B"), "Matrix")
W.scaled <- W / max(eigen(W)$values)

nc.sids$idx <- 1:nrow(nc.sids)
r.imp <- inla(NWPROP74M ~ 0 + f(idx, model = "generic1", Cmatrix = W.scaled),
  data = as.data.frame(nc.sids),
  control.predictor = list(compute = TRUE))

plot(r.imp$summary.random$idx[idx.na, "mean"], nc.sids$NWPROP74[idx.na])
abline(0, 1)


model = inla.rgeneric.define(inla.rgeneric.micar.model, debug = TRUE,
 n = nrow(nc.sids),
 x = nc.sids$NWPROP74M,
 idx.na = which(is.na(nc.sids$NWPROP74M)),
 W = W.scaled)

nc.sids$idxNA <- NA
formula = SID74 ~ 1 + f(idxNA, model = model) +
  f(idx, copy = "idxNA", fixed = FALSE,
   hyper = list(beta = list(prior = "normal", param = c(0, 0.001))))

r = inla(formula, data = as.data.frame(nc.sids),
 family = "poisson", E = EXP74,
 verbose = TRUE)
summary(r)

r.imp$summary.fitted.values[idx.na, "mean"]
r$summary.random$idx[idx.na, "mean"]
nc.sids$NWPROP74[idx.na]

becarioprecario/MIINLA documentation built on June 6, 2023, 12:45 a.m.