impacts() referring to wrong coefficients for SDEM and SLX with lagged intercept

[ruettenauer]

Dear Roger,

I just went through your teaching materials of ECS530. Thanks a lot for this valuable resource, it is really incredibly helpful!

However, I was a bit confused by seeing the indirect impacts after running a Durbin error model in ECS530_h19/ECS530_VII.Rmd, line 118. If I'm not confusing something here, the impacts() function after the Durbin error and the SLX seems to take up the wrong lagged coefficients from the model if a lagged intercept is included. See the example (taken from your materials) below:

library(spatialreg)
library(spdep)
library(sf)

nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
st_crs(nc) <- "+proj=longlat +datum=NAD27"
row.names(nc) <- as.character(nc$FIPSNO)
nc$ft.SID74 <- sqrt(1000)*(sqrt(nc$SID74/nc$BIR74) + sqrt((nc$SID74+1)/nc$BIR74))
nc$ft.NWBIR74 <- sqrt(1000)*(sqrt(nc$NWBIR74/nc$BIR74) + sqrt((nc$NWBIR74+1)/nc$BIR74))
gal_file <- system.file("weights/ncCR85.gal", package="spData")[1]
ncCR85 <- read.gal(gal_file, region.id=nc$FIPSNO)
closeAllConnections()
nc_lw <- nb2listw(ncCR85, style="B")
nc_lw_w <- nb2listw(ncCR85, style="W") # row-standardized to test without intercept

# SDEM with unstandardized W and intercept in lagX
m1b <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m1b)
summary(spatialreg::impacts(m1b))

# The indirect impact does not equal the lagged coef of ft.NWBIR74
all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.ft.NWBIR74"], 
          tolerance = 1e-5, check.attributes = FALSE)

# it equals the lagged intercept
all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.(Intercept)"], 
          tolerance = 1e-5, check.attributes = FALSE)


# Similarly, with more than one covariate, it starts with the intercept, and omits the last lag coefficient
m1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m1c)
summary(spatialreg::impacts(m1c))

# Indirect impact for "east" equals the lagged coefficient for "ft.NWBIR74"
all.equal(spatialreg::impacts(m1c)[[1]]$indirect["east"], m1c$coefficients["lag.ft.NWBIR74"], 
          tolerance = 1e-5, check.attributes = FALSE)


# Without intercept (row standardized W), everything works as expected
m1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE)
summary(m1c)
summary(spatialreg::impacts(m1c))



### Same problem occurs in impacts(lmSLX()) with lagged intercept
m2 <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE)
summary(m2)
summary(spatialreg::impacts(m2))

# but everything works as expected without intercept
m2b <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE)
summary(m2b)
summary(spatialreg::impacts(m2b))

To me, this looks like the lagged intercept is not taken into account when extracting the lagged coefficients for the impacts in SDEM and SLX (I did not test any other models so far)?

I'm happy to look into the code if you agree that there's something going wrong here.

Best wishes,
Tobias

[rsbivand]

Tobias, Thanks, it is likely that only row standardized weights were considered. Please investigate further; the involved code should have been modularized but most often was modified from other functions. The `Durbin=` formula is tricky, and that took attention away from other cases, like a lagged intercept. Roger Roger Bivand Falsensvei 32 5063 Bergen fre. 17. jul. 2020, 17:04 skrev Tobias Rüttenauer <notifications@github.com

…

: Dear Roger, I just went through your teaching materials of ECS530. Thanks a lot for this valuable resource, it is really incredibly helpful! However, I was a bit confused by seeing the indirect impacts after running a Durbin error model in ECS530_h19/ECS530_VII.Rmd, line 118. If I'm not confusing something here, the impacts() function after the Durbin error and the SLX seems to take up the wrong lagged coefficients from the model if a lagged intercept is included. See the example (taken from your materials) below: library(spatialreg) library(spdep) library(sf) nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE) st_crs(nc) <- "+proj=longlat +datum=NAD27" row.names(nc) <- as.character(nc$FIPSNO)nc$ft.SID74 <- sqrt(1000)*(sqrt(nc$SID74/nc$BIR74) + sqrt((nc$SID74+1)/nc$BIR74))nc$ft.NWBIR74 <- sqrt(1000)*(sqrt(nc$NWBIR74/nc$BIR74) + sqrt((nc$NWBIR74+1)/nc$BIR74))gal_file <- system.file("weights/ncCR85.gal", package="spData")[1]ncCR85 <- read.gal(gal_file, region.id=nc$FIPSNO) closeAllConnections()nc_lw <- nb2listw(ncCR85, style="B")nc_lw_w <- nb2listw(ncCR85, style="W") # row-standardized to test without intercept # SDEM with unstandardized W and intercept in lagXm1b <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE) summary(m1b) summary(spatialreg::impacts(m1b)) # The indirect impact does not equal the lagged coef of ft.NWBIR74 all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.ft.NWBIR74"], tolerance = 1e-5, check.attributes = FALSE) # it equals the lagged intercept all.equal(spatialreg::impacts(m1b)[[1]]$indirect, m1b$coefficients["lag.(Intercept)"], tolerance = 1e-5, check.attributes = FALSE) # Similarly, with more than one covariate, it starts with the intercept, and omits the last lag coefficientm1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74 + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE) summary(m1c) summary(spatialreg::impacts(m1c)) # Indirect impact for "east" equals the lagged coefficient for "ft.NWBIR74" all.equal(spatialreg::impacts(m1c)[[1]]$indirect["east"], m1c$coefficients["lag.ft.NWBIR74"], tolerance = 1e-5, check.attributes = FALSE) # Without intercept (row standardized W), everything works as expectedm1c <- spatialreg::errorsarlm(ft.SID74 ~ ft.NWBIR74 + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE) summary(m1c) summary(spatialreg::impacts(m1c)) ### Same problem occurs in impacts(lmSLX()) with lagged interceptm2 <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74 + east, weights=BIR74, data=nc, listw=nc_lw, Durbin=TRUE) summary(m2) summary(spatialreg::impacts(m2)) # but everything works as expected without interceptm2b <- spatialreg::lmSLX(ft.SID74 ~ ft.NWBIR74 + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE) summary(m2b) summary(spatialreg::impacts(m2b)) To me, this looks like the lagged intercept is not taken into account when extracting the lagged coefficients for the impacts in SDEM and SLX (I did not test any other models so far)? I'm happy to look into the code if you agree that there's something going wrong here. Best wishes, Tobias — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#13>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZ3BCFXHGEPDZ5XRD6AVDR4BR7HANCNFSM4O6MXRPA> .

