[R-sig-Geo] Spatial error model specification (spdep)

[Roger Bivand]

On Thu, 23 Feb 2012, Daniel Wiesmann wrote:

> Dear list
>
> I am trying to use spdep to estimate a spatial error model using data 
> aggregated by census divisions for an urban area.
>
> I have problems with the model specification with my data. When 
> estimating least squares (OLS) model, the reset-test indicates that the 
> correct specification includes quadratic terms.
>
> After estimating a model that passes the reset-test, I tested the 
> residuals for spatial autocorrelation using lm.morantest and found 
> evidence for it. Therefore I tried to run a spatial error model with the 
> same model formula.
>
> The problem I run into is that the function errorsarlm will not converge 
> when I include a quadratic term into the equation. My guess is that the 
> problem is related to the high correlation between x and x^2. The error 
> is the following:
>
>> Error in solve.default(asyvar, tol = tol.solve) :
>  system is computationally singular: reciprocal condition number = 1.92175e-14
>
> So my question is how to correctly specify the model? Or how to 
> successfully include quadratic terms?
>
> I have run into this problem with real data, but to illustrate the 
> problem, I added below an example that gives the same error with 
> artificial data.

Daniel:

First - this is a model posting in clarity and the provision of a working 
example. The only missing bit is the output of sessionInfo(), but that 
isn't important here. Congratulations! Helpers really appreciate postings 
like this!

The crucial information is given in ?errorsalm for the tol.solve= 
argument. If I rescale x by 0.1:

dat.errorsarlm <- errorsarlm(y ~ I(x/10) + I((x/10)^2), dat, nb)

works. In fact, the unscaled ML fit has converged, but the inversion 
needed to get the asymetric standard errors fails. So one can also try 
alternative fitting methods:

W <- as(as_dgRMatrix_listw(nb), "CsparseMatrix")
trMC <- trW(W, type="MC")
dat.errorsarlm <- errorsarlm(y ~ x + I(x^2), dat, nb, method="Matrix",
   trs=trMC)
summary(dat.errorsarlm)

gets near without the inversion problem, but arguably the SE of lambda is 
too large - there is a suspiciously large difference between the z-value 
test and the LR-test p-values.

Hope this helps

Roger

PS:

Being cheeky, this also works:

library(splines)
dat.errorsarlm <- errorsarlm(y ~ bs(x, 3), dat, nb)
summary(dat.errorsarlm)

and:

dat.errorsarlm <- errorsarlm(y ~ poly(x, 3), dat, nb)
summary(dat.errorsarlm)

See for instance recent papers by Daniel McMillen ...

>
> many thanks in advance,
>
> Daniel
>
> --------------------------------
>
> # Load libraries
> library(lmtest)
> library(spdep)
>
> # Create artificial data
> dat <- data.frame(x = 1:1000, x2 = (1:1000)^2)
> dat$y <- dat$x2 + rnorm(1000, 0, 10)
>
> # Tests suggest that quadratic model is the correct specification
> dat.lm <- lm(y ~ x, dat)
> reset(dat.lm)
>
> dat.lm <- lm(y ~ x + x2, dat)
> reset(dat.lm)
>
> # Create artificial neighbor matrix
> nb <- sample(c(1,0), 1e6, replace = T, prob = c(1, 100))
> nb <- matrix(nb, nrow = 1000)
> nb <- nb + t(nb)
> nb <- mat2listw(nb, style = 'W')
>
> # When I try to estimate spatial error model leads to an error
> dat.errorsarlm <- errorsarlm(y ~ x + x2, dat, nb)
>
>> Error in solve.default(asyvar, tol = tol.solve) :
>  system is computationally singular: reciprocal condition number = 1.92175e-14
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no