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.

Dear Prof. Gómez-Rubio, I found a minor error running the code in the JSS paper. I installed the INLAMSM package from github while the R code is the same as in the document that you sent me.

#Load SIDS data
library("rgdal")
nc.sids <- readOGR(system.file("shapes/sids.shp", package = "spData")[1],
                   verbose = FALSE)
proj4string(nc.sids) <- CRS("+proj=longlat +ellps=clrk66")

#Compute adjacency matrix, as nb object 'adj' and sparse matrix 'W'
library("spdep")
adj <- poly2nb(nc.sids)
W <- as(nb2mat(adj, style = "B"), "Matrix")

# First time period
r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
nc.sids$EXP74 <- r74 * nc.sids$BIR74
nc.sids$SMR74 <- nc.sids$SID74 / nc.sids$EXP74
nc.sids$NWPROP74 <- nc.sids$NWBIR74 / nc.sids$BIR74

# Second time period
r79 <- sum(nc.sids$SID79) / sum(nc.sids$BIR79)
nc.sids$EXP79 <- r79 * nc.sids$BIR79
nc.sids$SMR79 <- nc.sids$SID79 / nc.sids$EXP79
nc.sids$NWPROP79 <- nc.sids$NWBIR79 / nc.sids$BIR79

d <- data.frame(OBS = c(nc.sids$SID74, nc.sids$SID79),
                PERIOD = c(rep("74", 100), rep("79", 100)), 
                NWPROP = c(nc.sids$NWPROP74, nc.sids$NWPROP79),
                EXP = c(nc.sids$EXP74, nc.sids$EXP79))
# County-period index
d$idx <- 1:length(d$OBS)

library("INLAMSM")
#> Loading required package: Matrix
#> Loading required package: MCMCpack
#> Loading required package: coda
#> Loading required package: MASS
#> ##
#> ## Markov Chain Monte Carlo Package (MCMCpack)
#> ## Copyright (C) 2003-2020 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
#> ##
#> ## Support provided by the U.S. National Science Foundation
#> ## (Grants SES-0350646 and SES-0350613)
#> ##
library("INLA")
#> Loading required package: parallel
#> Loading required package: foreach
#> This is INLA_20.03.09 built 2020-03-09 09:12:35 UTC.
#> See www.r-inla.org/contact-us for how to get help.
packageVersion("INLAMSM")
#> [1] '0.2'

# Number of variables (i.e., periods)
k <- 2
# Define bivariate latent effect
model.indimcar <- inla.rgeneric.define(inla.rgeneric.indep.IMCAR.model,
                                       list(k = k, W = W))
A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
e  = rep(0, k)

IIMCAR <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.indimcar,
                                    extraconstr = list(A = as.matrix(A), e = e)),
               data = d,
               E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
               control.compute = list(dic = TRUE, waic = TRUE))
#> 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.
#> Error in inla.inlaprogram.has.crashed(): The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
#>   If this does not help, please contact the developers at <[email protected]>.

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

The same code runs without any issue if I specify model.indimcar <- inla.INDIMCAR.model(k = k, W = W) instead of the inla.rgeneric.define code.

obj$marginals.hyperpar[[1]] at Line 75 should be obj$marginals.hyperpar[[1 + spat.cor.param]]? Otherwise, it will refer to the gaussian observation precision for INDPMCAR and PMCAR.

https://github.com/becarioprecario/INLAMSM/blob/706d8340e973407b7e0eb18d582f59ef9b7b0a65/R/transform.R#L70-75

