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

[Adriele Giaretta Biase]

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.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151025/0d47ce79/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: gases_loog2pools_listaR.R
Tipo: application/octet-stream
Tamanho: 10165 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151025/0d47ce79/attachment.obj>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: gas.csv
Tipo: text/csv
Tamanho: 6651 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151025/0d47ce79/attachment.csv>