Olá Mauro, no dataset estão as variáveis: Time e y. Time é a variável independente que deve ser convertida em horas da seguinte forma: x=Time/60 # para deixar o tempo em horas E y é a variável resposta (gases) e dependente do modelo. Adriele. Em 26 de outubro de 2015 12:00, <r-br-request@listas.c3sl.ufpr.br> escreveu:
Enviar submissões para a lista de discussão R-br para r-br@listas.c3sl.ufpr.br
Para se cadastrar ou descadastrar via WWW, visite o endereço https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br ou, via email, envie uma mensagem com a palavra 'help' no assunto ou corpo da mensagem para r-br-request@listas.c3sl.ufpr.br
Você poderá entrar em contato com a pessoa que gerencia a lista pelo endereço r-br-owner@listas.c3sl.ufpr.br
Quando responder, por favor edite sua linha Assunto assim ela será mais específica que "Re: Contents of R-br digest..."
Tópicos de Hoje:
1. RES: Ajustes de modelos não-lineares em gases in vitro (Mauro Sznelwar)
----------------------------------------------------------------------
Message: 1 Date: Sun, 25 Oct 2015 22:29:57 -0200 From: "Mauro Sznelwar" <sznelwar@uol.com.br> To: <r-br@listas.c3sl.ufpr.br> Subject: [R-br] RES: Ajustes de modelos não-lineares em gases in vitro Message-ID: <000001d10f85$722d3f50$5687bdf0$@com.br> Content-Type: text/plain; charset="utf-8"
Eu tentei rodar e não reconheceu a variável ‘x’, não está no seu data set.
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.
--- Este email foi escaneado pelo Avast antivírus. https://www.avast.com/antivirus -------------- Próxima Parte ---------- Um anexo em HTML foi limpo... URL: < http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151025/d63a5a47/attac...
------------------------------
Subject: Legenda do Digest
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
------------------------------
Fim da Digest R-br, volume 58, assunto 21 *****************************************
-- 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.