Bayesian Regression using INLA

The INLA Approach to Bayesian models

The Integrated Nested Laplace Approximation, or INLA, approach is a recently developed, computationally simpler method for fitting Bayesian models [(Rue et al., 2009, compared to traditional Markov Chain Monte Carlo (MCMC) approaches. INLA fits models that are classified as latent Gaussian models, which are applicable in many settings (Martino & Rue, 2010. In general, INLA fits a general form of additive models such as:

\(\eta = \alpha + \sum_{j=1}^{nf} f^{(j)}(u_{ij}) + \sum_{k=1}^{n\beta}\beta_k z_{ki} + \epsilon_i\)

where \(\eta\) is the linear predictor for a generalized linear model formula, and is composed of a linear function of some variables u, \(\beta\) are the effects of covariates, z, and \(\epsilon\) is an unstructured residual (Rue et al., 2009). As this model is often parameterized as a Bayesian one, we are interested in the posterior marginal distributions of all the model parameters. Rue and Martino (2007) show that the posterior marginal for the random effects (x) in such models can be approximated as:

\(\tilde{p}(x_i|y) = \sum_k \tilde{p}(x_i|\theta_k, y) \tilde{p}(\theta_k|y) \Delta_k\)

via numerical integration (Rue & Martino, 2007; Schrodle & Held, 2011a, 2011b). The posterior distribution of the hyperparameters (\(\theta\)) of the model can also be approximated as:

\(\tilde{p}(\theta | y)) \propto \frac{p(x, \theta, y)}{\tilde{p}G(x| \theta,y)} \mid _{x} = x^*(\theta)\)

, where G is a Gaussian approximation of the posterior and \(x^*(\theta)\) is the mode of the conditional distribution of \(p(x|\theta,y)\). Thus, instead of using MCMC to find an iterative, sampling-based estimate of the posterior, it is arrived at numerically. This method of fitting the spatial models specified above has been presented by numerous authors (Blangiardo & Cameletti, 2015; Blangiardo et al., 2013; Lindgren & Rue, 2015; Martins et al., 2013; Schrodle & Held, 2011a, 2011b), with comparable results to MCMC.

Below, I show examples of using INLA to fit Bayesian regression models for data from US counties. One example will be a relatively small data set from Texas, while the other example will be all US counties. The US county example basically emulates the paper by Sparks et al 2012.

library(maptools)

## Loading required package: sp
## Checking rgeos availability: TRUE

library(spdep)

## Loading required package: Matrix

library(RColorBrewer)
library(lattice)
library(INLA)

## Loading required package: splines

#read in big US shapefile with lots of stuff in it
usdat<-readShapePoly("~/Google Drive//dem7263/data/output12711.shp",  proj4string=CRS("+proj=lcc"))

## Field name: 'SP_ID_1' changed to: 'SP_ID_1.1'

#get outcome, which are the numbers of infant deaths
usdat$infd03<-usdat$f1194803
#get denominator, which are the number of births  
usdat$births03<-usdat$f1254603
usdat$e<-usdat$births03*(sum(usdat$infd03)/sum(usdat$births03))
#subset to just have TX  
tx<-usdat[usdat$STATE==48,]

#form spatial neighbors
#In INLA, we don't need FIPS codes, we need a simple numeric index for our counties
usdat$struct<-1:nrow(usdat@data)
tx$struct<-1:nrow(tx@data)

#!!Be sure to set the row names of the neighbor list to the struct variable we just made!!
us.nb<-poly2nb(usdat, row.names=usdat$struct)
txnb<-poly2nb(tx, row.names = tx$struct)

#Save the neighbor list to an external file, INLA needs the list this way!!
nb2INLA(us.nb, file="~/Google Drive//dem7263/usa.gra")
nb2INLA(nb=txnb, file = "~/Google Drive//dem7263/txgraph.gra")

We can fit the OLS model using the bayesian framework as:

#Texas data
fit.l<-inla(I(infd03/births03)~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "gaussian", data=tx@data, control.compute = list(dic=T))
summary(fit.l)

## 
## Call:
## c("inla(formula = I(infd03/births03) ~ metro + giniz + hispz + blackz + ",  "    povz + ne_bad + isowbz, family = \"gaussian\", data = tx@data, ",  "    control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.6648          0.1372          0.1141          0.9160 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0072 0.0007     0.0058   0.0072     0.0086  0.0072   0
## metro       -0.0004 0.0009    -0.0022  -0.0004     0.0013 -0.0004   0
## giniz        0.0000 0.0005    -0.0010   0.0000     0.0010  0.0000   0
## hispz       -0.0003 0.0003    -0.0009  -0.0003     0.0003 -0.0003   0
## blackz       0.0021 0.0009     0.0004   0.0021     0.0039  0.0021   0
## povz        -0.0002 0.0006    -0.0013  -0.0002     0.0008 -0.0002   0
## ne_bad       0.0010 0.0013    -0.0015   0.0010     0.0034  0.0010   0
## isowbz       0.0009 0.0005    -0.0001   0.0009     0.0019  0.0009   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                              mean       sd 0.025quant
## Precision for the Gaussian observations 111795.21 24164.46   72738.65
##                                          0.5quant 0.975quant      mode
## Precision for the Gaussian observations 108882.54  167236.72 103102.02
## 
## Expected number of effective parameters(std dev): 165.24(9.519)
## Number of equivalent replicates : 1.531 
## 
## Deviance Information Criterion (DIC) ...: -2059.14
## Effective number of parameters .........: 163.74
## 
## Marginal log-Likelihood:  904.30 
## Posterior marginals for linear predictor and fitted values computed

#get a 95% credible interval for a regresion parameter:
inla.hpdmarginal(p = .95, fit.l$marginals.fixed$giniz)

##                      low       high
## level:0.95 -0.0009704781 0.00103873

inla.hpdmarginal(p = .95, fit.l$marginals.fixed$ne_bad)

##                     low        high
## level:0.95 -0.001519333 0.003441471

inla.hpdmarginal(p = .95, fit.l$marginals.fixed$isowbz)

##                     low        high
## level:0.95 -0.000118578 0.001907848

#US Data
fit.us<-inla(I(infd03/births03)~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "gaussian", data=usdat@data, control.compute = list(dic=T))
summary(fit.us)

## 
## Call:
## c("inla(formula = I(infd03/births03) ~ metro + giniz + hispz + blackz + ",  "    povz + ne_bad + isowbz, family = \"gaussian\", data = usdat@data, ",  "    control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5675          0.9649          0.2526          1.7851 
## 
## Fixed effects:
##                mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0066 1e-04     0.0063   0.0066     0.0068  0.0066   0
## metro        0.0007 2e-04     0.0003   0.0007     0.0011  0.0007   0
## giniz       -0.0001 1e-04    -0.0003  -0.0001     0.0001 -0.0001   0
## hispz       -0.0004 1e-04    -0.0006  -0.0004    -0.0002 -0.0004   0
## blackz       0.0016 1e-04     0.0014   0.0016     0.0018  0.0016   0
## povz         0.0010 1e-04     0.0007   0.0010     0.0013  0.0010   0
## ne_bad       0.0000 3e-04    -0.0005   0.0000     0.0006  0.0000   0
## isowbz       0.0006 1e-04     0.0004   0.0006     0.0008  0.0006   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                              mean       sd 0.025quant
## Precision for the Gaussian observations 193255.36 19374.34  158795.78
##                                          0.5quant 0.975quant      mode
## Precision for the Gaussian observations 191945.80  234876.96 189053.14
## 
## Expected number of effective parameters(std dev): 2332.88(44.31)
## Number of equivalent replicates : 1.313 
## 
## Deviance Information Criterion (DIC) ...: -26244.50
## Effective number of parameters .........: 2326.66
## 
## Marginal log-Likelihood:  12017.41 
## Posterior marginals for linear predictor and fitted values computed

inla.hpdmarginal(p = .95, fit.us$marginals.fixed$giniz)

##                      low         high
## level:0.95 -0.0003258645 0.0001277709

inla.hpdmarginal(p = .95, fit.us$marginals.fixed$ne_bad)

##                      low        high
## level:0.95 -0.0004676793 0.000555227

inla.hpdmarginal(p = .95, fit.us$marginals.fixed$isowbz)

##                     low         high
## level:0.95 0.0003972207 0.0008078396

OR, as a poisson model:

#Texas data
fit.p<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "poisson",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.p)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz, family = \"poisson\", data = tx@data, E = e + ",  "    1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5707          0.1054          0.0803          0.7565 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0200 0.0677    -0.1546  -0.0194     0.1113 -0.0183   0
## metro        0.0436 0.0830    -0.1187   0.0433     0.2072  0.0428   0
## giniz        0.0614 0.0489    -0.0343   0.0613     0.1578  0.0611   0
## hispz       -0.0048 0.0274    -0.0585  -0.0048     0.0490 -0.0049   0
## blackz       0.1316 0.0894    -0.0429   0.1313     0.3079  0.1305   0
## povz         0.0303 0.0433    -0.0550   0.0304     0.1149  0.0307   0
## ne_bad      -0.3154 0.0935    -0.5001  -0.3150    -0.1332 -0.3141   0
## isowbz       0.0039 0.0707    -0.1334   0.0034     0.1440  0.0023   0
## 
## The model has no random effects
## 
## The model has no hyperparameters
## 
## Expected number of effective parameters(std dev): 8.031(0.00)
## Number of equivalent replicates : 31.63 
## 
## Deviance Information Criterion (DIC) ...: 758.72
## Effective number of parameters .........: 8.029
## 
## Marginal log-Likelihood:  -420.33 
## Posterior marginals for linear predictor and fitted values computed

#or negative binomial
fit.nb<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "nbinomial",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.nb)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz, family = \"nbinomial\", data = tx@data, E = e + ",  "    1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5591          0.1981          0.0838          0.8410 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0011 0.0885    -0.1758  -0.0008     0.1720 -0.0002   0
## metro        0.0178 0.1044    -0.1871   0.0177     0.2227  0.0177   0
## giniz        0.0524 0.0752    -0.0953   0.0524     0.2000  0.0524   0
## hispz       -0.0194 0.0429    -0.1035  -0.0194     0.0648 -0.0195   0
## blackz       0.1411 0.1262    -0.1047   0.1404     0.3906  0.1390   0
## povz         0.0434 0.0708    -0.0965   0.0436     0.1816  0.0442   0
## ne_bad      -0.2211 0.1419    -0.4988  -0.2215     0.0585 -0.2222   0
## isowbz       0.0244 0.0961    -0.1605   0.0230     0.2170  0.0204   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                                       mean    sd
## size for the nbinomial observations (overdispersion) 15.31 3.038
##                                                      0.025quant 0.5quant
## size for the nbinomial observations (overdispersion)      10.15    15.04
##                                                      0.975quant  mode
## size for the nbinomial observations (overdispersion)      22.05 14.53
## 
## Expected number of effective parameters(std dev): 8.011(7e-04)
## Number of equivalent replicates : 31.71 
## 
## Deviance Information Criterion (DIC) ...: 797.46
## Effective number of parameters .........: 8.324
## 
## Marginal log-Likelihood:  -448.79 
## Posterior marginals for linear predictor and fitted values computed

