Spatial AR Maximum-Likelihood Estimation

Description

Estimates the model Y = ρ WY + X β + u by maximizing the log-likelihood function.

Usage

sarml(form,wmat=NULL,shpfile=NULL,wy=NULL,eigvar=NULL,startrho=NULL,
  print=TRUE,data=NULL)

Arguments

Details

The primary motivation for the sarml command is to provide a convenient way to estimate multiple spatial AR models without having to calculate the eigenvalues of W each time. Under the assumption that the errors, u, are independently and identically distributed normal, the log-likelihood function for the spatial AR model is

lnl = -n*log(pi)/2 - n*log(s2)/2 - ∑(u^2)/(2*s2) + ∑ log(1 - rho*eigvar)

where eigvar is the vector of eigenvalues of W. Though spdep provides a convenient and fast method for calculating the eigenvalues from a shape file, the calculation can nonetheless take a very long time for large data sets. The sarml command allows the user to input the vector of eigenvalues directly, which saves time when several models are estimated using the same W matrix. Unless a vector of eigenvalues is provided using the eigvar option, the eigenvalues are calculated from the shape file (provided using the shpfile option) using the spdep package.

Conditional on the value of ρ, the maximum likelihood estimate of β is simply the vector of coefficients from a regression of Y - ρ WY on X. The estimate of the error variance also has a closed form solution: σ^2 = ∑ u^2/n. Substituting these estimates into the log-likelihood functions leads to the following concentrated log-likelihood function:

lc = -n*(log(pi)+1)/2 - n*log(∑(u^2))/2 + ∑ log(1-rho*eigvar)

Working with the concentrated likelihood function reduces the optimization problem to a one-dimensional search for the value of ρ that maximizes lc. Unless a value is provided for startrho, the sarml procedure begins by using the optimize command to find the value of ρ that maximizes lc. This estimate of ρ (or the value provided by the startrho option) is then used to calculate the implied values of β and σ^2, and these values are used as starting values to maximize lnl using the nlm command.

The covariance matrix for the estimates of β and ρ, vmat, is the inverse of (1/σ^2)V. V has partitions V11 = X'X, V12 = X'WY, V21 = Y'W'X, and V22 = Y'W'WY + s2*∑ (eigvar/(1-rho*eigvar))^2.

Value

See Also

makew

qregspiv

qregbmat

qregsim1

qregsim2

qregcpar

qreglwr

Examples

library(spdep)

cmap <- readShapePoly(system.file("maps/CookCensusTracts.shp",
  package="McSpatial"))
cmap <- cmap[cmap$CHICAGO==1&cmap$CAREA!="O'Hare",]
samppop <- cmap$POPULATION>0&cmap$AREA>0
cmap <- cmap[samppop,]
cmap$lndens <- log(cmap$POPULATION/cmap$AREA)
lmat <- coordinates(cmap)
cmap$LONGITUDE <- lmat[,1]
cmap$LATITUDE  <- lmat[,2]
cmap$dcbd <- geodistance(longvar=cmap$LONGITUDE,latvar=cmap$LATITUDE,
  lotarget=-87.627800,latarget=41.881998)$dist

fit <- makew(shpfile=cmap,eigenvalues=TRUE)
wmat <- fit$wmat
eigvar <- fit$eigvar

# input w, calculate eigvar within sarml
fit <- sarml(lndens~dcbd,wmat=wmat,eigvar=eigvar,data=cmap)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(McSpatial)
Loading required package: lattice
Loading required package: locfit
locfit 1.5-9.1 	 2013-03-22
Loading required package: maptools
Loading required package: sp
Checking rgeos availability: TRUE
Loading required package: quantreg
Loading required package: SparseM

Attaching package: 'SparseM'

The following object is masked from 'package:base':

    backsolve

Loading required package: RANN
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/McSpatial/sarml.Rd_%03d_medium.png", width=480, height=480)
> ### Name: sarml
> ### Title: Spatial AR Maximum-Likelihood Estimation
> ### Aliases: sarml
> ### Keywords: Spatial AR Model Parametric Models
> 
> ### ** Examples
> 
> library(spdep)
Loading required package: Matrix
> 
> cmap <- readShapePoly(system.file("maps/CookCensusTracts.shp",
+   package="McSpatial"))
> cmap <- cmap[cmap$CHICAGO==1&cmap$CAREA!="O'Hare",]
> samppop <- cmap$POPULATION>0&cmap$AREA>0
> cmap <- cmap[samppop,]
> cmap$lndens <- log(cmap$POPULATION/cmap$AREA)
> lmat <- coordinates(cmap)
> cmap$LONGITUDE <- lmat[,1]
> cmap$LATITUDE  <- lmat[,2]
> cmap$dcbd <- geodistance(longvar=cmap$LONGITUDE,latvar=cmap$LATITUDE,
+   lotarget=-87.627800,latarget=41.881998)$dist
> 
> fit <- makew(shpfile=cmap,eigenvalues=TRUE)
> wmat <- fit$wmat
> eigvar <- fit$eigvar
> 
> # input w, calculate eigvar within sarml
> fit <- sarml(lndens~dcbd,wmat=wmat,eigvar=eigvar,data=cmap)
               Estimate Std. Error   z-value   Pr(>|z|)
(Intercept)  5.13908552 0.39796051 12.913557 0.00000000
dcbd        -0.01998899 0.01068662 -1.870468 0.06141883
rho          0.47351043 0.04051947 11.685997 0.00000000
sig2 =  0.9075362 
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>