testing nlme models for spatial correlation

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

--
Javier Moreira de Souza
Ingeniero Agr?nomo
099 406 006