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

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

> 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