lmSLX fits an lm model augmented with the spatially lagged RHS variables, including the lagged intercept when the spatial weights are not row-standardised. create_WX creates spatially lagged RHS variables, and is exposed for use in model fitting functions.

lmSLX(formula, data = list(), listw, na.action, weights=NULL,
 Durbin=TRUE, zero.policy=NULL)
create_WX(x, listw, zero.policy=NULL, prefix="")
# S3 method for SlX
impacts(obj, ...)
# S3 method for WXimpact
print(x, ...)
# S3 method for WXimpact
summary(object, ..., adjust_k=(attr(object, "type") == "SDEM"))
# S3 method for SlX
predict(object, newdata, listw, zero.policy=NULL, ...)

Arguments

formula

a symbolic description of the model to be fit. The details of model specification are given for lm()

data

an optional data frame containing the variables in the model. By default the variables are taken from the environment which the function is called.

listw

a listw object created for example by nb2listw

na.action

a function (default options("na.action")), can also be na.omit or na.exclude with consequences for residuals and fitted values - in these cases the spatial weights list will be subsetted to remove NAs in the data. It may be necessary to set zero.policy to TRUE because this subsetting may create no-neighbour observations. Note that only weights lists created without using the glist argument to nb2listw may be subsetted.

weights

an optional vector of weights to be used in the fitting process. Non-NULL weights can be used to indicate that different observations have different variances (with the values in weights being inversely proportional to the variances); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations (including the case that there are w_i observations equal to y_i and the data have been summarized) - lm

Durbin

default TRUE for lmSLX (Durbin model including WX); if TRUE, full spatial Durbin model; if a formula object, the subset of explanatory variables to lag

zero.policy

default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA

obj

A spatial regression object created by lmSLX

...

Arguments passed through

prefix

default empty string, may be “lag” in some cases

x, object

model matrix to be lagged; lagImpact objects created by impacts methods

adjust_k

default TRUE if SDEM else FALSE, adjust internal OLS SDEM standard errors by dividing by n rather than (n-k) (default changed and bug fixed after 0.7-8; standard errors now ML in SDEM summary and impacts summary and identical - for SLX use FALSE)

newdata

data frame in which to predict --- if NULL, predictions are for the data on which the model was fitted. Should have row names corresponding to region.id. If row names are exactly the same than the ones used for training, it uses in-sample predictors for forecast.

Value

The lmSLX function returns an “lm” object with a “mixedImps” list of three impact matrixes (impacts and standard errors) for direct, indirect and total impacts; total impacts calculated using gmodels::estimable().

See also

Author

Roger Bivand Roger.Bivand@nhh.no

Examples

