Julian Faraway 21 September 2020

In this example, we compare the Bayesian model output with the linear model fit. Although the interpretation of these two fits is different, we expect similar numerical results provided the priors are not informative.

See the introduction for an overview. Load the packages: (you may need to install them first)

```
library(ggplot2)
library(INLA)
```

```
Loading required package: Matrix
Loading required package: sp
Loading required package: parallel
Loading required package: foreach
This is INLA_20.03.17 built 2020-09-21 11:41:39 UTC.
See www.r-inla.org/contact-us for how to get help.
To enable PARDISO sparse library; see inla.pardiso()
```

```
library(faraway)
library(gridExtra)
library(brinla)
```

Load in the Chicago Insurance dataset as analyzed in Linear Models with R:

```
data(chredlin, package="faraway")
```

We take `involact`

as the response. Make some plots of the data:

```
p1=ggplot(chredlin,aes(race,involact)) + geom_point()
p2=ggplot(chredlin,aes(fire,involact)) + geom_point()
p3=ggplot(chredlin,aes(theft,involact)) + geom_point()
p4=ggplot(chredlin,aes(age,involact)) + geom_point()
p5=ggplot(chredlin,aes(income,involact)) + geom_point()
grid.arrange(p1,p2,p3,p4,p5)
```

Fit the standard linear model for reference purposes:

```
lmod <- lm(involact ~ race + fire + theft + age + log(income), chredlin)
sumary(lmod)
```

```
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.18554 1.10025 -1.08 0.28755
race 0.00950 0.00249 3.82 0.00045
fire 0.03986 0.00877 4.55 4.8e-05
theft -0.01029 0.00282 -3.65 0.00073
age 0.00834 0.00274 3.04 0.00413
log(income) 0.34576 0.40012 0.86 0.39254
n = 47, p = 6, Residual SE = 0.335, R-Squared = 0.75
```

Run the default INLA model and compare to the `lm`

output:

```
formula <- involact ~ race + fire + theft + age + log(income)
imod <- inla(formula, family="gaussian", data=chredlin)
cbind(imod$summary.fixed[,1:2],summary(lmod)$coef[,1:2])
```

```
mean sd Estimate Std. Error
(Intercept) -1.1853891 1.0979085 -1.1855396 1.1002549
race 0.0095020 0.0024844 0.0095022 0.0024896
fire 0.0398555 0.0087480 0.0398560 0.0087661
theft -0.0102944 0.0028121 -0.0102945 0.0028179
age 0.0083354 0.0027383 0.0083356 0.0027440
log(income) 0.3457063 0.3992695 0.3457615 0.4001234
```

The first two columns are the posterior means and SDs from the INLA fit while the third and fourth are the corresponding values from the lm fit. We can see they are almost identical. We expect this.

We can get the posterior for the SD of the error as:

```
bri.hyperpar.summary(imod)
```

```
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.33244 0.036546 0.27005 0.32914 0.41338 0.32274
```

which compares with the linear model estimate of:

```
summary(lmod)$sigma
```

```
[1] 0.33453
```

which gives a very similar result. But notice that the Bayesian output gives us more than just a point estimate.

The function `inla.hyperpar()`

is advertized to improve the estimates of
the hyperparameters (just sigma here) but makes very little difference
in this example:

```
bri.hyperpar.summary(inla.hyperpar(imod))
```

```
mean sd q0.025 q0.5 q0.975 mode
SD for the Gaussian observations 0.33245 0.036554 0.27005 0.32914 0.4134 0.32274
```

For the error SD

```
bri.hyperpar.plot(imod)
```

For the fixed parameters:

```
bri.fixed.plot(imod)
```

We can use these plots to judge which parameters may be different from zero. For age, fire, race and theft, we see the posterior density is well away from zero but for log(income) there is an overlap. This can be qualitatively compared to the p-values from the linear model output. The conclusions are compatible although the motivation and justifactions are different.

We need to set the control.predictor to compute the posterior means of the linear predictors:

```
result <- inla(formula, family="gaussian", control.predictor=list(compute=TRUE),data=chredlin)
ypostmean <- result$summary.linear.predictor
```

Compare these posterior means to the lm() fitted values:

```
plot(ypostmean[,1],lmod$fit, xlab="INLA posterior fitted means", ylab="LM fitted values")
abline(0,1)
```

Can see these are the same. Now compare the SDs from INLA with the SEs from the LM fit:

```
lmfit <- predict(lmod, se=TRUE)
plot(ypostmean[,2], lmfit$se.fit, xlab="INLA posterior SDs", ylab="LM fitted SEs")
abline(0,1)
```

We see that these are the same.

