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

Wed Aug 2 10:52:40 CEST 2017

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