[R-sig-Geo] Computational problems with errorsarlm

[Javier García]

Hello everybody:

 

I am trying to estimate a spatial error model, but I am facing to several
problems

 

1)  Running errorsarlm function the following message appears:

 

Warning messages:

1: In errorsarlm(y ~ z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 +  :

  inversion of asymptotic covariance matrix failed for tol.solve = 1e-10 

  número de condición recíproco = 3.80991e-16 - using numerical Hessian.

2: In sqrt(fdHess[1, 1]) : Se han producido NaNs

 

 

Getting the following results:

 

Approximate (numerical Hessian) standard error: NaN

z-value: NaN, p-value: NA

Wald statistic: NaN, p-value: NA

 

 

This can be easily “solved” changing tol.solve from 1.0e-10 to, for example,
1.0e-20. Doing this I get  the following results

 

Asymptotic standard error: 14.053

z-value: -44.177, p-value: < 2.22e-16

Wald statistic: 1951.6, p-value: < 2.22e-16

 

2)  However, I have a more serious problema: the estimate of lambda does not
make any sense

 

Lambda: -620.82, LR test value: 333.5, p-value: < 2.22e-16

 

Any idea about what it is happening? I am using a big dataset with 2800
observations (houses), 14 variables, and the spatial weight matrix has been
constructed “by hand” with the inverse of the inter-areas distances .
Moreover, several observations belong to the same area (in total we have
only 10 areas). As the intra-area distance is unknown but cannot be
considered zero, I calculate it as 1/(0.1*dist_min), being dist_min the
distance between the corresponding area and the nearest one (idea borrowed
from Pattanayak and Butry (2005) “Spatial complementarity of forest and
farms: accounting for ecosystem services”, American Journal of Agricultural
Economics). Could be due to my particular spatial weight matrix? Any
alternative? 

 

 

Cheers

Javi

 




JAVIER GARCÍA 

 

Departamento de Economía Aplicada III (Econometría y Estadística)

Facultad de Economía y Empresa (Sección Sarriko)
Avda. Lehendakari Aguirre 83

48015 BILBAO
T.: +34 601 7126 F.: +34 601 3754
  www.ehu.es 

http://www.unibertsitate-hedakuntza.ehu.es/p268-content/es/contenidos/inform
acion/manual_id_corp/es_manual/images/firma_email_upv_euskampus_bilingue.gif

 

 

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] average bearing of animal movement data

[Alice Domalik]

Hi there, 

I have seabird tracking data and I have used both the packages 'trip' and 
'move' to calculate the max distance travelled (using the function 'homedist' 
in 'trip', and 'distanceSummary' in 'move'). 
I would also like to describe the bearing of each animal when it is at its 
maximum displacement from the colony. I am wondering if anyone knows any 
packages that can calculate this. Alternatively, if someone knows how I can 
extract the coordinates of the location of maximum displacement. 

thanks so much! 

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] testing nlme models for spatial correlation

[Javier Moreira]

Thanks,
I try the random effects to the level of TRATAMIENTO, and works fine.
For the cross or nested, what hapens is, if you dont use all the 3 randoms
the degrees of freedom arent correct.
Thanks for your time.



El 2 ago. 2017 5:52 a. m., "Thierry Onkelinx" 
escribió:

