[ISLR] 4.분류(Classification) - 선형판별분석(LDA)

 

저번 포스팅에서는 로지스틱 회귀분석을 통한 분류를 알아보았는데요. 이번에는 또 다른 분류기법인 선형판별분석에 대해 알아보도록 하겠습니다.

 

 

 

1. 로지스틱 회귀분석과 선형판별분석(LDA)

 

선형판별분석은 로지스틱회귀분석과 접근하는 방법이 조금 다릅니다. 로지스틱 회귀분석은 일반적으로 2가지 범주에 대해서 적용하는 반면에 선형판별분석은 분류를 하기 위해서는 반응변수 Y가 2개 이상 유한 개의 범주로 구성되어 있어야 하는데요!

이때 반응변수 Y의 각 범주에 대해 설명변수 X의 분포를 모델링하고 베이즈 정리를 통해 $P(Y=k|X=x)$에 대한 추정치를 구합니다.

 

로지스틱 회귀분석을 이용하지 않고 선형판별분석을 이용해야 하는 경우는 아래와 같습니다.

1. 클래스가 잘 분리될 때 선형판별분석이 로지스틱 회귀분석보다 안정적이다.
   • 특정값을 기준으로 클래스가 명확하게 분리되는 경우 로지스틱 회귀분석의 계수는 매우 불안정해집니다.
2. 표본수가 적고 각 클래스에서 X의 분포가 근사적으로 정규분포를 따르면 선형판별분석이 안정적이다.
3. 범주의 수가 3개 이상일 때, 일반적으로 로지스틱 회귀분석보다 선형판별분석을 사용한다.

 

 

 

2. 기호와 베이즈정리

 

모델링에 앞서 아래와 같이 기호들을 정리하였습니다.

• K : 범주의 수 (K ≥ 1인 자연수)
• $\pi_k$ : 관측치가 k번째 클래스에서 나올 사전확률(prior probabilitay)로 데이터로부터 추정하거나 믿음(belief)를 반영하는 분포
• $f_k(X) = P(X=x|Y=k)$ : k번째 클래스에 속하는 관측치의 밀도함수(pmf, pdf)로 데이터로부터 추정
• $P_k(X) = P(Y=k|X)$ : $X=x$가 주어졌을 때, 해당 관측치의 반응변수가 k에 속하는 확률

 

우리의 목표는 $Y=k$인 관측치 $X=x$들에 대해 반응변수가 K에 속할 확률 $P_k(X) = P(Y=k|X=x)$의 값이 크게 나오도록 모델링을 해야합니다. 이때, 분포와 사전확률을 알고있으므로 베이즈 정리에 의해 아래와 같은 모델을 생각해 볼 수 있습니다. 

$$ Pr(Y=k | X=x) = \frac{\pi_kf_k(x)}{\sum_{i=1}^{K}\pi_if_i(x)}$$ 

베이즈 분류기는 $X=x$로 주어졌을때 $Pr(Y=k|X=x)$가 가장 큰 k번째 클래스로 분류하게 됩니다.

