Hierarchical CAR model in brms

Modeling
1

Hi,
I’m new into Bayesian analysis and I was trying to do a CAR model using brms. I was trying to replicate the results from the example with the Scotland Lip Cancer Dataset in https://mc-stan.org/users/documentation/case-studies/icar_stan.html, so I used these commands

library(brms)
library(spdep)

W.nb <- poly2nb(lipdata.sp, row.names = lipdata$name))
W <- nb2mat(W.nb, style="B")
lipdata$pcaff2<-pcaff*0.1
head(W)
          [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
Sutherland      0    0    0    0    0    0    0    0    0     0     0     0     0     0
Nairn           0    0    1    0    1    0    0    0    0     0     0     0     0     0
head(lipdata)
             observed  expected pcaff latitude longitude        name   pcaff2
Sutherland          9      1.4    16    57.29      5.50   Sutherland    1.6
Nairn              39      8.7    16    57.56      2.36        Nairn    1.6

fit<- brm(observed ~ pcaff2 + offset(log(expected)) + (1|name) + car(W.test2, gr=name, type="bym2"), data=lipdata, data2=list(W=W), family=poisson)

fit
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -0.44      0.13    -0.69    -0.18 1.01      684      911
pcaff2        0.63      0.13     0.36     0.89 1.01      632      651

but I’m not getting the same results. What am I doing wrong?

1 Like
2

Hi @BonScott,

The ICAR model with type = "BYM2" has the term (1|name) within it, so I’d get rid of that term in your model formula, like this:

fit<- brm(observed ~ pcaff2 + offset(log(expected)) + car(W.test2, gr=name, type="bym2"), data=lipdata, data2=list(W=W), family=poisson)

1 Like
3

Thanks. I’ve run the model without (1|name) and it gave me the same results, but if I run it without this term, I can’t get random effects estimates with ranef(fit)

4

I’m not a brms user so I can’t be much help with the best way to extract that term. I bet others who have used that could chime it with more help.

If you want that term specifically (rather than the sum of the spatially structured and unstructured parts, the entire BYM2 part) once you find the parameter you’re looking for you’ll probably also want to multiply it by its scale parameter and maybe use rho too, basically do some algebra to separate the spatial from the non-spatial term on the correct scale.

I’m working on an R package for spatial models that returns those parameters with the proper scaling etc as you’re looking for, though its a work in progress still.

1 Like
5

Hi @cmcd, would you not include (1|name) in any ICAR model, or just those with type = "BYM2?

6

@richard.turner If brms allows you to add an ICAR term to the model (not CAR and not BYM2), and then you decide to add the term (1|name) where name identifies each of your areal observations, then you have a common model specification (the BYM model). So in my opinion that’s not a bad choice (but just to be clear, I still don’t know much about brms).

If you’re using CAR or BYM2 instead of ICAR, then the term (1|name) is redundant.

1 Like