[ruettenauer]

I'll have a more detailed look into this and try to find a solution to the problem within the next couple of days. I'll let you know once I have any updates on the issue.

Best
Tobias

[ruettenauer]

So, if I got that right, the problem does not occur in the impacts() function, but in lmSLX() and errorsarlm() when the object "indirImps" is constructed. I made some very small changes in my fork, which solve the problem of the example I posted above: ruettenauer@0b5656f. It assumes that there will always be a lagged intercept if listw$style == "W". I hope I got his right.

However, I have not checked if this problem affects also other situations than the one in which I made these changes. Especially, I haven't checked if this also affects situations with Durbin= formula. I'll have a look into this, but will only be able to do this tomorrow or on Monday.

I also recognized that the lmSLX() and errorsarlm returns an error if the main formula has no intercept, e.g. in the example above:

m2c <- spatialreg::lmSLX(ft.SID74 ~ -1 + ft.NWBIR74  + east, weights=BIR74, data=nc, listw=nc_lw_w, Durbin=TRUE)

Is there a theoretical reason for this behaviour? I thought a model without intercept might be of interest e.g. if estimating a fixed effects model by using OLS on the pre-demeaned data.

Best
Tobias

[jtsayagog]

Hello, I think SLX and SDEM do not require a separate impact analysis. Elhorst has this table for the direct and indirect effects in these models, I attach it here. [image: image.png] However, I am trying to see the lag of the intercept as source of interaction. Best, Juan

…

On Fri, Jul 17, 2020 at 4:54 PM Tobias Rüttenauer ***@***.***> wrote: I'll have a more detailed look into this and try to find a solution to the problem within the next couple of days. I'll let you know once I have any updates on the issue. Best Tobias — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#13 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB5DBXUFAQ4GEGTK2NO5XPTR4DCARANCNFSM4O6MXRPA> .

[rsbivand]

The methods are needed to carry out linear combination of the standard errors. Roger Roger Bivand Falsensvei 32 5063 Bergen lør. 18. jul. 2020, 13:11 skrev Juan Tomás Sayago <notifications@github.com

…

: Hello, I think SLX and SDEM do not require a separate impact analysis. Elhorst has this table for the direct and indirect effects in these models, I attach it here. [image: image.png] However, I am trying to see the lag of the intercept as source of interaction. Best, Juan On Fri, Jul 17, 2020 at 4:54 PM Tobias Rüttenauer < ***@***.***> wrote: > I'll have a more detailed look into this and try to find a solution to the > problem within the next couple of days. I'll let you know once I have any > updates on the issue. > > Best > Tobias > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub > < #13 (comment)>, > or unsubscribe > < https://github.com/notifications/unsubscribe-auth/AB5DBXUFAQ4GEGTK2NO5XPTR4DCARANCNFSM4O6MXRPA > > . > — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#13 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZ3BGKYFMAMOAJRFDYHJ3R4F7NTANCNFSM4O6MXRPA> .

[ruettenauer]

Dear Roger,

I did some minor adaptions in the lmSLX() function to have a more flexible intercept handling but still receive the correct coefficients in the impacts method, which is in my fork of the package: ruettenauer@4336121. This allows to estimate models without a "direct intercept" in the main formula, and works no matter whether the lagged terms include an intercept or not. However, I'm not sure if I did something completely odd here, as for instance, I did not get the reason for the condition if (K == 1 && odd) warning("model configuration issue: no total impacts") else >>compute impacts<< (SLX_WX.R, old line 106, new line 115 in the commit). I skipped this condition in my recent commit.

I also uploaded some code with all the tests I performed: https://github.com/ruettenauer/mixed/blob/master/spatialreg_intercept-tests_lmslx_errorsarlm.R. All the included tests pass on my current fork version.

If you think all this makes sense, I'm happy to create a pull request into the spatialreg repository (I should probably turn my in-code comments into github comments or delete them before?)? However, let me know if I did something stupid here or if I did not consider important situations in the tests.

Best
Tobias

[ruettenauer]

Sorry, just recognized that I had an error in the earlier comment. Of course the first commit assumes that there is always an intercept if listw$style == "W" is not true. But the recent commit also allows to not have an intercept if listw$style != "W" and the Durbin formula includes a "-1".

[rsbivand]

@ruettenauer Tobias, may I please include your tests in https://github.com/ruettenauer/mixed/blob/master/spatialreg_intercept-tests_lmslx_errorsarlm.R in a tinytest test file in the package?

[ruettenauer]

Sure, you are more than welcome to use the tests, @rsbivand! I am glad if the tests are useful.

[rsbivand]

Definitely useful! I'm also looking at using multcomp::glht() directly internally, initially as an option.