Uma das maneiras mais comuns de se observar valores extremos (outliers) em nossos dados é usar o gráfico de boxplot em cada variável de interesse (veja no link como o boxplot é usado para classificar valores extremos).

Entretanto, quando temos dados multivariados é interessante observar os valores extremos de maneira conjunta, isto é, levando em consideração todas as variáveis dos dados. Neste caso, os boxplots univariados não nos servem mais¹. Detectar valores extremos em observações univariadas é diferente de detectá-los de maneira multivariada e o que era um outlier em uma variável pode não ser mais num contexo multivariado e vice-versa.

Neste roteiro, apresento os boxplots bivariados do pacote MVA, descritos em Everitt and Hothorn (2011) e o uso da Distância de Mahalanobis de cada observação ao centróide dos dados para diagnóstico de valores extremos.

Boxplots em um exemplo bivariado

Vamos usar um exemplo com apenas duas variáveis para facilitar a inspeção visual.

Escolhemos duas variáveis do conjunto de dados swiss do pacote datasets² e vamos primeiro inspecioná-las separadamente com os boxplots univariados.

# para determinar as observações outiers a serem nomeadas no plot
out <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | 
           x > quantile(x, 0.75) + 1.5 * IQR(x))
}
# OBS: é o critério usado na identificação dos oultiers nos boxplots

par(mfrow=c(1,2))
boxplot(swiss$Examination, main = "Examination")
boxplot(swiss$Education, main = "Education")
text(x=c(1,1,1,0.8,1.2), 
           y=swiss[out(swiss[,4]),4]-c(3,-2,2,1,-1), 
           label= rownames(swiss[out(swiss[,4]),]))
par(mfrow=c(1,1))

Para Examination não encontramos outliers segundo o método do boxplot, mas para Education temos 5 valores extremos.

Vamos agora desenhar um boxplot bivariado conforme Everitt and Hothorn (2011). Os boxplots bivariados são baseados nos cálculos “robustos” das medidas de locação, escala e correlação e consiste em um par de elipses concêntricas, uma das quais inclui 50% dos dados e a outra que delinea o limite dos dados e deixa de fora os possíveis outliers. Também são desenhadas linhas de regressão de x em y e de y em x, com sua interseção mostrando o parâmetro de locação (média das duas variáveis). O ângulo entre as linhas de regressão serão tão pequenos quanto maior for a correlação entre as variáveis.

Vamos usar a função bvbox do pacote MVA que acompanha o livro de Everitt and Hothorn (2011), para observar a relação entre as variáveis.

library(MVA)
bvbox(swiss[,c(3,4)], xlab = colnames(swiss)[3], 
      ylab = colnames(swiss)[4])
text(x = swiss[out(swiss[,4]),3], 
     y = swiss[out(swiss[,4]),4]-2.5, 
     rownames(swiss[out(swiss[,4]),]))

No boxplot bivariado, detectamos os mesmos outliers encontrados apenas em Education.

A desvantagem dos boxplots bivariados é que eles avaliam os outliers apenas para duas variáveis por vez e geralmente temos mais do que duas variáveis nos nossos dados. A seguir veremos uma forma mais geral de se avaliar os valores extremos multivariados.

Usando a distância de Mahalanobis

A distância de Mahalanobis (\(D^2_M\)) é uma medida de distância entre as observações (unidades amostrais) e pode também ser usada para medir a distância entre uma única observação multivariada e o centro da população (multivariado não se esqueça), da qual esta esta observação vem (Manly and Navarro Alberto (2017)). No nosso caso, vamos calcular a distância da observação ao centróide da distribuição dos dados.

O cálculo da distância de Mahalanobis leva em consideração a correlação entre as variáveis, e seu valor \(D^2_M\) pode ser pensado como um resíduo multivariado para a observação \(x\). \(D^2_M\) é uma medida de quão longe a obervação \(x\) está do centro da distribuição de todos os valores, levando em consideração todas as variáveis consideradas e suas covariâncias (Manly and Navarro Alberto (2017)). Um valor alto de \(D^2_M\) para a observação \(x\) pode indicar um possível outlier multivariado, se excluímos a possibilidade de erro de digitação/registro ou que a observação seja de outra distribuição.