#US data
fit.p_us<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "poisson",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.p_us)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz, family = \"poisson\", data = usdat@data, E = e + ",  "    1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5782          0.7740          0.2544          1.6067 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0348 0.0151     0.0050   0.0348     0.0643  0.0349   0
## metro        0.0649 0.0208     0.0241   0.0649     0.1058  0.0649   0
## giniz       -0.0016 0.0055    -0.0123  -0.0017     0.0092 -0.0018   0
## hispz       -0.0747 0.0075    -0.0893  -0.0747    -0.0600 -0.0747   0
## blackz       0.1389 0.0114     0.1165   0.1389     0.1614  0.1389   0
## povz         0.1448 0.0108     0.1237   0.1448     0.1659  0.1449   0
## ne_bad      -0.2301 0.0189    -0.2672  -0.2301    -0.1931 -0.2300   0
## isowbz      -0.0461 0.0126    -0.0708  -0.0461    -0.0213 -0.0461   0
## 
## The model has no random effects
## 
## The model has no hyperparameters
## 
## Expected number of effective parameters(std dev): 8.445(0.00)
## Number of equivalent replicates : 362.83 
## 
## Deviance Information Criterion (DIC) ...: 10303.27
## Effective number of parameters .........: 8.444
## 
## Marginal log-Likelihood:  -5205.15 
## Posterior marginals for linear predictor and fitted values computed

#or negative binusdat
fit.nb_us<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz, family = "nbinomial",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.nb_us)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz, family = \"nbinomial\", data = usdat@data, ",  "    E = e + 1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.5677          1.8257          0.2898          2.6833 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0333 0.0159     0.0020   0.0334     0.0644  0.0334   0
## metro        0.0601 0.0232     0.0146   0.0600     0.1055  0.0600   0
## giniz       -0.0046 0.0100    -0.0242  -0.0045     0.0149 -0.0045   0
## hispz       -0.0757 0.0101    -0.0955  -0.0757    -0.0559 -0.0757   0
## blackz       0.1178 0.0143     0.0897   0.1178     0.1460  0.1177   0
## povz         0.1574 0.0138     0.1303   0.1574     0.1844  0.1575   0
## ne_bad      -0.1969 0.0246    -0.2453  -0.1969    -0.1485 -0.1969   0
## isowbz      -0.0585 0.0164    -0.0907  -0.0585    -0.0262 -0.0586   0
## 
## The model has no random effects
## 
## Model hyperparameters:
##                                                       mean    sd
## size for the nbinomial observations (overdispersion) 52.45 4.545
##                                                      0.025quant 0.5quant
## size for the nbinomial observations (overdispersion)      44.10    52.25
##                                                      0.975quant  mode
## size for the nbinomial observations (overdispersion)      61.95 51.85
## 
## Expected number of effective parameters(std dev): 8.252(0.0051)
## Number of equivalent replicates : 371.32 
## 
## Deviance Information Criterion (DIC) ...: 10268.82
## Effective number of parameters .........: 8.635
## 
## Marginal log-Likelihood:  -5234.51 
## Posterior marginals for linear predictor and fitted values computed

