Skip to content

How to use create_WX() correctly

2 messages · Michael E. Rose, Roger Bivand

1 day later
#
On Wed, 20 May 2015, Michael E. Rose wrote:

            
The original use scenario for create_WX() was for internal calculations in 
GNS, SDM, SDEM and SLX models. It was expected that the first argument 
would be taken from model.matrix(), but partly ignored the possibility 
that the formula could be without an intercept. If the style of the 
spatial weights is W (row-standardised), the lagged intercept needed to be 
dropped, but other cases led to index mis-counting. So yes, this was a 
bug. From SVN revision 634 on R-forge, users can now do:

library(spdep)
data(oldcol)
lw <- nb2listw(COL.nb, style="W")
X <- COL.OLD[, c("INC", "HOVAL")]
WX <- create_WX(X, lw, prefix="my_lag")
head(WX)
X1 <- model.matrix(CRIME ~ INC + HOVAL - 1, data=COL.OLD)
WX1 <- create_WX(X1, lw, prefix="my_lag1")
head(WX1)
X2 <- model.matrix(CRIME ~ INC + HOVAL, data=COL.OLD)
WX2 <- create_WX(X2, lw, prefix="my_lag2")
head(WX2)

and get the same values in the WX* matrices.

With regard to impacts in the current model:

y = (I - \rho W)^{-1}(Xb + W*Xd)

where W != W*, I believe that a unit change in x_r (the r_th covariate) 
will be S(W)_r = (I - \rho W)^{-1}(Ib_r + W*d_r). If your number of 
observations (n) is moderate, you can calculate this as a dense matrix, 
taking the mean of the diagonal of S(W)_r as the direct impacts and the 
sum of all the elements in S(W)_r divided by n as the total impacts.

You could also do this by predicting with the original data as newdata, 
saving the prediction, incrementing x_r by 1 and replacing its lag in W*X 
by the lag of the incremented values, and predicing with the incremented 
newdata. The mean of the difference between the predictions is the total 
impact. See LeSage & Pace 2009 for details.

It will be wrong to see x_r and W*x_r as separate variables in terms of 
impacts, as the unit increment of x_r enters both terms.

Hope this helps,

Roger