Vamos calcular \(D^2_M\) para cada observação dos dados apresentados acima, e plotar estas distâncias. Um critério para o diagnóstico dos outliers a partir da distância de Mahalanobis é verificar se a distância calculada está acima de um limiar \(c^2\), que nada mais é do que a estatística do Qui-quadrado (\(\chi^2\)) calculada para um valor de probabilidade escolhido. No exemplo abaixo, eu escolhi um \(\alpha\) de 10% e fui na distribuição de (\(\chi^2\)) com 2 graus de liberdade (g.l. = número de variáveis sendo analisadas) saber qual é o quantil e desenhei a linha pontilhada no gráfico com este valor.

maha <- mahalanobis(swiss[,3:4], # dados
          center = colMeans(swiss[,3:4]), # médias das variáveis
          cov = cov(swiss[,3:4])) # matriz de covariâncias

#calculando o cˆ2 para delimitar o limiar
quant <- qchisq(0.1, 2, lower.tail = F)

plot(1:length(maha), maha, xlab = "Observações", 
     ylab= "Distância de Mahalanobis")
abline(h=quant, lty=2)

out <- maha[maha>quant]
text(c(42,44,46,48)-7, out, names(out))

Desta forma, para o critério usando a distância de Mahalanobis, 4 observações foram consideradas outliers bivariados. Este critério está dizendo que as 4 observações com maior distância tem pouca probabilidade de serem advindas da população com distribuição normal bivariada assumida nos dados.

Diagnosticando outliers multivariados

Agora vamos expandir nosso exemplo usando todas as variáveis existentes nos dados swiss. Aqui eu não vou fazer os boxplots bivariados (você já sabe como faz), mas quero apenas mostrar o plot da distância de Mahalanobis para ver se, com mais variáveis, temos diferentes outliers detectados.

maha2 <- mahalanobis(swiss, center = colMeans(swiss), 
                     cov = cov(swiss))

# lembre-se que os graus de liberdade agora são 6
quant2 <- qchisq(0.90, 6)

plot(1:length(maha2), maha2, xlab = "Observações", 
     ylab = "Distância de Mahalanobis")
abline(h = quant2, lty = 2)

out2 <- maha2[maha2 > quant2]
text(c(2.5,15,34,38,40,44), out2, names(out2))

Com mais variáveis, tivemos diferentes observações detectadas como outliers e por isso a importância de se analisar os valores extremos em dados multivariados de maneira conjunta!

Usando pacote `MVN`

O pacote MVN possui uma função que também detecta outliers multivariados a partir da distância de Mahalanobis. Esse pacote usa dois métodos robustos (ver detales aqui) para calcular a distância de Mahalanobis (a distância de Mahalanobis Robusta e a Ajustada), e portanto dão resultados um pouco diferentes do que fizemos anteriormente.

Outra diferença é que resultado é apresentado em um gráfico quantil-quantil (qqplot), com um dos eixos sendo a \(D^2_M\) e o outro os quantis de uma distribuição \(\chi^2\).

library(MVN)

outliers <- mvOutlier(swiss, alpha = 0.90, method = "quan") 
# tente method="adj.quan"

Observações finais

Vale lembrar que estes valores extremos não necessariamente são erros de digitação/medida, mas podem ser valores reais, porém pouco prováveis. A decisão sobre o que fazer com estes possíveis outliers detectados depende muito da natureza dos dados, do contexto de coleta dos dados e do tipo de análise a ser conduzida.

Uma sugestão feita por Manly and Navarro Alberto (2017) é a de conduzir análises com e sem os valores extremos e comparar os resultados. Se a conclusão é a mesma, então não há problemas. Mas se as conclusões dependem muito dos outliers, então é preciso muito mais cuidado com as decisões a serem tomadas.

Referências

Everitt, B, and T Hothorn. 2011. An Introduction to Applied Multivariate Analysis with R.

Manly, Bryan F. J., and Jorge A. Navarro Alberto. 2017. Multivariate Statistical Methods: A Primer. Fourth edition. Boca Raton: CRC Press.

eles também não servem se a quantidade de dados for muito grande, veja esse roteiro sobre gráficos alternativos aos boxplots nessa situação.↩
que não precisa ser carregado porque já vem junto como R base.↩

Diagnóstico de valores extremos em dados multivariados