## Model setup

• We have a count outcome (deaths and births), in counties over time, and a set of time-constant covariates.
• We have several options in the GLM framework with which to model these data, for example:

• Binomial – $y_{ij} \sim Bin(\pi_{ij}) \text{: } logit(\pi_{ij} ) = \beta_{0}+ x’\beta_k$

• Poisson – $y_{ij} \sim Pois(\lambda_{ij} E_{ij}) \text{: } log(\lambda_{ij} ) = log(E_{ij}) + \beta_{0}+ x’\beta_k$

• Negative Binomial – $y_{ij} \sim \text{Neg Bin} (\mu_{ij}, \alpha, E_{ij}) \text{: } log(\mu_{ij} ) = log(E_{ij}) + \beta_{0}+ x’\beta_k$

• In addition to various zero-inflated versions of these data.

## stat_bin() using bins = 30. Pick better value with binwidth.
## Warning: Removed 1849 rows containing non-finite values (stat_bin).  ## stat_bin() using bins = 30. Pick better value with binwidth.
## Warning: Removed 1849 rows containing non-finite values (stat_bin).  ## summarise() ungrouping output (override with .groups argument)  We can fit these model using the Bayesian framework with INLA.

First, we consider the basic GLM for the mortality outcome, with out any hierarchical structure. We can write this model as a Negative Binomial model, for instance as:

$\text{Deaths}_{ij} \sim NB(\mu_{ij}, \gamma)$ $\mu_{ij} = \text{log(E_d)}_{ij} + X’ \beta$

INLA will use vague Normal priors for the $$\beta$$’s, and we have other parameters in the model to specify priors for. INLA does not require you to specify all priors, as all parameters have a default prior specification. In this example, I will use a $$Gamma(1, .5)$$ prior for all hierarchical variance terms.

## ## Call:
## c("inla(formula = f1, family = \"nbinomial\", data = final.dat, E = ## E_d, ", " verbose = F, control.compute = list(waic = T), ## control.predictor = list(link = 1), ", " num.threads = 2)") ## Time used:
## Pre = 0.928, Running = 21.8, Post = 0.722, Total = 23.5 ## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -5.047 10.723 -26.102 -5.048 15.989 -5.047 0
## scale(pblack) 0.159 0.015 0.130 0.159 0.188 0.159 0
## scale(phisp) -0.025 0.013 -0.050 -0.025 0.001 -0.025 0
## scale(ppov) 0.041 0.015 0.012 0.041 0.070 0.041 0
## year 0.003 0.005 -0.008 0.003 0.013 0.003 0
## ## Model hyperparameters:
## mean sd 0.025quant
## size for the nbinomial observations (1/overdispersion) 0.624 0.009 0.608
## 0.5quant 0.975quant
## size for the nbinomial observations (1/overdispersion) 0.624 0.641
## mode
## size for the nbinomial observations (1/overdispersion) 0.624
## ## Expected number of effective parameters(stdev): 5.04(0.001)
## Number of equivalent replicates : 2124.92 ## ## Watanabe-Akaike information criterion (WAIC) ...: 114586.38
## Effective number of parameters .................: 10.27
## ## Marginal log-Likelihood: -57331.80 ## Posterior marginals for the linear predictor and
## the fitted values are computed

Plot our observed vs fitted values  ### Basic county level random intercept model

Now we add basic nesting of rates within counties, with a random intercept term for each county. This would allow there to be heterogeneity in the mortality rate for each county, over and above each county’s observed characteristics.

This model would be:

$\text{Deaths}_{ij} \sim NB(\mu_{ij}, \gamma)$ $\mu_{ij} = \text{log(E_d)}_{ij} + X’ \beta + u_j$ $u_j \sim \text{Normal} (0 , \tau_u)$

where $$\tau_u$$ here is the precision, not the variance and precision = 1/variance. INLA puts a log-gamma prior on the the precision by default.

