Skip to content

Problem with singular matrix errosarlm

3 messages · Roger Bivand, Defriez, Emma

Original
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
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
+  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
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))
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
Roger Bivand
#
On Fri, 20 Dec 2013, Defriez, Emma wrote:

            

      
      
      
    

        
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.

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Defriez, Emma
    

    

      
      #
    

  


  

      
Roger,

Thank you so much for your reply. The problem with the weights was that I was rescaling the variables before calculating the weights matrix. As lat and long are variables in my model they too are rescaled therefore having nonsense values when it came to calculating neighbours. It is now a much sparser matrix and the model takes only seconds to run.

Fixing this also seems to have fixed the errors I was getting.

Thank you again

Emma


PhD student
Imperial College London

-----Original Message-----
From: Roger Bivand [mailto:Roger.Bivand at nhh.no] 
Sent: 20 December 2013 21:36
To: Defriez, Emma
Cc: 'r-sig-geo at r-project.org'
Subject: Re: [R-sig-Geo] Problem with singular matrix errosarlm

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

        

            

      
      
      
    

        
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.

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        
--
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