Gráfico de barras com erro padrão

São muitas as opções para se fazer um gráfico de barras com a indicação do erro padrão (ou outra medida de dispersão) no R. Com uma rápida busca na rede pode-se ter dezenas de sugestões. Eu resolvi contribuir com mais uma.

Meu objetivo era representar a variação da densidade de algumas espécies medida em períodos noturnos e diurnos.

A partir de três tabelas, uma contendo as médias, outra o erro padrão e uma terceira com o número de observações, utilizo os comandos barplot, points, arrows e text para gerar o gráfico desejado.

Os dados das tabelas são copiados para o R pelo comando read.delim (clique aqui para mais detalhes) e transformados em matrizes pelo comando as.matriz.

Médias

	SP1	SP2	SP3	SP4
D	10	12	15	20
N	20	15	12	10

Erros padrões

	SP1	SP2	SP3	SP4
D	1	2	3	4
N	3	2	1	5

Número de observações

	SP1	SP2	SP3	SP4
D	30	40	35	10
N	20	50	20	27

# copia os dados das tabelas para o R
means <- as.matrix(read.delim("clipboard",row.names=1))
se <- as.matrix(read.delim("clipboard",row.names=1))
n <- as.matrix(read.delim("clipboard",row.names=1))
 
# calcula a amplitude da variação do erro
se.sup<-means+se
se.inf<-means-se

# desenha o gráfico
bp<-barplot(means,beside=T,ylim=c(0,max(se.sup*1.15)),
ylab="densidade",legend.text=c("Diurno","Noturno"),
args.legend=list(x = "topleft", bty="n"))
points(bp,means,pch=19)
arrows(bp,se.sup,bp,se.inf, code=3,angle=90,length=0.05)
text(bp,se.sup+1,n)

 

Análise de Variância e teste Kruskal–Wallis

É frequente termos a necessidade de verificar se as diferenças de uma variável medida em diversos grupos são significativas. Por exemplo, podemos querer testar a diferença dos comprimentos de peixes capturados em diferentes locais  ou da potência do motor de embarcações de diferentes grupos identificados em um cluster.

A seguir apresento uma sequencia básica para a realização da análise, incluindo testes à posteriori paramétricos e não paramétricos. Para a comparação múltipla não paramétrica é utilizado o pacote npmc. Os dados utilizados no exemplo são gerados no início da rotina.

# gera dados

GRP <- as.factor(rep(c("A","B","C","D"), c(30,50,40,20)))
UA <- round(rnorm(30,360,75),1)
UB <- round(rnorm(50,500,90),1)
UC <- round(rnorm(40,150,30),1)
UD <- round(rnorm(20,385,60),1)
U <- c(UA,UB,UC,UD)
dat.GRP <- data.frame(GRP,U)
rm(GRP,UA,UB,UC,UD,U)

# explora dados
attach(dat.GRP)
dat.GRP
summary(dat.GRP)
aggregate(U,by = list(GRP=GRP), FUN = "length")
aggregate(U,by = list(GRP=GRP), FUN = "min")
aggregate(U,by = list(GRP=GRP), FUN = "max")
aggregate(U,by = list(GRP=GRP), FUN = "mean")
aggregate(U,by = list(GRP=GRP), FUN = "sd")
aggregate(U,by = list(GRP=GRP), FUN = "median")
boxplot(U~GRP,notch=T,at=rank(tapply(U,GRP, median)))

# --- adaptado de Dalgaard 2008 Introductory Statistics with R
xbar <- tapply(U, GRP, mean)
s <- tapply(U, GRP, sd)
n <- tapply(U, GRP, length)
sem <- s/sqrt(n)
stripchart(U~GRP, method="jitter", jitter=0.05, pch=16, vert=T)
arrows(1:length(levels(GRP)),xbar+sem,1:length(levels(GRP)),xbar-sem,angle=90,code=3,length=.1)
lines(1:length(levels(GRP)),xbar,pch=4,type="b",cex=2)
rm(xbar,s,n,sem)
# ---

# ANOVA e teste de comparação múltipla de Tukey
aov.GRP<-aov(U~GRP)
summary(aov.GRP)
TukeyHSD(aov.GRP, ordered = TRUE)
plot(TukeyHSD(aov.GRP, ordered = TRUE))

# teste de Kruskal–Wallis e teste de comparação múltipla não paramétrica
kru.GRP<-kruskal.test(U~GRP)
kru.GRP

ori.names<-names(dat.GRP)
names(dat.GRP)<-c("class","var") # o pacote npmc requer que o nome da coluna que identifica o grupo seja class e a de valores var

library(npmc)
npmc(dat.GRP)

names(dat.GRP)<-ori.names
rm(ori.names)
detach(dat.GRP)