Or a simple, unstructured random county intercept model:

#Texas Model
fit.pre<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "iid"), family = "poisson",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.pre)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"iid\"), family = \"poisson\", ",  "    data = tx@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.6855          0.3068          0.0943          1.0866 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0202 0.0678    -0.1549  -0.0196     0.1113 -0.0184   0
## metro        0.0434 0.0831    -0.1190   0.0431     0.2071  0.0426   0
## giniz        0.0616 0.0490    -0.0343   0.0615     0.1581  0.0612   0
## hispz       -0.0049 0.0274    -0.0587  -0.0050     0.0490 -0.0050   0
## blackz       0.1316 0.0896    -0.0432   0.1312     0.3082  0.1305   0
## povz         0.0304 0.0434    -0.0551   0.0305     0.1151  0.0308   0
## ne_bad      -0.3148 0.0936    -0.4998  -0.3143    -0.1322 -0.3135   0
## isowbz       0.0036 0.0708    -0.1340   0.0031     0.1440  0.0020   0
## 
## Random effects:
## Name   Model
##  struct   IID model 
## 
## Model hyperparameters:
##                          mean       sd 0.025quant 0.5quant 0.975quant
## Precision for struct 18523.35 18316.96    1250.29 13124.05   66740.20
##                         mode
## Precision for struct 3389.29
## 
## Expected number of effective parameters(std dev): 8.184(0.1415)
## Number of equivalent replicates : 31.04 
## 
## Deviance Information Criterion (DIC) ...: 758.71
## Effective number of parameters .........: 8.183
## 
## Marginal log-Likelihood:  -420.41 
## Posterior marginals for linear predictor and fitted values computed

#US model
fit.pre_u<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "iid"), family = "poisson",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.pre_u)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"iid\"), family = \"poisson\", ",  "    data = usdat@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.6976          5.1539          0.3661          6.2175 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0303 0.0155    -0.0002   0.0304     0.0606  0.0304   0
## metro        0.0623 0.0219     0.0193   0.0623     0.1054  0.0623   0
## giniz       -0.0034 0.0078    -0.0187  -0.0034     0.0119 -0.0035   0
## hispz       -0.0753 0.0089    -0.0929  -0.0753    -0.0578 -0.0753   0
## blackz       0.1240 0.0130     0.0984   0.1240     0.1495  0.1240   0
## povz         0.1547 0.0123     0.1305   0.1547     0.1789  0.1547   0
## ne_bad      -0.2093 0.0220    -0.2523  -0.2093    -0.1662 -0.2093   0
## isowbz      -0.0554 0.0147    -0.0843  -0.0554    -0.0265 -0.0555   0
## 
## Random effects:
## Name   Model
##  struct   IID model 
## 
## Model hyperparameters:
##                        mean    sd 0.025quant 0.5quant 0.975quant   mode
## Precision for struct 141.38 26.09      98.76   138.36     200.59 132.13
## 
## Expected number of effective parameters(std dev): 149.09(16.08)
## Number of equivalent replicates : 20.55 
## 
## Deviance Information Criterion (DIC) ...: 10201.46
## Effective number of parameters .........: 149.41
## 
## Marginal log-Likelihood:  -5170.65 
## Posterior marginals for linear predictor and fitted values computed

