Re: [R-sig-Geo] 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" 
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 :
>
>> ​​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
>> 
>> ​
>> thanks!
>>
>> 2017-08-01 11:42 GMT-03:00 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 :
>>>
 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" 
 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
 

Re: [R-sig-Geo] 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 :

> ​​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
> 
> ​
> thanks!
>
> 2017-08-01 11:42 GMT-03:00 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 :
>>
>>> 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" 
>>> 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
>>>
>>> 

Re: [R-sig-Geo] 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

​
thanks!

2017-08-01 11:42 GMT-03:00 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 :
>
>> 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" 
>> 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 :
>>
>>> 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)),
>>>

Re: [R-sig-Geo] 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 :

> 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" 
> 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 :
>
>> 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
>> 

Re: [R-sig-Geo] 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" 
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 :

> 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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] 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 :

> 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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo