[R-br] Auto-correlação residual utilizando a gnls

[Walmes Zeviani]

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 em ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120911/a9321861/attachment.html>