마할라노비스 거리(Mahalanobis Distance)는 데이터 밀도(분산), 분포를 고려한 거리로 다변량 이상치를 판단하는 데 가장 대표적인 방법으로 정규성 조건등을 따지지 안아도 됨.

library(dplyr)
## 
## 다음의 패키지를 부착합니다: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(sp)

airquality
str(airquality);glimpse(airquality)
## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
## Rows: 153
## Columns: 6
## $ Ozone   <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 14, …
## $ Solar.R <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290, 27…
## $ Wind    <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9…
## $ Temp    <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58, 64…
## $ Month   <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ Day     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…
air = airquality %>% select(Ozone, Temp)
air = na.omit(air) #결측값 제거
air

다음 5단계를 순서대로 진행하여 마할라노비스 거리로 이상치를 찾음!

① Ozone와 Temp의 중심점을 구함.(평균) ② Ozone와 Temp의 공분산 행렬을 구함. ③ 중심점과 각 점의 마할라노비스 거리를 구함. ④ 카이제곱분포를 이용하여 Cut-Off 영역을 구함. ⑤ 거리가 Cut-Off 영역 밖이면 이상치로 판단함.

① Ozone와 Temp의 중심점을 구함.(평균)

(air.center = colMeans(air))
##    Ozone     Temp 
## 42.12931 77.87069

② Ozone와 Temp의 공분산 행렬을 구함.

(air.cov = cov(air))
##           Ozone      Temp
## Ozone 1088.2005 218.52121
## Temp   218.5212  89.97444

③ 중심점과 각 점의 마할라노비스 거리를 구함. 다음으로 Cut-Off를 구하기 위해 자유도가 2(변수 갯수)인 카이제곱분포의 95%분위를 구하고 그 분위수에 루트를 씌웁니다. 그런 다음 car 패키지의 ellipse() 함수를 이용하여 Cut-Off 영역을 구합니다.

#카이제곱분포의 95% 분위수
rad  = qchisq(p = 0.95 , df = ncol(air)) 
 #Cut-Off 영역의 반지름
rad  = sqrt(rad)
#Cut-Off 영역찾기
ellipse <- car::ellipse(center = air.center , shape = air.cov , radius = rad ,segments = 150 , draw = FALSE)

위의 코드처럼 car 패키지의 ellipse() 함수의 center 인자에는 데이터의 중심점, shape인자에는 공분산, radius 에는 Cut-off 영역의 반지름, segments 인자에는 Cut-Off 영역을 선분갯수를 입력합니다. 다음으로 Cut-Off 영역과 데이터의 산포도를 그려 이상치를 판별합니다.

ellipse <- as.data.frame(ellipse)
colnames(ellipse) <- colnames(air)

# 산포도 및 Cut-Off 영역 그리기
p<-ggplot(air , aes(x = Ozone , y = Temp)) +
    geom_point(size = 2) +
    geom_polygon(data = ellipse , fill = "orange" , color = "orange" , alpha = 0.2) +
    geom_point(aes(air.center[1] , air.center[2]) , size = 5 , color = "blue") +
    geom_text( aes(label = row.names(air)) , hjust = 1 , vjust = -1.5 ,size = 4 ) +
    ylab("Temp Values") + xlab("Ozone Values")
p

# Extract components
build <- ggplot_build(p)$data  #ggplot객체 p에서 모든 도형속성데이터 추출..
points <- build[[1]] #점좌표 DF생성
ell <- build[[2]] #타원(ellipse)의 원주좌표만 DF로 생성

# Find which points are inside the ellipse, and add this to the data
#point.in.polygon --> sp 패키지에서 제공하는 함수..타원내부에 있는 점엔 TRUE 바깥은 FALSE!
#도형점좌표 points[1:2] 에서 타원내부 점들만 in.ell로 로지컬컬럼(T/F)생성..
dat_ell <- data.frame(points[1:2], in.ell = as.logical(point.in.polygon(points$x, points$y, ell$x, ell$y)))

#타원도형 바깥에 존재하는 점들만 추출하여 DF생성--> 다른색(빨간색)으로 나타내기 위해..
out_dots <- dat_ell %>% filter(in.ell=="FALSE")

ggplot(air , aes(x = Ozone , y = Temp)) +
    geom_point(size = 2) +
    geom_polygon(data = ellipse , fill = "orange" , color = "orange" , alpha = 0.2) +
    geom_point(aes(air.center[1] , air.center[2]) , size = 5 , color = "blue") +
    geom_point(data=out_dots, aes(x=x, y=y) , size = 4 , color = "red") +
    geom_text( aes(label = row.names(air)) , hjust = 1 , vjust = -1.5 ,size = 4 ) +
    ylab("Temp Values") + xlab("Ozone Values")