[R-sig-ME] spatial correlation with nlme

[Thierry Onkelinx]

Dear Javier,

Please don't cross post unless you are asked to or in case your questions
isn't answered after a few days. I've answered in the original thread:
https://stat.ethz.ch/pipermail/r-sig-geo/2017-August/025877.html

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-02 3:43 GMT+02:00 Javier Moreira <javiermoreira op 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)
>
> 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
>
> ​
> 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")
>
>
> here its the data set
> ​
>  base_modelo3.csv
> <https://drive.google.com/file/d/0B6YImu-ZATe4bEFKTlAxci1PNW8/view?usp=
> drive_web>
> ​
> thanks!
>
> --
> Javier Moreira de Souza
> Ingeniero Agrónomo
> 099 406 006
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

	[[alternative HTML version deleted]]