Archive

Posts Tagged ‘as.layer’

Equação junto ao gráfico na lattice

Reposta à salinidade para quatro espécies de craca. A regulação osmótica está ligada à capacidade adaptativa do organismo ao ambiente. O gráfico mostra os valores observados e a curva de regressão estimada acompanhada de bandas de confiança (95%). As equações em cada gráfico representam o modelo ajustado que descreve a curva.

Os polinômios são os modelos de regressão mais empregados em algumas áreas do conhecimento. Isso não quer dizer que eles sejam modelos bons. Na verdade, isso depende do cada um define como um bom modelo. Os polinômios são uma aproximação local. Isso significa que eles não devem ser usados para predições fora da amplitude observada das covariáveis (região experimental). Em outras palavras, não se deve fazer extrapolações, mas pode-se, com certos cuidados, fazer interpolações. Os polinômios são desprovidos de interpretação a partir do segundo grau. A partir do segundo grau, apenas o intercepto tem interpretação pois é o valor da função quando a covariável (x) é zero. Por isso, interpretações em geral são perdas de tempo. O que pode-se fazer são julgamentos quando ao sinal do coeficiente de maior grau que por exemplo, num modelo quadrático, indica se a curvatura é para cima (ou não). Outro ponto é que os polinômios apresentam problema de colinearidade. Para superar essa deficiência pode-se centrar a covariável ou usar polinômios ortogonais. Apesar do polinômio já ser sem interpretação, o uso de polinômios ortogonais é menos interpretável ainda pois os coeficientes não são taxas em relação à x, x², x³…, mas para uma transformação linear de x, essa que induz a ortogonalidade. Resumindo, os polinômios ortogonais são úteis pois minimizam o problema de colineridade. Então, do ponto de vista de reportar uma equação de predição, aquela que descreve o valor esperado de y (resposta) em função de x, é mais útil usar o polinômio cru. Isso porque, mesmo sem ter interpretação, toda pessoa sabe traçar a curva de um polinômio, é algo familiar, visto desde o ensino médio, e novamente, não que isso seja um ponto positivo para um modelo. De qualquer forma, o post de hoje ensina a ajustar polinômios e representar no gráfico os valores observados, os valores ajustados com intervalo de confiança e a equação dentro do gráfico. Os procedimentos abaixo eu desenvolvi para atender demandas de consultoria. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# leitura

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

xyplot(Osm~Salinidade|Espécie, data=osmo, type=c("p","smooth"))

#-----------------------------------------------------------------------------
# vou desconsiderar o fato dos dados serem heterocedasticos e ajustar um lm
# em outra matéria trataremos os desvios de pressupostos
# ajuste usando polinômios ortogonais

m0 <- aov(Osm~Espécie*poly(Salinidade,2), data=osmo)
par(mfrow=c(2,2)); plot(m0); layout(1)

qqmath(~residuals(m0)|osmo$Espécie)
xyplot(sqrt(abs(residuals(m0)))~fitted(m0)|osmo$Espécie, type=c("p","smooth"))
# desvio de pressupostos, trataremos em outra matéria do blog

summary(m0) # inferências sem validade pois existe desvio de pressupostos

#-----------------------------------------------------------------------------
# partição das somas de quadrados, declara Salinidade dentro de Espécie

m1 <- aov(Osm~Espécie/poly(Salinidade,2), data=osmo)

names(coef(m1))
attr(terms(m1), "term.labels")
ncoef <- 2 # número de efeitos por expécie=2 (linear e quadrático)
ngrup <- 4 # número de níveis de espécie=4
split.list <- split(matrix(1:(ncoef*ngrup), ngrup, ncoef), f=1:ngrup)
names(split.list) <- levels(osmo$Espécie)
split.list <- lapply(split.list,
                     function(x){ names(x) <- c("lin","qua"); x})
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

split.list <- as.list(1:(ncoef*ngrup))
names(split.list) <- c(outer(levels(osmo$Espécie), c("lin","qua"),
                             paste, sep=": "))
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

#-----------------------------------------------------------------------------
# prepara as equações para colocar no gráfico, função para polinômiais

coefs <- coef(lm(Osm~-1+Espécie/(Salinidade+I(Salinidade^2)), data=osmo))
coefs <- as.numeric(formatC(coefs, format="f", digits=3, flag="-"))

get.expression <- function(coef, digits){
  coef <- as.numeric(formatC(coef, format="f", digits=digits, flag="-"))
  postext <- c("","*x", paste("*x^", 2:(length(coef)-1), sep=""))
  sig <- c(ifelse(coef[1]<0, "-", ""), ifelse(coef[-1]<0, "-", "+"))
  txt <- paste(paste(sig, abs(coef), postext, sep=""), collapse=" ")
  eval(parse(text=paste("expression(y==", txt, ")", sep="")))
}

coefs <- split(matrix(coefs, ncol=3), f=1:4)
eqn <- lapply(coefs, get.expression, digits=3)
eqn

#-----------------------------------------------------------------------------
# nomes científicos em itálico nas tarjas

fl <- as.expression(lapply(levels(osmo$Espécie),
    function(x){ bquote(italic(.(x))) }))
fl

#-----------------------------------------------------------------------------
# predição

pred <- expand.grid(Espécie=levels(osmo$Espécie), Salinidade=seq(0,40,by=2))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))

#-----------------------------------------------------------------------------
# carrega funções que serão necessárias para fazer as bandas de confiança

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
grep("ciH", ls(), value=TRUE) # nome das funções (painéis) que vamos usar

#-----------------------------------------------------------------------------
# gráfico final

#png("f035.png", width=600, height=500)
xyplot(Osm~Salinidade|Espécie, data=osmo, col=1,
       xlab=expression(Salinidade~(ppm)),
       ylab=expression(Osmolalidade~(mOsm~kg^{-1}~de~H[2]*O)),
       strip=strip.custom(bg="gray90", factor.levels=fl))+
  as.layer(xyplot(fit~Salinidade|Espécie, data=pred, type="l", col=1,
                  ly=pred$lwr, uy=pred$upr,
                  prepanel=prepanel.ciH, panel=panel.ciH))
dimLay <- trellis.currentLayout()
for(col in 1:ncol(dimLay)){
  for(row in 1:nrow(dimLay)){
    trellis.focus("panel", col, row, highlight=FALSE)
    grid.text(label=eqn[[dimLay[row,col]]],
              x=0.03, y=1-0.03, just=c("left","top"))
  }
}
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, ,

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

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