[R-sig-Geo] Problem with singular matrix errosarlm

[Roger Bivand]

On Fri, 20 Dec 2013, Defriez, Emma wrote:

> Hi,
>
> I am trying to fit a spatial regression model with errorsarlm in spdep, 
> but I have been running into some problems.
>
> I can get weights style W and B to work with method LU however not S or 
> C I get the following error (traceback below): Error in 
> solve.default(-(mat), tol.solve = tol.solve) : system is computationally 
> singular: reciprocal condition number = 0

My guess would be that the interval= argument should be given for S and C. 
In the LU case, the interval is set to -1,0.999 internally, and for each 
of these cases, the spatial coefficient may be being returned outside its 
domain, thus leading to failure in inverting the matrix returned by the 
numerical Hessian. Something similar may happen with Matrix - the Cholesky 
method.

>
> I have rescaled the RHS variables to mean 0 and standard deviation 1. I 
> have also tried removing all but two of the RHS variables which are 
> definitely not collinear to see if that was the problem but I get the 
> same error.
>
> Each data point represents a 1by1 degree cell for the world's oceans, 
> because I am using dnearneigh points at the equator have far fewer 
> neighbours than those nearer the poles. Could this unequal number of 
> neighbours be the reason for the error?
>
> Also this model has 4647 points and takes around 50minutes to run which 
> seems like a long time do you think this is due to the number of 
> neighbours? The full data set has 31609 so I imagine this is going to 
> take an extremely long time.
>

This is almost certainly caused by your choice of distance weights - you 
have 52% dense matrices, so lots of operations are being done on almost 
dense matrices. With 4647 points, you can use the eigen method too, which 
would give the domain interval exactly. You seem to have a mean number of 
links of 2461, so something rather odd is going on. The 1by1 degree cell 
for the world's land areas used in Bivand, Hauke & Kossowski (Geographical 
Analysis, 2013) with about 15K cells doesn't cause problems, but is very 
sparse.

Hope this clarifies,

Roger

> Code:
>
> chl<-(chl-mean(chl,na.rm=T))/sd(chl)
> sst<-(sst-mean(sst,na.rm=T))/sd(sst)
> chlm<-(chlm-mean(chlm,na.rm=T))/sd(chlm)
> sstm<-(sstm-mean(sstm,na.rm=T))/sd(sstm)
> depth<-(depth-mean(depth,na.rm=T))/sd(depth)
> lats<-(lats-mean(lats,na.rm=T))/sd(lats)
> longs<-(longs-mean(longs,na.rm=T))/sd(longs)
>
> ols<-lm(chl~sst+depth+chlm+sstm+lats+longs)
>
> coords<-cbind(longs,lats)
> coords<-as.matrix(coords)
>
> #Define neighbourhood here max is 200km
> nn<-dnearneigh(coords,0,200,longlat = TRUE)
>
> nn
> Neighbour list object:
> Number of regions: 4647
> Number of nonzero links: 11438816
> Percentage nonzero weights: 52.9707
> Average number of links: 2461.549
>
> #get inverse distances
> dsts <- nbdists(nn, coords)
> idw <- lapply(dsts, function(x) 1/(x))
>
> #Spatial weights
> nnweights<-nb2listw(nn, glist=idw, style='S')
>
>
> sem0_500<-errorsarlm(formula(ols), listw=nnweights,na.action=na.omit, method='LU',
> +               tol.solve=1e-16, control=list(returnHcov=FALSE))
> Error in solve.default(-(mat), tol.solve = tol.solve) :
>  system is computationally singular: reciprocal condition number = 0
>
>> traceback()
> 5: solve.default(-(mat), tol.solve = tol.solve)
> 4: solve(-(mat), tol.solve = tol.solve)
> 3: solve(-(mat), tol.solve = tol.solve)
> 2: getVmate(coefs, env, s2, trs, tol.solve = tol.solve, optim = con$optimHess,
>       optimM = con$optimHessMethod)
> 1: errorsarlm(formula(ols), listw = nnweights, na.action = na.omit,
>       method = "LU", tol.solve = 1e-16, control = list(returnHcov = FALSE))
>
> when I try to run it with different method I got a different error
>
>
>> sem0_500<-errorsarlm(formula(ols), listw=nnweights,na.action=na.omit, method='Matrix',
> +  tol.solve=1e-16, control=list(returnHcov=FALSE))
> Warning message: In .local(object, ...) :
>  Cholmod warning 'matrix not positive definite' at file ../Supernodal/t_cholmod_super_numeric.c, line 729
>
> Error in determinant(update(nChol, nW, 1/coef)) :
>  error in evaluating the argument 'x' in selecting a method for function 'determinant': Error in .local(object, ...) : CHOLMOD factorization was unsuccessful
>
>> traceback()
> 8: determinant(update(nChol, nW, 1/coef))
> 7: ifelse(coef > b, n * log(coef) + (.f * c(determinant(update(nChol,
>       nW, 1/coef))$modulus)), ifelse(coef < a, n * log(-(coef)) +
>       (.f * c(determinant(update(pChol, csrw, 1/(-coef)))$modulus)),
>       0))
> 6: Matrix_ldet(coef, env, which = which)
> 5: do_ldet(lambda, env)
> 4: f(arg, ...)
> 3: (function (arg)
>   -f(arg, ...))(0.930150734356483)
> 2: optimize(sar.error.f, interval = interval, maximum = TRUE, tol = con$tol.opt,
>       env = env)
> 1: errorsarlm(formula(ols), listw = nnweights, na.action = na.omit,
>       method = "Matrix", tol.solve = 1e-16, control = list(returnHcov = FALSE))
>
>> sessionInfo()
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C
> [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8
> [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8
> [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C
> [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] ncf_1.1-4    spdep_0.5-68 Matrix_1.1-0 sp_1.0-14
>
> loaded via a namespace (and not attached):
> [1] boot_1.3-9      coda_0.16-1     deldir_0.1-1    grid_3.0.2
> [5] lattice_0.20-24 LearnBayes_2.12 MASS_7.3-29     nlme_3.1-111
> [9] splines_3.0.2
>
> I feel I am doing something wrong when specifying the weights matrix but I can't work out what. Any help would be hugely appreciated.
>
>
> Thanks,
>
> Emma
>
> PhD student
> Imperial College London
>
> _______________________________________________
> 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, 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