A more complicated, Besag ICAR model:

fit.psp1<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "besag", graph = "~/Google Drive//dem7263/txgraph.gra"), family = "poisson",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.psp1)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"besag\", graph = \"~/Google Drive//dem7263/txgraph.gra\"), ",  "    family = \"poisson\", data = tx@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.3167          0.3921          0.0958          1.8046 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0200 0.0678    -0.1548  -0.0194     0.1114 -0.0183   0
## metro        0.0435 0.0831    -0.1189   0.0432     0.2071  0.0427   0
## giniz        0.0614 0.0490    -0.0344   0.0613     0.1578  0.0610   0
## hispz       -0.0048 0.0274    -0.0585  -0.0048     0.0491 -0.0049   0
## blackz       0.1319 0.0895    -0.0428   0.1315     0.3083  0.1308   0
## povz         0.0303 0.0433    -0.0551   0.0304     0.1149  0.0307   0
## ne_bad      -0.3151 0.0935    -0.5000  -0.3147    -0.1327 -0.3139   0
## isowbz       0.0038 0.0707    -0.1336   0.0033     0.1440  0.0022   0
## 
## Random effects:
## Name   Model
##  struct   Besags ICAR model 
## 
## Model hyperparameters:
##                          mean       sd 0.025quant 0.5quant 0.975quant
## Precision for struct 18085.21 18132.07    1131.83 12690.37   65785.32
##                         mode
## Precision for struct 3024.75
## 
## Expected number of effective parameters(std dev): 8.101(0.0666)
## Number of equivalent replicates : 31.35 
## 
## Deviance Information Criterion (DIC) ...: 758.64
## Effective number of parameters .........: 8.10
## 
## Marginal log-Likelihood:  -611.12 
## Posterior marginals for linear predictor and fitted values computed

#Change the default prior on the spatial variance:

fit.psp2<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "besag", graph = "~/Google Drive//dem7263/txgraph.gra", hyper = list(theta=list(prior="loggamma", param=c(.5, .005)))), family = "poisson",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.psp2)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"besag\", graph = \"~/Google Drive//dem7263/txgraph.gra\", ",  "    hyper = list(theta = list(prior = \"loggamma\", param = c(0.5, ",  "        0.005)))), family = \"poisson\", data = tx@data, E = e + ",  "    1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.7854          0.2865          0.1054          1.1772 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0138 0.0728    -0.1579  -0.0134     0.1281 -0.0127   0
## metro        0.0345 0.0849    -0.1316   0.0344     0.2015  0.0340   0
## giniz        0.0638 0.0512    -0.0363   0.0637     0.1647  0.0633   0
## hispz       -0.0041 0.0303    -0.0636  -0.0042     0.0555 -0.0042   0
## blackz       0.1588 0.0969    -0.0290   0.1580     0.3512  0.1564   0
## povz         0.0257 0.0479    -0.0693   0.0260     0.1191  0.0265   0
## ne_bad      -0.2899 0.1006    -0.4876  -0.2899    -0.0926 -0.2898   0
## isowbz       0.0074 0.0752    -0.1381   0.0066     0.1571  0.0052   0
## 
## Random effects:
## Name   Model
##  struct   Besags ICAR model 
## 
## Model hyperparameters:
##                        mean     sd 0.025quant 0.5quant 0.975quant  mode
## Precision for struct 172.04 166.09      28.19   123.40     609.94 67.99
## 
## Expected number of effective parameters(std dev): 13.30(2.418)
## Number of equivalent replicates : 19.09 
## 
## Deviance Information Criterion (DIC) ...: 757.46
## Effective number of parameters .........: 13.56
## 
## Marginal log-Likelihood:  -611.08 
## Posterior marginals for linear predictor and fitted values computed

1/fit.psp2$summary.hyperpar$mean

## [1] 0.005812638