f2~scale(pblack)+scale(phisp)+scale(ppov)+year+ #fixed effects
f(struct, model = "iid",param=c(1,.5)) #random effects mod2inla(formula = f2,data = final.dat,
family = "nbinomial", E = E_d,
control.compute = list(waic=T),
verbose = F) #total model summary
summary(mod2)
## ## Call:
## c("inla(formula = f2, family = \"nbinomial\", data = final.dat, E = ## E_d, ", " verbose = F, control.compute = list(waic = T), ## control.predictor = list(link = 1), ", " num.threads = 2)") ## Time used:
## Pre = 0.571, Running = 160, Post = 1.36, Total = 162 ## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -2.824 10.758 -23.945 -2.824 18.279 -2.824 0
## scale(pblack) 0.158 0.015 0.128 0.158 0.189 0.158 0
## scale(phisp) -0.041 0.014 -0.069 -0.041 -0.013 -0.041 0
## scale(ppov) 0.044 0.015 0.014 0.044 0.074 0.044 0
## year 0.001 0.005 -0.009 0.001 0.012 0.001 0
## ## Random effects:
## Name Model
## struct IID model
## ## Model hyperparameters:
## mean sd 0.025quant
## size for the nbinomial observations (1/overdispersion) 0.627 0.009 0.609
## Precision for struct 50.626 7.005 38.292
## 0.5quant 0.975quant
## size for the nbinomial observations (1/overdispersion) 0.627 0.644
## Precision for struct 50.138 65.780
## mode
## size for the nbinomial observations (1/overdispersion) 0.626
## Precision for struct 49.174
## ## Expected number of effective parameters(stdev): 125.34(15.33)
## Number of equivalent replicates : 85.47 ## ## Watanabe-Akaike information criterion (WAIC) ...: 114610.09
## Effective number of parameters .................: 66.26
## ## Marginal log-Likelihood: -57375.58 ## Posterior marginals for the linear predictor and
## the fitted values are computed

#### Marginal Distributions of hyperparameters

We can plot the posterior marginal of the hyperparameter in this model, in this case $$\sigma_u = 1/\tau_u$$

## low high
## level:0.95 0.01491462 0.02565338  final.dat$fitted_m2$summary.fitted.values$mean p1%>% filter(year%in%c(2000))%>% mutate(qrr=cut(fitted_m2, breaks = quantile(final.dat$fitted_m2, p=seq(0,1,length.out = 6)), include.lowest = T))%>%
ggplot()+geom_sf(aes(fill=qrr))+scale_colour_brewer(palette = "RdBu" )+scale_fill_brewer(palette = "RdBu", na.value="grey")+guides(fill=guide_legend(title="Relative Risk Quartile"))+ggtitle(label="Relative Risk Quartile - IID Model, 2000")+coord_sf(crs = 7603) p2%>%
filter(year%in%c(2007))%>%
mutate(qrr=cut(fitted_m2, breaks = quantile(final.dat$fitted_m2, p=seq(0,1,length.out = 6)), include.lowest = T))%>% ggplot()+geom_sf(aes(fill=qrr))+scale_colour_brewer(palette = "RdBu" )+scale_fill_brewer(palette = "RdBu", na.value="grey")+guides(fill=guide_legend(title="Relative Risk Quartile"))+ggtitle(label="Relative Risk Quartile - IID Model, 2007")+coord_sf(crs = 7603) library(gridExtra) ## ## Attaching package: 'gridExtra' ## The following object is masked from 'package:dplyr': ## ## combine  ## TableGrob (2 x 1) "arrange": 2 grobs ## z cells name grob ## 1 1 (1-1,1-1) arrange gtable[layout] ## 2 2 (2-2,1-1) arrange gtable[layout] ### BYM Model Model with spatial correlation – Besag, York, and Mollie (1991) model and temporal heterogeneity $\text{Deaths}_{ij} \sim NB(\mu_{ij}, \gamma)$ $\mu_{ij} = \text{log(E_d)}_{ij} + X’ \beta + u_j + v_j + \gamma_t$ Which has two random effects, one an IID random effect and the second a spatially correlated random effect, specified as a conditionally auto-regressive prior for the $$v_j$$’s. This is the Besag model: $v_j|v_{\neq j},\sim\text{Normal}(\frac{1}{n_i}\sum_{i\sim j}v_j,\frac{1}{n_i\tau})$ and $$u_j$$ is an IID normal random effect, $$\gamma_t$$ is also given an IID Normal random effect specification, and there are now three hyperparameters, $$\tau_u$$ and $$\tau_v$$ and $$\tau_{\gamma}$$ and each are given log-gamma priors. For the BYM model we must specify the spatial connectivity matrix in the random effect. #final.dat$year_c
f3~scale(pblack)+scale(phisp)+scale(ppov)+
f(struct, model = "bym", constr = T, scale.model = T, graph = H,param=c(1,.5))+
f(year, model="iid",param=c(1,.5)) #temporal random effect mod3inla(formula = f3,data = final.dat,
family = "nbinomial", E = E_d,
control.compute = list(waic=T),
verbose = F,
control.predictor = list(link=1)) #total model summary
summary(mod3)
## ## Call:
## c("inla(formula = f3, family = \"nbinomial\", data = final.dat, E = ## E_d, ", " verbose = F, control.compute = list(waic = T), ## control.predictor = list(link = 1), ", " num.threads = 2)") ## Time used:
## Pre = 0.737, Running = 138, Post = 1.26, Total = 140 ## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.115 0.129 -0.145 0.115 0.374 0.115 0
## scale(pblack) 0.157 0.016 0.126 0.158 0.189 0.158 0
## scale(phisp) -0.039 0.016 -0.069 -0.039 -0.007 -0.040 0
## scale(ppov) 0.043 0.016 0.012 0.043 0.075 0.043 0
## ## Random effects:
## Name Model
## struct BYM model
## year IID model
## ## Model hyperparameters:
## mean sd
## size for the nbinomial observations (1/overdispersion) 0.627 0.009
## Precision for struct (iid component) 51.094 7.099
## Precision for struct (spatial component) 1974.289 1903.577
## Precision for year 8.760 4.130
## 0.025quant 0.5quant
## size for the nbinomial observations (1/overdispersion) 0.609 0.627
## Precision for struct (iid component) 38.602 50.591
## Precision for struct (spatial component) 174.447 1425.658
## Precision for year 2.885 8.075
## 0.975quant mode
## size for the nbinomial observations (1/overdispersion) 0.644 0.628
## Precision for struct (iid component) 66.447 49.595
## Precision for struct (spatial component) 7055.730 496.592
## Precision for year 18.742 6.583
## ## Expected number of effective parameters(stdev): 133.75(15.30)
## Number of equivalent replicates : 80.09 ## ## Watanabe-Akaike information criterion (WAIC) ...: 114605.76
## Effective number of parameters .................: 69.81
## ## Marginal log-Likelihood: -56934.15 ## Posterior marginals for the linear predictor and
## the fitted values are computed  m3a inla.tmarginal(
function(x) (1/x),
mod3$marginals.hyperpar$Precision for struct (iid component))
m3b inla.tmarginal(
function(x) (1/x),
mod3$marginals.hyperpar$Precision for struct (spatial component))
m3c inla.tmarginal(
function(x) (1/x),
mod3$marginals.hyperpar$Precision for year) plot(m3a, type="l", main=c("Posterior distibution for between county variance", "- BYM model -"), xlim=c(0, .2), ylim=c(0, 300))
lines(m3b, col="red")
lines(m3c, col="green")
legend("topright", legend=c("BYM IID", "BYM Spatial", "Year"), col=c(1, "red", "green"), lty=c(1,1,1))  ## low high
## level:0.95 0.01475866 0.02544088
## low high
## level:0.95 0.00005416961 0.003970123
## low high
## level:0.95 0.03927999 0.2945931

