[R-sig-Geo] testing nlme models for spatial correlation
Javier Moreira
javiermoreira at gmail.com
Wed Aug 2 16:00:49 CEST 2017
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" <thierry.onkelinx at inbo.be>
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 <javiermoreira at gmail.com>:
>
>> 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
>> <https://drive.google.com/file/d/0B6YImu-ZATe4bEFKTlAxci1PNW8/view?usp=drive_web>
>>
>> thanks!
>>
>> 2017-08-01 11:42 GMT-03:00 Thierry Onkelinx <thierry.onkelinx at inbo.be>:
>>
>>> 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 <javiermoreira at gmail.com>:
>>>
>>>> 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" <thierry.onkelinx at inbo.be>
>>>> 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
>>>>
>>>> 2017-08-01 11:31 GMT+02:00 Javier Moreira <javiermoreira at gmail.com>:
>>>>
>>>>> hi,
>>>>> im trying to generate different models to account for spatial
>>>>> correlation.
>>>>> Im using nlme package, and mixed models, in order to compare two
>>>>> models,
>>>>> one that doesnt include the spatial correlation and one that does.
>>>>>
>>>>> Its a nested design, one that has 4 leves,
>>>>> BLOQUE/ AMBIENTE/ TRATAMIENTO/ SUBMUESTREO
>>>>> Its a harvest set data, with multiple point of data/ treatment, so the
>>>>> last
>>>>> level account for another term in the error for "sub-muestreo".
>>>>>
>>>>> My first problem its, when i try to add de correlation term to the
>>>>> model, i
>>>>> cant, when the random effects are taken to the level /SUBMUESTREO, and
>>>>> i
>>>>> have to leave it to the level of TRATAMIENTO.
>>>>> When i do that, i have 2 differences between models, the term
>>>>> accounting
>>>>> for sub-muestreo, and the spatial correlation.
>>>>>
>>>>> #MODELO 2##
>>>>> attach(base_modelo3)
>>>>> str(base_modelo3)
>>>>> data1=base_modelo3
>>>>> str(data1)
>>>>> data1=groupedData(REND.~1|BLOQUE/AMBIENTE/TRATAMIENTO/SUBMUESTREO,
>>>>> data=data1, units=list(y="(ton/ha)"))
>>>>> data1$TRATAMIENTO =factor(data1$TRATAMIENTO)
>>>>> data1$BLOQUE =factor(data1$BLOQUE)
>>>>> data1$AMBIENTE =factor(data1$AMBIENTE)
>>>>>
>>>>> modelo2_MM<-lme(REND.~1+TRATAMIENTO*AMBIENTE,
>>>>> random=~1|BLOQUE/AMBIENTE/TRATAMIENTO/SUBMUESTREO,
>>>>> weights=varComb(varIdent(form=~1|TRATAMIENTO)),
>>>>> data=data1,
>>>>> control=lmeControl(niterEM=150,msMaxIter=200))
>>>>> summary(modelo2_MM)
>>>>> anova(modelo2_MM)
>>>>>
>>>>> ##MODELO 4##
>>>>>
>>>>> modelo4_MM<-lme(REND.~1+TRATAMIENTO*AMBIENTE,
>>>>> random=~1|BLOQUE/AMBIENTE/TRATAMIENTO,
>>>>> weights=varComb(varIdent(form=~1|TRATAMIENTO)),
>>>>> correlation=corExp(form=~X+Y,nugget=T),
>>>>> data=data1,
>>>>> control=lmeControl(niterEM=150,msMaxIter=200))
>>>>> summary(modelo4_MM)
>>>>> anova(modelo4_MM)
>>>>>
>>>>> My second problem is, that i need to get the specific parameter for the
>>>>> error term that belongs to the spatial correlation, in order to map
>>>>> it. For
>>>>> what i watch, what lme does is send it to the general error, and so,
>>>>> what i
>>>>> could do is make the differences between the residuals of these two
>>>>> models.
>>>>> so, its essetial to get the exact same model, except for the
>>>>> correlation
>>>>> structure.
>>>>> If anybody knows how to get the specific term of error accounting for
>>>>> the
>>>>> correlation, it would be wonderful.
>>>>>
>>>>> E24=residuals(modelo24_3,type = "response")
>>>>> E40=residuals(modelo4_MM,type = "response")
>>>>> EE=E24-E40
>>>>> post=data.frame(E24,E40,EE,data1$X,data1$Y)
>>>>>
>>>>> coordinates(post)<-c("data1.X","data1.Y")
>>>>> coor_post<-coordinates(post)
>>>>>
>>>>> bubble(post,"E24",main="residuos modelo 2")
>>>>> bubble(post,"E40",main="residuos modelo 4")
>>>>> bubble(post,"EE",main="Est.espacial removida por Modelo 4")
>>>>>
>>>>> thanks!
>>>>>
>>>>> --
>>>>> Javier Moreira de Souza
>>>>> Ingeniero Agrónomo
>>>>> 099 406 006
>>>>>
>>>>> [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>>
>>>>
>>>
>>
>>
>> --
>> Javier Moreira de Souza
>> Ingeniero Agrónomo
>> 099 406 006
>>
>>
>>
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list