sarml: Spatial AR Maximum-Likelihood Estimation
In McSpatial: Nonparametric spatial data analysis

Description Usage Arguments Details Value See Also Examples

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

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

`form`	Model formula
`wmat`	The W matrix. If not specified, W will be calculated from the shape file. Default: W = NULL.
`shpfile`	Shape file. Needed unless (1) wy and eigvar are both provided, or (2) wmat and eigvar are provided
`wy`	The WY variable. Default: not specified; program attempts to calculate WY using wmat or shpfile.
`eigvar`	The vector of eigenvalues for W. Default: not provided. shpfile must be specified to calculate the eigenvalues within the sarml command.
`startrho`	A starting value for ρ. Default: startrho=0. Estimation will generally be faster if startrho = 0.
`print`	If print=F, no results are printed. Default: print=T.
`data`	A data frame containing the data. Default: use data in the current working directory

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.

`beta`	The estimated vector of coefficients, β.
`rho`	The estimated value of ρ.
`sig2`	The estimated error variance, σ^2.
`vmat`	The covariance matrix for (β, ρ^2).
`eigvar`	The vector of eigenvalues

makew

qreglwr

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)