Estou no celular e não tenho como avaliar de forma mais completa. Mas realizei análise dados semelhante s a um tempo atrás e a estrutura corAR1(form=~1|grupo) me atendeu, onde grupo é o fator de agrupamento ( acredito ser amostra...) uma vez que as medidas tomadas em tempos diferentes tem em comum a mesma amostra. Assim q puder detalho melhor o q fiz, se alguém não esclarecer antes Em 25/10/2015 9:32 PM, "Adriele Giaretta Biase" <adrielegbiase@gmail.com> escreveu:
Olá,
eu gostaria de saber se é possível usar estruturas de correlação nos resíduos quando ajusto modelos não-lineares em mensurações de dados de produção de gases *in vitro* (caracterizado por mensurações acumuladas ao longo do tempo) e sem repetição.
Breve descrição do banco de dados de gases *in vitro*: o experimento foi conduzido em frascos de (125 mL) que contem fruído ruminal com feno de alfalfa. Esses frascos são acoplados aos sensores que medem pressão dos gases no interior dos frascos com o conversor digital que registra automaticamente a cada 5 minutos, durante o período de 48 horas. As mensurações dos gases são, portanto acumuladas ao longo do tempo. Tenho quatro compostos (amostras) diferentes que chamei de (y, y2, y3 e y4), em que cada uma amostra contem 577 mensurações ao longo do tempo (x) que é dado em horas.
Eu preciso ajustar modelos não-lineares para cada amostra (exponencial, logístico e logístico bicompartimental). O modelo logístico bicompartimental foi o que melhor se ajustou sem usar estrutruras de correlação, porém, a análise residual possuí uma forte dependência (visto que os dados são acumulados, ou seja, a mensuração no tempo (t+1) depende da mensuração anterior do tempo t.
Então, eu teria que entrar com estruturas de correlações no modelo logístico bicompartimental para melhorar a análise residual?
Eu tentei algumas estruturas usando o argumento “correlation”, da função gnls. No quadro abaixo segue o seguinte resumo das estruturas inseridas no modelo, com os respectivos resultados gerados pela FAC ou os erros retornados pelo software R.
modelo
correlation
Erros
FAC - residual
M1
sem estrutrura de correlação
Sem erro
Tendência : possivelmente independência e heterogeneidade
M2
AR1: corAR1()
Erro: step halving factor reduced below minimum in NLS step
-
M3
AR1: corAR1(form=~x)
Sem erro
Tendência : possivelmente independência e heterogeneidade
M4
AR2: corAR1(p=2, q=0)
Coefficient matrix not invertible
-
M5
ARIMA: corARMA(p=2,q=1)
step halving factor reduced below minimum in NLS step
-
M6
ARIMA: corARMA(p=2,q=2)
Coefficient matrix not invertible
-
M7
ARIMA: corARMA(p=3,q=1)
Coefficient matrix not invertible
-
M8
AR(): corAR1(form=~1)
Erro: step halving factor reduced below minimum in NLS step
-
Observação: Eu fiz a análise usando a metodologia de séries temporais (modelo linear não estacionário – com uma diferença d=1) e a FAC e FACP mostraram que os resíduos eram ruído branco, como nos modelos ARIMA(p=1, d=1, q=0), ARIMA(p=2, d=1, q=0) e ARIMA(p=3, d=1, q=0). Seria possível ajustar modelos não-lineares (logístico bicompartimental) usando estruturas de correlações como o ARIMA?
Os dados e o script seguem em anexo com todos os modelos que tentei ajustar (somente para uma amostra, ou seja, correspondente as colunas chamadas de y e x). Agradeço se puderem me dar uma orientação.
rm(list=ls(all=TRUE))
#pacotes necessários para a análise.
require(manipulate)
require(nlstools)
require(car)
require(nlme)
require(lmtest)
# modelos estudados
#a * (1 - exp(-b * (t - c))) # modelo exponencial
#a /(1+exp(2+4*b*(c-t))) # modelo logistico
#a/(1+exp(2+4*b*(c-t)))+d/(1+exp(2+4*e*(c-t))) # modelo logistico 2 pool
#-----------------------------------------------------------------------
setwd('E:\\autoregressivo_modelo')
gas01<-read.csv("gas.csv",sep = ",", header = TRUE)
gas01
head(gas01)
attach(gas01)
t=Time/60 # para deixar o tempo em horas
dim(gas01)
gas=cbind(gas01,t)
names(gas) <- c("Time","y","y2","y3", "y4", "x") # variáveis que quero ajustar y (variavel dependente gases) e x o tempo em horas
head(gas)
dim(gas)
attach(gas)
#===================================================================
# Ajuste usando o gnls
#Ajuste do modelo exponencial -lote I considerando os erros independentes
mod1 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas, start=c(a=1.1179, b=0.0562, c=-9.5923, d=0.367, e=0.0226))
summary(mod1) # quadro de estimativas
R2(residuals(mod1), gas$y)
residuosnI<-summary(mod1)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod1)
hist(residuosnI)
# Ajuste do modelo Logístico 2 pool considerando - AR()
mod2 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corAR1(),
start=c(a=0.78, b=0.06, c=-5.5923, d=0.367, e=0.0226))
summary(mod2) # quadro de estimativas
R2(residuals(mod2), gas$y)
residuosnI<-summary(mod2)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod2)
hist(residuosnI)
plot(x,mod2$residuals)
plot(mod2)
# ------------------------------------------------------------------------
# Ajuste do modelo Logístico 2 pool considerando - AR()
mod3 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corAR1(form=~x),
start=c(a=0.78, b=0.06, c=-5.5923, d=0.367, e=0.0226))
summary(mod3) # quadro de estimativas
R2(residuals(mod3), gas$y)
residuosnI<-summary(mod3)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod3)
hist(residuosnI)
plot(x,mod3$residuals)
plot( ACF(mod3, maxLag =24), alpha = 0.05 )
# ------------------------------------------------------------------------
# Ajuste do modelo Logístico 2 pool considerando - AR2()
mod4 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corARMA(p=2,q=0),
start=c(a=0.78, b=0.06, c=-5.5923, d=0.367, e=0.0226))
summary(mod4) # quadro de estimativas
R2(residuals(mod4), gas$y)
residuosnI<-summary(mod4)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod4)
hist(residuosnI)
plot( ACF(mod4, maxLag =24), alpha = 0.05 )
# ------------------------------------------------------------------------
# Ajuste do modelo Logístico 2 pool considerando - ARIMA(p=2,q=1)
mod5 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corARMA(p=2,q=1),
start=c(a=1.1179, b=0.0562, c=-9.5923, d=0.367, e=0.0226))
summary(mod5) # quadro de estimativas
R2(residuals(mod5), gas$y)
residuosnI<-summary(mod5)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod5)
hist(residuosnI)
plot( ACF(mod5, maxLag =24), alpha = 0.05 )
# ------------------------------------------------------------------------
# Ajuste do modelo Logístico 2 pool considerando - ARIMA(p=2,q=2)
mod6 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corARMA(p=2,q=2),
start=c(a=1.1179, b=0.0562, c=-9.5923, d=0.367, e=0.0226))
summary(mod6) # quadro de estimativas
R2(residuals(mod6), gas$y)
residuosnI<-summary(mod6)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod6)
hist(residuosnI)
plot( ACF(mod6, maxLag =24), alpha = 0.05 )
# ------------------------------------------------------------------------
# Ajuste do modelo Logístico 2 pool considerando - ARIMA(p=3,q=1)
mod7 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corARMA(p=3,q=1),
start=c(a=1.1179, b=0.0562, c=-9.5923, d=0.367, e=0.0226))
summary(mod7) # quadro de estimativas
R2(residuals(mod7), gas$y)
residuosnI<-summary(mod7)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod6)
hist(residuosnI)
plot( ACF(mod7, maxLag =24), alpha = 0.05 )
Agradeço desde já,
Adriele.
-- Adriele Giaretta Biase. Mestre em Estatística e Experimentação Agropecuária - UFLA. Doutoranda em Estatística e Experimentação Agronômica - ESALQ/ USP Contato: (19) 8861-0619.
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.