Correcting spatial autocorrelation using SAR - different results with brms vs spatialreg packages?

Interfaces brms
1

I am using a simultaneous autoregressive lag model to account for spatial autocorrelation. When using a Moran I test to check model residuals for autocorrelation, it is still present with a brms SAR model. If I run the same model (I think it’s the same) in spatialreg and use a Moran I test on the residuals, the autocorrelation is gone. If I try it with a regular lm(), Moran’s I is lower than the one from the brms SAR. I am stumped. I am also new to brms. Is it not be appropriate to use Moran to test brms model residuals for spatial autocorrelation?

Thanks in advance for your thoughts!
C

## reproducible example
library(brms)
library(spatialreg)
library(spdep)
data(oldcol, package = "spdep")
listw = nb2listw(COL.nb)
# brms SAR
fit1 <- brm(CRIME ~ INC + HOVAL + sar(listw, type = "lag"), 
            data = COL.OLD, data2 = list(listw = listw),
            chains = 2, cores = 2)
# spatialreg SAR
fit2 <- spatialreg::lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw=listw)
# uncorrected model
fit3 <- lm(CRIME ~ INC + HOVAL, data=COL.OLD)
# extract residuals
bres <- residuals(fit1)[,"Estimate"]
sres <- residuals(fit2)
nres <- residuals(fit3)
# Moran test
moran.test(bres, listw) # brms SAR: I = 0.26
moran.test(sres, listw) # spatialreg SAR: I = 0.03
moran.test(nres, listw) # nothing: I = 0.23```



Operating System: Windows 10
brms 2.13.0
spatialreg 1.1.5
spdep 1.1.3
R 4.0.1
1 Like
2

By default, brms uses the mean predictions to compute residuals; predictions that do not include the residual autocorrelation. Hence the results you are seeing. I would try out residuals(method = "posterior_predict") which will include SAR terms in the predictions used to compute residuals.

3

Thanks so much for the quick reply! However this doesn’t seem to get us there. As above, the spatailreg SAR regression seems to deal with the autocorrelation well (I = 0.03) but the autocorrelation remains high in the brms posterior_predict residuals (I = 0.27). Any other thoughts?

Thanks in advance once again!
Colin

fit1 ← brm(CRIME ~ INC + HOVAL + sar(listw, type = “lag”),
data = COL.OLD, data2 = list(listw = listw),
chains = 2, cores = 2)
bres ← residuals(fit1, method = “posterior_predict”)[,“Estimate”]
moran.test(bres, listw) # brms post predect resid SAR: I = 0.27

4

You can see ?spatialreg::predict.sarlm for some information on how to calculate the de-trended residuals (i.e., residuals with the spatial trend cut out)

1 Like
5

Thanks. I’ll try to figure things out from there. Apologies, if I’m missing something. But based on Paul’s response above it seems the posterior_predict residuals from the brms SAR model should not be autocorrelated? The lagsarlm model fixes the autocorrelation well and I’m comfortable with. We’d like to replicate the results in brms and the continued autocorrelation in the brms SAR model that has me stuck.

C

6

I was hypothesising it fixes it. but I need to think of it in more detail and try it out myself. the problem is definitely not trivial and to me it isn’t even clear what the desired behavior should be (in brms or in other packages for that matter).

2 Likes
7

Thanks again! I’m quite sure I’m at the limit of my abilities but am sticking with it. I found my way to this discussion of different ways to calculate SAR residuals (Spatial econometrics -- computing residuals - Cross Validated) and thought I’d share in case it is useful.

Cheers,
Colin

8

I have checked again and brms incorporates the SAR parameters correclty in the model when using posterior_predict (you can double check by looking at the code in brms:::posterior_predict_gaussian_lagsar). Unfortunately, I am out of depth with respect to the moran.test and surrounding statistics so I don’t know how that should translate to the residuals. One potential problem I am seeing is that you use the Estimate column rather than the individual posterior draws, but I am spectulating here.

9

It looks like that method brms:::posterior_predict_gaussian_lagsar is taking draws from the SAR model, generically: y \sim Gau(\mu, \Sigma) with the SAR specification of the covariance matrix. That will generate spatial autocorrelation (SA) in the residuals, y - \mu (as it should).

The model implicitly has a spatial trend. What spdep is doing is calculating the implicit spatial trend embedded in the covariance matrix, then it calculates residuals as basically y - mu - SA (just the general idea). Those ‘detrended’ residuals are what @ColinKG is getting from spdep.

When you remove autocorrelation from the residuals, then the Moran coefficient becomes a small negative number (like -0.05).

2 Likes
10

Thanks @cmcd ! That’s exactly the information we were looking for. Very clear now.

1 Like