Dados

Como exemplo, vamos utilizar os dados de medidas morfológicas de pardais que sobreviveram e que não sobreviveram à uma tempestado disponibilizados no material que acompanha o livro de Manly and Navarro Alberto (2017) (baixe os dados aqui).

pardais <- read.csv("Bumpus sparrows.csv",header=TRUE)

Vamos primeiro observar as matrizes de correlação (que são as matrizes de covariância padronizadas pelos desvios padrões) para cada grupo (sobreviventes e não sobreviventes), pois são mais fáceis de se comparar:

cov(pardais[pardais$Survivorship == "S", -1])

##                Total_length Alar_extent L_beak.head L_humerous L_keel_sternum
## Total_length      11.047619        9.10   1.5566667     0.8700      1.2861905
## Alar_extent        9.100000       17.50   1.9100000     1.3100      0.8800000
## L_beak.head        1.556667        1.91   0.5313333     0.1890      0.2396667
## L_humerous         0.870000        1.31   0.1890000     0.1760      0.1325000
## L_keel_sternum     1.286190        0.88   0.2396667     0.1325      0.5749048

cor(pardais[pardais$Survivorship == "NS", -1])

##                Total_length Alar_extent L_beak.head L_humerous L_keel_sternum
## Total_length      1.0000000   0.7761963   0.6769768  0.6824212      0.6568714
## Alar_extent       0.7761963   1.0000000   0.6978185  0.7845546      0.6200093
## L_beak.head       0.6769768   0.6978185   1.0000000  0.8347046      0.5698878
## L_humerous        0.6824212   0.7845546   0.8347046  1.0000000      0.6677936
## L_keel_sternum    0.6568714   0.6200093   0.5698878  0.6677936      1.0000000

Para facilitar a inspeção visual, eu fiz a diferença entre as correlações do grupo não-sobrevivente com o sobrevivente.

cor(pardais[pardais$Survivorship == "NS", -1]) - cor(pardais[pardais$Survivorship == "S", -1])

##                Total_length Alar_extent L_beak.head L_humerous L_keel_sternum
## Total_length     0.00000000  0.12172887  0.03446994 0.05850172      0.1465157
## Alar_extent      0.12172887  0.00000000  0.07144870 0.03811285      0.3425714
## L_beak.head      0.03446994  0.07144870  0.00000000 0.21665701      0.1362510
## L_humerous       0.05850172  0.03811285  0.21665701 0.00000000      0.2512489
## L_keel_sternum   0.14651569  0.34257144  0.13625102 0.25124892      0.0000000

Aqui, podemos ver que, apesar de muitas correlações terem uma diferença pequena, todas as diferenças são positivas indicando que o grupo não-sobrevivente possui maior variabilidade nas medidas morfológicas do que o grupo sobrevivente. Seriam estas diferenças significativas?

Teste M de Box (Box M-test)

Este teste compara a variação em duas ou mais amostras multivariadas. É considerado um testes sensível à premissa de normalidade multivariada, e portanto é possivel que um resultado significativo seja devido à não normalidade dos dados do que à matrizes de covariância não homogêneas (Manly and Navarro Alberto (2017)).

Vamos primeiro observar a normalidade multivariada dos dados (roteiro útil aqui).

library(MVN)
mardiaTest(pardais[,-1], qqplot=T)

##    Mardia's Multivariate Normality Test 
## --------------------------------------- 
##    data : pardais[, -1] 
## 
##    g1p            : 5.139891 
##    chi.skew       : 41.97578 
##    p.value.skew   : 0.1941829 
## 
##    g2p            : 35.93328 
##    z.kurtosis     : 0.3904202 
##    p.value.kurt   : 0.6962258 
## 
##    chi.small.skew : 45.45482 
##    p.value.small  : 0.1110146 
## 
##    Result         : Data are multivariate normal. 
## ---------------------------------------

Agora vamos apliar o teste M de box que está disponível no pacote biotools com a função boxM. Nesta função, damos a matriz de variáveis quantitativas que queremos comparar e o vetor que identifica os grupos a serem comparados.

library(biotools)

## ---
## biotools version 3.1

boxM(pardais[,-1], pardais[,1])

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  pardais[, -1]
## Chi-Sq (approx.) = 10.408, df = 15, p-value = 0.7933

Com este teste, concluimos que as diferenças que observamos nas matrizes de correlação não são signficativas. Então, podemos dizer que há homocedasticidade multivariada nos dados.

Método PERMDIST: uma generalização do teste de Levene univariado

Podemos também comparar a homocedasticidade multivariada com uma generalização do teste de homogeneidade de variâncias de Levene proposto por Anderson (2006). Esta proposta está baseada na distância das observações aos centróides dos grupos no espaço multivariado (distância de mahalanobis). Para testar se as variâncias dos grupos são diferentes, as distâncias entre as observações dos grupos ao centróide do grupo são comparadas por ANOVA. Quando usamos a distância euclidiana, temos o análogo multivariado do teste de Levene, mas é possível usar outras distâncias.

A função betadisper() do pacote vegan faz o teste e o help(betadisper) provém bastante informação interessante sobre essa abordagem, vale a pena conferir!

library(vegan)
pardist <- dist(pardais[,-1], "euclidian")
mod <- betadisper(pardist, pardais[,1])
anova(mod)

## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value Pr(>F)
## Groups     1   19.6 19.6021  2.3021 0.1359
## Residuals 47  400.2  8.5149

O resultado deste método também indica homocedasticidade entre as matrizes de covariância para os pardais sobreviventes e não sobreviventes.

O legal desta abordagem é que conseguimos desenhar um gráfico dos eixos de uma Análise de Coordenadas Principais ( PCoA), comparando os grupos.

plot(mod)

Neste gráfico vemos que os pardais que não-sobrevieram tem maior variabilidade (as distâncias de suas observações são mais espalhadas), como havíamos suposto comaprando as matrizes de correlação. E também vemos que não parece haver diferença na variabilidade entre os pardais que sobreviveram ou não.

Referências

Anderson, M.J. (2006) Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62, 245–253.