Archive

Posts Tagged ‘barchart’

Média e desvio-padrão de muitas variáveis separado por grupos

Gráfico de barras com as médias e desvios-padrões para variáveis resposta de um experimento avaliando cultivares de banana e meios de cultura em cultura de tecidos.

Em estatística descritiva é comum a síntese dos dados em medidas de posição e dispersão. As medidas mais frequentemente usadas são a média e o desvio-padrão. Não raramente, existe a necessidade de obter essas medidas separando as observações por alguma variável categórica do conjunto de dados. Nesse post apresento algumas formas de obter isso. Estou certo de que as soluções que apresento abaixo não cobrem 5% das disponíveis considerando a extensa quantidade de pacotes/funções do R ligadas à análise descritiva, portanto, considere isso apenas como um ponto de partida. Por fim, sintetizo em um gráfico os resultamos aplicados em um conjunto de dados reais. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando um conjunto de dados artificial

da <- expand.grid(A=1:3, B=1:2, r=1:15)
da <- cbind(da, as.data.frame(matrix(rnorm(nrow(da)*3), nrow(da))))
str(da)

#-----------------------------------------------------------------------------
# usando funções do reshape e plyr

require(reshape)
require(plyr)

ddply(da[,-c(1:3)], .(da$A, da$B), mean) # as médias
ddply(da[,-c(1:3)], .(da$A, da$B), sd)   # os desvios-padrões

db <- melt(da, id.vars=1:3, id.measures=4:6)
str(db)

#-----------------------------------------------------------------------------
# fazendo empilhamento das observações

ddply(db[,-3], .(A, B, variable), summarise, m=mean(value), sd=sd(value))

#-----------------------------------------------------------------------------
# usando funções do pacote doBy

library(doBy)

summaryBy(V1+V2+V3~A+B, data=da,
          FUN=function(x){ c(m=mean(x), sd=sd(x)) })

#-----------------------------------------------------------------------------
# usando funções do Hmisc

require(Hmisc)

with(da, bystats(cbind(V1, V2, V3), A, B,
                 fun=function(x){ c(m=colMeans(x), s=sd(x)) }))

#-----------------------------------------------------------------------------
# dados reais

bnn <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/banana.txt",
                  header=TRUE, sep="\t")
str(bnn)
bnn <- melt(bnn, id.measures=4:12)
bnn <- ddply(bnn, .(cult, meio, variable),
             summarise, m=mean(value), sd=sd(value))
str(bnn)

#-----------------------------------------------------------------------------
# gráfico de barras com desvio-padrão

require(lattice)
require(RColorBrewer)
display.brewer.all()
cols <- brewer.pal(nlevels(bnn$meio), "Set3")

png("f028.png", 500, 500)
trellis.par.set(superpose.polygon=list(col=cols),
                strip.background=list(col="gray90"))
barchart(m~cult|variable, groups=meio, data=bnn,
         sdv=bnn$sd, scales="free", layout=c(3,3),
         xlab="Cultivar", ylab=expression(bar(x)+sd),
         key=list(title="Meio de cultura", cex.title=1.2,
           columns=3, text=list(levels(bnn$meio)),
           rectangle=list(col=trellis.par.get("superpose.polygon")$col)),
         prepanel=function(y, sdv, subscripts=subscripts, ...){
           ly <- as.numeric(y+sdv[subscripts])
           list(ylim=range(y, ly, finite=TRUE))
         },
         panel=function(x, y, subscripts, groups, sdv, box.ratio, ...){
           panel.barchart(x, y, subscripts=subscripts,
                          groups=groups, box.ratio=box.ratio, ...)
           d <- 1/(nlevels(groups)+nlevels(groups)/box.ratio)
           g <- (as.numeric(groups[subscripts])-1); g <- (g-median(g))*d
           panel.arrows(as.numeric(x)+g, y, as.numeric(x)+g, y+sdv[subscripts],
                        angle=90, length=0.025)
         })
dev.off()

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

Como adicionar legendas em gráficos da lattice

Gráfico com bandas de confiança, legenda e assinatura de autor e versão.

No post Como fazer legendas em gráficos eu mostrei procedimentos para gráficos do pacote graphics. Nesse post vou apresentar procedimento para fazer legendas em gráficos da lattice, desenvolvido por Deepayan Sarkar. Iremos empregar as opções auto.key=, key=, colorkey=, e as funções draw.key() e draw.colorkey(). O pacote grid, desenvolvido por Paul Murrell será necessário para posicionarmos as legendas na página gráfica. Nos exemplos temos legendas para pontos, linhas, pontos e linhas, barras, escala de cores para gráficos de três dimensões, gráfico com duas legendas, posicionamento da legenda com o mouse e gráfico com legenda para representar o ajuste de um modelo de regressão com valores observados, preditos e intervalo de confiança.

#-----------------------------------------------------------------------------
# pacotes

require(lattice) # para fazer os gráficos
require(grid)    # para posicionar as legendas

#-----------------------------------------------------------------------------
# dados de indice agronômico em função da cultivar de milho e dose de N

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

#-----------------------------------------------------------------------------
# legenda para gráfico de dispersão

# gráfico com legenda padrão
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE, auto.key=TRUE)

# gráfico com legenda em colunas e alteração dos padrões dos pontos
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(columns=3),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# muda a posição da legenda
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(space="right", columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# coloca a legenda dentro do gráfico
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(corner=c(0.95,0.05), columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas

db <- with(da, aggregate(indice, list(cultivar=cultivar, dose=dose), mean))
col <- trellis.par.get("superpose.line")$col[1:nlevels(db$cultivar)]

# gráfico com legendas de linhas, nomes *antes* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         text=list(levels(db$cultivar)),
         lines=list(col=col)))

# gráfico com legendas de linhas, nomes *depois* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         lines=list(col=col),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas e pontos

xyplot(x~dose, groups=cultivar, data=db, type="o", pch=19,
       key=list(corner=c(0.95,0.05), columns=1,
         type="o", divide=1,
         lines=list(col=col, pch=19),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para pontos e reta de regressão

# ajuste de modelo e criação de objeto com valores observados e preditos
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
dc$indice <- predict(m0, newdata=dc)
dc$ind <- "predito"
da$ind <- "obsevado"
dc <- rbind(da[,-3], dc)
str(dc)

# legenda fora do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(space="top", columns=2, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

# legenda dentro do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(corner=c(0.95,0.05), columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

#-----------------------------------------------------------------------------
# dados de mortalidade de aves em função do tipo de resfriamento do aviário

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

#-----------------------------------------------------------------------------
# legenda para gráficos de barras

mb <- with(ma, aggregate(mortes, list(galpao=galpao, sistema=asper), sum))
str(mb)

col <- c("forestgreen", "palegreen")

# legenda fora da região gráfica
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda fora da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda dentro da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           corner=c(0.05,0.95), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com aspersão","sem aspersão"))))

#-----------------------------------------------------------------------------
# escolhendo a posição da legenda com o mouse

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col)
trellis.focus(name="page")
loc <- grid.locator(unit="npc") # clicar no gráfico
trellis.unfocus()

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de\nresfriamento", cex.title=1.1,
           x=loc$x, y=loc$y, corner=c(0.5,0.5), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com\naspersão","sem\naspersão"))))

#-----------------------------------------------------------------------------
# nota de rodapé no gráfico com assinatura, data e versão

# ajuste de modelo e criação de objeto com valores observados, preditos e IC
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
indice <- predict(m0, newdata=dc, interval="confidence")
dc <- cbind(dc, stack(data.frame(indice)))
names(dc)[3] <- "indice"
dc <- rbind(da[,-3], dc)

#png("f017.png", w=500, h=350)
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(3,1),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
       key=list(space="top", columns=3, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))))
