[R-br] RES: Ajustes de modelos não-lineares em gases in vitro

Domingo Outubro 25 22:29:57 BRST 2015

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/attachment.html>