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)

Gráfico de dipersão com boxplots

O gráfico de dispersão e uma forma bastante ilustrativa de representar a correlação entre duas variáveis. Informações sobre a distribuição das variáveis da ordenada e da abscissa podem ser representadas por boxplots colocados na margem do gráfico de dispersão.

Para o exemplo abaixo utilizei os dados e as curvas de regressão do tópico Relação comprimento-peso de peixes

# Plotagem do gráfico de dispersão

zones=matrix(c(2,0,1,3),ncol=2,byrow=TRUE)
layout(zones,widths=c(20,5),heights=c(5,20))
par(mar=c(8,8,0,0))
plot(P~C,subset=G=="M",pch=20,xlab="Ct (mm)",ylab="Pt (g)",cex.axis=1.3,cex.lab=1.3)

# para plotagem da regressão linear
curve(exp(coef(lm.CPm)[1])*x^(coef(lm.CPm)[2]),col="blue",add=T)
posicao<-locator(1) # clique no gráfico para indicar onde colocar a legenda
legend(posicao,bty="n",cex=1.5,text.col="blue",legend=substitute(Pt==a%*%Ct^b,list(a=exp(coef(lm.CPm)[1]),b=round(coef(lm.CPm)[2],3))))
posicao$y <- posicao$y-500
legend(posicao,bty="n",cex=1.5,text.col="blue",legend=substitute(R^2==r,list(r=round(summary(lm.CPm)$adj.r.squared,3))))

# para plotagem da regressão não linear
curve(coef(nls.CPm)[1]*x^coef(nls.CPm)[2], col="darkgreen",add=T)
posicao<-locator(1) # clique no gráfico para indicar onde colocar a legenda
legend(posicao,bty="n",cex=1.5,text.col="darkgreen",legend=substitute(Pt==a%*%Ct^b,list(a=coef(nls.CPm)[1],b=round(coef(nls.CPm)[2],3))))
posicao$y <- posicao$y-500
legend(posicao,bty="n",cex=1.5,text.col="darkgreen",legend=substitute(R^2==r,list(r=round(Rsq.ad(nls.CPm),3))))

# para plotagem dos boxplots
par(mar=c(0,8,0,0))
boxplot(C,subset=G=="M",axes=FALSE, space=0, horizontal=TRUE)
par(mar=c(8,0,0,0))
boxplot(P,subset=G=="M",axes=FALSE, space=0)
 

Como ordenar as categorias de um boxplot por suas medianas

A técnica gráfica de boxplot é muito utilizada na análise exploratória de dados, etapa essencial em um trabalho.

No entanto o R, por padrão, coloca os fatores em ordem alfabética. Uma forma de reordenar estes fatores é através da função ordered, como descrito neste tópico.

Outra forma muito útil é através do argumento at e das funções rank e tapply, exemplificado a baixo.

categ <- rep(c("c","b","a"),c(10,10,10))

valor<-c(rnorm(10,5,2),rnorm(10,12,4),rnorm(10,8,3))

dados1<-data.frame(categ,valor)

# boxplot normal, não ordenado

boxplot(dados1$valor~dados1$categ)

# boxplot ordenado

boxplot(valor~categ,

at=rank(tapply(dados1$valor,dados1$categ, median)))

A forma de ordenamento acima é a mais simples, mas se houver valores de mediana iguais as caixas sairão sobrepostas. Uma forma alternativa de ordenamento é:


categ2 <- with(dados1, factor(categ,
levels=levels(categ)[order(tapply(valor,categ,median))]))
boxplot(valor~categ2,dados1)
rm(categ2)


PS: encontrei esta dica no seguinte link, onde outras possibilidades também são apresentadas.