DEM 7473 - Spatial Hierarchical GLMM(s) using the INLA Approximation
Corey S. Sparks, Ph.D.
November 12, 2018
Multi - Level Models
We can also do hierarchical models with spatial effects in INLA. This example looks at correlation among MSAs in obesity.
#load brfss
library(car)
library(knitr)
library(dplyr)
library(INLA)
brfssurl<-"https://github.com/coreysparks/data/blob/master/brfss_14.Rdata?raw=true"
load(url(brfssurl))
set.seed(12345)
#samps<-sample(1:nrow(brfss_14), size = 90000, replace=F)
#brfss_14<-brfss_14[samps,]
#The names in the data are very ugly, so I make them less ugly
nams<-names(brfss_14)
#we see some names are lower case, some are upper and some have a little _ in the first position. This is a nightmare.
newnames<-gsub(pattern = "x_",replacement = "",x = nams)
names(brfss_14)<-tolower(newnames)
#BMI
brfss_14$bmi<-ifelse(is.na(brfss_14$bmi5)==T, NA, brfss_14$bmi5/100)
brfss_14$obese<-ifelse(brfss_14$bmi>30,1,0)
#Poor or fair self rated health
#brfss_14$badhealth<-ifelse(brfss_14$genhlth %in% c(4,5),1,0)
brfss_14$badhealth<-Recode(brfss_14$genhlth, recodes="4:5=1; 1:3=0; else=NA")
#race/ethnicity
brfss_14$black<-Recode(brfss_14$racegr3, recodes="2=1; 9=NA; else=0")
brfss_14$white<-Recode(brfss_14$racegr3, recodes="1=1; 9=NA; else=0")
brfss_14$other<-Recode(brfss_14$racegr3, recodes="3:4=1; 9=NA; else=0")
brfss_14$hispanic<-Recode(brfss_14$racegr3, recodes="5=1; 9=NA; else=0")
brfss_14$race_eth<-Recode(brfss_14$racegr3, recodes="1='nhwhite'; 2='nh black'; 3='nh other';
4='nh multirace'; 5='hispanic'; else=NA", as.factor = T)
brfss_14$race_eth<-relevel(brfss_14$race_eth, ref = "nhwhite")
#insurance
brfss_14$ins<-ifelse(brfss_14$hlthpln1==1,1,0)
#income grouping
brfss_14$inc<-ifelse(brfss_14$incomg==9, NA, brfss_14$incomg)
#education level
brfss_14$educ<-Recode(brfss_14$educa, recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad';
5='3somecol'; 6='4colgrad';9=NA", as.factor=T)
#brfss_14$educ<-relevel(brfss_14$educ, ref='0Prim')
#employment
brfss_14$employ<-Recode(brfss_14$employ, recodes="1:2='Employed'; 2:6='nilf';
7='retired'; 8='unable'; else=NA", as.factor=T)
brfss_14$employ<-relevel(brfss_14$employ, ref='Employed')
#marital status
brfss_14$marst<-Recode(brfss_14$marital, recodes="1='married'; 2='divorced'; 3='widowed';
4='separated'; 5='nm';6='cohab'; else=NA", as.factor=T)
brfss_14$marst<-relevel(brfss_14$marst, ref='married')
#Age cut into intervals
brfss_14$agec<-cut(brfss_14$age80, breaks=c(0,24,39,59,79,99), include.lowest = T)I want to see how many people we have in each MSA in the data:
#Now we will begin fitting the multilevel regression model with the msa
#that the person lives in being the higher level
head(data.frame(name=table(brfss_14$mmsaname),id=unique(brfss_14$mmsa)))## name.Var1
## 1 Aberdeen, SD, Micropolitan Statistical Area
## 2 Aguadilla-Isabela, PR, Metropolitan Statistical Area
## 3 Albuquerque, NM, Metropolitan Statistical Area
## 4 Allentown-Bethlehem-Easton, PA-NJ, Metropolitan Statistical Area
## 5 Anchorage, AK, Metropolitan Statistical Area
## 6 Atlanta-Sandy Springs-Roswell, GA, Metropolitan Statistical Area
## name.Freq id
## 1 620 10100
## 2 544 10380
## 3 1789 10740
## 4 1095 10900
## 5 1785 11260
## 6 2776 12060#people within each msa
#How many total MSAs are in the data?
length(table(brfss_14$mmsa))## [1] 132#countiesHigher level predictors
We will often be interested in factors at both the individual AND contextual levels. To illustrate this, I will use data from the American Community Survey measured at the MSA level. Specifically, I use the DP3 table, which provides economic characteristics of places, from the 2010 5 year ACS Link.
library(acs)
#Get 2010 ACS median household incomes for tracts in Texas
msaacs<-geo.make(msa="*")
acsecon<-acs.fetch( endyear=2010, span=5, geography=msaacs, variable = c("B19083_001","B17001_001","B17001_002", "B03002_001","B03002_004", "B03002_012" ))
colnames(acsecon@estimate)## [1] "B19083_001" "B17001_001" "B17001_002" "B03002_001" "B03002_004"
## [6] "B03002_012"msaecon<-data.frame(gini=acsecon@estimate[, "B19083_001"],
ppoverty=acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"],
pblack=acsecon@estimate[,"B03002_004"]/acsecon@estimate[, "B03002_001"],
phisp=acsecon@estimate[,"B03002_012"]/acsecon@estimate[, "B03002_001"],
giniz=scale(acsecon@estimate[, "B19083_001"]),
ppovertyz=scale(acsecon@estimate[, "B17001_002"]/acsecon@estimate[, "B17001_001"]))
msaecon$zpblack<-scale(msaecon$pblack)
msaecon$zphisp<-scale(msaecon$phisp)
msaecon$ids<-paste(acsecon@geography$metropolitanstatisticalareamicropolitanstatisticalarea)Let’s see the geographic variation in these economic indicators:
library(tigris)## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.##
## Attaching package: 'tigris'## The following object is masked from 'package:graphics':
##
## plotmsa<-core_based_statistical_areas(cb=T, year = 2014)## Warning in readOGR(dsn = tmp, layer = shape, encoding = "UTF-8", verbose =
## FALSE, : Z-dimension discardedmsa_ec<-geo_join(msa, msaecon, "CBSAFP", "ids", how="inner")
tx_ec<-msa_ec[grep(msa_ec$NAME, pattern = "TX"), ]
library(RColorBrewer)
library(sp)
spplot(tx_ec, "gini", at=quantile(tx_ec$gini), col.regions=brewer.pal(n=6, "Reds"), col="transparent", main="Gini Coefficient")
spplot(tx_ec, "phisp", at=quantile(tx_ec$phisp), col.regions=brewer.pal(n=6, "Reds"), col="transparent", main="Percent Hispanic")
Create spatial information for higher level units
#See what counties are in the brfss data
tx_ec$struct<-1:dim(tx_ec)[1]
city.est.dat<-tx_ec@data[,c( "giniz","ppovertyz", "zpblack", "zphisp", "struct")]
city.est.dat$obese<-NA
head(city.est.dat)## giniz ppovertyz zpblack zphisp struct obese
## 18 0.5617655 0.11998153 1.0868041 -0.35892680 1 NA
## 21 0.6583485 -0.43879447 0.4031748 0.83881520 2 NA
## 25 0.2398224 0.04468834 0.2847100 -0.30408008 3 NA
## 34 0.3042110 -0.18356177 0.9809344 -0.08645997 4 NA
## 56 1.6563723 1.20997059 -0.4320280 3.18954054 5 NA
## 72 0.6261542 0.05361338 -0.2224854 -0.07489272 6 NAbrfss_14$cbsa<-as.character(brfss_14$mmsa)
indat<-merge(brfss_14, tx_ec, by.x="cbsa", by.y="CBSAFP", all.x=T)
brf.est<-indat[, c("giniz","ppovertyz", "zpblack", "zphisp", "struct", "obese")]
brf.est<-brf.est[order(brf.est$struct),]
head(brf.est)## giniz ppovertyz zpblack zphisp struct obese
## 65291 0.9480973 1.349245 -0.5109484 3.903845 7 0
## 65292 0.9480973 1.349245 -0.5109484 3.903845 7 0
## 65293 0.9480973 1.349245 -0.5109484 3.903845 7 1
## 65294 0.9480973 1.349245 -0.5109484 3.903845 7 0
## 65295 0.9480973 1.349245 -0.5109484 3.903845 7 1
## 65296 0.9480973 1.349245 -0.5109484 3.903845 7 0##Here is where I add the cities that need to be estimated to the rest of the data
m.est<-rbind(city.est.dat, brf.est)
struct.in<-unique(brf.est$struct)
m.est$comp<-ifelse(m.est$struct%in%struct.in ,1,0)
m.est$rm<-ifelse(m.est$comp==1&is.na(m.est$obese)==T,1,0)
m.est<-m.est[-which(m.est$rm==1),]
m.est<-m.est[is.na(m.est$struct)==F,]
m.est<-m.est[order(m.est$struct),]
#
# fake_dat<-expand.grid(race_eth=levels(brfss_14$race_eth), agec=levels(brfss_14$agec), CBSAFP=levels(as.factor(tx_ec$CBSAFP) ))
# fake_dat<-merge(fake_dat, tx_ec, by="CBSAFP")
library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`nbs<-knearneigh(coordinates(tx_ec), longlat = T, k = 5)
nbs<-knn2nb(nbs, row.names = tx_ec$struct, sym = T)
mat <- nb2mat(nbs, style="B",zero.policy=TRUE)
colnames(mat) <- rownames(mat)
mat <- as.matrix(mat[1:dim(mat)[1], 1:dim(mat)[1]])fit_est<- inla(obese~ giniz+zpblack+zphisp+
f(struct, model="bym", graph=mat,constr=TRUE, scale.model=TRUE),
family = "binomial",Ntrials = 1,
data=m.est, num.threads = 2,
control.predictor = list(link=1))
# control.inla=list(strategy='gaussian'))
#
summary(fit_est)##
## Call:
## c("inla(formula = obese ~ giniz + zpblack + zphisp + f(struct, model = \"bym\", ", " graph = mat, constr = TRUE, scale.model = TRUE), family = \"binomial\", ", " data = m.est, Ntrials = 1, control.predictor = list(link = 1), ", " num.threads = 2)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.7965 210.9759 0.7979 214.5703
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -1.0170 0.1465 -1.2830 -1.0262 -0.6987 -1.0422 0.0009
## giniz 0.0803 0.0906 -0.1091 0.0833 0.2519 0.0882 0.0010
## zpblack -0.0315 0.1912 -0.4368 -0.0252 0.3395 -0.0149 0.0012
## zphisp 0.0848 0.0587 -0.0397 0.0871 0.1959 0.0910 0.0009
##
## Random effects:
## Name Model
## struct BYM model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for struct (iid component) 896.83 4749.38 17.18
## Precision for struct (spatial component) 1952.66 2101.58 116.53
## 0.5quant 0.975quant mode
## Precision for struct (iid component) 197.82 5872.76 38.85
## Precision for struct (spatial component) 1313.13 7579.98 307.20
##
## Expected number of effective parameters(std dev): 5.815(0.6675)
## Number of equivalent replicates : 1443.93
##
## Marginal log-Likelihood: -5120.29
## Posterior marginals for linear predictor and fitted values computedm.est$est.obese<-fit_est$summary.fitted.values$mean
m.est.est<-tapply(m.est$est.obese, m.est$struct, mean, na.rm=T)
m.est.est<-data.frame(struct=names(unlist(m.est.est)), obeseest=unlist(m.est.est))
#m.est<-m.est[is.na(m.est$bmi)==T,]
msa.est<-merge(tx_ec, m.est.est, by.y="struct", by.x="struct", all.x=T, sort=F)
head(msa.est@data)## struct CSAFP CBSAFP AFFGEOID GEOID NAME
## 1 1 <NA> 45500 310M200US45500 45500 Texarkana, TX-AR
## 2 2 206 19100 310M200US19100 19100 Dallas-Fort Worth-Arlington, TX
## 3 3 <NA> 37580 310M200US37580 37580 Paris, TX
## 4 4 346 32220 310M200US32220 32220 Marshall, TX
## 5 5 204 28780 310M200US28780 28780 Kingsville, TX
## 6 6 206 11980 310M200US11980 11980 Athens, TX
## LSAD ALAND AWATER gini ppoverty pblack phisp
## 1 M1 5290957359 213758310 0.458 0.1722829 0.23590709 0.05030979
## 2 M1 24034136107 923118950 0.461 0.1341090 0.14606033 0.26626250
## 3 M2 2349622248 67103114 0.448 0.1671391 0.13049096 0.06019864
## 4 M2 2330820963 40969626 0.450 0.1515457 0.22199305 0.09943552
## 5 M2 6060005913 1803211452 0.492 0.2467476 0.03629285 0.69009795
## 6 M2 2262941841 193237933 0.460 0.1677488 0.06383223 0.10152109
## giniz ppovertyz zpblack zphisp ids obeseest
## 1 0.5617655 0.11998153 1.0868041 -0.35892680 45500 0.2651257
## 2 0.6583485 -0.43879447 0.4031748 0.83881520 19100 0.2901915
## 3 0.2398224 0.04468834 0.2847100 -0.30408008 37580 0.2655554
## 4 0.3042110 -0.18356177 0.9809344 -0.08645997 32220 0.2658052
## 5 1.6563723 1.20997059 -0.4320280 3.18954054 28780 0.3600754
## 6 0.6261542 0.05361338 -0.2224854 -0.07489272 11980 0.2773943library(mapview)
clrs <- colorRampPalette(brewer.pal(8, "Blues"))
mapview(msa.est,"obeseest", col.regions=clrs, map.types="OpenStreetMap")Multi-level model in INLA
INLA is very good at doing multi-level models, and much faster than other modeling strategies. Plus we can build in the correlated random effects, such as the Besag model at the higher level of analysis. Below, we fit three models. A basic random intercept model: This corresponds to a multilevel logistic model with a higher level variables as predictors and can be written: \[logit(y_{ij}) = \beta_{0j} + \sum {\beta x_i}\]
\[\beta_{0j} = \beta_0 + u_j\]
with \[u_j \sim N(0, \tau_u)\]
library(INLA)
indat<-indat[is.na(indat$struct)==F,]
fit_in1<- inla(obese~ race_eth+agec+f(struct, model="iid"),
family = "binomial", Ntrials =1,
data=indat, num.threads = 2)
summary(fit_in1)##
## Call:
## c("inla(formula = obese ~ race_eth + agec + f(struct, model = \"iid\"), ", " family = \"binomial\", data = indat, Ntrials = 1, num.threads = 2)" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.1792 4.2317 0.1187 6.5295
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -2.0848 0.1450 -2.3747 -2.0830 -1.8048 -2.0795
## race_ethhispanic 0.5356 0.0577 0.4220 0.5357 0.6486 0.5359
## race_ethnh black 0.8601 0.0897 0.6838 0.8601 1.0358 0.8603
## race_ethnh multirace 0.2368 0.2268 -0.2189 0.2404 0.6719 0.2477
## race_ethnh other -0.6386 0.1808 -1.0047 -0.6347 -0.2947 -0.6268
## agec(24,39] 0.8669 0.1512 0.5761 0.8649 1.1692 0.8610
## agec(39,59] 1.2867 0.1435 1.0119 1.2845 1.5749 1.2800
## agec(59,79] 1.1911 0.1437 0.9158 1.1889 1.4797 1.1845
## agec(79,99] 0.3096 0.1701 -0.0209 0.3085 0.6466 0.3061
## kld
## (Intercept) 0
## race_ethhispanic 0
## race_ethnh black 0
## race_ethnh multirace 0
## race_ethnh other 0
## agec(24,39] 0
## agec(39,59] 0
## agec(59,79] 0
## agec(79,99] 0
##
## Random effects:
## Name Model
## struct IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for struct 10010.16 15928.03 41.07 2144.24 56482.98
## mode
## Precision for struct 69.52
##
## Expected number of effective parameters(std dev): 10.15(1.52)
## Number of equivalent replicates : 827.20
##
## Marginal log-Likelihood: -5012.441/fit_in1$summary.hyperpar## mean sd 0.025quant 0.5quant
## Precision for struct 9.98985e-05 6.278241e-05 0.02434962 0.0004663662
## 0.975quant mode
## Precision for struct 1.770445e-05 0.01438416m<- inla.tmarginal(
function(x) (1/x),
fit_in1$marginals.hyperpar$`Precision for struct`)
inla.hpdmarginal(.95, marginal=m)## low high
## level:0.95 5.533341e-06 0.02030636plot(m, type="l", main=c("Posterior distibution for between MSA variance", "- Random Intercept model -"))
That is our individual level, random intercept model. Now I will fit a model that includes MSA level demographic characteristics, the Gini index, %black and %hispanic. This corresponds to a multi-level model with higher level predictors:
\[y_{ij} = \beta_{0j} + \sum {\beta x_i} + \sum {\gamma z_j}\]
\[\beta_{0j} = \beta_0 + \sum {\gamma z_j}+ u_j\]
with \[u_j \sim N(0, \tau_u)\]
fit_in2<- inla(obese~ race_eth+agec+giniz+zpblack+zphisp+f(struct, model="iid"),
data=indat,
family="binomial", Ntrials = 1)
summary(fit_in2)##
## Call:
## c("inla(formula = obese ~ race_eth + agec + giniz + zpblack + zphisp + ", " f(struct, model = \"iid\"), family = \"binomial\", data = indat, ", " Ntrials = 1)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.2699 5.2320 0.1576 7.6595
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -2.1674 0.1864 -2.5166 -2.1738 -1.7750 -2.1819
## race_ethhispanic 0.5260 0.0609 0.4065 0.5260 0.6454 0.5260
## race_ethnh black 0.8586 0.0900 0.6817 0.8587 1.0349 0.8588
## race_ethnh multirace 0.2382 0.2269 -0.2176 0.2419 0.6735 0.2491
## race_ethnh other -0.6368 0.1808 -1.0030 -0.6329 -0.2929 -0.6250
## agec(24,39] 0.8705 0.1512 0.5796 0.8685 1.1729 0.8646
## agec(39,59] 1.2852 0.1435 1.0102 1.2830 1.5734 1.2785
## agec(59,79] 1.1869 0.1439 0.9112 1.1847 1.4758 1.1802
## agec(79,99] 0.2994 0.1704 -0.0315 0.2982 0.6368 0.2959
## giniz 0.0898 0.0795 -0.0863 0.0943 0.2321 0.0994
## zpblack -0.0632 0.1498 -0.4174 -0.0514 0.1934 -0.0411
## zphisp 0.0039 0.0486 -0.1093 0.0079 0.0869 0.0121
## kld
## (Intercept) 0.0008
## race_ethhispanic 0.0000
## race_ethnh black 0.0000
## race_ethnh multirace 0.0000
## race_ethnh other 0.0000
## agec(24,39] 0.0000
## agec(39,59] 0.0000
## agec(59,79] 0.0000
## agec(79,99] 0.0000
## giniz 0.0035
## zpblack 0.0066
## zphisp 0.0039
##
## Random effects:
## Name Model
## struct IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for struct 8888.77 15329.01 16.48 834.40 54400.43 32.91
##
## Expected number of effective parameters(std dev): 12.77(0.96)
## Number of equivalent replicates : 657.78
##
## Marginal log-Likelihood: -5029.871/fit_in2$summary.hyperpar## mean sd 0.025quant 0.5quant
## Precision for struct 0.0001125015 6.523578e-05 0.06067202 0.001198462
## 0.975quant mode
## Precision for struct 1.838221e-05 0.03038934m2<- inla.tmarginal(
function(x) (1/x),
fit_in2$marginals.hyperpar$`Precision for struct`)
inla.hpdmarginal(.95, marginal=m2)## low high
## level:0.95 7.739928e-06 0.04089993plot(m2, type="l", main=c("Posterior distibution for between MSA variance", "- Multi-level model -"))
Finally, we model the county level means using the Besag-York and Mollie model that we had used on the areal data last week. In this example, we are modeling the correlation between MSAs at the higher level, instead of like last week, where we were modeling the correlation among neighboring counties.
This is a multi level model with correlation between neighboring MSAs, assuming a Besag, York and Mollie, 1991 convolution prior for the random effects \[y_{ij} = \beta_{0} +\sum {\beta x_i} + \sum {\gamma z_j} + u_j + v_j\]
\[\beta_{0j} = \beta_0 + \sum {\gamma z_j}+ u_j+v_j\]
with
\[u_j \sim N(0, \tau_u)\]
\[v_j|v_{\neq j},\sim\text{Normal}(\frac{1}{n_i}\sum_{i\sim j}v_j,\frac{1}{n_i\tau})\]
fit_in3<- inla(obese~ race_eth+agec+giniz+zpblack+zphisp+
f(struct, model="bym", graph=mat,constr=TRUE, scale.model=TRUE),
family = "binomial", Ntrials = 1,
data=indat,
control.predictor = list(link=1))
summary(fit_in3)##
## Call:
## c("inla(formula = obese ~ race_eth + agec + giniz + zpblack + zphisp + ", " f(struct, model = \"bym\", graph = mat, constr = TRUE, scale.model = TRUE), ", " family = \"binomial\", data = indat, Ntrials = 1, control.predictor = list(link = 1))" )
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 2.5203 141.4598 0.6014 144.5814
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -2.1043 0.2128 -2.5134 -2.1083 -1.6714 -2.1150
## race_ethhispanic 0.5267 0.0609 0.4070 0.5267 0.6462 0.5267
## race_ethnh black 0.8590 0.0900 0.6820 0.8591 1.0354 0.8592
## race_ethnh multirace 0.2386 0.2269 -0.2174 0.2422 0.6740 0.2495
## race_ethnh other -0.6372 0.1809 -1.0035 -0.6332 -0.2931 -0.6254
## agec(24,39] 0.8720 0.1512 0.5811 0.8700 1.1744 0.8661
## agec(39,59] 1.2801 0.1435 1.0051 1.2779 1.5683 1.2734
## agec(59,79] 1.1771 0.1438 0.9016 1.1749 1.4659 1.1704
## agec(79,99] 0.2836 0.1702 -0.0471 0.2824 0.6208 0.2800
## giniz 0.0735 0.0991 -0.1330 0.0764 0.2628 0.0812
## zpblack -0.0978 0.2141 -0.5515 -0.0906 0.3200 -0.0782
## zphisp -0.0108 0.0657 -0.1496 -0.0083 0.1146 -0.0040
## kld
## (Intercept) 0.0004
## race_ethhispanic 0.0000
## race_ethnh black 0.0000
## race_ethnh multirace 0.0000
## race_ethnh other 0.0000
## agec(24,39] 0.0000
## agec(39,59] 0.0000
## agec(59,79] 0.0000
## agec(79,99] 0.0000
## giniz 0.0011
## zpblack 0.0016
## zphisp 0.0012
##
## Random effects:
## Name Model
## struct BYM model
##
## Model hyperparameters:
## mean sd 0.025quant
## Precision for struct (iid component) 196.08 345.00 13.34
## Precision for struct (spatial component) 2313.66 2704.55 144.71
## 0.5quant 0.975quant mode
## Precision for struct (iid component) 98.66 987.31 33.28
## Precision for struct (spatial component) 1486.66 9489.25 385.43
##
## Expected number of effective parameters(std dev): 13.97(0.6252)
## Number of equivalent replicates : 600.95
##
## Marginal log-Likelihood: -4988.31
## Posterior marginals for linear predictor and fitted values computed1/fit_in3$summary.hyperpar## mean sd
## Precision for struct (iid component) 0.0051000060 0.0028985706
## Precision for struct (spatial component) 0.0004322156 0.0003697472
## 0.025quant 0.5quant
## Precision for struct (iid component) 0.074961576 0.0101356012
## Precision for struct (spatial component) 0.006910347 0.0006726505
## 0.975quant mode
## Precision for struct (iid component) 0.0010128487 0.030050173
## Precision for struct (spatial component) 0.0001053824 0.002594482m3_sp<- inla.tmarginal(
function(x) (1/x),
fit_in3$marginals.hyperpar$`Precision for struct (spatial component)`)
m3_iid<- inla.tmarginal(
function(x) (1/x),
fit_in3$marginals.hyperpar$`Precision for struct (iid component)`)
inla.hpdmarginal(.95, marginal=m3_sp)## low high
## level:0.95 3.330795e-05 0.004717095inla.hpdmarginal(.95, marginal=m3_iid)## low high
## level:0.95 0.0002015436 0.05420964plot(m3_sp, type="l", main=c("Posterior distibution for between Spatial MSA variance", "- Multi-level model -"), xlim=c(0, .015))
lines(m3_iid, col=2,lty=2)
legend("topright", legend=c("Spatial Variance", "IID Variance"), col=c(1,2), lty=c(1,2))