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 |
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 |
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 |
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 |
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 |
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 |
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 |
Free forum by Nabble | Edit this page |