수학적으로 좀 더 자세하게 알고 싶은 분들은 베이지안 분포 포스팅(https://moogie.tistory.com/131)을 참고해주세요.

 

 

 

3. 선형판별분석(LDA) 모델

 

• 가정(Assumption)
  1. $f_k(X)$가 다변량정규분포라고 가정
  2. 모든 K개 클라스에서 등분산성을 만족한다고 가정 (=각 클라스의 분산은 공분산행렬을 가짐)
  3. $\pi_k$, $u_k$, $\sum^2$를 추정해야함 (for i=1, ... ,K)
     • 추정은 아래와 같이 합니다. (일변량 정규분포의 경우)
     • $\hat{\pi_k} = \frac{n_k}{\sum_{i=1}^{K}n_i}$
     • $\hat{u_k} = \frac{\sum_{i:y_i=k}x_i}{n_k} $
     • $\hat{\sigma^2} = \frac{ \sum_{k=1}^{K} \sum_{i:y_i=k}^{} (x_i-\hat{u_k})^2 } {n-K}$

• 예측(Prediction)
  1. 가정한 모수들의 추정치를 가지고 베이즈 정리에 의해 유도된 식에 대입 후 정리
  2. p=1인 일변량 정규분포의 경우 아래식을 만족하는 k를 예측값으로 할당 
     • $Argmax_k(x \frac{\hat{u_k}}{\hat{\sigma^2}} - \frac{\hat{u_k^2}}{2\sigma^2} + log(\hat{\pi_k}))$
  3. p>1인 다변량 정규분포의 경우 아래식을 만족하는 k를 예측값으로 할당
     
     • $Argmax_k(x^T\sum^{-1}u_k - \frac{1}{2}u_k^T\sum^{-1}u_k + log(\pi_k))$

 

참고로, 일변량 및 다변량 정규분포는 다음과 같이 표기합니다.

$X \sim N(u, \sigma^2) \quad f(x) : \frac{1}{\sqrt{2\pi\sigma^2}}exp(-\frac{(x-u)^2}{2\sigma^2})$

$X \sim N(u, \sum) \quad f(x) : \frac{1}{(2\pi)^{p/2}|\sum|^{1/2}}exp(-\frac{1}{2}(x-u)^T\sum^{-1}(x-u))$

 

 

4. 요약(Summary)

 

 

5. 예제(Example)

로지스틱 회귀분석과 비교하기 위해서 이전 포스팅에서 사용한 대학원 진학과 관련된 데이터를 이용해 선형판별분석을 진행하였습니다.

사용한 데이터 셋의 반응변수(Y)는 admit로 대학원진학 여부를 나타내는 이진 변수로 1이면 대학원에 진학했음을 나타냅니다. 또한 모델에 사용한 설명변수(X)의 이름과 뜻은 다음과 같습니다. (GRE : 시험점수, GPA : 학점, rank : 학부의 명성(1~4, 작을수록 명성이 좋음))

data <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
data$rank <- factor(data$rank)

 

LDA모델에 적합하기에 앞서, 범주별 설명변수의 정규성을 검정해야 하는데요. 범주별로 설명변수의 분포가 정규분포를 따르지 않는다면 Box-Cox변환과 같은 변환을 고려해야 합니다. 이번 예시에서는 KS검정을 사용한 ad.test 적용한 결과 정규분포를 따른다는 것을 알 수 있습니다. 

data %>% group_by(admit) %>% group_map(~nortest::ad.test(.$gre)$p.value)

또한, 아래 시각화에서 어느정도의 skewness가 있지만 정규분포라고 생각할 수 있습니다.

data %>% 
  mutate(admit = as_factor(admit)) %>%
  ggplot(mapping=aes(x=gre, fill=admit)) +
  geom_histogram(aes(y = after_stat(density)), color = "white", alpha=0.5,
                 position = "identity", bins=15) +
  geom_density(position="identity", alpha=0.5)

 

R에서는 MASS::lda( ) 또는 Tidymodeling에서 parsnip::discrim_linear( ) 함수를 사용해서 적합할 수 있습니다. 

lda.fit <- MASS::lda(admit ~ gre + gpa + rank, data = data)
print(lda.fit)

• Prior Probability 부분에서 사전확률로 0과 1의 해당되는 값은 68.25%와, 31.75%입니다.
• Group means는 반응변수의 범주별로 설명변수의 평균값을 나타낸 부분입니다. 범주별로 상당히 차이가 있는 변수들이 존재하네요.
• Coefficienmts of linear discriminants는 결정경계를 생성할 때 사용합니다.

$$LD1 = 0.003gre + 1.1gpa - 1.15rank2 - 2.07rank3 - 2.3rank4$$

 

아래와 같이 결정경계를 만들었다는 것을 확인할 수 있습니다. 

data %>% tidytable::get_dummies(cols = rank) %>% dplyr::select(-rank_1) %>% 
  mutate(LD1 = ld[1]*gre + ld[2]*gpa + ld[3]*rank_2 + ld[4]*rank_3 + ld[5]*rank_4)  %>% 
  cbind(., predict(lda.fit, .)$posterior, ld1 = predict(lda.fit, .)$x, class=predict(lda.fit,. )$class) %>% 
  ggplot(mapping=aes(x=ld1.LD1, fill = factor(class))) + geom_histogram()

 

마지막으로 gre값과 rank에 따른 모델의 예측값(admit가 1일 확률)을 시각적으로 표현하였습니다.

predict_grid <- data %>% modelr::data_grid(gre = seq_range(gre, 500, expand = 0.1), rank, gpa = median(gpa)) 
pred <- predict(lda.fit, newdata = predict_grid)
prob <- pred$posterior[,2]
predict_grid %>% mutate(prob) %>% 
  ggplot(mapping=aes(x=gre, y=prob)) + geom_line(mapping=aes(color = rank))