[R-sig-Geo] spatial durbin model with large dataset

[Roger Bivand]

On Wed, 23 Jun 2010, Mi Diao wrote:

> Thanks a lot for your help, Roger!
>
> I ran into another problem when I try to calibrate spatial lag and
> spatial durbin model with this dataset using spdep. The error messages
> that I got are:
>
> Error in solve.default(-(mat), tol.solve = tol.solve) :
>  Lapack routine dgesv: system is exactly singular

After an error, always report the output of traceback(), and always show 
the command in full that led to the error. In addition, always show the 
output of sessionInfo(). On examining the entrails, it looks (guesswork) 
as though the matrix output of the adjusted finite difference Hessian 
of the coefficients seems to be singular, which suggests that perhaps the 
X variables and the WX are strongly correlated. Try:

SD_fit <- lagsarlm(y ~ x, data=df, listw=lw, type="mixed",
   method="Matrix", tr=trMat, control=list(compiled_sse=FALSE))
SD_fit <- lagsarlm(y ~ x, data=df, listw=lw, type="mixed",
   method="Matrix", control=list(compiled_sse=FALSE))
SD_fit <- lagsarlm(y ~ x, data=df, listw=lw, type="mixed",
   method="Matrix", control=list(compiled_sse=TRUE))
SD_fit <- lagsarlm(y ~ x, data=df, listw=lw, type="mixed",
   method="Matrix", control=list(fdHess=FALSE))

and if any of them succeed, do

summary(SD_fit)

and see if any of the coefficients are aliased. The other warnings about 
the line optimisation reseting on NA suggests that there are issues with 
your spatial weights too - the fact that you can fit the model in GeoDa 
doesn't mean that the fit is OK - with odd weights, you can get odd 
results. Does the lagsarlm() fit with type="lag" give you the same results 
as GeoDa? Are your weights in spdep the same as in GeoDa? Are you maybe 
using inverse distance or similar odd weights?

If this doesn't make sense, please put your code, a copy of the GeoDa lag 
model fit, and save()-ed input data objects into a zip archive and post a 
link to them.

Hope this helps,

Roger

>
> In addition: There were 38 warnings (use warnings() to see them)
>> warnings()
> Warning messages:
> 1: In determinant(x, TRUE) : This version of the Matrix package returns
> |determinant(L)| instead of determinant(A), i.e., a
> *DIFFERENT* value.
> If still necessary, do change your code, following
> http://matrix.r-forge.r-project.org
>
> 2: In optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE,  ... :
>  NA/Inf replaced by maximum positive value
> ......
>
> I have been able to calibrate the spatial lag model with the same
> dataset using Geoda. Could anyone help me figure this out?
>
> Thanks a lot!
>
> Mi
>
>
> On Mon, Jun 21, 2010 at 5:13 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> On Sun, 20 Jun 2010, Mi Diao wrote:
>>
>>> Dear all,
>>>
>>> I am trying to calibrate a Spatial Durbin Model for a dataset with
>>> 52000 observations. Geoda can handle a dataset of this size, but it
>>> can only estimate spatial lag and spatial error model. Can the spdep
>>> package of R handle my dataset?
>>
>> Yes. At a recent demo, I ran the 25k ?house example on a small memory older
>> laptop without trouble (http://spatial.nhh.no/R/sea10.zip). With 50k and
>> many RHS variables, you may find a busy Windows 1Gb machine gasps a bit - if
>> so, shut down other programs (IE etc.) before starting R. Look at ?lagsarlm,
>> then ?impacts. Typically, you'll need something like:
>>
>> # lw is your weights object, probably from nb2listw()
>> # df is the data.frame with your data
>>
>> # first generate sparse matrix product traces
>> W <- as(as_dgRMatrix_lw(lw), "CsparseMatrix")
>> trMat <- trW(W, type="MC")
>>
>> SD_fit <- lagsarlm(y ~ x, data=df, listw=lw, type="mixed",
>>  method="Matrix", tr=trMat, control=list(compiled_sse=TRUE))
>> # type="mixed" fits a Spatial Durbin model;
>> # method could also be "Chebyshev" or "MC" for approximations,
>> # "Matrix" uses updating sparse Cholesky Jacobians, see ?do_ldet
>> # for details. "Matrix" needs symmetric neighbours, some other methods
>> # may not. Beware of k-nearest neighbour schemes, for the same reasons
>> # as in GeoDa
>>
>> summary(SD_fit)
>>
>> # now generate samples from the fitted model to assess the
>> # significance of the impacts
>> imp_SD_fit <- impacts(SD_fit, tr=trMat, R=2000)
>> summary(imp_SD_fit)
>> plot(imp_SD_fit)
>> summary(imp_SD_fit, zstats=TRUE, short=TRUE)
>>
>> Hope this helps,
>>
>> Roger
>>
>>>
>>> Thanks a lot for your kind help!
>>>
>>> Mi
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> Roger Bivand
>> Economic Geography Section, Department of Economics, Norwegian School of
>> Economics and Business Administration, 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
>>
>>
>
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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