[R-sig-Geo] ML estimate of neighbors radius

[Roger Bivand]

On Thu, 13 Dec 2012, Gregory Macfarlane wrote:

> I am trying to identify the maximum likelihood estimates in an SDM model 
> (a hedonic home price model, with observations being 5,000 individual 
> homes), $$y=\rho W y + X\beta + WX\lambda + \epsilon$$
>
> I have been able to successfully estimate parameters $\rho, \beta,\lambda$ using the `lagsarlm` function in the `spdep` package. In my case, $W$ is an inverse distance matrix limited to a radius $r$, or
> $$W_{ij} = \begin{cases}
>  1/d_{ij} & \text{for $d_{ij} \leq r$ miles}\\
>  0  & \text{for $d_{ij} > r$ miles}
>  \end{cases}$$
> I would like to identify $r$ simultaneously with $\rho, \beta,\lambda$. 
> I have even written a likelihood function to do this in R (making use of 
> the neighbors and lag functions in `spdep`)
>
> sdm.lik <- function(params, y, X, GeoPoints){
>  n <- length(y)
>  rho <- params[1]
>  radius <- params[2]
>  sigma2 <- params[3]
>  delta <- params[-c(1,2,3)]
>  # make W subject to radius
>  points.dnn <- dnearneigh(GeoPoints, 0, radius*5280)
>  dists   	 <- nbdists(points.dnn, GeoPoints)
>  dists.inv  <- lapply(dists, function(x) 1/x)
>  Wlist <- nb2listw(points.dnn, glist=dists.inv, style="W", zero.policy=TRUE)
>  W <- listw2mat(Wlist)
>  # Lag X variables
>  Z <- cbind(X[,1], lag.listw(Wlist, X[,1], zero.policy=TRUE))
>  for(i in 2:ncol(X)){
>    Z <- cbind(Z, X[,i], lag.listw(Wlist, X[,1], zero.policy=TRUE))
>  }
>  # Calculate log likelihood
>  e <- y - ((rho * W) %*% y) - (Z %*% delta)
>  lagmatrix <- diag(1,nrow=length(y), ncol=length(y)) - rho*W
>  logL <- ( -0.5*n*log(pi * sigma2) + log( det(lagmatrix) )
>                                           -(t(e) %*% e) / (2*sigma2) )
>  return(-logL)
> }
>
> Optimizing this function is far too slow to be practical, however. 
> Unsurprisingly, the highest cost operation is the log-determinant.

Yes, all of the components are present. See ?do_ldet in spdep. Also note 
that log(det()) will underflow as det() approaches zero. Do also note that 
W should be sparse - see the same help page and the code of the setup 
functions for insight. You certainly do not want to be handling dense 
matrices. Note that by endogenising W jointly, you make life very hard for 
yourself - you could take a coarse grid of radius values instead and 
optimize rho and beta for fixed radius, which would run very much faster. 
If you change W at each pass, none of the available methods for speeding 
up the logdet will benifit from running the setup function only once - you 
need now to set up the framework, such as solving the eigenproblem, each 
time. Probably "Matrix_J" will work best in your current setting.

Hope this helps,

Roger

>   1. Do you have tips for computing this determinant more quickly?
>   2. Does `lagsarlm` already have functionality to do this?
>   2.a). Can I hack it to make it do this?
>
> Greg Macfarlane
> Georgia Tech | Civil Engineering
>
>
> 	[[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