Hi folks,
I've been running some multilevel (random effects) models where I have 1049
residents nested within 52 neighborhood clusters. My outcome variable is a
survey-based measure of residents' perceived levels of neighborhood
problems. It can reasonably be modeled as a continuous, normally distributed
variable. My request for assistance is related to setting up a conditional
autoregressive model for the neighborhood-level intercepts in WinBUGS, so
below I show the model as it stands before adding the CAR structure and how
I figured out that spatial autocorrelation should probably be incorporated
into the model. Then I ask my specific questions about how to get CAR
working correctly in WinBUGS.
I'm using R2WinBUGS to send the data out to WinBUGS and coda to read the
MCMC results back into R. I've tested the model shown below, which
essentially assumes that the level 2 neighborhood residuals (denoted u.j)
are independent of each other and follow a normal distribution.
# WinBUGS Code
model {
# LEVEL 1 MODEL for Residents, i = 1, ..., n = 1049:
# y.ij = B.0j + B.1*age.36.55 + B.2*age.56up + B.3*female +
# B.4*race.B + B.5*race.H + B.6*race.O + B.7*marital.M +
# B.8*marital.D + B.9*marital.W + B.10*educ.NoHS +
# B.11*educ.UGD + B.12*educ.PGD + B.13*employed +
# B.14*inc.0.15 + B.15*inc.15.25 + B.16*inc.25.45 +
# B.17*homeowner + B.18*kids + r.ij
# y.ij = Outcome score for each resident
# B.0j = Random intercept for each cluster
# B.1 to B.18 = Regression coefficients for level 1 predictors
# r.ij = Level 1 residuals, ~ Normal(mean = 0, var = s.y^2)
for (i in 1:n) {
y[i] ~ dnorm(y.hat[i], tau.y) # Distribution of the outcome scores
# Level 1 Fixed Effects (omit the r.ij term)
y.hat[i] <- B.0j[clusterid[i]] + B.1*age.36.55[i] + B.2*age.56up[i] +
B.3*female[i] + B.4*race.B[i] + B.5*race.H[i] +
B.6*race.O[i] +
B.7*marital.M[i] + B.8*marital.D[i] + B.9*marital.W[i] +
B.10*educ.NoHS[i] + B.11*educ.UGD[i] + B.12*educ.PGD[i] +
B.13*employed[i] + B.14*inc.0.15[i] + B.15*inc.15.25[i] +
B.16*inc.25.45[i] + B.17*homeowner[i] + B.18*kids[i]
}
#---------------------------------------------------------------------------
-
# LEVEL 1 FIXED EFFECTS
# Non-informative priors for B.1 to B.18: Normal(mean = 0, var = 100^2)
B.1 ~ dnorm(0, .0001)
B.2 ~ dnorm(0, .0001)
B.3 ~ dnorm(0, .0001)
B.4 ~ dnorm(0, .0001)
B.5 ~ dnorm(0, .0001)
B.6 ~ dnorm(0, .0001)
B.7 ~ dnorm(0, .0001)
B.8 ~ dnorm(0, .0001)
B.9 ~ dnorm(0, .0001)
B.10 ~ dnorm(0, .0001)
B.11 ~ dnorm(0, .0001)
B.12 ~ dnorm(0, .0001)
B.13 ~ dnorm(0, .0001)
B.14 ~ dnorm(0, .0001)
B.15 ~ dnorm(0, .0001)
B.16 ~ dnorm(0, .0001)
B.17 ~ dnorm(0, .0001)
B.18 ~ dnorm(0, .0001)
# LEVEL 1 RANDOM EFFECT (Variance component based on r.ij residual term)
tau.y ~ dgamma(2, 1) # tau.y = precision for y = 1/(s.y^2)
#---------------------------------------------------------------------------
-
# LEVEL 2 MODEL for Clusters, j = 1, ..., J = 52:
# B.0j = G.00 + G.2*cvperson.cl + G.3*hv.cl + u.j
# B.0j = Random intercept for each cluster
# G.00 = Mean of cluster-level intercepts (grand mean of the outcome)
# G.2 & G.3 = Regression coefficients for level 2 predictors
# u.j = Level 2 residuals, ~ Normal(mean = 0, var = s.B0j^2)
for (j in 1:J) {
B.0j[j] ~ dnorm(B.0j.hat[j], tau.B0j)
# Level 2 Fixed Effects
B.0j.hat[j] <- G.00 + G.2*cvperson.cl[j] + G.3*hv.cl[j]
}
# Non-informative priors for G.00, G.2, & G.3: Normal(mean = 0, var =
100^2)
G.00 ~ dnorm(0, .0001)
G.2 ~ dnorm(0, .0001)
G.3 ~ dnorm(0, .0001)
#---------------------------------------------------------------------------
-
# LEVEL 2 RANDOM EFFECT (Variance component based on u.j term)
tau.B0j ~ dgamma(2, 1) # tau.B0j = precision for B.0j = 1/(s.B0j^2)
#---------------------------------------------------------------------------
-
# CALCULATE THE ICC & OTHER PARAMETERS TO SAVE.
s2.y <- 1/tau.y # Level 1 variance
s.y <- sqrt(s2.y) # Level 1 SD
s2.B0j <- 1/tau.B0j # Level 2 variance (intercept)
s.B0j <- sqrt(s2.B0j) # Level 2 SD (intercept)
tot.var <- s2.y + s2.B0j # Total variance
ICC <- s2.B0j/tot.var # Intraclass correlation (ICC)
}
The model runs fine, but I had reasons to expect spatial autocorrelation
between the neighborhoods that decreases with increasing distance between
them. Since I had the boundaries and attributes of the clusters in a
SpatialPolygonsDataFrame named clusters, I extracted the u.j residuals and
added them as an additional variable to the clusters object and ran a
Moran's I test as shown in the R code below. The inverse-distance weights
are applied to other neighborhoods within half the diagonal distance across
the study area. Clusters further apart than that get weights of zero.
># Create distance-based neighbor list for the clusters
>clusters.nb <- dnearneigh(x = coordinates(clusters), d1 = 0, d2 =
clusters.diag/2,
+ row.names = row.names(clusters), longlat = F)
>
># Create a row-standarized, inverse-distance spatial weights list.
> dsts <- nbdists(clusters.nb, coords = coordinates(clusters), longlat = F)
> idw <- lapply(dsts, function(x) 1/x)
> clusters.idw <- nb2listw(clusters.nb, glist = idw, style = "W")
>
># Copy level 2 residuals from HLM model into the clusters data frame
>clusters$hlm.7.u.j <- as.vector(HLM.L2$hlm.7.u.j) # HLM Model 7
>
># Run Moran's I test (see Bivand, Pebesma, & Gomez-Rubio, 2008, pp.
258-264)
># on level 2 residuals from HLM Model 7
>
lm.morantest.exact(lm(hlm.7.u.j ~ 1, data = clusters), listw = clusters.idw)
Global Moran's I statistic with exact p-value
data:
model:lm(formula = hlm.7.u.j ~ 1, data = clusters)
weights: clusters.idw
Exact standard deviate = 1.7477, p-value = 0.04026
alternative hypothesis: greater
sample estimates:
[1] 0.04342825
So, as expected, there is some evidence of spatial autocorrelation in those
neighborhood-level residuals. Great. So now I want to integrate that CAR
structure back into the WinBUGS model so I can have the autocorrelation
parameter estimated there and also see how it would affect the level 2
regression coefficients, the DIC, and so on. Here's where I run into
difficulties. I think the way to do this is to replace the level 2 model
with something like this:
B.0j[1:J] ~ car.proper(B.0j.hat[], C[], adj[], num[], M[], tau.B0j, gamma)
for (j in 1:J) {
# Level 2 Fixed Effects
B.0j.hat[j] <- G.00 + G.2*cvperson.cl[j] + G.3*hv.cl[j]
}
# Non-informative priors for G.00, G.2, & G.3: Normal(mean = 0, var = 100^2)
G.00 ~ dflat()
G.2 ~ dnorm(0, .0001)
G.3 ~ dnorm(0, .0001)
gamma.min <- min.bound(C[], adj[], num[], M[])
gamma.max <- max.bound(C[], adj[], num[], M[])
gamma ~ dunif(gamma.min, gamma.max)
tau.B0j ~ dgamma(2, 1) # tau.B0j = precision for B.0j = 1/(s.B0j^2)
Is that the right way to do this? If so, then I need to set the arguments of
car.proper (particularly C[] and M[]). My best guess so far is that C[]
should correspond to the my.weights object created in R code below. Is that
correct? I'd like to get some help either confirming or correcting that
guess. Finally, how do I set M[] to get this to work correctly?
All the examples I've seen are actually for Poisson models of the outcome,
whereas I'm working with a normal distribution on the outcome. So I'm not
sure how to adapt the method for setting the M[] values in those examples to
apply to my case.
# Create the adjacency, weights, and number of neighbors vectors WinBUGS
needs
# for a CAR spatial random effect
my.adj <- listw2WB(clusters.idw)$adj
my.weights <- listw2WB(clusters.idw)$weights
my.num <- listw2WB(clusters.idw)$num
Any help in getting this CAR structure properly specified in the WinBUGS
model would be deeply appreciated.
Steven J. Pierce, M.S.
Doctoral Student in Ecological/Community Psychology
Department of Psychology
Michigan State University
Web: http://www.psychology.msu.edu/eco/
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