I'm running the example script from your repository on Windows 10
`R version 4.1.1 (2021-08-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] grid parallel stats graphics grDevices utils datasets
[8] methods base

other attached packages:
[1] rgdal_1.5-23 dplyr_1.0.7 INLAMSM_0.2-3 MCMCpack_1.5-0
[5] MASS_7.3-54 coda_0.19-4 survey_4.1-1 survival_3.2-11
[9] INLA_21.02.23 foreach_1.5.1 Matrix_1.3-4 tigris_1.4.1
[13] spdep_1.1-8 spData_0.3.10 sp_1.4-5 haven_2.4.3
[17] mapview_2.10.0 sf_1.0-2

loaded via a namespace (and not attached):
[1] mcmc_0.9-7 nlme_3.1-152
[3] matrixStats_0.60.1 satellite_1.0.2
[5] webshot_0.5.2 gmodels_2.18.1
[7] httr_1.4.2 tools_4.1.1
[9] utf8_1.2.2 R6_2.5.1
[11] KernSmooth_2.23-20 DBI_1.1.1
[13] colorspace_2.0-2 raster_3.4-13
[15] tidyselect_1.1.1 leaflet_2.0.4.1
[17] compiler_4.1.1 quantreg_5.86
[19] cli_3.0.1 leafem_0.1.6
[21] SparseM_1.81 expm_0.999-6
[23] scales_1.1.1 classInt_0.4-3
[25] proxy_0.4-26 rappdirs_0.3.3
[27] systemfonts_1.0.2 stringr_1.4.0
[29] digest_0.6.27 foreign_0.8-81
[31] rmarkdown_2.10 svglite_2.0.0
[33] base64enc_0.1-3 pkgconfig_2.0.3
[35] htmltools_0.5.2 fastmap_1.1.0
[37] htmlwidgets_1.5.3 rlang_0.4.11
[39] rstudioapi_0.13 farver_2.1.0
[41] generics_0.1.0 jsonlite_1.7.2
[43] crosstalk_1.1.1 gtools_3.9.2
[45] magrittr_2.0.1 s2_1.0.6
[47] Rcpp_1.0.7 munsell_0.5.0
[49] fansi_0.5.0 lifecycle_1.0.0
[51] stringi_1.7.4 yaml_2.2.1
[53] maptools_1.1-1 gdata_2.18.0
[55] forcats_0.5.1 crayon_1.4.1
[57] deldir_0.2-10 lattice_0.20-44
[59] splines_4.1.1 hms_1.1.0
[61] leafpop_0.1.0 knitr_1.33
[63] pillar_1.6.2 uuid_0.1-4
[65] boot_1.3-28 codetools_0.2-18
[67] stats4_4.1.1 LearnBayes_2.15.1
[69] wk_0.5.0 glue_1.4.2
[71] evaluate_0.14 mitools_2.4
[73] leaflet.providers_1.9.0 png_0.1-7
[75] vctrs_0.3.8 MatrixModels_0.5-0
[77] purrr_0.3.4 assertthat_0.2.1
[79] xfun_0.25 e1071_1.7-8
[81] class_7.3-19 conquer_1.0.2
[83] tibble_3.1.4 iterators_1.0.13
[85] units_0.7-2 ellipsis_0.3.2
[87] brew_1.0-6 `

When I get to run the model:
r <- inla(OBS ~ 1 + f(idx, model = model), data = d, E = EXP, family = "poisson", control.compute = list(config = TRUE), control.predictor = list(compute = TRUE))

Rstudio throws the error
Fatal error: Unable to initialize the JIT

Any suggestions on a fix?

Dear all, thank you very much for working on this R package. It helped me a lot.

I created this github issue since I want to implement a (slightly different) multivariate CAR effect than what you defined in the package and, while I was studying your code, I faced some doubts/issues so I hope you can help me. I created a github issue instead of an email since here it's easier to show my doubts with the code, sorry if that's a problem.

More precisely I don't understand the definition of the log.prior function in the MCAR and IMCAR generic cases:

I read from here and here that the priors should be specified on theta (the internal representation of the parameters) and, if I define a prior using the natural scale of the parameters, then the internal prior on theta is obtained using the transfomation of r.v. theorem. For example, when you defined a uniform prior on the autocorrelation parameter, you defined that the prior on theta[1] (the internal representation of the autocorrelation parameter) is given by the derivative of the inverse of the transformation (i.e. -theta[1] - 2 * log(1 + exp(-theta[1]))) times a constant (i.e. the pdf of the uniform r.v.). I don't understand how to generalize and apply these ideas to the multivariate content (i.e. the Wishart distribution). Could you please explain the code from lines 307 to 313?

Moreover, I think I found a small typo here:

since when you define the internal prior on the log-precisions I think you should consider theta[-1] or theta[2:(k + 1)], i.e. exclude theta[1], which is the autocorrelation parameter, from theta.

Dear Prof. Gómez-Rubio, this is again a minor issue but I was thinking that maybe it would be useful implementing a new function for transforming the hyperparameters of an INDIMCAR model into the original scale (i.e. the log precisions into variances).

I know that the transformations are really easy (just function(x) exp(-x)) but I think that the function could still be useful for an easy comparison of the variances of the four types of model (INDIMCAR, IMCAR, PCAR and MMODEL).