#US Models
fit.psp1_u<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "besag", graph = "~/Google Drive//dem7263/usa.gra"), family = "poisson",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.psp1_u)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"besag\", graph = \"~/Google Drive//dem7263/usa.gra\"), ",  "    family = \"poisson\", data = usdat@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.5401          8.2291          0.3298         10.0990 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0279 0.0156    -0.0028   0.0279     0.0583  0.0280   0
## metro        0.0581 0.0221     0.0147   0.0581     0.1015  0.0581   0
## giniz        0.0008 0.0080    -0.0148   0.0008     0.0164  0.0007   0
## hispz       -0.0565 0.0117    -0.0793  -0.0565    -0.0334 -0.0566   0
## blackz       0.1073 0.0148     0.0782   0.1073     0.1363  0.1073   0
## povz         0.1475 0.0131     0.1217   0.1476     0.1733  0.1476   0
## ne_bad      -0.1634 0.0243    -0.2110  -0.1634    -0.1157 -0.1634   0
## isowbz      -0.0658 0.0155    -0.0962  -0.0659    -0.0354 -0.0659   0
## 
## Random effects:
## Name   Model
##  struct   Besags ICAR model 
## 
## Model hyperparameters:
##                       mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for struct 70.71 14.69      46.84    68.98     104.29 65.51
## 
## Expected number of effective parameters(std dev): 110.90(11.68)
## Number of equivalent replicates : 27.63 
## 
## Deviance Information Criterion (DIC) ...: 10099.51
## Effective number of parameters .........: 111.24
## 
## Marginal log-Likelihood:  -7506.17 
## Posterior marginals for linear predictor and fitted values computed

1/fit.psp1_u$summary.hyperpar$mean

## [1] 0.01414282

#Change the default prior on the spatial variance:

fit.psp2_u<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "besag", graph = "~/Google Drive//dem7263/usa.gra", hyper = list(theta=list(prior="loggamma", param=c(.5, .005)))), family = "poisson",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.psp2_u)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"besag\", graph = \"~/Google Drive//dem7263/usa.gra\", ",  "    hyper = list(theta = list(prior = \"loggamma\", param = c(0.5, ",  "        0.005)))), family = \"poisson\", data = usdat@data, E = e + ",  "    1e-04, control.compute = list(dic = T))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.3087          8.2034          0.3854          9.8975 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0278 0.0156    -0.0030   0.0278     0.0582  0.0279   0
## metro        0.0580 0.0221     0.0145   0.0580     0.1014  0.0580   0
## giniz        0.0008 0.0080    -0.0149   0.0008     0.0166  0.0008   0
## hispz       -0.0562 0.0118    -0.0792  -0.0562    -0.0331 -0.0563   0
## blackz       0.1072 0.0149     0.0780   0.1072     0.1363  0.1072   0
## povz         0.1475 0.0132     0.1216   0.1475     0.1734  0.1475   0
## ne_bad      -0.1630 0.0244    -0.2108  -0.1630    -0.1152 -0.1630   0
## isowbz      -0.0659 0.0155    -0.0963  -0.0659    -0.0354 -0.0660   0
## 
## Random effects:
## Name   Model
##  struct   Besags ICAR model 
## 
## Model hyperparameters:
##                       mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for struct 68.21 13.86      45.57    66.61      99.80 63.41
## 
## Expected number of effective parameters(std dev): 113.35(11.69)
## Number of equivalent replicates : 27.03 
## 
## Deviance Information Criterion (DIC) ...: 10100.05
## Effective number of parameters .........: 113.67
## 
## Marginal log-Likelihood:  -7501.95 
## Posterior marginals for linear predictor and fitted values computed

1/fit.psp2_u$summary.hyperpar$mean

## [1] 0.0146609

You notice less variation in the US hyperparameter than in the Texas hyperparameter, that’s becuase we have more data in the US example.

Compare the marginal distributions for the precision parameter from each model. The precision = 1/\(\sigma^2\), so we need to transform it to a variance scale so we can interpret it as we normally would!

1/fit.psp1$summary.hyperpar$mean

## [1] 5.52938e-05

m<- inla.tmarginal(
        function(x) 1/x,
        fit.psp1$marginals.hyperpar$'Precision for struct')

plot(m, type="l", main="Posterior distibution for model hyperparameter - Texas data", xlim=c(0, .05))

m2<- inla.tmarginal(
        function(x) 1/x,
        fit.psp2$marginals.hyperpar$'Precision for struct')
lines(m2,lty=3, lwd=3    )

Quite a different story, depending on the prior for the Texas example

