MLE estimation of spatial effects radius 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$$
The likelihood function of this model is
$$\ln(\mathcal{L}) = \dfrac{n\ln(\pi\sigma^2)}{2} + \ln|I - \rho W| + \dfrac{e'e}{2\sigma^2}
$$
with $e= y - \rho W y - X\beta - WX\lambda$  
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. 


*

*Do you have tips for computing this determinant more quickly?

*Does lagsarlm already have functionality to do this?

*Can I hack it to make it do this?

 A: This turns out to be a two-part answer
1. Log-Determinant
The spdep package does contain tools to efficiently calculate the log-determinant of a sparse weights matrix. The following lines of code drop the calculation time from 11 seconds to about 2. The code is not very clean because the methods are not intended for use by typical statistical consumers.
# Calculate log-determinant of (I- rho*W) - using methods in the spdep package
env <- new.env(parent=globalenv())
assign("listw", Wlist, envir=env)
assign("can.sim", can.sim, envir=env)
assign("similar", FALSE, envir=env)
assign("family", "SAR", envir=env)
assign("n", n, envir=env)
Matrix_setup(env, Imult=2, super=FALSE)
get("similar", envir=env)
logdet <- do_ldet(rho, env)

2. Maximum Likelihood Radius
Roger Bivand, the package author, explained to me that such a method is not really valid, and I discovered this myself when I couldn't find any traction on the likelihood function. What he instead recommended was to simply estimate the same model specification over a grid of radii and select the radius that provided the maximum point. The figure below shows how the likelihood value changes for higher sample sizes and different radii on my dataset.

