sacsarlm: Spatial simultaneous autoregressive SAC model estimation

Description Usage Arguments Details Value Control arguments Author(s) References See Also Examples

View source: R/sacsarlm.R

Description

Maximum likelihood estimation of spatial simultaneous autoregressive “SAC/SARAR” models of the form:

y = rho W1 y + X beta + u, u = lambda W2 u + e

where rho and lambda are found by nlminb or optim() first, and beta and other parameters by generalized least squares subsequently

Usage

1
2
3
4
sacsarlm(formula, data = list(), listw, listw2 = NULL, na.action, type="sac",
 method = "eigen", quiet = NULL, zero.policy = NULL, tol.solve = 1e-10,
 llprof=NULL, interval1=NULL, interval2=NULL, trs1=NULL, trs2=NULL,
 control = list())

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

listw2

a listw object created for example by nb2listw, if not given, set to the same spatial weights as the listw argument

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 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.

type

default "sac", may be set to "sacmixed" for the Manski model to include the spatially lagged independent variables added to X using listw; when "sacmixed", the lagged intercept is dropped for spatial weights style "W", that is row-standardised weights, but otherwise included

method

"eigen" (default) - the Jacobian is computed as the product of (1 - rho*eigenvalue) using eigenw, and "spam" or "Matrix" for strictly symmetric weights lists of styles "B" and "C", or made symmetric by similarity (Ord, 1975, Appendix C) if possible for styles "W" and "S", using code from the spam or Matrix packages to calculate the determinant; "LU" provides an alternative sparse matrix decomposition approach. In addition, there are "Chebyshev" and Monte Carlo "MC" approximate log-determinant methods.

quiet

default NULL, use !verbose global option value; if FALSE, reports function values during optimization.

zero.policy

default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE (default) assign NA - causing sacsarlm() to terminate with an error

tol.solve

the tolerance for detecting linear dependencies in the columns of matrices to be inverted - passed to solve() (default=1.0e-10). This may be used if necessary to extract coefficient standard errors (for instance lowering to 1e-12), but errors in solve() may constitute indications of poorly scaled variables: if the variables have scales differing much from the autoregressive coefficient, the values in this matrix may be very different in scale, and inverting such a matrix is analytically possible by definition, but numerically unstable; rescaling the RHS variables alleviates this better than setting tol.solve to a very small value

llprof

default NULL, can either be an integer, to divide the feasible ranges into a grid of points, or a two-column matrix of spatial coefficient values, at which to evaluate the likelihood function

trs1, trs2

default NULL, if given, vectors for each weights object of powered spatial weights matrix traces output by trW; when given, used in some Jacobian methods

interval1, interval2

default is NULL, search intervals for each weights object for autoregressive parameters

control

list of extra control arguments - see section below

Details

Because numerical optimisation is used to find the values of lambda and rho, care needs to be shown. It has been found that the surface of the 2D likelihood function often forms a “banana trench” from (low rho, high lambda) through (high rho, high lambda) to (high rho, low lambda) values. In addition, sometimes the banana has optima towards both ends, one local, the other global, and conseqently the choice of the starting point for the final optimization becomes crucial. The default approach is not to use just (0, 0) as a starting point, nor the (rho, lambda) values from gstsls, which lie in a central part of the “trench”, but either four values at (low rho, high lambda), (0, 0), (high rho, high lambda), and (high rho, low lambda), and to use the best of these start points for the final optimization. Optionally, nine points can be used spanning the whole (lower, upper) space.

Value

A list object of class sarlm

type

“sac”

rho

lag simultaneous autoregressive lag coefficient

lambda

error simultaneous autoregressive error coefficient

coefficients

GLS coefficient estimates

rest.se

asymptotic standard errors if ase=TRUE, otherwise approximate numeriacal Hessian-based values

ase

TRUE if method=eigen

LL

log likelihood value at computed optimum

s2

GLS residual variance

SSE

sum of squared GLS errors

parameters

number of parameters estimated

logLik_lm.model

Log likelihood of the non-spatial linear model

