I've run into an issue where R INLA isn't computing the fitted marginal values. I first had it with my own dataset, and have been able to reproduce it following an example from this book. I suspect there must be some configuration I need to change, or maybe INLA isn't working well with something under the hood? Anyways here is the code:
library("rgdal")
boston.tr <- readOGR(system.file("shapes/boston_tracts.shp",
package="spData")[1])
#create adjacency matrices
boston.adj <- poly2nb(boston.tr)
W.boston <- nb2mat(boston.adj, style = "B")
W.boston.rs <- nb2mat(boston.adj, style = "W")
boston.tr$CMEDV2 <- boston.tr$CMEDV
boston.tr$CMEDV2 [boston.tr$CMEDV2 == 50.0] <- NA
#define formula
boston.form <- log(CMEDV2) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) +
AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT)
boston.tr$ID <- 1:length(boston.tr)
#run model
boston.iid <- inla(update(boston.form, . ~. + f(ID, model = "iid")),
data = as.data.frame(boston.tr),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
control.predictor = list(compute = TRUE)
)
When I look at the output of this model, it species that the fitted values were computed:
summary(boston.iid)
Call:
c("inla(formula = update(boston.form, . ~ . + f(ID, model = \"iid\")), ", " data = as.data.frame(boston.tr),
control.compute = list(dic = TRUE, ", " waic = TRUE, cpo = TRUE), control.predictor = list(compute = TRUE))"
)
Time used:
Pre = 0.981, Running = 0.481, Post = 0.0337, Total = 1.5
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 4.376 0.151 4.080 4.376 4.672 4.376 0
CRIM -0.011 0.001 -0.013 -0.011 -0.009 -0.011 0
ZN 0.000 0.000 -0.001 0.000 0.001 0.000 0
INDUS 0.001 0.002 -0.003 0.001 0.006 0.001 0
CHAS1 0.056 0.034 -0.010 0.056 0.123 0.056 0
I(NOX^2) -0.540 0.107 -0.751 -0.540 -0.329 -0.540 0
I(RM^2) 0.007 0.001 0.005 0.007 0.010 0.007 0
AGE 0.000 0.001 -0.001 0.000 0.001 0.000 0
log(DIS) -0.143 0.032 -0.206 -0.143 -0.080 -0.143 0
log(RAD) 0.082 0.018 0.047 0.082 0.118 0.082 0
TAX 0.000 0.000 -0.001 0.000 0.000 0.000 0
PTRATIO -0.031 0.005 -0.040 -0.031 -0.021 -0.031 0
B 0.000 0.000 0.000 0.000 0.001 0.000 0
log(LSTAT) -0.329 0.027 -0.382 -0.329 -0.277 -0.329 0
Random effects:
Name Model
ID IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 169.24 46.04 99.07 160.46 299.72 141.30
Precision for ID 42.84 3.40 35.40 43.02 49.58 43.80
Deviance Information Criterion (DIC) ...............: -996.85
Deviance Information Criterion (DIC, saturated) ....: 1948.94
Effective number of parameters .....................: 202.49
Watanabe-Akaike information criterion (WAIC) ...: -759.57
Effective number of parameters .................: 337.73
Marginal log-Likelihood: 39.74
CPO and PIT are computed
Posterior marginals for the linear predictor and
the fitted values are computed
However, when I try to inspect those fitted marginal values, there is nothing there:
> boston.iid$marginals.fitted.values
NULL
Interestingly enough, I do get a summary of the posteriors, so they must be getting computed somehow?
> boston.iid$summary.fitted.values
mean sd 0.025quant 0.5quant 0.975quant mode
fitted.Predictor.001 2.834677 0.07604927 2.655321 2.844934 2.959994 2.858717
fitted.Predictor.002 3.020424 0.08220780 2.824525 3.034319 3.149766 3.052558
fitted.Predictor.003 3.053759 0.08883760 2.841738 3.071530 3.188051 3.094010
fitted.Predictor.004 3.032981 0.09846662 2.801099 3.056692 3.175215 3.084842
Any ideas on what I'm mis-specifying in the call. I have set compute = T which is what I had seen causing issues on the R-INLA forums.
The developers intentionally disabled computing the marginals to make the model faster.
To enable it, you can add these to the inla arguments:
control.predictor=list(compute=TRUE)
control.compute=list(return.marginals.predictor=TRUE)
So it looks something like this:
boston.form <- log(CMEDV2) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) +
AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT)
boston.tr$ID <- 1:length(boston.tr)
#run model
boston.iid <- inla(update(boston.form, . ~. + f(ID, model = "iid")),
data = as.data.frame(boston.tr),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, return.marginals.predictor=TRUE),
control.predictor = list(compute = TRUE)
)
boston.iid$summary.fitted.values
boston.iid$marginals.fitted.values