R-sig-geo

Question about the functions calculating SAR and CAR models in spdep package

Classic

List

Threaded

1 message

Rachel Chocron

Question about the functions calculating SAR and CAR models in spdep package

Hello!

I sent the question below to Mr. Roger Bivand, and I want to share his
answer with the users of the R-sig-geo list. See his answers colored redinline:

I am using the spdep R package for spatial statistical analysis of
ecological data. At first I would like to thank you for implementing and
maintaining the package which helped me a lot in my work.
I would like to consult you with the following question:
I tried using the functions: errorsarlm and spautolm for calculation of SAR
and CAR models.
When the model formula I specified contained explaining variables with
different scales I received the error:* system is computationally singular*

The solution I found is setting tol.solve to a very low value: 2e-28.
After doing that, the functions worked and gave results.

Here are the commands I used:
*
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

remove(list=ls())
library(spdep)

BBS.data <-
read.table("J:/BBS/SpeciesResults/BBSHeterogeneityClimateSpeciesNicheWidth.csv",header=TRUE,sep=",")
attach(BBS.data)

pointsDF=SpatialPointsDataFrame(coords=cbind(POINT_X,POINT_Y),data=BBS.data)

neighbors=dnearneigh(pointsDF,0,200000)

weightsb<-nb2listw(neighbors,style="B",zero.policy=TRUE)

lm1<-lm(Average_species_num1996.2006~
sqr_X+mean_400m+sqr_mean_400m+summer_temperature+sqr_temperature+year_precipitation+sqr_precipitation+RANGE_400m+sqr_RANGE_400m,data=pointsDF)
summary(lm1)
AIC(lm1)

#When tol.solve = 2e-27, the function returns an error: system is
computationally singular
errorsar1<-errorsarlm(Average_species_num1996.2006~
sqr_X+mean_400m+sqr_mean_400m+summer_temperature+sqr_temperature+year_precipitation+sqr_precipitation+RANGE_400m+sqr_RANGE_400m,data=pointsDF,weightsb,na.action=na.omit,zero.policy=TRUE,tol.solve=2e-28)

#When tol.solve = 2e-28, the function returns an error: system is
computationally singular
errorsar1b<-spautolm(Average_species_num1996.2006~
sqr_X+mean_400m+sqr_mean_400m+summer_temperature+sqr_temperature+year_precipitation+sqr_precipitation+RANGE_400m+sqr_RANGE_400m,data=pointsDF,
listw=weightsb, family="SAR", method="full",
verbose=TRUE,zero.policy=TRUE,tol.solve=2e-29)

#When tol.solve = 2e-27, the function returns an error: system is
computationally singular
errorcar1 <- spautolm(Average_species_num1996.2006~
sqr_X+mean_400m+sqr_mean_400m+summer_temperature+sqr_temperature+year_precipitation+sqr_precipitation+RANGE_400m+sqr_RANGE_400m,data=pointsDF,
listw=weightsb, family="CAR", method="full",
verbose=TRUE,zero.policy=TRUE,tol.solve=2e-28)
summary(errorcar1)

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
*

Here is the output of summary() on the variables in the modes:

> summary(sqr_X) Min. 1st Qu. Median Mean 3rd Qu. Max.
7.089e+06 2.821e+11 1.095e+12 1.497e+12 2.339e+12 9.188e+12 >
summary(mean_400m) Min. 1st Qu. Median Mean 3rd Qu.
Max. NA's
3.354 194.400 340.400 603.200 780.500 3109.000 1 >
summary(sqr_mean_400m) Min. 1st Qu. Median Mean 3rd Qu. Max.
NA's
11 37780 115900 747400 609200 9668000 1 >
summary(summer_temperature) Min. 1st Qu. Median Mean 3rd Qu.
Max.
6.03 15.71 18.67 19.01 22.41 30.21 >
summary(sqr_temperature) Min. 1st Qu. Median Mean 3rd Qu.
Max.
36.36 246.80 348.70 379.70 502.10 912.70 >
summary(sqr_precipitation) Min. 1st Qu. Median Mean 3rd Qu.
Max.
8818 233500 822000 856700 1267000 4429000 > summary(RANGE_400m)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
6.0 65.0 123.5 214.0 268.8 1516.0 1 >
summary(sqr_RANGE_400m) Min. 1st Qu. Median Mean 3rd Qu. Max.
NA's
36 4225 15250 101600 72230 2298000 1 >
summary(year_precipitation) Min. 1st Qu. Median Mean 3rd Qu.
Max.
93.9 483.2 906.7 846.9 1126.0 2104.0

-------------------------------------------------------------------------------------------------------------------------------------

Here is the output of sessionInfo():

R version 2.15.1 (2012-06-22)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=Hebrew_Israel.1255 LC_CTYPE=Hebrew_Israel.1255
LC_MONETARY=Hebrew_Israel.1255 LC_NUMERIC=C
[5] LC_TIME=Hebrew_Israel.1255

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] spdep_0.5-46 coda_0.15-2 deldir_0.0-19 maptools_0.8-16
foreign_0.8-50 nlme_3.1-104 MASS_7.3-18 Matrix_1.0-6
[9] lattice_0.20-6 boot_1.3-4 sp_0.9-99

loaded via a namespace (and not attached):
[1] grid_2.15.1 tools_2.15.1

-----------------------------------------------------------------------------------------------------------------------------------------------

My question is: Are the resulting coefficients reliable?
One of the reasons my explaining variables are with different scales is
that I am inserting a linear and a square coefficient of the same
explanatory variable.
I didn't want to rescale the variables because I thought it would make the
interpretation of the results harder.

The help page explains that you should rescale. Rescale the base variable,
then the square will not be so problematic. If the variable is in metres,
use kilometres - there is no difference in interpretation.

The coefficients in errorsarlm() are not affected by the near-sigularity of
the matrix being inverted, but the standard error may be. The coefficients
in spautolm() may be affected.

> I compared the resulting coefficients to the coefficients calculated by
> regular OLS of the same models and they are identical in sign and similar
> in value.
>

This is what you would expect, because the spatial covariance only affects
the disturbance term in these models. The summary() of errorsarlm output
includes the option to show the results of a Hausman test, which checks
that the SAR error model and OLS model coefficients are similar in value -
if the test reports a significant result, the SAR error model is
misspecified. There is no similar test for CAR models, but the
interpretation is the same.

Hope this helps,

Roger

[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo