Auto-correlação residual utilizando a gnls
Olá Pessoal, gostaria de contar com a experiência de vocês. Estou com problemas no momento de estimar o parâmetro de auto-correlação residual utilizando a gnls. Fiz o teste Durbin-Watson e como conclusão os resíduos são dependentes. Mas quando estimo os parâmetros de acordo com a fórmula abaixo a função me retorna que a estimativa do parâmetro de autocorrelação, "Phi1", é zero. Minha variável independente é daf e a variável dependente é mfs2. reg3=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf),correlation=corAR1(form=~daf)) Li a documentação da "nlme" mas não consigo compreender a diferença de declarar a correlação de outras formas, como por exemplo: reg2=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf),correlation=corAR1(form=~1|daf)) reg3=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf),correlation=corAR1()) Se possível alguém poderia esclarecer qual eu deveria utilizar e qual a diferença entre estas formas. DETALHE: Apenas na reg3 o parâmetro Phi1 é diferente de zero, mas utilizando este método ele dá diferente de zero mesmo quando os resíduos são independentes. Grato. Tales Jesus Fernandes Doutorando em Estatística UFLA Universidade Federal de Lavras
No primeiro (~daf) você diz que toda observação tem correlação com as demais e isso é uma função de daf. No segundo (~1|daf) não sei como ele interpretou mas se daf é uma variável não categórica isso não faz sentido, no caso de ser categórica, indivíduo por exemplo, e usar (~daf|individuo) a correlação depende da daf mas só dentro do mesmo indivíduo, observações entre indivíduos são independentes. Na terceira está sendo usado aquilo assumido na formula do seu groupedData. Se nada foi declarado, possivelmente a ordem dos registros tá sendo usada. Dificilmente você vai obter estimativas iguais aos valores usados na simulação ou valores sob sua hipótese. Mesmo simulando sob independência os seus dados raramente você terá estimativa zero. Nesse caso, o que deve acontecer com boa frequência é ser não significativo. Como exercício vale a pena você simular e estimar desses modelos para entender o que tá sendo assumido com cada formula. À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Prof. Walmes, muito obrigado. Eu não estava compreendendo a diferença entre as situações. Lendo a sua descrição acredito que a forma mais correta seria a reg3 (eu já desconfiava disso). Mas o problema é que quando utilizo esta função a assíntota superior do modelo fica sub-estimada. Dessa forma o ajuste fica muito ruim. Veja as figuras que são geradas no CRM abaixo. Saberia me dizer o que está acontecendo ? Ou me indicar a o melhor caminho para me livrar deste problema ? Segue o CMR: rm(list=ls()) library(car) library(nlme) library(qpcR) # para calcular o R2 daf=c(96,111,126,141,158,174,189,204,219,235,250,265,279,293) mfs2=c(0.014,0.027,0.033,0.189,0.439,0.602,0.653,0.671,0.775,0.814,0.911,1.012,1.059,1.005) # Logístico reg1=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf),correlation=corAR1(form=~daf)) reg3=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf),correlation=corAR1()) Rsq(reg1) Rsq(reg3) plot(daf,mfs2, main="Modelo Logístico") lines(daf,fitted(reg1),col="green") lines(daf,fitted(reg3),col="red") # Gompertz mod1=gnls(mfs2~(alfa*exp(-exp(beta-(gama*daf)))), start=c(alfa=1.5,beta=3,gama=0.02),weights=varExp(form=~daf),correlation=corAR1(form=~daf)) mod3=gnls(mfs2~(alfa*exp(-exp(beta-(gama*daf)))), start=c(alfa=1.5,beta=3,gama=0.02),weights=varExp(form=~daf),correlation=corAR1()) Rsq(mod1) Rsq(mod3) plot(daf,mfs2, main="Modelo Gompertz") lines(daf,fitted(mod1),col="green") lines(daf,fitted(mod3),col="red") Atenciosamente. Tales Jesus Fernandes Doutorando em Estatística UFLA Universidade Federal de Lavras
Tales, Tome cuido com o exagero. O R e outros aplicativos possuem muitas funções e flexibilidade na hora de criar modelos mas nós temos que ser cautelosos, entender o modelo que estamos propondo e críticos o suficiente para aceitar esses modelos. Você não tem tantos dados assim e está se propondo a modelar variância e correlação simultaneamente. Não tô dizendo que você esteja errado. Você vai ver abaixo que existe problemas de convergência, ou talvez de identificabilidade. Você tá propondo um modelo grande demais. No final eu coloco um duplo exponencial só pra você ver que existe mais de uma solução possível. rm(list=ls()) library(car) library(nlme) library(qpcR) # para calcular o R2 require(nls2) daf=c(96,111,126,141,158,174,189,204,219,235,250,265,279,293) diff(daf) # espaçamento não exatamente iguais, AR(1) assume equidistância mfs2=c(0.014,0.027,0.033,0.189,0.439,0.602,0.653,0.671,0.775,0.814,0.911,1.012,1.059,1.005) ordem <- seq_along(daf) plot(mfs2~daf) args(SSlogis) n0 <- nls(mfs2~SSlogis(daf, A, x0, S)) par(mfrow=c(2,2)); plot(as.lm(n0)); layout(1) # na minha opinião, esse tamanho de amostra é insuficiente para afirmar # qualquer padrão de variância e correlação... acf(residuals(n0)) # indicaria fracamente médias móveis pacf(residuals(n0)) # indica AR(2) # Logístico reg1=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(form=~daf)) summary(reg1) anova(reg1, n0) ## não diferiu no meu modelo mais simples (homocedástico e independente) reg3=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(form=~1)) summary(reg3) anova(reg3, n0) reg4 <- gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(value=0.01, form=~ordem)) summary(reg4) # veja que modelo 3 e 4 são iguais, pois corAR1() assume a ordem das observações no vetor anova(reg4, n0) # essa significancia é problema de mínimo local, veja que dando um chute mais # adequado para assintota, a estimativa de phi é outra reg5 <- gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1,beta=3,gama=0.015), # chute para alfa=1 weights=varExp(form=~daf), correlation=corAR1(form=~ordem)) summary(reg4) anova(reg5, n0) logLik(reg4) logLik(reg5) plot(daf,mfs2, main="Modelo Logístico") lines(daf,fitted(n0),col=1) # sem hetero e correl, not bad lines(daf,fitted(reg1),col=2) # usando AR1(~daf), not bad lines(daf,fitted(reg3),col=3) # usando AR1(~1), completamente sem sentido lines(daf,fitted(reg4),col=4) # usando AR1(~ordem), mesma coisa lines(daf,fitted(reg5),col=5) # usando AR1(~ordem), mínimo global # o padrão dos dados sugere duas logísticas, como se fosse crescimento de insetos # cada logística tem relação com as ecdises reg6 <- nls(mfs2~A1/(1+exp(-(daf-d1)/S1))+(A2-A1)/(1+exp(-(daf-d2)/S2)), start=list(A1=0.6, d1=170, S1=20, A2=1, d2=250, S2=15)) summary(reg6) confint.default(reg6) # aponta que S1 e S2 podem ser iguais, faz até sentido par(mfrow=c(2,2)); plot(as.lm(reg6)); layout(1) pred <- data.frame(daf=seq(90, 300, length=100)) pred$y <- predict(reg6, newdata=pred) plot(daf,mfs2, main="Modelo Logístico") lines(y~daf, data=pred) abline(h=coef(reg6)[c(1,4)], lty=3) abline(v=coef(reg6)[c(2,5)], lty=3) lapply(list(n0, reg3, reg4, reg5, reg6), logLik) # ultimo modelo com a maior verossimilhança À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Obrigado pelas colaborações, foram muito esclarecedoras! Walmes, no fim o que você chama de dupla exponencial, é apenas uma reparametrização no logístico duplo não é? Como ficaria para esta situação um gompertz duplo, ou dupla exponencial reparametrizando o modelo de gompertz? Att; Tales Jesus Fernandes Doutorando em Estatística UFLA Universidade Federal de Lavras ________________________________ De: Walmes Zeviani <walmeszeviani@gmail.com> Para: Tales Fernandes <talesest@yahoo.com.br> Cc: "r-br@listas.c3sl.ufpr.br" <r-br@listas.c3sl.ufpr.br> Enviadas: Terça-feira, 11 de Setembro de 2012 11:28 Assunto: Re: [R-br] Auto-correlação residual utilizando a gnls Tales, Tome cuido com o exagero. O R e outros aplicativos possuem muitas funções e flexibilidade na hora de criar modelos mas nós temos que ser cautelosos, entender o modelo que estamos propondo e críticos o suficiente para aceitar esses modelos. Você não tem tantos dados assim e está se propondo a modelar variância e correlação simultaneamente. Não tô dizendo que você esteja errado. Você vai ver abaixo que existe problemas de convergência, ou talvez de identificabilidade. Você tá propondo um modelo grande demais. No final eu coloco um duplo exponencial só pra você ver que existe mais de uma solução possível. rm(list=ls()) library(car) library(nlme) library(qpcR) # para calcular o R2 require(nls2) daf=c(96,111,126,141,158,174,189,204,219,235,250,265,279,293) diff(daf) # espaçamento não exatamente iguais, AR(1) assume equidistância mfs2=c(0.014,0.027,0.033,0.189,0.439,0.602,0.653,0.671,0.775,0.814,0.911,1.012,1.059,1.005) ordem <- seq_along(daf) plot(mfs2~daf) args(SSlogis) n0 <- nls(mfs2~SSlogis(daf, A, x0, S)) par(mfrow=c(2,2)); plot(as.lm(n0)); layout(1) # na minha opinião, esse tamanho de amostra é insuficiente para afirmar # qualquer padrão de variância e correlação... acf(residuals(n0)) # indicaria fracamente médias móveis pacf(residuals(n0)) # indica AR(2) # Logístico reg1=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(form=~daf)) summary(reg1) anova(reg1, n0) ## não diferiu no meu modelo mais simples (homocedástico e independente) reg3=gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(form=~1)) summary(reg3) anova(reg3, n0) reg4 <- gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1.5,beta=3,gama=0.015), weights=varExp(form=~daf), correlation=corAR1(value=0.01, form=~ordem)) summary(reg4) # veja que modelo 3 e 4 são iguais, pois corAR1() assume a ordem das observações no vetor anova(reg4, n0) # essa significancia é problema de mínimo local, veja que dando um chute mais # adequado para assintota, a estimativa de phi é outra reg5 <- gnls(mfs2~((alfa)/(1 + exp(beta - gama*daf))), start=c(alfa=1,beta=3,gama=0.015), # chute para alfa=1 weights=varExp(form=~daf), correlation=corAR1(form=~ordem)) summary(reg4) anova(reg5, n0) logLik(reg4) logLik(reg5) plot(daf,mfs2, main="Modelo Logístico") lines(daf,fitted(n0),col=1) # sem hetero e correl, not bad lines(daf,fitted(reg1),col=2) # usando AR1(~daf), not bad lines(daf,fitted(reg3),col=3) # usando AR1(~1), completamente sem sentido lines(daf,fitted(reg4),col=4) # usando AR1(~ordem), mesma coisa lines(daf,fitted(reg5),col=5) # usando AR1(~ordem), mínimo global # o padrão dos dados sugere duas logísticas, como se fosse crescimento de insetos # cada logística tem relação com as ecdises reg6 <- nls(mfs2~A1/(1+exp(-(daf-d1)/S1))+(A2-A1)/(1+exp(-(daf-d2)/S2)), start=list(A1=0.6, d1=170, S1=20, A2=1, d2=250, S2=15)) summary(reg6) confint.default(reg6) # aponta que S1 e S2 podem ser iguais, faz até sentido par(mfrow=c(2,2)); plot(as.lm(reg6)); layout(1) pred <- data.frame(daf=seq(90, 300, length=100)) pred$y <- predict(reg6, newdata=pred) plot(daf,mfs2, main="Modelo Logístico") lines(y~daf, data=pred) abline(h=coef(reg6)[c(1,4)], lty=3) abline(v=coef(reg6)[c(2,5)], lty=3) lapply(list(n0, reg3, reg4, reg5, reg6), logLik) # ultimo modelo com a maior verossimilhança À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
participantes (3)
-
andrebvs -
Tales Fernandes -
Walmes Zeviani