False Discovery Rate
07 May 2017 | Feature Selection
이번 글에서는 변수선택(Feature Selection) 기법으로 널리 쓰이고 있는 False Discovery Rate(FDR)에 대해 살펴보도록 하겠습니다. 이번 글 역시 고려대 김성범 교수님 강의를 정리하였음을 먼저 밝힙니다. 그럼 시작하겠습니다.
FDR의 문제 의식
제품의 ‘정상’, ‘불량’을 판별하는 이진분류 문제를 푼다고 칩시다. (정상, 불량이든 긍정, 부정이든 어떤 것이든 상관없습니다) 우리가 가진 데이터의 독립변수는 $x_1, x_2, x_3$ 세 개이고 하나의 행(레코드)은 관측치 하나에 대응한다고 가정하겠습니다. 그렇다면 세 변수 가운데 정상, 불량을 판별하는 데 가장 중요한 변수는 무엇일까요?
위 예시에서 $x_1$은 정상, 불량 범주를 판별하는 데 전혀 도움이 안 될 것입니다. 반면 $x_2$는 상대적으로 중요할 겁니다. 정상 관측치에선 작은 값(평균 : 10)을, 불량 관측치에선 큰 값(200)을 갖기 때문입니다. $x_3$는 $x_1$보다는 중요한 변수지만 $x_3$보다는 덜 중요한 변수일 겁니다.
지금까지 말씀드린 내용을 이를 수식으로 나타내면 아래와 같습니다. 이 식에 따라 좌변을 계산하면 각각 0, 190, 1.6이 됩니다.
\[\left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| \ge \tau\]
그렇다면 $τ$가 얼마나 커야 유의미(=$x_i$가 중요변수)할까요? 이를 엄밀히 따져 보려면 통계학 기법을 써야 합니다. 다시 말해 귀무가설(Null Hypothesis)과 대립가설(Alternative Hypothesis)을 아래와 같이 설정하고, 변수마다 각각 가설검정(Hypothesis Test)을 실시해야 한다는 것이지요.
\[{ H }_{ 0 }\quad :\quad \left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| =0\\ { H }_{ 1 }\quad :\quad \left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| >0\]
그런데 여기서 문제가 있습니다. 우리는 대개 변수가 100개를 훌쩍 넘는 다변량 데이터를 취급한다는 점입니다. 모든 가설을 동시에 검정하면 $x_i$가 중요하지 않은 변수($H_0$가 참)인데도 중요변수($H_0$를 기각)라는 결과가 나오는 경향이 있습니다. 이를 다중성(Multiplicity) 문제라고 합니다.
개별 가설검정의 유의수준(significance level)이 0.01, 다시 말해 개별 가설검정의 신뢰도가 0.99인 상황에서 가설검정을 실시한다고 합시다. 유의수준이란 귀무가설 $H_0$가 참인데 잘못해서 기각할 확률(제1종 오류)의 최대값을 뜻합니다.
그럼 변수가 100개인 다변량 데이터에 대해 변수선택을 하기 위해 100번의 가설검정을 실시한다면 100회 가설검정 전체의 신뢰도는 크게 낮아집니다. 다시 말해 중요하지 않은 변수인데도 중요변수라는 결과가 나올 가능성이 높아진다는 얘기입니다. 아래 표를 볼까요?
검정 횟수
검정 전체의 신뢰도
1
$0.99^1=0.99$
2
$0.99^2=0.98$
10
$0.99^{10}=0.90$
100
$0.99^{100}=0.37$
이를 그래프로 나타내면 아래와 같습니다.
FDR은 이러한 다중가설검정 문제를 해결하기 위해 제안된 방법입니다. 물론 개별 유의수준($α$)을 ‘목표값/변수개수’로 정하고 기존대로 가설검정을 수행하는 기법(Family-wise Test)이 제시되기도 했습니다. 예컨대 전체 목표 유의수준이 0.01이고 변수개수가 10개라면 개별 유의수준을 0.001로 놓고 가설검정을 각각 실시하는 방식이죠. 그러나 이 방식은 중요한 변수($H_0$가 거짓)인데도 중요하지 않다($H_0$가 참)는 결론을 내는 경우가 많은 단점이 있다고 합니다.
FDR 절차
FDR은 제1종 오류를 최대한 적게 범하면서도 중요한 가설을 최대한 많이 찾을 수 있도록($H_0$ 기각) 한다는 점이 강점입니다. 검정해야할 다중가설(데이터의 변수 개수)이 $m$개일 때 FDR은 아래와 같이 정의됩니다.
구분
$H_0$ 채택
$H_0$ 기각
Total
$H_0$가 참
$U$
$V$
$m_0$
$H_0$가 거짓
$T$
$S$
$m_1$
Total
$W$
$R$
$m$
\[FDR=E\left[ \frac { V }{ V+S } \right] =E\left[ \frac { V }{ R } \right]\]
위 수식에 정의된 FDR을 해석하면 이렇습니다. FDR은 기각된 귀무가설($H_0$) 가운데 잘못 기각된 가설이 차지하는 비율의 평균으로 사용자가 지정해주는 값입니다. 바꿔 말해 $m$회 가설검정 가운데 각 $H_0$이 참인데 기각된 비율, 즉 검정 전체 유의수준에 가까운 의미라는 뜻입니다. 실제로 FDR 제안자는 FDR을 $α$라고 표기하고 있습니다.
그럼 FDR 절차를 살펴볼까요? 검정해야할 다중가설, 즉 데이터 변수 개수가 15개라고 합시다. 목표로 하는 FDR 값을 0.05로 정해보겠습니다. 우선 각각의 가설에 대해 $p-value$를 구한 뒤 이를 오름차순으로 정렬하고, 순위(rank)를 매깁니다. 아래 표와 같습니다.
$p-value$는 $H_0$이 참일 확률을 의미하므로, $α$를 0.05로 놓고 기존대로 단일 가설검정을 실시했다면 1~9번 변수 가 중요변수라는 결과를 얻었을 겁니다. Family-wise Test 방식대로 ‘목표 유의수준/변수 개수=0.05/15’로 가설검정을 실시한다면, 그 값이 0.0033이므로 1~3번 변수만 중요하다는 결론을 내렸을 겁니다.
FDR 절차에서 $i$는 rank, $m$은 변수 개수, $α$는 사용자가 정한 FDR 값입니다. 예컨대 1번 변수는 $1/15$에 0.05를 곱해 계산합니다. 마찬가지로 2번 변수는 $2/15$에 0.05를 곱해 구합니다.
이 값이 $p-value$보다 큰 변수 랭크의 최대값이 FDR 절차가 뽑은 중요변수의 개수가 됩니다. 위 표 기준으로 설명하자면 이렇습니다. 표 아래쪽부터 시작해 위쪽으로 훑어 가면서 $i/m*α$ 값이 $p-value$보다 큰 지 여부를 검토합니다. 처음으로 $p-value$ 값을 넘어서는 지점은 바로 4번 변수이네요. 따라서 1~4번 변수를 중요한 변수라고 결론을 내게 됩니다.
예제에서도 알 수 있듯 FDR은 단일 가설검정에 비해 제1종 오류가 낮고(비중요 변수 배제), Family-wise Test 방식에 비해 중요한 가설을 많이 찾아낸다는 장점이 있습니다.
FDR 예제
그럼 실제 데이터를 가지고 문제를 풀어보겠습니다. UCI data 가운데 하나를 가져왔습니다. 개요는 아래와 같습니다.
Ionosphere Data Set
데이터 설명 : 전리층(ionosphere)에 대한 레이더 관측치를 모은 데이터
관측치 개수 : 351개
독립변수 : 34개
종속변수 : 전리층의 패턴 유무 (good : 전리층에 특정 패턴 존재, bad : 패턴 없음)
모든 독립변수에 대해 T-test를 수행하였습니다. ($H_0$=두 집단의 평균이 같다, $H_1$=두 집단 평균이 다르다) 변수별 p-value는 아래 표와 같습니다.
변수
P-value
변수
P-value
변수
P-value
변수
P-value
1
2.303168e-11
11
3.060316e-03
21
6.897335e-05
31
7.025167e-07
2
NaN
12
8.189137e-03
22
5.051969e-02
32
5.547445e-01
3
2.339468e-15
13
9.549388e-04
23
3.350101e-04
33
3.748258e-06
4
5.941660e-02
14
1.341207e-03
24
9.191420e-01
34
2.747245e-01
5
3.153540e-15
15
1.017555e-04
25
1.162854e-03
6
1.988959e-02
16
1.507825e-02
26
9.795861e-01
7
7.577585e-13
17
1.036061e-01
27
5.055064e-02
8
6.995369e-04
18
4.976238e-02
28
4.980595e-01
9
5.545817e-07
19
3.255247e-02
29
2.570992e-05
10
4.515815e-02
20
5.604890e-01
30
9.475607e-01
Individual Test($α=0.05$) 결과 선택된 중요 변수는 22개입니다. (p-value가 0.05보다 작은 변수 선택)
1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 28, 30, 32
Family-Wise Test($α=0.05$, $α/m$, $m$=전체 변수 개수) 결과 선택된 중요 변수는 15개입니다. (p-value가 0.0015보다 작은 변수 선택)
1, 3, 5, 7, 8, 9, 13, 14, 15, 21, 23, 25, 29, 31, 33
위 계산결과를 토대로 $α$값을 0.05로 두고 FDR을 수행한 결과는 아래 표와 같습니다.
변수 ID
rank
P-value
$i/m*α$
3
1
2.339468e-15
0.001470588
5
2
3.153540e-15
0.002941176
7
3
7.577585e-13
0.004411765
1
4
2.303168e-11
0.005882353
9
5
5.545817e-07
0.007352941
31
6
7.025167e-07
0.008823529
33
7
3.748258e-06
0.010294118
29
8
2.570992e-05
0.011764706
21
9
6.897335e-05
0.013235294
15
10
1.017555e-04
0.014705882
23
11
3.350101e-04
0.016176471
8
12
6.995369e-04
0.017647059
13
13
9.549388e-04
0.019117647
25
14
1.162854e-03
0.020588235
14
15
1.341207e-03
0.022058824
11
16
3.060316e-03
0.023529412
12
17
8.189137e-03
0.025000000
16
18
1.507825e-02
0.026470588
6
19
1.988959e-02
0.027941176
19
20
3.255247e-02
0.029411765
10
21
4.515815e-02
0.030882353
18
22
4.976238e-02
0.032352941
22
23
5.051969e-02
0.033823529
27
24
5.055064e-02
0.035294118
4
25
5.941660e-02
0.036764706
17
26
1.036061e-01
0.038235294
34
27
2.747245e-01
0.039705882
28
28
4.980595e-01
0.041176471
32
29
5.547445e-01
0.042647059
20
30
5.604890e-01
0.044117647
24
31
9.191420e-01
0.045588235
30
32
9.475607e-01
0.047058824
26
33
9.795861e-01
0.048529412
2
34
NaN
0.050000000
위 표에서 $i/m*0.05$가 $p-value$보다 크다는 조건을 만족하는 rank의 최대값은 19입니다. 그 결과 19개 변수가 선택됐습니다.
3, 5, 7, 1, 9, 31, 33, 29, 21, 15, 23, 8, 13, 25, 14, 11, 12, 16, 6
FDR level 변화에 따른 중요 변수 선택 결과는 아래와 같습니다. FDR level이 클 수록 선택되는 변수 개수가 많아지는 경향을 보였습니다.
FDR level
선택변수(개)
0.01
16
0.02
17
0.03
18
0.04
19
0.06
20
0.08
24
0.09
25
0.14
26
0.35
27
0.61
28
0.64
30
Ionosphere Data Set을 7:3 비율로 학습셋과 검증셋으로 나누고, 학습셋에 대해 LDA 모델을 적합한 뒤 검증셋을 대상으로 단순정확도를 측정하였습니다. 통계적으로 유의한 일부 변수만으로 구축한 모델은 전체 변수를 모두 사용한 모델의 성능에 크게 뒤지지 않음을 확인할 수 있었습니다.
구분
전체
Individual
Family-wise
FDR
학습에 사용된 변수(개)
33
19
9
14
단순정확도
0.8762
0.8476
0.8286
0.8190
위 문제를 푸는 데 사용한 R 코드는 아래와 같습니다.
# load & split data
data <- read.csv('data.csv',header=F)
good <- data[which(data[,35] == 'g'),-35]
bad <- data[which(data[,35] == 'b'),-35]
# T-test
Ttest <- c()
for (i in 1:dim(good)[2]) {
Ttest[i] <- t.test(good[,i], bad[,i])$p.value
}
# individual test
which(Ttest<0.05)
# family-wise test (패키지 사용)
# bonferroni : p-values are multiplied by the number of comparisons
family.wise.test=p.adjust(Ttest, method = "bonferroni")
which(family.wise.test<0.05)
# family-wise test (자작)
which(Ttest < 0.05/length(Ttest))
# FDR (패키지 사용)
library(fdrtool)
fdr <- fdrtool(Ttest, statistic='pvalue', cutoff.method = 'fndr')
fdr$param
# FDR (자작)
fdr.level <- 0.05
var_table <- matrix(0, nrow=length(Ttest),ncol=4)
colnames(var_table) <- c('id', 'rank', 'P-value', 'i/m*FDR')
var_table[,1] <- c(1:length(Ttest))
var_table[,3] <- Ttest
var_table <- var_table[order(var_table[,3]),]
var_table[,2] <- c(1:length(Ttest))
var_table[,4] <- var_table[,2] / length(Ttest) * fdr.level
max(which(var_table[,3] < var_table[,4])) # 중요 변수 개수
var_table[1:max(which(var_table[,3] < var_table[,4])),1] # 중요 변수 추출
# all variables model
library(MASS)
data[,35] <- as.numeric(data[,35])
colnames(data)[35] <- 'Target'
data <- data[,-2]
trn_idx <- sample(1:dim(data)[1], round(0.7*dim(data)[1]))
train_data <- data[trn_idx,]
test_data <- data[-trn_idx,]
allmodel <- lda(Target~.,train_data)
# var selection model
train_data_good <- train_data[which(train_data$Target == 1),-train_data$Target]
train_data_bad <- train_data[which(train_data$Target == 2),-train_data$Target]
Ttest <- c()
for (i in 1:dim(train_data_good)[2]) {
Ttest[i] <- t.test(train_data_good[,i], train_data_bad[,i])$p.value
}
family.wise.test=p.adjust(Ttest, method = "bonferroni")
individual_model <- lda(Target~.,train_data[,c(which(Ttest<0.05),34)])
familywise_model <- lda(Target~.,train_data[,c(which(family.wise.test<0.05),34)])
# fdr model
fdr.level <- 0.05
var_table <- matrix(0, nrow=length(Ttest),ncol=4)
colnames(var_table) <- c('id', 'rank', 'P-value', 'i/m*FDR')
var_table[,1] <- c(1:length(Ttest))
var_table[,3] <- Ttest
var_table <- var_table[order(var_table[,3]),]
var_table[,2] <- c(1:length(Ttest))
var_table[,4] <- var_table[,2] / length(Ttest) * fdr.level
important_vars <- var_table[1:max(which(var_table[,3] < var_table[,4])),1]
fdr_model <- lda(Target~.,train_data[,c(important_vars,34)])
# prediction
all <- predict(allmodel, test_data)
individual <- predict(individual_model, test_data)
familywise <- predict(familywise_model, test_data)
fdr <- predict(fdr_model, test_data)
# calculate ACC
ACC_all <- sum(diag(table(test_data$Target, all$class)))/length(all$class)
ACC_indi <- sum(diag(table(test_data$Target, individual$class)))/length(individual$class)
ACC_fdr <- sum(diag(table(test_data$Target, fdr$class)))/length(fdr$class)
ACC_family <- sum(diag(table(test_data$Target, familywise$class)))/length(familywise$class)
반복샘플링 기반 FDR
지금까지 설명해드린 FDR 기법은 우리가 가진 데이터가 특정 확률분포를 따르고, 이로부터 적당한 방법을 취해 $p-value$를 계산할 수 있어야만 적용이 가능합니다. 하지만 데이터가 특정 분포를 따르지 않거나, 그 분포를 모른다고 할 때는 $p-value$ 계산이 어렵습니다. 반복샘플링 기반 FDR은 말 그대로 반복샘플링을 통해 분포에 대한 가정없이 $p-value$를 계산할 수 있는 기법입니다.
가정과 절차는 이렇습니다. 특정 변수를 기준으로 관측치를 랜덤하게 섞습니다. 만약 해당 변수가 범주 분류에 의미 있는 변수라면 실제 데이터의 통계량이 관측치를 랜덤하게 섞은 데이터의 통계량보다 훨씬 클 겁니다. 아래 그림처럼요.
통계적 유의성을 확보하기 위해서는 이러한 과정을 충분히 많이 반복해야 합니다. 만약 아래처럼 100번을 반복 수행했다고 가정해보겠습니다. p-value는 아래와 같이 구할 수 있습니다.
\({ p-value }_{ { X }_{ 1 } }=\frac { { T }_{ { x }_{ 1 } }^{ (i) }\quad 중에 \quad{ T }_{ { x }_{ 1 } } 보다\quad 큰\quad값의\quad개수}{ 반복횟수\times변수의\quad수 } =\frac { { T }_{ { x }_{ 1 } }^{ (i) }\quad중에 \quad{ T }_{ { x }_{ 1 } } 보다\quad 큰\quad값의\quad개수 }{ 100\times 2 } \\ { p-value }_{ { X }_{ 2 } }=\frac { { T }_{ { x }_{ 2 } }^{ (i) }\quad 중에\quad{ T }_{ { x }_{ 2 } } 보다\quad 큰\quad값의\quad개수 }{ 반복횟수\times변수의\quad수} =\frac { { T }_{ { x }_{ 2 } }^{ (i) } \quad중에\quad{ T }_{ { x }_{ 2 } } 보다\quad 큰\quad값의\quad개수 }{ 100\times 2 }\)
이렇게 p-value를 구하게 되면 이후 절차는 지금까지 설명해드린 FDR 절차 그대로 따르면 됩니다.
이진분류 이외 문제에 FDR 적용하기
p-value만 구하면 이후엔 FDR 기법을 적용할 수 있습니다. 정규분포를 따르는 데이터의 이진분류 문제의 경우 T-test를 수행하면 p-value를 구할 수 있었고, 특정 데이터 분포를 가정할 수 없는 상황이라면 반복 샘플링 기반으로 p-value를 구할 수 있습니다.
그렇다면 $y$가 연속형 숫자인 회귀 문제에서 중요 변수를 찾아내기 위해 FDR을 적용해야 한다면 p-value를 어떻게 구해야 할까요? 귀무가설과 대립가설을 아래와 같이 설정하고 모든 변수에 대해 가설검정을 실시합니다.
\[{ H }_{ 0 }\quad :\quad { \beta }_{ i }=0\\ { H }_{ 1 }\quad :\quad { \beta }_{ i }\neq 0\]
p-value는 귀무가설이 맞을 확률을 의미하므로 위 가설검정에서 p-value가 작게 나올 경우 $H_0$는 기각하게 됩니다. 바꿔 말해 $β_i$가 중요하다는 이야기이죠. 선형회귀 모델의 계수들은 t분포를 따른다고 합니다. 검정통계량 $t$는 아래와 같이 구합니다.
\(t=\frac { \hat { { \beta }_{ i } } -0 }{ sd\{ \hat { { \beta }_{ i } } \} }\)
3개 이상의 범주 분류 문제의 경우 p-value를 어떻게 구할까요? 두 가지 방법이 있습니다. 첫번째는 범주 두 개를 한 세트로 묶어서 이진 분류 문제처럼 p-value를 구하는 겁니다. 나머지는 분산분석(ANOVA) 기법을 적용하는 것입니다.
이번 글에서는 변수선택(Feature Selection) 기법으로 널리 쓰이고 있는 False Discovery Rate(FDR)에 대해 살펴보도록 하겠습니다. 이번 글 역시 고려대 김성범 교수님 강의를 정리하였음을 먼저 밝힙니다. 그럼 시작하겠습니다.
FDR의 문제 의식
제품의 ‘정상’, ‘불량’을 판별하는 이진분류 문제를 푼다고 칩시다. (정상, 불량이든 긍정, 부정이든 어떤 것이든 상관없습니다) 우리가 가진 데이터의 독립변수는 $x_1, x_2, x_3$ 세 개이고 하나의 행(레코드)은 관측치 하나에 대응한다고 가정하겠습니다. 그렇다면 세 변수 가운데 정상, 불량을 판별하는 데 가장 중요한 변수는 무엇일까요?
위 예시에서 $x_1$은 정상, 불량 범주를 판별하는 데 전혀 도움이 안 될 것입니다. 반면 $x_2$는 상대적으로 중요할 겁니다. 정상 관측치에선 작은 값(평균 : 10)을, 불량 관측치에선 큰 값(200)을 갖기 때문입니다. $x_3$는 $x_1$보다는 중요한 변수지만 $x_3$보다는 덜 중요한 변수일 겁니다.
지금까지 말씀드린 내용을 이를 수식으로 나타내면 아래와 같습니다. 이 식에 따라 좌변을 계산하면 각각 0, 190, 1.6이 됩니다.
\[\left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| \ge \tau\]그렇다면 $τ$가 얼마나 커야 유의미(=$x_i$가 중요변수)할까요? 이를 엄밀히 따져 보려면 통계학 기법을 써야 합니다. 다시 말해 귀무가설(Null Hypothesis)과 대립가설(Alternative Hypothesis)을 아래와 같이 설정하고, 변수마다 각각 가설검정(Hypothesis Test)을 실시해야 한다는 것이지요.
\[{ H }_{ 0 }\quad :\quad \left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| =0\\ { H }_{ 1 }\quad :\quad \left| \overline { { x }_{ i,normal } } -\overline { { x }_{ i,abnormal } } \right| >0\]그런데 여기서 문제가 있습니다. 우리는 대개 변수가 100개를 훌쩍 넘는 다변량 데이터를 취급한다는 점입니다. 모든 가설을 동시에 검정하면 $x_i$가 중요하지 않은 변수($H_0$가 참)인데도 중요변수($H_0$를 기각)라는 결과가 나오는 경향이 있습니다. 이를 다중성(Multiplicity) 문제라고 합니다.
개별 가설검정의 유의수준(significance level)이 0.01, 다시 말해 개별 가설검정의 신뢰도가 0.99인 상황에서 가설검정을 실시한다고 합시다. 유의수준이란 귀무가설 $H_0$가 참인데 잘못해서 기각할 확률(제1종 오류)의 최대값을 뜻합니다.
그럼 변수가 100개인 다변량 데이터에 대해 변수선택을 하기 위해 100번의 가설검정을 실시한다면 100회 가설검정 전체의 신뢰도는 크게 낮아집니다. 다시 말해 중요하지 않은 변수인데도 중요변수라는 결과가 나올 가능성이 높아진다는 얘기입니다. 아래 표를 볼까요?
| 검정 횟수 | 검정 전체의 신뢰도 |
|---|---|
| 1 | $0.99^1=0.99$ |
| 2 | $0.99^2=0.98$ |
| 10 | $0.99^{10}=0.90$ |
| 100 | $0.99^{100}=0.37$ |
이를 그래프로 나타내면 아래와 같습니다.
FDR은 이러한 다중가설검정 문제를 해결하기 위해 제안된 방법입니다. 물론 개별 유의수준($α$)을 ‘목표값/변수개수’로 정하고 기존대로 가설검정을 수행하는 기법(Family-wise Test)이 제시되기도 했습니다. 예컨대 전체 목표 유의수준이 0.01이고 변수개수가 10개라면 개별 유의수준을 0.001로 놓고 가설검정을 각각 실시하는 방식이죠. 그러나 이 방식은 중요한 변수($H_0$가 거짓)인데도 중요하지 않다($H_0$가 참)는 결론을 내는 경우가 많은 단점이 있다고 합니다.
FDR 절차
FDR은 제1종 오류를 최대한 적게 범하면서도 중요한 가설을 최대한 많이 찾을 수 있도록($H_0$ 기각) 한다는 점이 강점입니다. 검정해야할 다중가설(데이터의 변수 개수)이 $m$개일 때 FDR은 아래와 같이 정의됩니다.
| 구분 | $H_0$ 채택 | $H_0$ 기각 | Total |
|---|---|---|---|
| $H_0$가 참 | $U$ | $V$ | $m_0$ |
| $H_0$가 거짓 | $T$ | $S$ | $m_1$ |
| Total | $W$ | $R$ | $m$ |
위 수식에 정의된 FDR을 해석하면 이렇습니다. FDR은 기각된 귀무가설($H_0$) 가운데 잘못 기각된 가설이 차지하는 비율의 평균으로 사용자가 지정해주는 값입니다. 바꿔 말해 $m$회 가설검정 가운데 각 $H_0$이 참인데 기각된 비율, 즉 검정 전체 유의수준에 가까운 의미라는 뜻입니다. 실제로 FDR 제안자는 FDR을 $α$라고 표기하고 있습니다.
그럼 FDR 절차를 살펴볼까요? 검정해야할 다중가설, 즉 데이터 변수 개수가 15개라고 합시다. 목표로 하는 FDR 값을 0.05로 정해보겠습니다. 우선 각각의 가설에 대해 $p-value$를 구한 뒤 이를 오름차순으로 정렬하고, 순위(rank)를 매깁니다. 아래 표와 같습니다.
$p-value$는 $H_0$이 참일 확률을 의미하므로, $α$를 0.05로 놓고 기존대로 단일 가설검정을 실시했다면 1~9번 변수 가 중요변수라는 결과를 얻었을 겁니다. Family-wise Test 방식대로 ‘목표 유의수준/변수 개수=0.05/15’로 가설검정을 실시한다면, 그 값이 0.0033이므로 1~3번 변수만 중요하다는 결론을 내렸을 겁니다.
FDR 절차에서 $i$는 rank, $m$은 변수 개수, $α$는 사용자가 정한 FDR 값입니다. 예컨대 1번 변수는 $1/15$에 0.05를 곱해 계산합니다. 마찬가지로 2번 변수는 $2/15$에 0.05를 곱해 구합니다.
이 값이 $p-value$보다 큰 변수 랭크의 최대값이 FDR 절차가 뽑은 중요변수의 개수가 됩니다. 위 표 기준으로 설명하자면 이렇습니다. 표 아래쪽부터 시작해 위쪽으로 훑어 가면서 $i/m*α$ 값이 $p-value$보다 큰 지 여부를 검토합니다. 처음으로 $p-value$ 값을 넘어서는 지점은 바로 4번 변수이네요. 따라서 1~4번 변수를 중요한 변수라고 결론을 내게 됩니다.
예제에서도 알 수 있듯 FDR은 단일 가설검정에 비해 제1종 오류가 낮고(비중요 변수 배제), Family-wise Test 방식에 비해 중요한 가설을 많이 찾아낸다는 장점이 있습니다.
FDR 예제
그럼 실제 데이터를 가지고 문제를 풀어보겠습니다. UCI data 가운데 하나를 가져왔습니다. 개요는 아래와 같습니다.
Ionosphere Data Set
데이터 설명 : 전리층(ionosphere)에 대한 레이더 관측치를 모은 데이터
관측치 개수 : 351개
독립변수 : 34개
종속변수 : 전리층의 패턴 유무 (good : 전리층에 특정 패턴 존재, bad : 패턴 없음)
모든 독립변수에 대해 T-test를 수행하였습니다. ($H_0$=두 집단의 평균이 같다, $H_1$=두 집단 평균이 다르다) 변수별 p-value는 아래 표와 같습니다.
| 변수 | P-value | 변수 | P-value | 변수 | P-value | 변수 | P-value |
|---|---|---|---|---|---|---|---|
| 1 | 2.303168e-11 | 11 | 3.060316e-03 | 21 | 6.897335e-05 | 31 | 7.025167e-07 |
| 2 | NaN | 12 | 8.189137e-03 | 22 | 5.051969e-02 | 32 | 5.547445e-01 |
| 3 | 2.339468e-15 | 13 | 9.549388e-04 | 23 | 3.350101e-04 | 33 | 3.748258e-06 |
| 4 | 5.941660e-02 | 14 | 1.341207e-03 | 24 | 9.191420e-01 | 34 | 2.747245e-01 |
| 5 | 3.153540e-15 | 15 | 1.017555e-04 | 25 | 1.162854e-03 | ||
| 6 | 1.988959e-02 | 16 | 1.507825e-02 | 26 | 9.795861e-01 | ||
| 7 | 7.577585e-13 | 17 | 1.036061e-01 | 27 | 5.055064e-02 | ||
| 8 | 6.995369e-04 | 18 | 4.976238e-02 | 28 | 4.980595e-01 | ||
| 9 | 5.545817e-07 | 19 | 3.255247e-02 | 29 | 2.570992e-05 | ||
| 10 | 4.515815e-02 | 20 | 5.604890e-01 | 30 | 9.475607e-01 |
Individual Test($α=0.05$) 결과 선택된 중요 변수는 22개입니다. (p-value가 0.05보다 작은 변수 선택)
1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 28, 30, 32
Family-Wise Test($α=0.05$, $α/m$, $m$=전체 변수 개수) 결과 선택된 중요 변수는 15개입니다. (p-value가 0.0015보다 작은 변수 선택)
1, 3, 5, 7, 8, 9, 13, 14, 15, 21, 23, 25, 29, 31, 33
위 계산결과를 토대로 $α$값을 0.05로 두고 FDR을 수행한 결과는 아래 표와 같습니다.
| 변수 ID | rank | P-value | $i/m*α$ |
|---|---|---|---|
| 3 | 1 | 2.339468e-15 | 0.001470588 |
| 5 | 2 | 3.153540e-15 | 0.002941176 |
| 7 | 3 | 7.577585e-13 | 0.004411765 |
| 1 | 4 | 2.303168e-11 | 0.005882353 |
| 9 | 5 | 5.545817e-07 | 0.007352941 |
| 31 | 6 | 7.025167e-07 | 0.008823529 |
| 33 | 7 | 3.748258e-06 | 0.010294118 |
| 29 | 8 | 2.570992e-05 | 0.011764706 |
| 21 | 9 | 6.897335e-05 | 0.013235294 |
| 15 | 10 | 1.017555e-04 | 0.014705882 |
| 23 | 11 | 3.350101e-04 | 0.016176471 |
| 8 | 12 | 6.995369e-04 | 0.017647059 |
| 13 | 13 | 9.549388e-04 | 0.019117647 |
| 25 | 14 | 1.162854e-03 | 0.020588235 |
| 14 | 15 | 1.341207e-03 | 0.022058824 |
| 11 | 16 | 3.060316e-03 | 0.023529412 |
| 12 | 17 | 8.189137e-03 | 0.025000000 |
| 16 | 18 | 1.507825e-02 | 0.026470588 |
| 6 | 19 | 1.988959e-02 | 0.027941176 |
| 19 | 20 | 3.255247e-02 | 0.029411765 |
| 10 | 21 | 4.515815e-02 | 0.030882353 |
| 18 | 22 | 4.976238e-02 | 0.032352941 |
| 22 | 23 | 5.051969e-02 | 0.033823529 |
| 27 | 24 | 5.055064e-02 | 0.035294118 |
| 4 | 25 | 5.941660e-02 | 0.036764706 |
| 17 | 26 | 1.036061e-01 | 0.038235294 |
| 34 | 27 | 2.747245e-01 | 0.039705882 |
| 28 | 28 | 4.980595e-01 | 0.041176471 |
| 32 | 29 | 5.547445e-01 | 0.042647059 |
| 20 | 30 | 5.604890e-01 | 0.044117647 |
| 24 | 31 | 9.191420e-01 | 0.045588235 |
| 30 | 32 | 9.475607e-01 | 0.047058824 |
| 26 | 33 | 9.795861e-01 | 0.048529412 |
| 2 | 34 | NaN | 0.050000000 |
위 표에서 $i/m*0.05$가 $p-value$보다 크다는 조건을 만족하는 rank의 최대값은 19입니다. 그 결과 19개 변수가 선택됐습니다.
3, 5, 7, 1, 9, 31, 33, 29, 21, 15, 23, 8, 13, 25, 14, 11, 12, 16, 6
FDR level 변화에 따른 중요 변수 선택 결과는 아래와 같습니다. FDR level이 클 수록 선택되는 변수 개수가 많아지는 경향을 보였습니다.
| FDR level | 선택변수(개) |
|---|---|
| 0.01 | 16 |
| 0.02 | 17 |
| 0.03 | 18 |
| 0.04 | 19 |
| 0.06 | 20 |
| 0.08 | 24 |
| 0.09 | 25 |
| 0.14 | 26 |
| 0.35 | 27 |
| 0.61 | 28 |
| 0.64 | 30 |
Ionosphere Data Set을 7:3 비율로 학습셋과 검증셋으로 나누고, 학습셋에 대해 LDA 모델을 적합한 뒤 검증셋을 대상으로 단순정확도를 측정하였습니다. 통계적으로 유의한 일부 변수만으로 구축한 모델은 전체 변수를 모두 사용한 모델의 성능에 크게 뒤지지 않음을 확인할 수 있었습니다.
| 구분 | 전체 | Individual | Family-wise | FDR |
|---|---|---|---|---|
| 학습에 사용된 변수(개) | 33 | 19 | 9 | 14 |
| 단순정확도 | 0.8762 | 0.8476 | 0.8286 | 0.8190 |
위 문제를 푸는 데 사용한 R 코드는 아래와 같습니다.
# load & split data
data <- read.csv('data.csv',header=F)
good <- data[which(data[,35] == 'g'),-35]
bad <- data[which(data[,35] == 'b'),-35]
# T-test
Ttest <- c()
for (i in 1:dim(good)[2]) {
Ttest[i] <- t.test(good[,i], bad[,i])$p.value
}
# individual test
which(Ttest<0.05)
# family-wise test (패키지 사용)
# bonferroni : p-values are multiplied by the number of comparisons
family.wise.test=p.adjust(Ttest, method = "bonferroni")
which(family.wise.test<0.05)
# family-wise test (자작)
which(Ttest < 0.05/length(Ttest))
# FDR (패키지 사용)
library(fdrtool)
fdr <- fdrtool(Ttest, statistic='pvalue', cutoff.method = 'fndr')
fdr$param
# FDR (자작)
fdr.level <- 0.05
var_table <- matrix(0, nrow=length(Ttest),ncol=4)
colnames(var_table) <- c('id', 'rank', 'P-value', 'i/m*FDR')
var_table[,1] <- c(1:length(Ttest))
var_table[,3] <- Ttest
var_table <- var_table[order(var_table[,3]),]
var_table[,2] <- c(1:length(Ttest))
var_table[,4] <- var_table[,2] / length(Ttest) * fdr.level
max(which(var_table[,3] < var_table[,4])) # 중요 변수 개수
var_table[1:max(which(var_table[,3] < var_table[,4])),1] # 중요 변수 추출
# all variables model
library(MASS)
data[,35] <- as.numeric(data[,35])
colnames(data)[35] <- 'Target'
data <- data[,-2]
trn_idx <- sample(1:dim(data)[1], round(0.7*dim(data)[1]))
train_data <- data[trn_idx,]
test_data <- data[-trn_idx,]
allmodel <- lda(Target~.,train_data)
# var selection model
train_data_good <- train_data[which(train_data$Target == 1),-train_data$Target]
train_data_bad <- train_data[which(train_data$Target == 2),-train_data$Target]
Ttest <- c()
for (i in 1:dim(train_data_good)[2]) {
Ttest[i] <- t.test(train_data_good[,i], train_data_bad[,i])$p.value
}
family.wise.test=p.adjust(Ttest, method = "bonferroni")
individual_model <- lda(Target~.,train_data[,c(which(Ttest<0.05),34)])
familywise_model <- lda(Target~.,train_data[,c(which(family.wise.test<0.05),34)])
# fdr model
fdr.level <- 0.05
var_table <- matrix(0, nrow=length(Ttest),ncol=4)
colnames(var_table) <- c('id', 'rank', 'P-value', 'i/m*FDR')
var_table[,1] <- c(1:length(Ttest))
var_table[,3] <- Ttest
var_table <- var_table[order(var_table[,3]),]
var_table[,2] <- c(1:length(Ttest))
var_table[,4] <- var_table[,2] / length(Ttest) * fdr.level
important_vars <- var_table[1:max(which(var_table[,3] < var_table[,4])),1]
fdr_model <- lda(Target~.,train_data[,c(important_vars,34)])
# prediction
all <- predict(allmodel, test_data)
individual <- predict(individual_model, test_data)
familywise <- predict(familywise_model, test_data)
fdr <- predict(fdr_model, test_data)
# calculate ACC
ACC_all <- sum(diag(table(test_data$Target, all$class)))/length(all$class)
ACC_indi <- sum(diag(table(test_data$Target, individual$class)))/length(individual$class)
ACC_fdr <- sum(diag(table(test_data$Target, fdr$class)))/length(fdr$class)
ACC_family <- sum(diag(table(test_data$Target, familywise$class)))/length(familywise$class)
반복샘플링 기반 FDR
지금까지 설명해드린 FDR 기법은 우리가 가진 데이터가 특정 확률분포를 따르고, 이로부터 적당한 방법을 취해 $p-value$를 계산할 수 있어야만 적용이 가능합니다. 하지만 데이터가 특정 분포를 따르지 않거나, 그 분포를 모른다고 할 때는 $p-value$ 계산이 어렵습니다. 반복샘플링 기반 FDR은 말 그대로 반복샘플링을 통해 분포에 대한 가정없이 $p-value$를 계산할 수 있는 기법입니다.
가정과 절차는 이렇습니다. 특정 변수를 기준으로 관측치를 랜덤하게 섞습니다. 만약 해당 변수가 범주 분류에 의미 있는 변수라면 실제 데이터의 통계량이 관측치를 랜덤하게 섞은 데이터의 통계량보다 훨씬 클 겁니다. 아래 그림처럼요.
통계적 유의성을 확보하기 위해서는 이러한 과정을 충분히 많이 반복해야 합니다. 만약 아래처럼 100번을 반복 수행했다고 가정해보겠습니다. p-value는 아래와 같이 구할 수 있습니다.
\({ p-value }_{ { X }_{ 1 } }=\frac { { T }_{ { x }_{ 1 } }^{ (i) }\quad 중에 \quad{ T }_{ { x }_{ 1 } } 보다\quad 큰\quad값의\quad개수}{ 반복횟수\times변수의\quad수 } =\frac { { T }_{ { x }_{ 1 } }^{ (i) }\quad중에 \quad{ T }_{ { x }_{ 1 } } 보다\quad 큰\quad값의\quad개수 }{ 100\times 2 } \\ { p-value }_{ { X }_{ 2 } }=\frac { { T }_{ { x }_{ 2 } }^{ (i) }\quad 중에\quad{ T }_{ { x }_{ 2 } } 보다\quad 큰\quad값의\quad개수 }{ 반복횟수\times변수의\quad수} =\frac { { T }_{ { x }_{ 2 } }^{ (i) } \quad중에\quad{ T }_{ { x }_{ 2 } } 보다\quad 큰\quad값의\quad개수 }{ 100\times 2 }\)
이렇게 p-value를 구하게 되면 이후 절차는 지금까지 설명해드린 FDR 절차 그대로 따르면 됩니다.
이진분류 이외 문제에 FDR 적용하기
p-value만 구하면 이후엔 FDR 기법을 적용할 수 있습니다. 정규분포를 따르는 데이터의 이진분류 문제의 경우 T-test를 수행하면 p-value를 구할 수 있었고, 특정 데이터 분포를 가정할 수 없는 상황이라면 반복 샘플링 기반으로 p-value를 구할 수 있습니다.
그렇다면 $y$가 연속형 숫자인 회귀 문제에서 중요 변수를 찾아내기 위해 FDR을 적용해야 한다면 p-value를 어떻게 구해야 할까요? 귀무가설과 대립가설을 아래와 같이 설정하고 모든 변수에 대해 가설검정을 실시합니다.
\[{ H }_{ 0 }\quad :\quad { \beta }_{ i }=0\\ { H }_{ 1 }\quad :\quad { \beta }_{ i }\neq 0\]p-value는 귀무가설이 맞을 확률을 의미하므로 위 가설검정에서 p-value가 작게 나올 경우 $H_0$는 기각하게 됩니다. 바꿔 말해 $β_i$가 중요하다는 이야기이죠. 선형회귀 모델의 계수들은 t분포를 따른다고 합니다. 검정통계량 $t$는 아래와 같이 구합니다. \(t=\frac { \hat { { \beta }_{ i } } -0 }{ sd\{ \hat { { \beta }_{ i } } \} }\) 3개 이상의 범주 분류 문제의 경우 p-value를 어떻게 구할까요? 두 가지 방법이 있습니다. 첫번째는 범주 두 개를 한 세트로 묶어서 이진 분류 문제처럼 p-value를 구하는 겁니다. 나머지는 분산분석(ANOVA) 기법을 적용하는 것입니다.
Comments