AIC_lm.model

AIC of the non-spatial linear model

method

the method used to calculate the Jacobian

call

the call used to create this object

residuals

GLS residuals

tarX

model matrix of the GLS model

tary

response of the GLS model

y

response of the linear model for rho=0

X

model matrix of the linear model for rho=0

opt

object returned from numerical optimisation

pars

starting parameter values for final optimization, either given or found by trial point evaluation

mxs

if default input pars, optimal objective function values at trial points

fitted.values

Difference between residuals and response variable

se.fit

Not used yet

rho.se

if ase=TRUE, the asymptotic standard error of rho, otherwise approximate numeriacal Hessian-based value

lambda.se

if ase=TRUE, the asymptotic standard error of lambda

resvar

the asymptotic coefficient covariance matrix for (s2, rho, lambda, B)

zero.policy

zero.policy for this model

aliased

the aliased explanatory variables (if any)

LLNullLlm

Log-likelihood of the null linear model

fdHess

the numerical Hessian-based coefficient covariance matrix for (rho, lambda, B) if computed

resvar

asymptotic coefficient covariance matrix

optimHess

FALSE

timings

processing timings

na.action

(possibly) named vector of excluded or omitted observations if non-default na.action argument used

Control arguments

fdHess:

default NULL, then set to (method != "eigen") internally; use fdHess to compute an approximate Hessian using finite differences when using sparse matrix methods with fdHess from nlme; used to make a coefficient covariance matrix when the number of observations is large; may be turned off to save resources if need be

LAPACK:

default FALSE; logical value passed to qr in the SSE log likelihood function

Imult:

default 2; used for preparing the Cholesky decompositions for updating in the Jacobian function

cheb_q:

default 5; highest power of the approximating polynomial for the Chebyshev approximation

MC_p:

default 16; number of random variates

MC_m:

default 30; number of products of random variates matrix and spatial weights matrix

super:

default FALSE using a simplicial decomposition for the sparse Cholesky decomposition, if TRUE, use a supernodal decomposition

opt_method:

default “nlminb”, may be set to “L-BFGS-B” to use box-constrained optimisation in optim

opt_control:

default list(), a control list to pass to nlminb or optim

pars:

default NULL, for which five trial starting values spanning the lower/upper range are tried and the best selected, starting values of rho and lambda

npars

default integer 4L, four trial points; if not default value, nine trial points

pre_eig1, pre_eig2

default NULL; may be used to pass pre-computed vectors of eigenvalues

Author(s)

Roger Bivand [email protected]

References

Anselin, L. 1988 Spatial econometrics: methods and models. (Dordrecht: Kluwer); LeSage J and RK Pace (2009) Introduction to Spatial Econometrics. CRC Press, Boca Raton.

Roger Bivand, Gianfranco Piras (2015). Comparing Implementations of Estimation Methods for Spatial Econometrics. Journal of Statistical Software, 63(18), 1-36. https://www.jstatsoft.org/v63/i18/.

Bivand, R. S., Hauke, J., and Kossowski, T. (2013). Computing the Jacobian in Gaussian spatial autoregressive models: An illustrated comparison of available methods. Geographical Analysis, 45(2), 150-179.

See Also

lm, lagsarlm, errorsarlm, summary.sarlm, eigenw, impacts.sarlm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
data(oldcol)
COL.sacW.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, 
 nb2listw(COL.nb, style="W"))
summary(COL.sacW.eig, correlation=TRUE)
W <- as(nb2listw(COL.nb, style="W"), "CsparseMatrix")
trMatc <- trW(W, type="mult")
summary(impacts(COL.sacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE)
COL.msacW.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, 
 nb2listw(COL.nb, style="W"), type="sacmixed")
summary(COL.msacW.eig, correlation=TRUE)
summary(impacts(COL.msacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE)

Example output

Loading required package: sp
Loading required package: Matrix

Call:
sacsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw(COL.nb, 
    style = "W"))

Residuals:
      Min        1Q    Median        3Q       Max 
-37.32081  -5.33662  -0.20219   6.59672  23.25604 

