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

[Adriele Giaretta Biase]

Olá pessoal,


eu enviei uma dúvida sobre qual melhor estrutura de correlação residual
devo considerar num modelo não linear (logístico bicompartimental) para
ajuste de dados com produção de gases in vitro, que possui uma alta
dependência residual.

Eu organizei melhor o script e os dados, em anexo (novo anexo), em que os
dados que estão na planilha, a primeira coluna é o x variável independente
representando o tempo. E o y são as mensurações de gases.
Então a nova planilha está exatamente as informações no qual preciso obter
o ajuste do conjunto.

Tentei inserir algumas estruturas de correlação, mas não deu certo. Alguém
sabe como devo ajustar esse conjunto de dados (qual melhor estrutura de
correlação residual que posso considerar) usando modelos não lineares na
função gnls?



Adriele



#-----------------------------------------------------------------------
setwd('C:\\Users\\Adriele\\Desktop')
gas<-read.csv("gases_lista_R.csv",sep = ",", header = TRUE)
head(gas)
attach(gas)
#-----------------------------------------------------------------------

# Ajuste usando o gnls
#Ajuste do modelo exponencial -lote I considerando os erros independentes (
o que não ocorre!!)
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

# função para o cálculo do R2
R2 <- function(residuals, observed){
1-(sum(residuals)^2)/sum((observed-mean(observed))^2)
}

R2(residuals(mod1), gas$y)
residuosnI<-summary(mod1)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod1)
hist(residuosnI)
plot( ACF(mod1, maxLag =24), alpha = 0.05 )

# ------------------------------------------------------------------------

# 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)
plot( ACF(mod2, maxLag =24), alpha = 0.05 )
# ------------------------------------------------------------------------


# 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 )
# ------------------------------------------------------------------------

# Ajuste do modelo Logístico 2 pool considerando - AR() a ordem das
obsevações
mod8 <- gnls(y ~ a/(1+exp(2+4*b*(c-x)))+d/(1+exp(2+4*e*(c-x))),
data=gas,correlation=corAR1(form=~1),
start=c(a=0.78, b=0.06, c=-5.5923, d=0.367, e=0.0226))
summary(mod8) # quadro de estimativas
R2(residuals(mod8), gas$y)
residuosnI<-summary(mod8)$resid
shapiro.test(residuosnI)
ks.test(rnorm(577),residuosnI)
plot(mod6)
hist(residuosnI)

-- 
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/20151028/353d93de/attachment.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: gases_lista_R.csv
Tipo: text/csv
Tamanho: 8481 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151028/353d93de/attachment.csv>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: gases_loog2pools_listaR.R
Tipo: application/octet-stream
Tamanho: 8476 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20151028/353d93de/attachment.obj>