Error including categorical covariates in the model

[agila5]

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

^{Created on 2020-05-07 by the reprex package (v0.3.0)}

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.

[agila5]

Actually the solution is really simple, I just had to transform the list structure into a data.frame 😅

[becarioprecario]

Hi,

Actually the solution is really simple, I just had to transform the list structure into a data.frame 😅

Data should be a list because you can vectors of variables and the matrix with the categorical covariate. I am not sure why it doesn’t work. In the code of the paper we include the intercepts as categorical covariates.

Best,

Virgilio