testing nlme models for spatial correlation

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