#US example
1/fit.psp1_u$summary.hyperpar$mean

## [1] 0.01414282

um<- inla.tmarginal(
        function(x) 1/x,
        fit.psp1_u$marginals.hyperpar$'Precision for struct')

plot(um, type="l", main="Posterior distibution for model hyperparameter - US data")

1/fit.psp2_u$summary.hyperpar$mean

## [1] 0.0146609

um2<- inla.tmarginal(
        function(x) 1/x,
        fit.psp2_u$marginals.hyperpar$'Precision for struct')
lines(um2,lty=3, lwd=3    )

The two different priors show very similar results in terms of the posterior distribution of the county level variance.

Let’s have a look at the raw vs fitted SMR for infant mortality:

tx$smr<-tx$infd03/tx$e
tx$fittedsmr<-fit.psp2$summary.fitted.values$mean
plot(tx$smr, tx$fittedsmr)

brks<-quantile(c(tx$smr, tx$fittedsmr), p=c(0, .25, .5, .75, 1), na.rm=T)

brks

##        0%       25%       50%       75%      100% 
## 0.0000000 0.8090127 1.0081150 1.1274526 4.7424671

spplot(tx, c("smr", "fittedsmr"), col.regions=brewer.pal(n=5, "RdBu"), at=brks , main="Raw and Fitted SMR")

spplot(tx, "fittedsmr", col.regions=brewer.pal(n=5, "RdBu"), at=brks, main="Fitted SMR")

tx$spre<-fit.psp2$summary.random$struct$mean
spplot(tx, "spre", col.regions=brewer.pal(n=5, "RdBu"), at=quantile(tx$spre, p=c(0, .25, .5, .75, 1)), main="Fitted Spatial Random Effect")

##US example
usdat$smr<-usdat$infd03/usdat$e
usdat$fittedsmr<-fit.psp2_u$summary.fitted.values$mean
plot(usdat$smr, usdat$fittedsmr)

brks<-quantile(c(usdat$smr, usdat$fittedsmr), p=seq(0, 1,length.out = 9), na.rm=T)

brks

##        0%     12.5%       25%     37.5%       50%     62.5%       75% 
## 0.0000000 0.6108312 0.8172425 0.9131458 1.0001002 1.1017731 1.2424354 
##     87.5%      100% 
## 1.5388164 8.6480283

spplot(usdat, c("smr", "fittedsmr"), col.regions=brewer.pal(n=9, "Greys"), at=brks , main="Raw and Fitted SMR")

spplot(usdat, "fittedsmr", col.regions=brewer.pal(n=9, "Greys"), at=brks, main="US County Fitted SMR")

usdat$spre<-fit.psp2_u$summary.random$struct$mean
spplot(usdat, "spre", col.regions=brewer.pal(n=9, "Greys"), at=quantile(usdat$spre, p=seq(0, 1,length.out = 9)), main="US County Fitted Spatial Random Effect")

The Besag, York and Mollie model. This model has two random effects, one which is spatially correlated, and the other which is iid normal. In INLA, and when using the BYM model, the first 1:n rows of the random effect are the combined spatial and “nonspatial” effects.

#Texas data
fit.tx_bym<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "bym", graph = "~/Google Drive//dem7263/txgraph.gra"), family = "poisson",E=e+.0001, data=tx@data, control.compute = list(dic=T))
summary(fit.tx_bym)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"bym\", graph = \"~/Google Drive//dem7263/txgraph.gra\"), ",  "    family = \"poisson\", data = tx@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          0.8384          0.7886          0.1069          1.7339 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.0207 0.0690    -0.1578  -0.0201     0.1133 -0.0190   0
## metro        0.0408 0.0837    -0.1228   0.0406     0.2057  0.0401   0
## giniz        0.0627 0.0501    -0.0352   0.0626     0.1614  0.0623   0
## hispz       -0.0055 0.0284    -0.0611  -0.0055     0.0503 -0.0055   0
## blackz       0.1359 0.0920    -0.0434   0.1355     0.3176  0.1346   0
## povz         0.0305 0.0448    -0.0579   0.0306     0.1181  0.0309   0
## ne_bad      -0.3066 0.0959    -0.4959  -0.3063    -0.1191 -0.3057   0
## isowbz       0.0021 0.0726    -0.1388   0.0015     0.1460  0.0003   0
## 
## Random effects:
## Name   Model
##  struct   BYM model 
## 
## Model hyperparameters:
##                                             mean      sd 0.025quant
## Precision for struct (iid component)     2072.21 1911.92     200.33
## Precision for struct (spatial component) 1572.40 1682.88      71.62
##                                          0.5quant 0.975quant   mode
## Precision for struct (iid component)      1534.28    7137.04 573.04
## Precision for struct (spatial component)  1047.72    6071.74 169.48
## 
## Expected number of effective parameters(std dev): 10.21(1.505)
## Number of equivalent replicates : 24.87 
## 
## Deviance Information Criterion (DIC) ...: 758.45
## Effective number of parameters .........: 10.37
## 
## Marginal log-Likelihood:  -378.10 
## Posterior marginals for linear predictor and fitted values computed