pushViewport(viewport())
grid.text(label=paste("Walmes Zeviani --",
            format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
            R.version.string), rot=90,
          x=unit(1, "npc")-unit(2, "mm"),
          y=unit(0, "npc")+unit(2, "mm"),
          just=c("left", "bottom"),
          gp=gpar(cex=0.8, col="gray50"))
popViewport()
#dev.off()

#-----------------------------------------------------------------------------
# legenda para gráficos de 3 dimensões

x <- seq(-10,10,l=50)
dd <- expand.grid(x=x, y=x)
dd$z <- with(dd, 500+x+2*y-0.5*x*y-x^2-y^2)

# gráfico sem cores e sem legenda
wireframe(z~x+y, data=dd)

# gráfico com cores e legenda
wireframe(z~x+y, data=dd, drape=TRUE)

# remove a legenda
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE)

# coloca a legenda embaixo
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=list(space="bottom"))

# redimensionando a legenda
wireframe(z~x+y, data=dd, drape=TRUE,
          colorkey=list(width=1, height=0.7))

# gráfico de níveis de cor
levelplot(z~x+y, data=dd)

# com contornos e redimensiomento da legenda
levelplot(z~x+y, data=dd, contour=TRUE,
          colorkey=list(width=1, height=0.7))

#-----------------------------------------------------------------------------
# usando a draw.key(), usar fora da função gráfica ou dentro da panel()

xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(2,2),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico")
draw.key(list(space="top", columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))),
         draw=TRUE,
         vp=viewport(x=unit(0.75, "npc"), y=unit(0.65, "npc")))

#-----------------------------------------------------------------------------
# usando a draw.colorkey(), usar fora da função gráfica ou dentro da panel()

wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE,
          panel=function(...){
            panel.wireframe(...)
            draw.colorkey(list(space="bottom", at=seq(250,500,25),
                               height=0.8, width=1), draw=TRUE,
                          vp=viewport(x=unit(0.5,"npc"), y=unit(0.07,"npc")))
          })

#-----------------------------------------------------------------------------
# usando duas legendas

de <- expand.grid(w=gl(3,5,la="w"), z=gl(4,3,la="z"))
de$x <- with(de, rnorm(w, as.numeric(w), 0.5))
de$y <- with(de, rnorm(z, as.numeric(z), 0.5))
pch <- c(1,2,5,7); col=1:3

xy <- xyplot(y~x, data=de, groups=z:w,
             pch=rep(pch, e=3), col=rep(col, 4))

print(xy, position=c(0,0,0.8,1), more=FALSE)
draw.key(list(text=list(levels(de$z)), points=list(pch=pch, col=1),
         space="top", columns=1, title="Níveis de Z", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.6, "npc")))
draw.key(list(text=list(levels(de$w)), points=list(pch=15, col=col),
         space="top", columns=1, title="Níveis de W", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.4, "npc")))

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

Fatorial com desdobramento da interação, gráfico e tabela de médias

Rendimento de grãos em função dos níveis de potássio e umidade do solo. Letras minúsculas comparam níveis de potassio em um mesmo nível de umidade, as letras maiúsculas fazem o inverso.

O R não é um aplicativo específico de análise de experimentos. Mesmo assim, ele oferece funções que permitem a análise dos mais diversos tipos de experimentos, ficando a cargo do usuário saber usar corretamente as funções/pacotes que o R oferece.

Pode-se dizer que os experimentos fatoriais são os mais comuns na ciência. Nas ciências agrárias são amplamente usados. Têm-se uma tradição de representar a interação entre os fatores (qualitativos) empregando testes de médias para níveis de um fator fixando níveis dos demais.

Nessa ridícula, apresento procedimentos para obter o desdobramento de interação em um fatorial duplo. Os dados são do rendimento de grãos de soja em função dos fatores umidade do solo e dose de fertilizante potássico. O experimento foi instalado em casa de vegetação com as unidades experimentais (vasos) agrupadas em blocos que controlaram variações de luminosidade/umidade/ventilação.