> Dear Javier,
>
> It looks like you have only one observation for each combination of BLOQUE/
> AMBIENTE/ TRATAMIENTO/ SUBMUESTREO. That is an observation level random
> effect (OLRE) which doesn't make sense with a Gaussian distribution. The
> ORLE and the residual would model the exact same thing. Model2 doesn't make
> sense.
>
> Two other problems: BLOQUE has only 3 levels and AMBIENTE and TRAITAMIENTO
> are both used in the fixed and the random part. See
> http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#
> should-i-treat-factor-xxx-as-fixed-or-random. AMBIENTE and TRAITAMIENTO
> are rather crossed than nested. See http://bbolker.github.io/
> mixedmodels-misc/glmmFAQ.html#nested-or-crossed and
> https://www.muscardinus.be/2017/07/lme4-random-effects/
>
> I can't reproduce the error you get on model4. The output seems
> reasonable, although it has the same problems as model2.
>
> Your dataset is suitable for the SPDE approach in INLA.
>
> I would recommend that you consult a (local) statistician.
>
> Best regards,
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to say
> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
> 2017-08-01 17:03 GMT+02:00 Javier Moreira :
>
>> ​​Sorry, the error is
>>
>> Error in corFactor.corSpatial(object) :
>>   NA/NaN/Inf in foreign function call (arg 1)
>> In addition: Warning messages:
>> 1: In min(unlist(attr(object, "covariate"))) :
>>   no non-missing arguments to min; returning Inf
>> 2: In min(unlist(attr(object, "covariate"))) :
>>   no non-missing arguments to min; returning Inf
>>
>> and i attach the data to this mail.
>> ​
>>  base_modelo3.csv
>> 
>> ​
>> thanks!
>>
>> 2017-08-01 11:42 GMT-03:00 Thierry Onkelinx :
>>
>>> Dear Javier,
>>>
>>> Your problem is hard to understand without a reproducible example. You
>>> only gives the code, not the data nor the error message.
>>>
>>> Best regards,
>>>
>>> ir. Thierry Onkelinx
>>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>>> and Forest
>>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>>> Kliniekstraat 25
>>> 1070 Anderlecht
>>> Belgium
>>>
>>> To call in the statistician after the experiment is done may be no more
>>> than asking him to perform a post-mortem examination: he may be able to say
>>> what the experiment died of. ~ Sir Ronald Aylmer Fisher
>>> The plural of anecdote is not data. ~ Roger Brinner
>>> The combination of some data and an aching desire for an answer does not
>>> ensure that a reasonable answer can be extracted from a given body of data.
>>> ~ John Tukey
>>>
>>> 2017-08-01 16:36 GMT+02:00 Javier Moreira :
>>>
 Thanks Thierry,
 I would check that info.
 Any ideas why if i choose the finest level of random effects as
 /SUBMUESTREO and run model 4 (with correlation) wont let me?
 If i undesrtand you wel, you adress more the second question i made, im
 all right?
 Thanks!


 El 1 ago. 2017 8:10 a. m., "Thierry Onkelinx" 
 escribió:

 Dear Javier,

 The correlation structure in nlme only works on the residuals within
 the finest level of random effect. Observations in different random effect
 are independent.

 Have a look at the INLA package (http://www.r-inla.org). INLA allows
 for correlated random effects. So you have spatial correlation among the
 random effects (instead of among residuals). INLA has options for
 correlations along a 2D regular grid, a neighbourhood graph or a mesh. If
 you want an book on this, I recommend Zuur et al (2017) Beginner's Guide to
 Spatial, Temporal and Spatial-Temporal Ecological Data Analysis with 
 R-INLA.

 Best regards,


 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for Nature
 and Forest
 team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
 Kliniekstraat 25
 1070 Anderlecht
 

Re: [R-sig-Geo] bivariate spatial correlation in R

[Rafael Pereira]

Roger,


Thank you for your response.


I recognize the data is not ideal and the analysis has limitations because
of the lack of information on population displacements that might have
occurred over the years. Nonetheless, there are plenty of data + literature
showing how the spatial distribution of income classes and land use
patterns is fairly stable over time in this city, particularly for a short
timescales like in this analysis. Having said this, I believe these two
questions (1) what socioeconomic classes have gained more accessibility?
and (2) “were wealthier areas in 2010 able to attract more changes to
accessibility?” in the end ask the same thing but with different phrasings,
though your phrasing (2) is more precise/correct.



On the more technical discussion, I see your point that spatial AND
temporal correlation in my data would make permutation bootstrap
inappropriate to generate significance levels, thus making bivariate
Moran’s I biased. Thank you very much for the clarifications! This has been
very helpful and I will have a closer look at which spatial regression
models are more appropriate for my analysis.


On a side note, do you think the function to calculate bivariate Moran’s I
is correct?  And could it be incorporated in the next version of spdep? If
so, please give credit to Rogério Barbosa, the researcher who proposed the
code in Stack Overflow.

best,

Rafael HM Pereira
http://urbandemographics.blogspot.com


On Mon, Jul 31, 2017 at 10:52 PM, Roger Bivand  wrote:

> Rafael,
>
> I'm sorry, but there is no way you can logically "analyze who benefits the
> recent changes in the transport system in terms of access to jobs" from the
> data you have.
>
> Even if you had aggregate household income data for 2014 and 2017 (not for
> 2010 only), you would not know whether wealthier families had not dispaced
> poorer families as accessibility improved. You need individual data, either
> survey or register, preferably panel, to show that changes in accessibility
> change the economic welfare of households controlling for movement of
> households. The timestamps on the data make any attempt to do this very
> risky; the real findings from a hypothetical surevey-based panel might be
> completely different, especially if poorer households were displaced (also
> indirectly, through rising house prices driven by improved accessibility).
> Gauging the welfare effects of transport investments is very hard to
> instrument.
>
> The closest I could get was to map deciles of the change in access (more
> negatives than positives) and compare the aspatial income distributions:
>
> library(spdep)
> library(rgdal)
> map <- readOGR(dsn=".", layer="test_map")
> library(classInt)
> cI <- classIntervals(map$diffaccess, n=10, style="quantile")
> library(RColorBrewer)
> ybrpal <- brewer.pal(6, "YlOrBr")
> fC <- findColours(cI, ybrpal)
> qtm(map, fill="diffaccess", fill.breaks=cI$brks, format="Europe2")
> map$faccess <- factor(findInterval(map$diffaccess, cI$brks,
>   all.inside=TRUE), labels=names(attr(fC, "table")))
> qtm(map, fill="diffaccess", fill.breaks=cI$brks, format="Europe2")
> acc_income <- split(map$income, map$faccess)
> do.call("rbind", lapply(acc_income, summary))
> dens <- lapply(acc_income, density)
> plot(1, ylab="", xlab="", type="n", xlim=c(-2000, 11000), ylim=c(0,
>   0.002))
> for (i in seq(along=dens)) lines(dens[[i]], col=i)
> legend("topright", legend=names(dens), col=1:length(dens), lty=1, bty="n")
>
> These density curves really do not suggest any clear relationship, other
> than that some areas with increased accessibility had higher incomes in
> 2010.
>
> You can examine the reverse relationship - were aggregate areas that were
> more wealthy in 2010 able to attract more changes to accessibility? The
> answer seems to be yes, they were able to do this:
>
> nb <- poly2nb(map)
> lw <- nb2listw(nb, style = "W", zero.policy = T)
> lm.morantest(lm(diffaccess ~ I(income/1000), map), lw)
> # SLX model
> summary(lmSLX(diffaccess ~ I(income/1000), map, lw))
> lm.morantest(lmSLX(diffaccess ~ I(income/1000), map, lw), lw)
> # Spatial Durbin error model - SDEM
> obj <- errorsarlm(diffaccess ~ I(income/1000), map, lw, etype="emixed")
> summary(impacts(obj))
> summary(impacts(lmSLX(diffaccess ~ I(income/1000), map, lw)))
> LR.sarlm(lmSLX(diffaccess ~ I(income/1000), map, lw), obj)
>
> It would be possible to run lm.morantest.sad() on the output of the SDEM
> model taking global spatial autocorrelation into account. If you need that,
> follow up in this thread.
>
> Bivariate Moran's I should not be used in this case, but could be used in
> other cases - use in change over time is troubling because randomisation
> will not be a good guide as t=1 and t=2 are subject to temporal as well as
> spatial autocorrelation, so you cannot use permutation bootstrap to find a
> usable measure of significance.
>
> Hope this clarifies, and thanks for the code.
>
> Roger
>
> On Sun, 30 Jul 2017, 

Re: [R-sig-Geo] testing nlme models for spatial correlation

[Thierry Onkelinx]

Dear Javier,

It looks like you have only one observation for each combination of BLOQUE/
AMBIENTE/ TRATAMIENTO/ SUBMUESTREO. That is an observation level random
effect (OLRE) which doesn't make sense with a Gaussian distribution. The
ORLE and the residual would model the exact same thing. Model2 doesn't make
sense.

Two other problems: BLOQUE has only 3 levels and AMBIENTE and TRAITAMIENTO
are both used in the fixed and the random part. See
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#should-i-treat-factor-xxx-as-fixed-or-random.
AMBIENTE and TRAITAMIENTO are rather crossed than nested. See
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#nested-or-crossed
and https://www.muscardinus.be/2017/07/lme4-random-effects/

I can't reproduce the error you get on model4. The output seems reasonable,
although it has the same problems as model2.

Your dataset is suitable for the SPDE approach in INLA.

I would recommend that you consult a (local) statistician.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2017-08-01 17:03 GMT+02:00 Javier Moreira :

> ​​Sorry, the error is
>
> Error in corFactor.corSpatial(object) :
>   NA/NaN/Inf in foreign function call (arg 1)
> In addition: Warning messages:
> 1: In min(unlist(attr(object, "covariate"))) :
>   no non-missing arguments to min; returning Inf
> 2: In min(unlist(attr(object, "covariate"))) :
>   no non-missing arguments to min; returning Inf
>
> and i attach the data to this mail.
> ​
>  base_modelo3.csv
> 
> ​
> thanks!
>
> 2017-08-01 11:42 GMT-03:00 Thierry Onkelinx :
>
>> Dear Javier,
>>
>> Your problem is hard to understand without a reproducible example. You
>> only gives the code, not the data nor the error message.
>>
>> Best regards,
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>> and Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>>
>> To call in the statistician after the experiment is done may be no more
>> than asking him to perform a post-mortem examination: he may be able to say
>> what the experiment died of. ~ Sir Ronald Aylmer Fisher
>> The plural of anecdote is not data. ~ Roger Brinner
>> The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of data.
>> ~ John Tukey
>>
>> 2017-08-01 16:36 GMT+02:00 Javier Moreira :
>>
>>> Thanks Thierry,
>>> I would check that info.
>>> Any ideas why if i choose the finest level of random effects as
>>> /SUBMUESTREO and run model 4 (with correlation) wont let me?
>>> If i undesrtand you wel, you adress more the second question i made, im
>>> all right?
>>> Thanks!
>>>
>>>
>>> El 1 ago. 2017 8:10 a. m., "Thierry Onkelinx" 
>>> escribió:
>>>
>>> Dear Javier,
>>>
>>> The correlation structure in nlme only works on the residuals within the
>>> finest level of random effect. Observations in different random effect are
>>> independent.
>>>
>>> Have a look at the INLA package (http://www.r-inla.org). INLA allows
>>> for correlated random effects. So you have spatial correlation among the
>>> random effects (instead of among residuals). INLA has options for
>>> correlations along a 2D regular grid, a neighbourhood graph or a mesh. If
>>> you want an book on this, I recommend Zuur et al (2017) Beginner's Guide to
>>> Spatial, Temporal and Spatial-Temporal Ecological Data Analysis with R-INLA.
>>>
>>> Best regards,
>>>
>>>
>>> ir. Thierry Onkelinx
>>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>>> and Forest
>>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>>> Kliniekstraat 25
>>> 1070 Anderlecht
>>> Belgium
>>>
>>> To call in the statistician after the experiment is done may be no more
>>> than asking him to perform a post-mortem examination: he may be able to say
>>> what the experiment died of. ~ Sir Ronald Aylmer Fisher
>>> The plural of anecdote is not data. ~ Roger Brinner
>>> The combination of some data and an aching desire for an answer does not
>>> ensure that a reasonable answer can be extracted from a given body of data.
>>> ~ John Tukey
>>>
>>>