1/fit.tx_bym$summary.hyperpar$mean

## [1] 0.0004825768 0.0006359723

#US data
fit.us_bym<-inla(infd03~metro+giniz+hispz+blackz+povz+ne_bad+isowbz + f(struct, model = "bym", graph = "~/Google Drive//dem7263/usa.gra"), family = "poisson",E=e+.0001, data=usdat@data, control.compute = list(dic=T))
summary(fit.us_bym)

## 
## Call:
## c("inla(formula = infd03 ~ metro + giniz + hispz + blackz + povz + ",  "    ne_bad + isowbz + f(struct, model = \"bym\", graph = \"~/Google Drive//dem7263/usa.gra\"), ",  "    family = \"poisson\", data = usdat@data, E = e + 1e-04, control.compute = list(dic = T))" )
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##          1.3513         42.4625          0.7703         44.5840 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept)  0.0278 0.0156    -0.0029   0.0279     0.0583  0.0280   0
## metro        0.0579 0.0221     0.0145   0.0579     0.1014  0.0579   0
## giniz        0.0009 0.0080    -0.0149   0.0008     0.0166  0.0008   0
## hispz       -0.0569 0.0117    -0.0798  -0.0569    -0.0337 -0.0571   0
## blackz       0.1071 0.0149     0.0779   0.1071     0.1362  0.1071   0
## povz         0.1476 0.0132     0.1217   0.1476     0.1734  0.1477   0
## ne_bad      -0.1635 0.0244    -0.2115  -0.1636    -0.1156 -0.1636   0
## isowbz      -0.0659 0.0155    -0.0964  -0.0660    -0.0354 -0.0660   0
## 
## Random effects:
## Name   Model
##  struct   BYM model 
## 
## Model hyperparameters:
##                                             mean      sd 0.025quant
## Precision for struct (iid component)     3288.92 2278.51     697.90
## Precision for struct (spatial component)   73.44   15.73      48.08
##                                          0.5quant 0.975quant    mode
## Precision for struct (iid component)      2729.63    9217.78 1761.33
## Precision for struct (spatial component)    71.52     109.59   67.68
## 
## Expected number of effective parameters(std dev): 115.96(14.49)
## Number of equivalent replicates : 26.42 
## 
## Deviance Information Criterion (DIC) ...: 10102.88
## Effective number of parameters .........: 116.53
## 
## Marginal log-Likelihood:  -4689.67 
## Posterior marginals for linear predictor and fitted values computed

1/fit.us_bym$summary.hyperpar$mean

## [1] 0.0003040515 0.0136171355

usdat$bymre<-fit.us_bym$summary.random$struct$mean[1:nrow(usdat@data)]
spplot(usdat, "bymre", col.regions=brewer.pal(n=9, "Greys"), at=quantile(usdat$bymre, p=seq(0, 1,length.out = 9)), main="US County Fitted Spatial Random Effect")

Exceedence probabilities

Here we calculate our exceedence probabilities, in this case, I calculate the probability that a county has a SMR greater than 1.25.

a2<-1.25
usdat$inlaprob<-unlist(lapply(fit.us_bym$marginals.fitted.values, function(X){
  1-inla.pmarginal(a2, X)
}))
hist(usdat$inlaprob)

usdat$exceed.cut<-cut(usdat$inlaprob,breaks =  c(0,0.2,0.5,0.8,.95,1), include.lowest = T)

spplot(usdat,"exceed.cut" , col.regions=brewer.pal(n=6, name="Greys"), key.space="bottom", main= expression(paste("Exceedence Probability Infant SMR ","Pr( ",theta," >1.25"," )") ))