Suppose we wish to predict the response at a new set of inputs. We add a case for the new inputs and set the response to missing:

```
newobs = data.frame(race=25,fire=10,theft=30,age=65,involact=NA,income=11,side="s")
chpred = rbind(chredlin,newobs)
```

We fit the model as before:

```
respred <- inla(formula, family="gaussian", control.predictor=list(compute=TRUE), data=chpred)
respred$summary.fixed
```

```
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -1.1853891 1.0979084 -3.3512580 -1.1854262 0.9787808 -1.1853992 2.5514e-05
race 0.0095020 0.0024844 0.0046010 0.0095019 0.0143991 0.0095020 2.5522e-05
fire 0.0398555 0.0087480 0.0225979 0.0398553 0.0570993 0.0398556 2.5526e-05
theft -0.0102944 0.0028121 -0.0158419 -0.0102945 -0.0047513 -0.0102944 2.5527e-05
age 0.0083354 0.0027383 0.0029334 0.0083353 0.0137331 0.0083354 2.5525e-05
log(income) 0.3457063 0.3992695 -0.4419491 0.3456954 1.1327297 0.3457100 2.5513e-05
```

The output is the same as before as expected. The fitted value for the last observation (the new 48th observation) provides the posterior distribution for the linear predictor corresponding to the new set of inputs.

```
respred$summary.fitted.values[48,]
```

```
mean sd 0.025quant 0.5quant 0.975quant mode
fitted.Predictor.48 0.51265 0.054281 0.4057 0.51265 0.6196 0.51265
```

We might compare this to the standard linear model prediction:

```
predict(lmod, newdata=newobs, interval = "confidence", se=TRUE)
```

```
$fit
fit lwr upr
1 0.51266 0.40292 0.62239
$se.fit
[1] 0.054337
$df
[1] 41
$residual.scale
[1] 0.33453
```

We get the essentially the same outcome. These answers only reflect variability in the linear predictor.

Usually we want to include the variation due to a new error. This requires more effort. One solution is based on simulation. We can sample from the posterior we have just calculated and sample from the posterior for the error distribution. Adding these together, we will get a sample from the posterior predictive distribution. Although this can be done quickly, it goes against the spirit of INLA by using simulation.

We can compute an approximation for the posterior predictive
distribution using some trickery. We set up an index variable `idx`

. We
add a *random effect* term using this index with the `f()`

function.
This is a mean zero Gaussian random variable whose variance is a
hyperparameter. This term functions exactly the same as the usual error
term. Since the model will also add an error term, we would be stuck
with two essentially identical terms. But we are able to knock out the
usual error term by specifying that it have a very high precision
(i.e. it will vary just a little bit). In this way, our random effect
term will take almost all the burden of representing this variation.

```
chpred$idx <- 1:48
eformula <- involact ~ race + fire + theft + age + log(income) + f(idx, model="iid", hyper = list(prec = list(param=c(1, 0.01))))
respred <- inla(eformula, data = chpred, control.predictor = list(compute=TRUE), control.family=list(initial=12,fixed=TRUE))
bri.hyperpar.summary(respred)
```

```
mean sd q0.025 q0.5 q0.975 mode
SD for idx 0.33315 0.036626 0.27063 0.32984 0.41427 0.32342
```

We see that the SD for the index variable is about the same as the previous SD for the error. Look at the predictive distribution:

```
respred$summary.fitted.values[48,]
```

```
mean sd 0.025quant 0.5quant 0.975quant mode
fitted.Predictor.48 0.51266 0.33893 -0.15516 0.51265 1.1804 0.51265
```

We can compare this to the linear model prediction interval:

```
predict(lmod,newdata=newobs,interval="prediction")
```

```
fit lwr upr
1 0.51266 -0.17179 1.1971
```

The results are similar but not the same. Notice that the lower end of the intervals are both negative which should not happen with a strictly positive response. But that’s not the fault of this process but of our choice of model.

Here is a posterior predictive density for the prediction:

```
ggplot(data.frame(respred$marginals.fitted.values[[48]]),aes(x,y))+geom_line()+ylab("Density")+xlab("involact")
```

We see the problem of having density for values less than zero.

For more on this, see the Google groups INLA mailing list discussion.

Bayesian diagnostic methods can use these statistics. We need to ask for them at the time we fit the model:

```
imod <- inla(formula, family="gaussian", data=chredlin, control.compute=list(cpo=TRUE))
```

The CPO is the probability density of an observed response based on the model fit to the rest of the data. We are interested in small values of the CPO since these represent unexpected response values. Plot the CPOs:

```
n = nrow(chredlin)
plot(1:n,imod$cpo$cpo, ylab="CPO",type="n")
text(1:n,imod$cpo$cpo, 1:n)
```

