Dear Prof. G贸mez-Rubio, I have yet another question 馃槄 I checked the code included in the paper on how to include a covariate in the model estimating a coefficient for each level of the random effect, and I think that the idea is to create a matrix all NA and fill the matrix in the appropriate rows and columns with the values of the covariate.
I tried to apply the same idea to a character/factor covariate and it doesn't work:
# packages
options(width = 120)
suppressPackageStartupMessages({
library(INLA)
library(INLAMSM)
library(rgdal)
library(spdep)
})
#Load SIDS data
nc.sids <- readOGR(system.file("shapes/sids.shp", package="spData")[1])
#> OGR data source with driver: ESRI Shapefile
#> Source: "C:\Users\Utente\Documents\R\win-library\3.6\spData\shapes\sids.shp", layer: "sids"
#> with 100 features
#> It has 22 fields
proj4string(nc.sids) <- CRS("+proj=longlat +ellps=clrk66")
#Compute adjacency matrix, as nb object 'adj' and sparse matrix 'W'
adj <- poly2nb(nc.sids)
W <- as(nb2mat(adj, style = "B"), "Matrix")
# Compute expected cases
r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
nc.sids$EXP74 <- r74 * nc.sids$BIR74
r79 <- sum(nc.sids$SID79) / sum(nc.sids$BIR79)
nc.sids$EXP79 <- r79 * nc.sids$BIR79
# Format data
d <- list(
OBS = c(nc.sids$SID74, nc.sids$SID79),
intercept = rep(c("1974", "1979"), each = nrow(W)),
EXP = c(nc.sids$EXP74, nc.sids$EXP79)
)
# County-period index
d$idx <- 1:length(d$OBS)
# Add covariate as a character vector
nc.sids$NWPROP74 <- nc.sids$NWBIR74 / nc.sids$BIR74
nc.sids$NWPROP79 <- nc.sids$NWBIR79 / nc.sids$BIR79
n <- nrow(W)
NWPROP_cut <- matrix(NA_character_, ncol = 2, nrow = 2 * n)
NWPROP_cut[1:n, 1] <- as.character(cut(nc.sids$NWPROP74, 3))
NWPROP_cut[n + 1:n, 2] <- as.character(cut(nc.sids$NWPROP79, 3))
head(NWPROP_cut)
#> [,1] [,2]
#> [1,] "(0.000719,0.259]" NA
#> [2,] "(0.000719,0.259]" NA
#> [3,] "(0.000719,0.259]" NA
#> [4,] "(0.000719,0.259]" NA
#> [5,] "(0.516,0.773]" NA
#> [6,] "(0.516,0.773]" NA
tail(NWPROP_cut)
#> [,1] [,2]
#> [195,] NA "(0.00248,0.255]"
#> [196,] NA "(0.255,0.507]"
#> [197,] NA "(0.255,0.507]"
#> [198,] NA "(0.255,0.507]"
#> [199,] NA "(0.255,0.507]"
#> [200,] NA "(0.255,0.507]"
d$NWPROP_cut <- NWPROP_cut
# Define the independent IMCAR effect
k <- 2
model <- inla.INDIMCAR.model(k = k, W = W)
A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
e <- rep(0, k)
# Fit the model without that covariate
mod <- inla(
formula = OBS ~ 0 + intercept +
f(idx, model = model, extraconstr = list(A = as.matrix(A), e = e)),
family = "poisson",
data = d,
E = EXP
)
#> Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'rgeneric' in section 'latent' is marked as 'experimental'; changes may appear at any time.
#> Use this model with extra care!!! Further warnings are disabled.
# Add the covariate
mod <- inla(
formula = OBS ~ 0 + intercept + NWPROP_cut +
f(idx, model = model, extraconstr = list(A = as.matrix(A), e = e)),
family = "poisson",
data = d,
E = EXP
)
#> Error in `[[<-.data.frame`(`*tmp*`, i, value = structure(c(1L, 1L, 1L, : replacement has 400 rows, data has 200
Did I miss something here or do I have to define the character matrix in some other way?
Thank you very much for your help.