testing nlme models for spatial correlation

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

testing nlme models for spatial correlation

Javier Moreira
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Thierry Onkelinx
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 <[hidden email]>:

> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Javier Moreira
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" <[hidden email]>
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 <[hidden email]>:

> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Thierry Onkelinx
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 <[hidden email]>:

> 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" <[hidden email]>
> 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 <[hidden email]>:
>
>> 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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Javier Moreira
​​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 <[hidden email]>:

> 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 <[hidden email]>:
>
>> 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" <[hidden email]>
>> 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 <[hidden email]>:
>>
>>> 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
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>


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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Thierry Onkelinx
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 <[hidden email]>:

> ​​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 <[hidden email]>:
>
>> 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 <[hidden email]>:
>>
>>> 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" <[hidden email]>
>>> 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 <[hidden email]>:
>>>
>>>> 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
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>>
>>>
>>
>
>
> --
> Javier Moreira de Souza
> Ingeniero Agrónomo
> 099 406 006
>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: testing nlme models for spatial correlation

Javier Moreira
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" <[hidden email]>
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 <[hidden email]>:
>
>> ​​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 <[hidden email]>:
>>
>>> 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 <[hidden email]>:
>>>
>>>> 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" <[hidden email]>
>>>> 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 <[hidden email]>:
>>>>
>>>>> 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
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>>
>>>>
>>>
>>
>>
>> --
>> Javier Moreira de Souza
>> Ingeniero Agrónomo
>> 099 406 006
>>
>>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...