Observations 3, 6 and 35 have particularly low values. We might compare these to the jacknife residuals (an analogous quantity) from the standard linear model:

```
plot(imod$cpo$cpo, abs(rstudent(lmod)), xlab="CPO",ylab="Absolute residuals", type="n")
text(imod$cpo$cpo, abs(rstudent(lmod)), 1:n)
```

We can see the two statistics are related. They pick out the same three observations.

The PIT is the probability of a new response less than the observed response using a model based on the rest of the data. We’d expect the PIT values to be uniformly distributed if the model assumptions are correct.

```
pit <- imod$cpo$pit
uniquant <- (1:n)/(n+1)
plot(uniquant, sort(pit), xlab="uniform quantiles", ylab="Sorted PIT values")
abline(0,1)
```

Looks good - no indication of problems. Might be better to plot these on a logit scale because probabilities very close to 0 or 1 won’t be obvious from the above plot.

```
plot(logit(uniquant), logit(sort(pit)), xlab="uniform quantiles", ylab="Sorted PIT values", main="Logit scale")
abline(0,1)
```

This marks out three points as being somewhat unusual:

```
which(abs(logit(pit)) > 4)
```

```
[1] 3 6 35
```

The same three points picked out by the CPO statistic.

We can compare the PIT values and posterior means for the fitted values (which are the same as the fitted values from the linear model fit).

```
plot(predict(lmod), pit, xlab="Posterior fitted means", ylab="PIT")
```

This indicates some problems for low values. Since the response is strictly positive and this is not reflected in our models, this does spot a real flaw in our model.

Incidentally, a similar conclusion can be made from the standard linear model diagnostic:

```
plot(lmod,1)
```

which looks fairly similar. In this example, we might try something other than a Gaussian response which reflects the zero response values.

We can also plot the PITs against the predictors - for example:

```
plot(chredlin$race, pit, xlab="Minority percentage", ylab="PIT")
```

Which shows no problems with the fit.

We need to ask to get the Watanabe AIC (WAIC) and deviance information criterion DIC calculated:

```
imod <- inla(formula, family="gaussian", data=chredlin, control.compute=list(dic=TRUE, waic=TRUE))
```

which produces:

```
c(WAIC=imod$waic$waic, DIC=imod$dic$dic)
```

```
WAIC DIC
40.151 38.611
```

We do not have many predictors here so exhaustive search is feasible. We
fit all models that include race but with some combination of the other
four predictors. Use `lm()`

to get the AIC but INLA to get DIC and WAIC

```
listcombo <- unlist(sapply(0:4,function(x) combn(4, x, simplify=FALSE)),recursive=FALSE)
predterms <- lapply(listcombo, function(x) paste(c("race",c("fire","theft","age","log(income)")[x]),collapse="+"))
coefm <- matrix(NA,16,3)
for(i in 1:16){
formula <- as.formula(paste("involact ~ ",predterms[[i]]))
lmi <- lm(formula, data=chredlin)
result <- inla(formula, family="gaussian", data=chredlin, control.compute=list(dic=TRUE, waic=TRUE))
coefm[i,1] <- AIC(lmi)
coefm[i,2] <- result$dic$dic
coefm[i,3] <- result$waic$waic
}
rownames(coefm) <- predterms
colnames(coefm) <- c("AIC","DIC","WAIC")
round(coefm,4)
```

```
AIC DIC WAIC
race 62.032 62.216 62.416
race+fire 50.744 50.990 54.014
race+theft 63.925 64.171 64.073
race+age 54.026 54.273 54.304
race+log(income) 57.883 58.129 58.572
race+fire+theft 43.970 44.303 46.623
race+fire+age 47.306 47.639 50.031
race+fire+log(income) 51.412 51.744 54.445
race+theft+age 54.252 54.585 54.597
race+theft+log(income) 59.591 59.923 60.360
race+age+log(income) 54.804 55.137 55.152
race+fire+theft+age 36.876 37.321 39.117
race+fire+theft+log(income) 45.569 46.013 48.202
race+fire+age+log(income) 49.272 49.717 51.793
race+theft+age+log(income) 55.215 55.660 55.810
race+fire+theft+age+log(income) 38.028 38.611 40.151
```

The model corresponding to the minimum value of the criterion is the
same for all three. We would choose `involact ~ race+fire+theft+age`

.

We see that the standard linear modeling and the Bayesian approach using INLA produce very similar results. On one hand, one might decide to stick with the standard approach since it is easier. But on the other hand, this does give the user confidence to use the Bayesian approach with INLA because it opens up other possibilities with the model fitting and interpretation. The consistency between the approaches in the default settings is reassuring.