Uma alternativa para apresentar bandas de confiança em regressão

August 19, 2011
By

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

Valores preditos para índice agronômico de milho em função da dose de nitrogênio e cultivar. Regiões em cinza são os intervalos de confiança (95%) para os valores preditos.

Neste post novamente vou apresentar a análise dos dados de índice agronômico de milho que iniciei em outros posts (Predição em modelos de regressão com blocos, Como fazer regressão polinomial com diferentes graus). A única novidade será a representação gráfica, que ao invés de separar as curvas por nível de cultivar, vou apresentar no mesmo gráfico usando uma região escura para representar os intervalos de confiança para os valores preditos. Esse gráfico facilita perceber se os intervalos têm sobreposição ou não. No entanto, ter (ou não) sobreposição não aponta aceitação (ou rejetição) da hipótese de igualdade dos valores preditos. Hipóteses devem ser definidas antes da realização do experimento e avaliadas por testes adequados. Foi necessário desenvolver uma função complemento para os gráficos da lattice para produzir as regiões de confiança. No script abaixo eu reproduzo as análises e faço o gráfico que mencionei. Até a próxima ridícula!

#-----------------------------------------------------------------------------
# dados de índice agronômico de milho

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# ajuste do modelo. veja os posts anteriores a esse mencionados no texto.

levels(da$cultivar)
da$Pi.C <- with(da, ifelse(da$cultivar=="Pioneer-B815", dose^3, 0))
m3 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+Pi.C, data=da,
         contrast=list(bloco=contr.sum))
summary(m3)

#-----------------------------------------------------------------------------
# faz o gráfico final de predição pelo método matricial

# data.frame da predição
pred <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,10))
pred$Pi.C <- with(pred, ifelse(cultivar=="Pioneer-B815", dose^3, 0))

# matriz de predição sem colunas para os efeitos de blocos
X <- model.matrix(~cultivar*(dose+I(dose^2))+Pi.C, data=pred)

# vetor de estimativas sem os efeitos de blocos
b <- coef(m3)[-grep("bloco", names(coef(m3)))]

# margem de erro, produto do erro padrão e quantil
me <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
              se.fit=TRUE)$se*qt(0.975, df=df.residual(m3))

# data.frame com valores preditos e IC
pred$fit <- X%*%b
pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me
str(pred)

#-----------------------------------------------------------------------------
# paineis feitos para lattice para apresentar os intervalos de confiança como
# regiões de outra cor, úteis para sobrepor ICs

prepanel.ciH <- function(x, y, ly, uy, subscripts, ...){
  x <- as.numeric(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  list(ylim=range(uy, ly, finite=TRUE), xlim=range(x))
}

panel.ciH <- function(x, y, ly, uy, subscripts, ...){
  y <- as.numeric(y)
  x <- as.numeric(x)
  or <- order(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  panel.polygon(c(x[or], x[rev(or)]),
                c(ly[or], uy[rev(or)]),
                col=1, alpha=0.03, border=NA)
  panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
  panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
  panel.xyplot(x, y, subscripts=subscripts, ...)
}

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC pelo método matricial

require(lattice)
require(latticeExtra) # as.layer()
require(grid)

#png("f020.png", w=500, h=350)
with(pred,
     xyplot(fit~dose, groups=cultivar, type="l", lwd=2,
            strip=strip.custom(bg="gray90"),
            xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
            sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
            ly=lwr, uy=upr,
            prepanel=prepanel.ciH,
            panel=panel.superpose,
            panel.groups=panel.ciH,
            key=list(space="top", columns=1,
              points=list(pch=1), text=list("valores observados"),
              lines=list(col=1), text=list("valores preditos"),
              rectangle=list(col="gray90", border=NA), text=list("IC 95%"))
            ))+
  as.layer(xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE))
draw.key(list(columns=1,
              lines=list(lty=1, lwd=2,
                col=trellis.par.get("superpose.symbol")$col[1:3]),
              text=list(levels(da$cultivar))),
         draw=TRUE, vp=viewport(x=unit(0.75, "npc"), y=unit(0.25, "npc")))
#dev.off()

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

Tags: , , , , ,

Comments are closed.