This indicates very low spatially correlated variance in these data.

## Exceedence probabilities

In Bayesian spatial models that are centered on an epidemiological type of outcome, it is common to examine the data for spatial clustering. One way to do this is to examine the clustering in the relative risk from one of these GLMM models. For instance if $$\theta$$ is the relative risk $\theta = exp(\beta_0 + \beta_1*x_1 + u_j)$ from one of our Negative binomial models above. We can use the posterior marginals of the relative risk to ask $$\theta \gt \theta^*$$ where $$\theta^*$$ is a specific level of excess risk, say 50% extra or $$\theta > 1.25$$. If the density, or $$\text{Pr}(\theta \gt \theta^*)$$ is high, then there is evidence that the excess risk is not only high, but significantly high.

To get the exceedence probabilities from one of our models, we can use the inla.pmarginal() function to ask if $$\text{Pr}(\theta \gt \theta^*)$$  So, we see lots of occasions where the exceedence probability is greater than .9. We can visualize these in a map.

final.dat\$exceedprob

final.dat%>%
filter(year%in%c(2007))%>%
mutate(qrr=cut(exceedprob, breaks = c(0, .5, .9, .95, .99, 1), include.lowest = T))%>%
ggplot()+geom_sf(aes(fill=qrr))+scale_colour_brewer(palette = "Blues" )+scale_fill_brewer(palette = "Blues", na.value="grey")+guides(fill=guide_legend(title=""))+ggtitle(label=expression(paste("Exceedence Probability Relative Risk ","Pr( ",theta," >1.25"," ) - 2007") ))+coord_sf(crs = 7603)  #library(mapview) #map1%
# filter(year%in%c(2007))%>%
# mutate(qrr=cut(exceedprob, breaks = c(0, .5, .9, .95, .99, 1), include.lowest = T)) #clrs
#mapView(as(map1, "Spatial"), zcol="qrr", legend=T, col.regions=clrs, map.types="OpenStreetMap")

Which shows several areas of the south where risk the infant mortality rate is signficantly higher than the national rate, with high posterior probability.