data(oldcol, package="spdep")
lw <- spdep::nb2listw(COL.nb, style="W")
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
summary(COL.SLX)
#> 
#> Call:
#> lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
#>     data = as.data.frame(x), weights = weights)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -36.536  -7.835   0.474   8.349  25.594 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  75.0287     6.6260  11.323 1.26e-14 ***
#> INC          -1.1089     0.3738  -2.967  0.00485 ** 
#> HOVAL        -0.2897     0.1014  -2.858  0.00649 ** 
#> lag.INC      -1.3710     0.5613  -2.443  0.01867 *  
#> lag.HOVAL     0.1918     0.2003   0.957  0.34369    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 10.93 on 44 degrees of freedom
#> Multiple R-squared:  0.6088,	Adjusted R-squared:  0.5732 
#> F-statistic: 17.12 on 4 and 44 DF,  p-value: 1.553e-08
#> 
summary(impacts(COL.SLX))
#> Impact measures (SlX, estimable, n-k):
#>           Direct   Indirect       Total
#> INC   -1.1089293 -1.3709725 -2.47990173
#> HOVAL -0.2897283  0.1917608 -0.09796753
#> ========================================================
#> Standard errors:
#>          Direct  Indirect     Total
#> INC   0.3738129 0.5612771 0.4965456
#> HOVAL 0.1013673 0.2003335 0.2028016
#> ========================================================
#> Z-values:
#>          Direct   Indirect      Total
#> INC   -2.966535 -2.4425945 -4.9943086
#> HOVAL -2.858202  0.9572079 -0.4830709
#> 
#> p-values:
#>       Direct    Indirect Total     
#> INC   0.0030118 0.014582 5.9047e-07
#> HOVAL 0.0042605 0.338462 0.62905   
#> 
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL + I(HOVAL^2), data=COL.OLD, listw=lw, Durbin=TRUE)
summary(impacts(COL.SLX))
#> Impact measures (SlX, estimable, n-k):
#>                  Direct     Indirect       Total
#> INC        -0.947594274 -1.275338647 -2.22293292
#> HOVAL      -0.777427839 -0.355048446 -1.13247628
#> I(HOVAL^2)  0.004639919  0.005608104  0.01024802
#> ========================================================
#> Standard errors:
#>                 Direct   Indirect      Total
#> INC        0.398832844 0.59687399 0.56325440
#> HOVAL      0.464456540 0.98213200 1.07348981
#> I(HOVAL^2) 0.004226259 0.00908385 0.01025212
#> ========================================================
#> Z-values:
#>               Direct   Indirect      Total
#> INC        -2.375918 -2.1366966 -3.9465877
#> HOVAL      -1.673844 -0.3615079 -1.0549483
#> I(HOVAL^2)  1.097879  0.6173709  0.9996008
#> 
#> p-values:
#>            Direct   Indirect Total     
#> INC        0.017505 0.032623 7.9273e-05
#> HOVAL      0.094161 0.717720 0.29145   
#> I(HOVAL^2) 0.272258 0.536990 0.31750   
#> 
summary(COL.SLX)
#> 
#> Call:
#> lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
#>     data = as.data.frame(x), weights = weights)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -33.974  -7.764   0.907   6.820  24.395 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)    92.459835  19.278187   4.796 2.06e-05 ***
#> INC            -0.947594   0.398833  -2.376   0.0221 *  
#> HOVAL          -0.777428   0.464457  -1.674   0.1016    
#> I.HOVAL.2.      0.004640   0.004226   1.098   0.2785    
#> lag.INC        -1.275339   0.596874  -2.137   0.0385 *  
#> lag.HOVAL      -0.355048   0.982132  -0.362   0.7195    
#> lag.I.HOVAL.2.  0.005608   0.009084   0.617   0.5403    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 10.99 on 42 degrees of freedom
#> Multiple R-squared:  0.6224,	Adjusted R-squared:  0.5684 
#> F-statistic: 11.54 on 6 and 42 DF,  p-value: 1.361e-07
#> 
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL + I(HOVAL^2), data=COL.OLD, listw=lw, Durbin=~INC)
summary(impacts(COL.SLX))
#> Impact measures (SlX, estimable, n-k):
#>                  Direct  Indirect        Total
#> INC        -1.079064628 -1.010896 -2.089960575
#> HOVAL      -0.634518755        NA -0.634518755
#> I(HOVAL^2)  0.003455273        NA  0.003455273
#> ========================================================
#> Standard errors:
#>                Direct  Indirect      Total
#> INC        0.38471071 0.4552667 0.44650193
#> HOVAL      0.44760078        NA 0.44760078
#> I(HOVAL^2) 0.00411036        NA 0.00411036
#> ========================================================
#> Z-values:
#>                Direct  Indirect      Total
#> INC        -2.8048728 -2.220448 -4.6807425
#> HOVAL      -1.4175997        NA -1.4175997
#> I(HOVAL^2)  0.8406254        NA  0.8406254
#> 
#> p-values:
#>            Direct    Indirect Total     
#> INC        0.0050336 0.026388 2.8584e-06
#> HOVAL      0.1563077 NA       0.15631   
#> I(HOVAL^2) 0.4005579 NA       0.40056   
#> 
summary(COL.SLX)
#> 
#> Call:
#> lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
#>     data = as.data.frame(x), weights = weights)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -35.908  -7.085   0.770   7.147  24.684 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 83.679968   9.264846   9.032  1.4e-11 ***
#> INC         -1.079065   0.384711  -2.805  0.00747 ** 
#> HOVAL       -0.634519   0.447601  -1.418  0.16335    
#> I.HOVAL.2.   0.003455   0.004110   0.841  0.40510    
#> lag.INC     -1.010896   0.455267  -2.220  0.03159 *  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 10.96 on 44 degrees of freedom
#> Multiple R-squared:  0.607,	Adjusted R-squared:  0.5712 
#> F-statistic: 16.99 on 4 and 44 DF,  p-value: 1.716e-08
#> 
COL.SLX <- lmSLX(CRIME ~ INC, data=COL.OLD, listw=lw)
summary(COL.SLX)
#> 
#> Call:
#> lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
#>     data = as.data.frame(x), weights = weights)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -49.297  -6.658   0.669   7.315  29.826 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  73.9849     6.2081  11.918 1.16e-15 ***
#> INC          -1.5889     0.3564  -4.458 5.28e-05 ***
#> lag.INC      -1.0859     0.4812  -2.257   0.0288 *  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Residual standard error: 11.65 on 46 degrees of freedom
#> Multiple R-squared:  0.5353,	Adjusted R-squared:  0.5151 
#> F-statistic: 26.49 on 2 and 46 DF,  p-value: 2.214e-08
#> 
summary(impacts(COL.SLX))
#> Impact measures (SlX, estimable, n-k):
#>        Direct  Indirect     Total
#> INC -1.588901 -1.085867 -2.674768
#> ========================================================
#> Standard errors:
#>        Direct  Indirect    Total
#> INC 0.3564039 0.4811809 0.407313
#> ========================================================
#> Z-values:
#>        Direct  Indirect     Total
#> INC -4.458147 -2.256671 -6.566861
#> 
#> p-values:
#>     Direct     Indirect Total     
#> INC 8.2671e-06 0.024029 5.1387e-11
#> 
if (FALSE) {
crds <- cbind(COL.OLD$X, COL.OLD$Y)
mdist <- sqrt(sum(diff(apply(crds, 2, range))^2))
dnb <- spdep::dnearneigh(crds, 0, mdist)
dists <- spdep::nbdists(dnb, crds)
f <- function(x, form, data, dnb, dists, verbose) {
  glst <- lapply(dists, function(d) 1/(d^x))
  lw <- spdep::nb2listw(dnb, glist=glst, style="B")
  res <- logLik(lmSLX(form=form, data=data, listw=lw))
  if (verbose) cat("power:", x, "logLik:", res, "\n")
  res
}
opt <- optimize(f, interval=c(0.1, 4), form=CRIME ~ INC + HOVAL,
 data=COL.OLD, dnb=dnb, dists=dists, verbose=TRUE, maximum=TRUE)
glst <- lapply(dists, function(d) 1/(d^opt$maximum))
lw <- spdep::nb2listw(dnb, glist=glst, style="B")
SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
summary(SLX)
summary(impacts(SLX))
}
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
pslx0 <- predict(COL.SLX)
pslx1 <- predict(COL.SLX, newdata=COL.OLD, listw=lw)
all.equal(pslx0, pslx1)
#> [1] TRUE
COL.OLD1 <- COL.OLD
COL.OLD1$INC <- COL.OLD1$INC + 1
pslx2 <- predict(COL.SLX, newdata=COL.OLD1, listw=lw)
sum(coef(COL.SLX)[c(2,4)])
#> [1] -2.479902
mean(pslx2-pslx1)
#> [1] -2.479902