Os procedimentos incluem o ajuste do modelo, o teste sobre os efeitos dos fatores pela análise de variância, a checagem das pressuposições do modelo, o fatiamento das somas de quadrados, o teste de médias. No final são apresentados procedimentos para representar os resultados em gráfico e tabela de dupla entrada. É importante dizer que os fatores deveriam ser analisados como quantitativos, onde espera-se uma interpretação mais interessante sobre o fenômeno. Entretanto, o método abaixo é ilustrativo do procedimento para fatores qualitativos. Em uma futura oportunidade analisarei dos dados considerando os fatores como quantitativos.

#-----------------------------------------------------------------------------
# importa e prepara os dados
 
rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt",
                   header=TRUE, sep="\t")
rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
str(rend)
 
#-----------------------------------------------------------------------------
# análise gráfica
 
require(lattice)
xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
 
#-----------------------------------------------------------------------------
# ajuste do modelo e faz checagem
 
m0 <- aov(rg~bloc+A*K, data=rend)
summary(m0)
par(mfrow=c(2,2)); plot(m0); layout(1)
 
#-----------------------------------------------------------------------------
# desdobrando somas de quadrados para a variação de K dentro de A e vice-versa

m1 <- aov(rg~bloc+K/A, data=rend)
desAinK <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
                  grep, x=names(coef(m1)[m1$assign==3]))
summary(m1, split=list("K:A"=desAinK))

m2 <- aov(rg~bloc+A/K, data=rend)
desKinA <- sapply(paste("A", levels(rend$A), sep=""), simplify=FALSE,
                  grep, x=names(coef(m2)[m2$assign==3]))
summary(m2, split=list("A:K"=desKinA))
 
#-----------------------------------------------------------------------------
# desdobrando a interação em testes de Tukey
 
require(agricolae)
require(plyr)

KinA <- sapply(levels(rend$A), simplify=FALSE,
               function(a){
                 with(subset(rend, A==a),
                      HSD.test(rg, K,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

KinA <- ldply(KinA, NULL)
KinA$M <- gsub(" ", "", KinA$M, fixed=TRUE)
KinA$trt <- as.factor(as.numeric(as.character(KinA$trt)))
str(KinA)

AinK <- sapply(levels(rend$K), simplify=FALSE,
               function(k){
                 with(subset(rend, K==k),
                      HSD.test(rg, A,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

AinK <- ldply(AinK, NULL)
AinK$M <- toupper(gsub(" ", "", AinK$M, fixed=TRUE))
AinK$trt <- as.factor(as.numeric(as.character(AinK$trt)))
str(AinK)

#-----------------------------------------------------------------------------
# combinando os data.frames para obter o gráfico e a tabela

KinA$comb <- paste(KinA$trt, KinA$.id)
AinK$comb <- paste(AinK$.id, AinK$trt)
desd <- merge(KinA, AinK, by=c("comb","comb"))
desd$let <- with(desd, paste(M.y, M.x, sep=""))
desd <- desd[,c("trt.x","trt.y","means.x","let")]
names(desd) <- c("K","A","means","let")
str(desd)

#-----------------------------------------------------------------------------
# gráfico das médias com as letras sobre as barras

#png("f005.png", w=500, h=250)
barchart(means~K|A, data=desd, horiz=FALSE, layout=c(3,1),
         ylab="Rendimento de grãos (g/vaso)", ylim=c(NA,45),
         xlab="Dose de potássio (mg/dm³)",
         key=list(text=list("Níveis de umidade do solo (%)")),
         sub=list("*letras minúsculas comparam médias em um mesmo nível de umidade.", cex=0.8, font=3),
         strip=strip.custom(bg="gray90"), col=colors()[51],
         panel=function(x, y, subscripts, ...){
           panel.barchart(x, y, subscripts=subscripts, ...)
           panel.text(x, y, label=desd[subscripts,"let"], pos=3)
         })
#dev.off()

#-----------------------------------------------------------------------------
# tabela de dupla entrada com as médias das combinações A e K

require(reshape)
desd$means.let <- with(desd, paste(means, let))
cast(desd, K~A, value="means.let")

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