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
W.nb <- poly2nb(lipdata.sp, row.names = lipdata$name))
W <- nb2mat(W.nb, style="B")
[,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
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)
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?
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)
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
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.
Hi @cmcd, would you not include
(1|name) in any ICAR model, or just those with
type = "BYM2?
@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
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.