Type: sac 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept) 47.783766   9.902659  4.8253 1.398e-06
INC         -1.025894   0.326326 -3.1438  0.001668
HOVAL       -0.281651   0.090033 -3.1283  0.001758

Rho: 0.36807
Asymptotic standard error: 0.19668
    z-value: 1.8714, p-value: 0.061285
Lambda: 0.16668
Asymptotic standard error: 0.29661
    z-value: 0.56196, p-value: 0.57415

LR test value: 10.285, p-value: 0.0058432

Log likelihood: -182.2348 for sac model
ML residual variance (sigma squared): 95.604, (sigma: 9.7777)
Number of observations: 49 
Number of parameters estimated: 6 
AIC: 376.47, (AIC for lm: 382.75)

 Correlation of coefficients 
            sigma rho   lambda (Intercept) INC  
rho         -0.10                               
lambda       0.03 -0.76                         
(Intercept)  0.09 -0.90  0.68                   
INC         -0.04  0.38 -0.29  -0.59            
HOVAL       -0.01  0.09 -0.07  -0.22       -0.41

Impact measures (sac, trace):
          Direct   Indirect      Total
INC   -1.0632722 -0.5601501 -1.6234223
HOVAL -0.2919129 -0.1537847 -0.4456977
========================================================
Simulation results (asymptotic variance matrix):
========================================================
Simulated z-values:
         Direct   Indirect     Total
INC   -3.161211 -0.8247916 -1.807086
HOVAL -3.108806 -0.6425311 -1.329634

Simulated p-values:
      Direct    Indirect Total   
INC   0.0015711 0.40949  0.070749
HOVAL 0.0018784 0.52053  0.183639

Call:
sacsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw(COL.nb, 
    style = "W"), type = "sacmixed")

Residuals:
     Min       1Q   Median       3Q      Max 
-37.8045  -6.5244  -0.2207   5.9944  22.8691 

Type: sacmixed 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept) 50.92026   68.25721  0.7460 0.455664
INC         -0.95072    0.44033 -2.1591 0.030841
HOVAL       -0.28650    0.09994 -2.8667 0.004148
lag.INC     -0.69261    1.69113 -0.4096 0.682132
lag.HOVAL    0.20852    0.28702  0.7265 0.467546

Rho: 0.31557
Asymptotic standard error: 0.9458
    z-value: 0.33365, p-value: 0.73864
Lambda: 0.15415
Asymptotic standard error: 1.0643
    z-value: 0.14484, p-value: 0.88484

LR test value: 12.07, p-value: 0.016837

Log likelihood: -181.3422 for sacmixed model
ML residual variance (sigma squared): 93.149, (sigma: 9.6514)
Number of observations: 49 
Number of parameters estimated: 8 
AIC: 378.68, (AIC for lm: 382.75)

 Correlation of coefficients 
            sigma rho   lambda (Intercept) INC   HOVAL lag.INC
rho         -0.36                                             
lambda       0.34 -0.98                                       
(Intercept)  0.36 -1.00  0.98                                 
INC         -0.25  0.68 -0.67  -0.69                          
HOVAL        0.17 -0.46  0.45   0.44       -0.60              
lag.INC     -0.34  0.95 -0.94  -0.96        0.57 -0.43        
lag.HOVAL   -0.28  0.77 -0.76  -0.79        0.57 -0.40  0.62  

Impact measures (sacmixed, trace):
          Direct   Indirect      Total
INC   -1.0317003 -1.3693141 -2.4010144
HOVAL -0.2768608  0.1629265 -0.1139344
========================================================
Simulation results (asymptotic variance matrix):
========================================================
Simulated z-values:
         Direct   Indirect      Total
INC   -2.699179 -0.6204524 -0.9801006
HOVAL -2.434386  0.2417365 -0.1046863

Simulated p-values:
      Direct    Indirect Total  
INC   0.0069511 0.53496  0.32704
HOVAL 0.0149171 0.80898  0.91662

spdep documentation built on April 4, 2018, 1:03 a.m.