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)
 

Ajuste do modelo de von Bertalanffy

O modelo de crescimento de von Bertalanffy é muito utilizado para descrever a variação de comprimento de peixes, moluscos e crustáceos ao longo do tempo.

A seguir apresento um passo-a-passo para ajustar a curva aos dados de comprimento (Lt) na idade (t), analisar os parâmetros e fazer o gráfico.
O ajuste é feito de forma não linear pela função nls,  os intervalos de confiança das estimativas são calculados com a função confint e o coeficiente de determinação (R2) pela função Rsq do pacote qpcR. As funções expression e substitute são utilizadas para escrever as equações no gráfico.

t	Lt
1	102,0
2	167,0
3	219,4
4	260,7
5	294,9
6	323,2
7	343,0
8	369,5
9	401,7
10	410,0

 

# carrega pacote para cálculo do R2

library("qpcR")

# importa dados da área de transferência

dat.tL <- read.delim("clipboard",dec=",")

attach(dat.tL)

# ajuste do modelo

vb.pargo <- nls(Lt~Linf*(1-exp(-k*(t-t0))),start=list(Linf=500,k=0.2,t0=0))
summary(vb.pargo)

Formula: Lt ~ Linf * (1 - exp(-k * (t - t0)))

Parameters:
      Estimate Std. Error t value Pr(>|t|)   
Linf 501.51567   19.81444  25.311 3.84e-08 ***
k      0.16185    0.01541  10.504 1.55e-05 ***
t0    -0.46264    0.13937  -3.319   0.0128 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.587 on 7 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 5.881e-07

# calcula intervalo de confiança dos parâmetros

confint(vb.pargo)

            2.5%       97.5%
Linf 461.8889677 562.3393409
k      0.1250699   0.1996879
t0    -0.8411141  -0.1621223


1-(deviance(vb.pargo)/((length(Lt)-1)*var(Lt))) # R2 "na mão"

[1] 0.9976615

Rsq(vb.pargo) # R2 pelo pacote qpcR

[1] 0.9976615

Rsq.ad(vb.pargo) # R2 ajustado pelo pacote qpcR
[1] 0.9969934 

# desenha o gráfico. Em Windows substituir ["\U221E"] por [infinity]

plot(Lt~t,xlab="idade (anos)",ylab="comprimento total (mm)",
xlim=range(0,10),ylim=range(0,500),cex.lab=1.2,
main=expression(L[i]==L["\U221E"]*"["*1-e^{-k(t-t[0])}*"]"),cex.main=1.5)
 

#desenha a curva

curve(coef(vb.pargo)[1]*(1-exp(-coef(vb.pargo)[2]*(x-coef(vb.pargo)[3]))),add=T,col="tomato1")
 

# coloca a legenda, deve-se clicar no gráfico para indicar o local da legenda

legend(locator(1),bty="n",legend=substitute(L[i]==Linf%*%"["*1-e^{-k%*%(t-t0)}*"]", list(Linf=round(coef(vb.pargo)[1],1),k=round(coef(vb.pargo)[2],2),t0=-round(coef(vb.pargo)[3],2))),cex=1.5)

detach(dat.tL)