Análise de resíduos em regressão não linear

May 29, 2011
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

Conjunto de gráficos para análise dos resíduos de um modelo de regressão não linear.

A análise dos resíduos de um modelo é feita para verificar a plausividade das pressuposições envolvidas. Os modelos linerares de regressão clássico, ou seja, aqueles em que as observações são realizações independentes (independência) e apresentam a mesma dispersão (homogeneidade de variância), podem ser ajustados no R com a função lm() e aov(). Para objetos dessas classes existe um método plot.lm() que apresenta os gráficos de análise de resíduo. Porém, modelos de regressão não linear podem ser ajustados com a função nls() que não possui um método para análise de resíduos.

Para aplicação de inferência (testes de hipótese, intervalos de confiança, etc), o modelo não linear é aproximado linearmente. Com isso eu quero dizer que podemos usar as mesmas técnicas de análise de resíduos para modelos de regressão não linear. A matriz gradiente do modelo (de derivadas da função em relação ao vetor de parâmetros) pode ser usada dentro da função lm() e com isso podemos obter facilmente os gráficos de análise de resíduos.

Com isso ainda podemos obter a análise de variância para o modelo de regressão não linear, fazer a partição ortogonal da soma de quadrados total em devido ao modelo de regressão e devido aos desvios de regressão. No então, esse não deveria ser o quadro apresentado se compararmos o que obtemos com modelo de regressão linear. No modelo linear particionamos a soma de quadrados total corrigida para o intercepto. O nosso modelo de regressão não linear se torna um modelo de intercepto se V e D forem zero. Portando, o modelo nulo seria aquele em que estimamos apenas A e que podemos comparar com o modelo completo por meio da função/método anova.nls(). Esse é o teste de hipótese para o modelo de regressão não linear. Vale lembrar que alguns modelos não lineares não podem ser reduzidos a modelos de intercepto via restrição paramétrica. Para esses modelos deve-se ter cuidados ao usar o R², inclusive.

Abaixo apresento o procedimento para ajuste de um modelo não linear, extração da matriz gradiente e obtenção dos gráficos de resíduos. No final eu fiz uma função (R2()) para facilitar a obtenção do quadro de análise de variãncia e R². Não é porque você tem uma função para isso que vai sair aplicando ela em qualquer modelo não linear, ok? Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de postássio liberado em função do tempo

klib <- data.frame(k=c(51.03, 57.76, 26.60, 60.65, 87.07, 64.67,
                     91.28, 105.22, 72.74, 81.88, 97.62, 90.14,
                     89.88, 113.22, 90.91, 115.39, 112.63, 87.51,
                     104.69, 120.58, 114.32, 130.07, 117.65, 111.69,
                     128.54, 126.88, 127.00, 134.17, 149.66, 118.25,
                     132.67, 154.48, 129.11, 151.83, 147.66, 127.30),
                   t=rep(c(15, 30, 45, 60, 75, 90,
                     120, 150, 180, 210, 240, 270), each=3))

#-----------------------------------------------------------------------------
# ajustando o modelo de regressão não linear aos dados

n0 <- nls(k~A*t/(V+t)+D*t, data=klib, start=list(A=90, V=15, D=0.21))
summary(n0)

#-----------------------------------------------------------------------------
# extraindo a matriz gradiente avaliada nas estimativas dos parâmetros

F <- attr(n0$m$fitted(), "gradient")
F

#-----------------------------------------------------------------------------
# passando a matriz gradiente para a lm(), importante remover intercepto

m0 <- lm(k~-1+F, data=klib)

#-----------------------------------------------------------------------------
# gráfico de análise dos resíduos

#png("f007.png", w=500, h=500);
par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1))
plot(m0)
mtext("Análise de resíduos para modelo de regressão não linear",
      outer=TRUE, line=-2, cex=1.4)
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# veja que o ajuste é o mesmo pelas medidas abaixo

cbind(fitted(m0), fitted(n0))       # valores ajustados
cbind(residuals(m0), residuals(n0)) # valores preditos
c(summary(m0)$sig, summary(n0)$sig) # estimativa do desvio padrão residual
vcov(m0); vcov(n0)                  # matriz de covariância das estimativas

#-----------------------------------------------------------------------------
# quadro de anova com SQ de regressão e SQ de resíduo

anova(m0) # partição da soma de quadrados total
anova(n0, lm(k~1, klib)) # SQ do modelo não linear corrigido para intercepto

#-----------------------------------------------------------------------------
# função que retorna a anova e R2 para modelos de regressão não linear

R2 <- function(nls.obj){
  da <- eval(nls.obj$data)
  resp.name <- all.vars(summary(nls.obj)$formula)[1]
  form <- paste(resp.name, "~1", sep="")
  m0 <- lm(form, da)
  an <- anova(nls.obj, m0)
  sqn <- deviance(nls.obj)
  sqe <- deviance(m0)
  r2 <- 1-(sqn/sqe)
  aov <- data.frame(fv=c("regression","residuals"),
                    gl=c(-an$Df[2],an$Res.Df[1]),
                    sq=c(-an$Sum[2],an$Res.Sum[1]))
  aov$qm <- aov$sq/aov$gl
  aov$F <- c(aov$qm[1]/aov$qm[2], NA)
  aov$"Pr(>F)" <- c(1-pf(aov$F[1], df1=aov$gl[1], df2=aov$gl[2]), NA)
  names(aov) <- c(" ","Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
  return(list(anova=aov, R2=r2))
}

R2(n0)

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

Tags: , , , , , , , , , ,

Comments are closed.