Model setup
 We have a count outcome (deaths and births), in counties over time, and a set of timeconstant 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 zeroinflated versions of these data.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1849 rows containing nonfinite values (stat_bin).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1849 rows containing nonfinite 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 ## ## WatanabeAkaike information criterion (WAIC) ...: 114586.38
## Effective number of parameters .................: 10.27
## ## Marginal logLikelihood: 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 loggamma 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),
control.predictor = list(link=1),
num.threads = 2,
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 ## ## WatanabeAkaike information criterion (WAIC) ...: 114610.09
## Effective number of parameters .................: 66.26
## ## Marginal logLikelihood: 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 (11,11) arrange gtable[layout]
## 2 2 (22,11) 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 autoregressive prior for the \(v_j\)’s. This is the Besag model:
\[v_jv_{\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 loggamma 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),
num.threads = 2,
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 ## ## WatanabeAkaike information criterion (WAIC) ...: 114605.76
## Effective number of parameters .................: 69.81
## ## Marginal logLikelihood: 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.
Like!! Great article post.Really thank you! Really Cool.
These are actually great ideas in concerning blogging.
Thank you ever so for you article post.
I am regular visitor, how are you everybody? This article posted at this web site is in fact pleasant.
bookmarked!!, I like your blog!