Arquivo

Posts Tagged ‘latticeExtra’

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

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

#-----------------------------------------------------------------------------
Anúncios

Gráfico com escala logarítmica em notação de potência

Gráfico (verde) em escala logarítmica com notação de potência.

O pacote latticeExtra traz funções que permitem aprimorar as funções gráficas disponíveis no pacote lattice, ambos desenvolvidos por Deepayan Sarkar. Recentemente precisei confeccionar gráficos em escala logarítmica para os quais eu tinha um enorme roteiro de funções construídas à custo de muita dedicação e tempo. Felizmente, hoje não preciso mais reaprender as implementações que fiz porque a latticeExtra fornece funções para escrever as escalas dos eixos em notação de potência.

As funções disponibilizadas são: xscale.components.logpower, xscale.components.fractions, xscale.components.log10ticks, xscale.components.log10.3 e
xscale.components.subticks, todas elas na versão para o eixo y também. Essas funções permitem, por exemplo, eixos escritos com potência e fração e os traços marcas (tick) para datas considerando o diferente número de dias de cada mês. Tudo isso você pode conferir rodando os exemplos das funções. Até a próxima ridícula.

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

require(lattice)

data(Earthquake, package="MEMSS")

p1 <-
  xyplot(accel~distance, data=Earthquake, pch=19, lwd=2,
         prepanel=prepanel.loess, col="blue",
         type=c("p", "g", "smooth"), scales=list(log=2),
         xlab="Distância do epicentro (km)",
         ylab="Aceleração máxima horizontal (g)")

require(latticeExtra)

p2 <-
  xyplot(accel~distance, data=Earthquake, pch=19, lwd=2,
         prepanel=prepanel.loess, col="forestgreen",
         type=c("p", "g", "smooth"), scales=list(log=2),
         xlab="Distância do epicentro (km)",
         ylab="Aceleração máxima horizontal (g)",
         xscale.components=xscale.components.logpower,
         yscale.components=yscale.components.logpower)

#png("f011.png", w=500, h=300)
print(p1, split=c(1,1,2,1), more=TRUE)
print(p2, split=c(2,1,2,1))
#dev.off()

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