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

[Defriez, Emma]

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

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.

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