Hard skills/Statistic (43)

층화 무작위배정을 실제로 하는 방법


본 포스팅에서는 임상시험에서 층화 무작위배정 을 실제로 어떻게 하는지에 대해 알아보겠습니다.층화 무작위배정 은 공변량에 대해 층화하여 무작위 배정을 하는 방법입니다. 층화 무작위배정 을 하는 이유는 1) 해당 공변량 대해 Randomization 되어있도록 해서 bias 를 줄이고자 하는 것입니다. 이는 treatment/control group 간에 공변량의 분포가 다른 경우, bias 가 생기기 때문입니다. 하지만 공변량에 대해 treatment/control group 이 randomization 되어있다면, 해당 이론적으로는 공변량에 대해서는 분포가 다를 수 없어 이로 인한 bias 가 생기지 않게 됩니다. 2) 또한 공변량을 보정하면 공변량의 effect 를 제외하여 treatment effect size 의 참값을 더욱 잘 추정할 수 있습니다. 이를 통해 검정력을 증가시킬 수 있습니다. 


다음과 같은 상황을 예로 들어 봅시다. 


1. 총 4개의 기관이 임상 시험에 참여

2. 기관당 임상시험 참여자는 20명

3. 나이를 공변량으로 잡고 Block size 는 6,6,6,2 로 설정한다. 

구체적으로 기관 별로 아래와 같은 

- 20~30세 : 6명

- 30~40세 : 6명

- 40~50세 : 6명

- 50~60세 : 2명


4. treatment/control group 의 배정 비율은 1:1 이다.

즉, Block size가 6이면 treatment 3명, control 3명을 무작위 배정한다.


전략


1. 기관별로 층화 -> 똑같은 코드를 4번 돌리는 것으로 구현 

2. 나이를 사이즈 6,6,6,2 로 층화 무작위 배정 -> for 문으로 구현

3. 블록별 무작위 배정을 위해 블록 별로 0~1사이의 난수를 발생하고 (이중 for 문으로 구현 ), 이 값을 기준으로 sorting 해서 중위수를 기준으로 두 그룹으로 나눈 후, A, B 약을 배정한다. 

blockrand <- function(R, N, blocksize){ result <- data.frame(block=numeric(N), item=numeric(N), drug=character(N), rand=numeric(N), random_no=character(N), stringsAsFactors=FALSE) index <- 1 for (i in c(1:ceiling(N/blocksize))){ if (i == (ceiling(N/blocksize))){ blocksize <- N - (i-1)*blocksize # 2 } for (j in c(1:blocksize)){ result[index, 'block'] <- i result[index, 'item'] <- j if (j <= blocksize/2){ result[index, 'drug'] <- 'A' } else { result[index, 'drug'] <- 'B' } result[index, 'rand'] <- runif(1) index = index + 1 } } result <- result[order(result$block, result$rand), ] result$random_no <- c(paste0(R,'-R', as.character(c(1:(ceiling(N/blocksize)*blocksize))))) return (result) } inst1 <- blockrand('01', 20, 6) inst2 <- blockrand('02', 20, 6) inst3 <- blockrand('03', 20, 6) inst4 <- blockrand('04', 20, 6)

이를 위하여 blockrand 라는 함수를 만듭니다. 이 함수는 기관명 R, 기관별 대상자 수 N, 블록 사이즈 blocksize 를 parameter 로 받아 결과 테이블을 dataframe 형식으로 생성하여 반환하는 함수입니다. 


결과 예시

 

block item         drug       rand random_no

3      1    3 A 0.03323463     01-R1

5      1    5   B 0.23377875     01-R2

1      1    1 A 0.58872123     01-R3

4      1    4   B 0.64914300     01-R4

2      1    2 A 0.86294140     01-R5

6      1    6   B 0.86923215     01-R6

8      2    2 A  0.01529889     01-R7

9      2    3 A 0.11754757     01-R8

11     2    5   B 0.13073496     01-R9

10     2    4   B 0.46934007    01-R10

7      2    1 A  0.56171635    01-R11

12     2    6   B 0.61339224    01-R12

14     3    2 A 0.11352796    01-R13

13     3    1 A 0.62667742    01-R14

15     3    3 A 0.82257959    01-R15

18     3    6   B 0.85910418    01-R16

17     3    5   B 0.87995944    01-R17

16     3    4   B 0.92778095    01-R18

20     4    2   B 0.70824728    01-R19

19     4    1 A 0.76212643    01-R20


이러한 층화 무작위배정은 머신러닝에서 특정 변수의 분포를 train/valid/test set 에서 맞추어줘서 모델의 일반화 성능을 보다 정확히 평가하고 싶은 경우에 사용되기도 합니다. 


  • 2021.05.25 11:17

    비밀댓글입니다

다중 비교 (Multiple comparison) 문제와 보정 방법 


본 포스팅에서는 다중 비교가 무엇인지, 어떤 방법으로 보정할 수 있는지를 알아보겠습니다. 다중 분석이 문제가 된다는 점을 이해하기 위해 먼저 예를 설명하겠습니다.폐암과 연관성이 있는 요인을 찾기 위한 연구를 수행한다고 해봅시다. 예를 들어, 나이, 체질량 지수, 식습관, 음주습관, 운동량 등 100 가지의 요인에 대해 폐암과의 연관성을 유의 수준 0.05 하에서 비교하고자 합니다. 이를 위해 폐암 환자와 정상인을 수집했고, 연속형 변수에 대해서는. T-test, 범주형 변수는 카이제곱 검정을 수행하여 최종적으로, 운동량, BMI, 학력, 우유 섭취량, 경제력이 폐암과 연관성이 있다고 나왔다고 합시다. 이 때, 이 변수들이 폐암과의 연관성이 있다고 결론 내릴 수 있을까요? 물론 각각의 검정은 유의 수준 0.05 하에서 검정한 것이기 때문에, 연관성이 실제로 있을 가능성이 높은 것이 사실입니다. 하지만 100가지 요인에 대한 검정을 하나의 연구로 볼 때, 모든 변수가 연관성이 실제로 없더라도 평균적으로 5개의 잘못된 결과를 얻게 됩니다. 만약에 이러한 잘못된 결론이 학계에 보고되면 잘못된 과학적 결론을 내릴 수 있습니다. 따라서, 개별 검정 뿐 아니라 한 연구의 에러를 줄이는 방법이 필요합니다. 즉, 한 연구에서 나올 수 있는 잘못된 결과를 줄이자 라는 것이 다중검정 보정의 핵심입니다. 


여기서 family-wise type 1 error 의 개념이 등장합니다. Family-wise type 1 error rate (FWER) 란 한 연구에서 적어도 한 개의 잘못된 결론 (false positive)이 나올 수 있는 확률을 의미합니다. 만약 이것이 0.05 라면 한 연구에서 적어도 1개의 잘못된 결론이 나올 확률이 0.05 라는 것입니다. FWER 를 통제하는 방법으로 잘 알려진 것이 본페로니 보정 방법입니다. 본페로니 보정은  test 의 수가 n 이고, FWER 를 0.05 로 통제할 때, 개별 테스트의 유의수준을 alpha/m 으로 설정시킵니다. 모든 검정이 실제로 연관성이 없는 경우 (null 인 경우), 아래 식이 m 이 클 수록 대략적으로 만족됩니다. 


$$ FWER= 1-(1-\alpha/m)^m = \alpha  $$


<본페로니 방법을 통한 FWER 통제>


하지만 본페로니 방법의 문제는 너무 보수적 (strict/conservative) 이라는 것입니다. 보수적이라는 뜻은 귀무가설을 웬만해선 기각하지 않는다는 것으로, false positive 는 상당히 줄일 수 있지만, false negative 는 많아지게 됩니다. 예를 들어 test 의 개수가 1000개인데, 연관성이 없는 것이 900개, 연관성이 있는 것이 100개 일 때, 900개 중에 1개라도 잘못나올 확률은 1-(1-0.05/1000)^900 < 0.05 입니다. 물론 false positive 가 적어져서 좋긴 하지만, false negative 가 많아진다는 문제가 생깁니다. 따라서, 다중 비교 검정의 핵심은 어떻게 false positive 를 줄이면서, false negative 도 줄일 수 있는가? 입니다.


다중 검정 보정 방법


그러면 어떻게 false positive, false-negative 의 타협을 보는가 (FWER 은 본페로니 수준이면서, false negative 를 줄일 수 있는가)에 문제에 있어서 도입되는 한 가지 방법이 바로 multi-step 보정 방법입니다. Multi-step 보정법은 step-down 방법, step-up 방법으로 나뉘는데, 예를 들어, step-down 방법은 Holm's 방법, step-up 방법은 Hochberg 방법이 있습니다. 이러한 multi-step 보정 방법은 fwe 는 그대로 두면서 false negative 는 줄이는 방법으로 알려져 있습니다. 이를 수리적으로 보이는 것은 다소 어렵지만 simulation 을 통해서 bonferonni 방법보다 효율적이라는 것이 많이 알려져 있습니다. Multi-step 보정 방법은 모든 검정에서 나온 p-value 를 정렬 (sorting) 한 후, 각 검정마다 각기 다른 p-value cutoff 를 적용시키는 방법입니다. Step-down 방법은 p-value 가 가장 작은 검정부터, step-up 방법은 p-value 가 큰 검정부터 귀무가설 기각 여부를 보게됩니다. 


FWER 이 아닌 False Discovery Rate (FDR) 을 줄이는 방법이 최근 많이 사용되고 있습니다. FDR 을 다중검정에서 사용할 때의 의미는 false postive, false negative 를 줄이는데 집중할 것이 아니라, 내가 귀무가설을 기각한 검정 중 틀린 것 (false positive/true positive+false positive = discoveries) 의 비율을 줄이자는 것입니다. 최근 FDR 통제에 많이 쓰이는 방법 중 하나가 Benjamin-Hochberg 방법입니다. 이 방법은 FDR 를 잘 통제한다고 알려진 방법입니다. 


$$ FDR = \frac{false_{positive}}{true_{positive}+false_{positive}} $$


Holm's 방법, Hochberg 방법, Benjamin-Hochberg 방법 모두 간단하며, 검색해보면 쉽게 알 수 있기 때문에 본 포스팅에서는 따로 설명하지 않고, 아래 시각적으로 설명한 그림을 첨부하였습니다. Benjamin-Hochberg 방법의 경우 이전 제 포스팅에서 설명했습니다.



 

<Holm's, Hochberg, Benjamin-Hochberg Procedure 의 비교>



다중 검정 보정이 필요한 상황


1. 수많은 요인들에 대한 연관성 분석을 수행할 때


이 예로는 유전체학 분야를 들 수 있습니다. genomics 분야에서 microarray 와 같은 high-throughput 기술의 발달됨에 따라, 수많은 유전자 변이 마커, 유전자 발현 (gene expression) 과 표현형의 연관성을 보는 연구가 수행되었습니다. 이러한 Genomics 분야의 발전은 다중 검정 보정이 생겨난 이유이기도 합니다. 일반적으로 한 사람에서 300만개의 단일염기다형성 (single-nucleotide polymorphism) 이 있습니다. 이러한 마커들과 질병의 연관성을 보는 연구에서 multiple comparison 문제가 생길 수 밖에 없고, 수 많은 false positive 가 생기게 됩니다. 이 분야의 default 라고 얘기할 수 있는 보정 방법은 FDR 을 통제하는 Benjamin-Hochberg 방법입니다. 


2. 임상 시험 


임상 시험의 수많은 상황에서 다중 비교 문제가 생기게 됩니다. 신약 개발을 예를 든다면, 아래와 같은 상황이 발생할 수 있습니다. 


1) 비교 그룹 수가 3개 이상인 경우 :ANOVA 검정 후 유의해서, 사후검정을 할 때, 다중 검정이 발생하게 됩니다. 예를 들어, 대조약, 저투여군, 고투여군 3그룹을 비교할 때 생길 수 있습니다. 이러한 임상시험 세팅에서는 다중성을 반드시 보정해주어야합니다. 


2) 하위군 분석 : 임상 시험 대상자를 어떤 공변량을 기준으로 하위군으로 나누어 신약의 효과를 평가할 때 발생합니다. 대표적으로 다중 검정이 발생하는 사례로, 하위군 분석 결과를 확증적으로 인정받기 위해서는 임상시험계획서에 이를 명기하고 다중검정 보정법을 제시해야합니다. 


3) 중간 분석 : 임상 시험에서 시험 기간이 종료되기 전에 중간에 데이터를 오픈해서 테스트를 하는 경우도 있습니다. 이 경우도 대표적으로 다중 검정이 발생하는 사례로, 각 절차에서 O-brien fleming 방법을 사용한 다중검정 보정법이 종종 사용됩니다.

변수 종류별 시각화 및 검정 방법


얼마전, 변수 종류별로 사용할 수 있는 시각화 및 검정 방법을 간단하게 요약한 표를 발견해 공유합니다. 출처는 이곳 입니다. 표를 그린 기본 아이디어를 보면, 종속변수를 Binary, Nominal, Ordinal, Interval, Normal, Censored Interval 으로 나누었습니다. 독립변수의 경우, Binary Categorical (Paired/Unpaired), Categorical (Paired/Unpared), Normal, Multivariate 으로 나누었습니다. 이 때, ordinal 의 경우, 연속형 범주로 또는 범주형 변수 상황에 따라 둘 다 가능합니다.


그래프의 경우 요약하면, Barplot, Boxplot, Scatter plot 이 기본입니다. 카테고리x카테고리 = Barplot, 카테고리x연속형 = Box plot, 연속형x연속형 = Scatter plot 으로 기본적으로 그리면 됩니다. 여기서 Box plot 의 경우 전통적인 통계에서 많이 쓰이는 그래프이지만, Histogram with different colors or side by side 또는 violin plot 도 선호됩니다. 


대부분의 기초 통계에서 배우는 검정은 대부분 이 표안에 속해있습니다. 검정에서 하고자 하는 것은 변수 X와 변수 Y 가 연관성이 있는가? 입니다. 검정도 마찬가지로, X와 Y의 종류에 따라 다양하게 존재합니다. 아래표의 대부분의 검정은 기초 통계에서 배우는 내용입니다. Paired data 에 적용하는 맥니마 검정은 설문지/심리학 연구 등에 자주 사용되는 통계적 검정 방법인데 이전 제 포스팅에서 다루었습니다. 코크란의 Q test는 맥니마 검정의 확장으로 종속변수가 두 개 이상일 때 사용하는 통계적 검정 방법입니다. 


출처

https://www.r-bloggers.com/overview-of-statistical-tests/

통계적 검정의 종류


통계적 검정에 있어 p-value 를 보고 기각 여부를 판단하는 방법도 있지만, 소위 '신뢰구간이 0 을 포함하는지' 를 보고 가설검정을 할 수도 있다. 왜 신뢰구간이 0을 포함하지 않으면 두 그룹의 차이가 유의하다고 볼 수 있을까? 본 포스팅에서는 우월성 검정과 동등성 검정에 대하여 신뢰구간과 양측 가설검정이 어떤 관계를 갖고 있는지 설명한다. 


1. 우월성 검정 (Superiority Test)


우월성 검정은 A가 B보다 우월함을 보이는 검정으로 일반적인 통계적 검정 (two-sample t-test) 이다. two-sample t-test 에 대한 기초내용은 이전 포스트에서 정리하였다. 


$$ H_0 : \mu_1 = \mu_2 $$

$$ H_0 : \mu_1 \ne \mu_2 $$


우월성 검정을 단측검정이 아니라 양측검정으로 하는 이유는 신뢰구간과 검정의 결과를 동일하게 맞추기 위해서이다. 예를 들어, A 약이 B 약보다 우월함을 보이고 싶을 때 검정을 통해서 구한 p-value 만 제시하는 것보다 구간 추정값을 제시하는 것이 더욱 많은 정보를 제공한다. 따라서 양측검정을 하면, 구간 추정을 통한 검정과 같은 판단을 내릴 수 있기 때문에 양측검정이 더 선호된다. 일반적으로 우월성 검정의 결과로 양측검정에 대한 p-value 와 함께 95 % 신뢰 구간을 제시된다. 


양측검정과 95 % 신뢰구간이 0을 포함하는지 여부가 같은 판단을 내리는 이유 


양측검정의 p-value


합동 분산을 이용한 two-sample t-test 를 예로 들어보자. p-value는 귀무가설하에서 더 극단적인 데이터를 관찰할 확률이다. 여기서 극단적인 확률이란 표본의 통계량이 대립가설을 지지하는 쪽으로 나오는 것을 말한다. 따라서 부등호의 방향이 [X > 통계량] 임을 알 수 있다. (단 통계량이 양수인 경우)


$$ P= 2 P[X > \frac{\tilde X_1 - \tilde X_2 - (\mu_1-\mu_2)}{\sqrt{s^2_p(1/n_1+1/n_2)}}] $$

$$ X \sim t(n_1+n_2-2) $$


위 식에서 p-value는 귀무가설 하에서라는 정보를 이용하면 mu_1 - mu_2 = 0 이고, 이를 통해 p-value 를 계산할 수 있다.  


95% 신뢰구간


$$ (\tilde X_1 - \tilde X_2 - t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)} ,  \tilde X_1-\tilde X_2 + t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)})$$


이 때, 통계량의 차이가 양수이고, 0을 포함하지 않는 경우는 구간 추정값의 left bound 가 0 보다 큰 경우인데, p-value 는 0.05 보다 작음을 알 수 있다. 통계량의 차이가 음수인 경우도 마찬가지로 right bound 가 0보다 작을 때의 p-value 가 0.05 보다 작기 때문에 동일하다. 이를 식으로 설명하면 아래와 같다. 


$$ If, \tilde X_1 - \tilde X_2 - t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)} > 0 $$

$$ P < 2P[X > t_{0.025}] = 0.05$$ 


2. 동등성 검정 (equivalence test) 


동등성검정은 A와 B가 의미있는 차이가 없음을 보이는 목표를 갖는 검정이다. 동등성 검정에서는 delta 라는 미리 정의된 값을 이용한다. 


$$ H_0 : \mu_1-\mu_2 \ge \delta $$

$$ H_1 : \mu_1-\mu_2 < \delta $$


동등성 검정 결과는 95 % 신뢰구간이 (-delta, delta) 안에 들어오는지 여부와 같다. 신뢰구간이 delta 사이에 들어오면 양측검정의 귀무가설을 기각할 수 있다. 왜 이렇게 되는지 아래에서 설명한다. 


양측검정의 p-value


우월성검정과 비교해 부등호의 방향을 주의해야 한다.


$$ P= 2 P[X < \frac{\tilde X_1 - \tilde X_2 - (\delta)}{\sqrt{s^2_p(1/n_1+1/n_2)}}]  $$

$$ X \sim t(n_1+n_2-2) $$


95 % 신뢰구간


95 % 신뢰구간은 검정과 독립적으로 나오는 값이기 때문에 1. 과 같다. 


$$ (\tilde X_1 - \tilde X_2 - t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)} ,  \tilde X_1-\tilde X_2 + t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)})$$


만약, 위 신뢰구간의 right bound 가 delta 보다 작다면, left bound 도 자연스레 -delta 보다 클 것이다. (신뢰구간은 대칭적이기 때문에) 따라서 아래식을 가정했을 때, p-value가 0.05 보다 작음을 보이면 된다.


$$ If, \tilde X_1 - \tilde X_2 + t_{0.025}\sqrt{s^2_p(1/n_1+1/n_2)} < \delta $$

$$ P < 2P[X < -t_{0.025}] = 2P[X > t_{0.025}] = 0.05 $$


따라서 동등성 검정도 따로 통계량을 구하지 않고도, 신뢰구간을 통해 바로 검정할 수 있다는 것을 알 수 있다. 

모분산의 추정과 검정


어떤 공장에서 일정 크기의 볼트를 만든다고 해보자. 제조 공정에서 볼트의 지름이 일정해야 품질 기준을 만족할 수 있다. 볼트의 지름이 10mm 여야하는데, 평균적으로는 10mm 이지만, 분산이 존재해  불량품이 일정 부분 존재한다. 이런 상황에서 볼트 지름의 분산을 확인하는 방법으로 모분산에 대한 추론을 할 수 있다. 


모분산의 추정에는 카이제곱 분포를 이용한다. 카이제곱분포는 표준정규분포에서 n개의 샘플을 뽑아서 그것들의 제곱을 모두 합한 것을 하나의 확률 변수 Q 로 정의하자. 카이제곱분포는 Q 가 따르는 분포이다. 카이제곱 분포는 하나의 모수 n 을 갖고, 이를 자유도라 한다.


$$ Q = \sum^{n}_{i=1} Z_i $$

$$ Q \sim \chi (n) $$


이 때,  아래와 같은 확률 변수를 정의하자. 


$$ X \sim N(\mu, \sigma^2) $$


모분산은 아래 통계량이 자유도가 n-1 인 통계량을 따른다는 것으로부터 추정된다. 


$$ \frac{(n-1)s^2}{\sigma^2} \sim \chi(n-1) $$


증명은 아래와 같다. 


$$ \sum^{n}_{i=1} \frac{(X_i - \mu)^2}{\sigma^2} = \frac{(n-1)s^2}{\sigma^2} + \frac{(\tilde X-\mu)^2}{\sigma^2/n} $$

$$ \sum^{n}_{i=1} \frac{(X_i - \mu)^2}{\sigma^2} \sim \chi(n) $$

$$ \frac{(\tilde X-\mu)^2}{\sigma^2/n} \sim \chi(1) $$


따라서 chi-square 분포의 additivity 에 의해 아래와 같다. 


$$ \frac{(n-1)s^2}{\sigma^2} \sim \chi(n-1) $$ 


한 가지 유의할점은 모집단 X 가 정규분포를 따라야 이러한 분포가 타당하다는 것이다. 


예제 


1) 공장에서 볼트의 분산이 1mm 보다 작아야 현재 생산 공정을 유지하기로 했다. n=100 개의 샘플을 고른 후, 표본 분산을 추정했을 때, 1.2mm 가 나왔다면, 모분산이 1mm 보다 크다고 볼 수 있는지 유의수준 0.05 하에서 검정하라  


검정 통계량은 아래와 같이 계산된다.


$$ \frac{(99)1.2}{1} = 111.8 $$

$$ Q \sim \chi(99) $$

$$ P (Q > 111.8) = 0.178 $$


p-value 가 0.178로 0.05 보다 크므로 모분산이 1mm 보다 크다고 할 수 없다. 


2) 모분산의 95% 신뢰구간을 구하라


모분산의 신뢰구간은 아래와 같은 식으로 구할 수 있다. 아래 신뢰구간은 통계량이 카이제곱 분포를 따른다는 사실을 통해 쉽게 확인할 수 있다. 


$$ [\frac{(n-1)s^2}{\chi_{0.025}(n-1)}, \frac{(n-1)s^2}{\chi_{0.975}(n-1)}] $$


이 때, n=100, s^2 = 1.2, chi_0.975 = 73.36, chi_0.025 = 128.422 를 대입하면 신뢰구간은 아래와 같다.


[0.92, 1.61] 


이 신뢰구간은 양측 검정을 할 때에도 활용할 수 있다. 


참고

http://www.uniwise.co.kr/pdfupload/lecture_upload/R201500241/stlsm_30.pdf


두 모집단의 모비율의 차이에 대한 검정


앞선 포스트에서 중심극한정리에 의한 정규 근사를 통해 표본 비율의 분포를 구하고 검정하는 방법을 다루었다. 이번에는 두 모집단에 대해서 비율에 차이가 있는지를 검정하는 방법에 대해서 정리해보려고 한다. 앞선 one sample 비율 검정과 다른점은 one sample 검정에서는 가설검정을 할 때, 귀무가설 하에서의 모비율을 검정 통계량 계산에 활용할 수 있었던 반면, two sample 비율 검정에서는 귀무가설이 두 모비율이 같다고 정의되는 경우가 일반적이기 때문에 모비율이 정의되지 않아 모비율을 다른 방법을 통해 추정해야 한다는 것이다. 


모델 


각각 n_1, n_2 명 존재하는 집단 A, B 에서 어떤 법안에 찬성하는지 여부를 투표를 통해 결정하는 상황을 예로 들어보자. 집단 A, B 에서 사람의 찬성 여부를 각각 확률변수 X, Y 라 할 때 다음과 같이 모델링할 수 있다. 


$$ X \sim Be(p_1) $$

$$ Y \sim Be(p_2) $$ 

 

표본 비율은 아래와 같이 정의된다. 


$$ \tilde p_1 = \frac{\sum^{n_1}_{i=1} X_i}{n_1} $$

$$ \tilde p_2 = \frac{\sum^{n_2}_{i=1} Y_i}{n_2}$$


표본비율은 베르누이 분포의 표본 평균이라고 볼 수 있기 때문에 중심극한 정리에 의해 표본비율은 아래 분포로 근사된다. 


$$ \tilde p_1 \sim N(p_1, p_1(1-p_1)/n_1) $$

$$ \tilde p_2 \sim N(p_2, p_2(1-p_2)/n_2) $$


가설 검정


$$ H_0 : p_1 = p_2 $$

$$ H_1 : p_1 \ne p_2 $$


검정 통계량 


아래 분포를 통해 검정한다. 

$$ \tilde p_1 - \tilde p_2 \sim N(p_1 - p_2, p_1(1-p_1)/n_1 + p_2(1-p_2)/n_2)  $$


검정 통계량은 다음과 같다. 

$$ \frac{(\tilde p_2 - \tilde p_2) - (p_1-p_2)}{\sqrt{p_1(1-p_1)/n_1 + p_2(1-p_2)/n_2}} $$


이 때, p_1-p_2 는 귀무가설 하에서 0이 되는데, 분모의 p_1, p_2는 사라지지 않는다. 이 모비율을 추정하는 방법에 따라 검정통계량은 아래 2가지 정도로 계산해볼 수 있다. 


1. 표본비율을 통해 모비율을 추정 


$$ \frac{(\tilde p_1 - \tilde p_2)}{\sqrt{\tilde p_1(1-\tilde p_1)/n_1 + \tilde p_2(1-\tilde p_2)/n_2}} $$


2. 합동 비율을 추정 


$$ \tilde p = \frac{\sum^{n_1}_{i=1} X_i+ \sum^{n_2}_{i=1} Y_i}{n_1+n_2} $$

$$ \frac{(\tilde p_1 - \tilde p_2)}{\sqrt{\tilde p(1-\tilde p)/n_1 + \tilde p(1-\tilde p)/n_2}} = \frac{(\tilde p_1 - \tilde p_2)}{\sqrt{(1/n_1+1/n_2)\tilde p(1-\tilde p)}} $$


예제


문제: 실험 참여자 (암 환자) 에 대해 대조약 (기존약) 과 시험약을 투약하고, 12개월 간 관찰아혀 사망한 환자의 수를 관찰한다. 이 때, 처리집단이 대조집단에 비해 더 사망자가 적다는 것을 비율의 차이를 검정하는 방법을 통해 확인해보자.  


 

 표본 크기 

사망자 환자 수 

 대조집단

 300  

50 

 처리집단

 200 

20  


$$ n_1 = 300, \tilde p_1 = 1/6 $$

$$ n_2 = 200, \tilde p_2 = 1/10 $$


검정 통계량은 아래와 같다. 


$$ \frac{(\tilde p_1 - \tilde p_2)}{\sqrt{\tilde p_1(1-\tilde p_1)/n_1 + \tilde p_2(1-\tilde p_2)/n_2}} $$

$$ \frac{(1/6 - 1/10)}{\sqrt{1/6(1-1/6)/300 + 1/10(1-1/10)/200}} =  0.0666/0.0302 = 2.205 $$

$$ P(Z > 2.205) = 0.013 $$ 


양측검정을 하기위해 곱하기 2를 해주면 p-value 가 0.026 으로 0.05 보다 작으므로 두 모비율은 다르다고 결론 내릴 수 있다. 마지막으로 R 의 stat 패키지를 통해 구한 p-value 가 손으로 구한 값과 같은지를 확인해보자.


prop.test(n=c(300,200), x=c(50,20), correct=F)


correct option 은 R 의 기본옵션인 연속성 보정을 사용하지 않는다는 것을 의미한다. 


> prop.test(n=c(300,200), x=c(50,20), correct=F)


        2-sample test for equality of proportions without continuity

        correction


data:  c(50, 20) out of c(300, 200)

X-squared = 4.4297, df = 1, p-value = 0.03532

alternative hypothesis: two.sided

95 percent confidence interval:

 0.007445812 0.125887521

sample estimates:

   prop 1    prop 2

0.1666667 0.1000000


이 때, X-squared 값에 root 를 씌우면 Z-score 를 계산할 수 있다. Z-score는 4.4297에 sqrt를 하면, 2.104685 이다. 1. 표본비율을 통해 모비율을 추정 하는 방법을 통해 구한 Z-score 와 다르다는 것을 알 수 있다. 2. 합동 비율을 추정 방법을 통해 Z-score 를 구하면, Z-score 가 2.10468 로 R 의 결과와 같다는 것을 확인할 수 있다. 


(1/6-1/10)/sqrt((1/300+1/200)*0.14*0.86) = 2.10468


R 에서는 합동 비율을 추정하는 방법을 통해 비율의 차이를 검정한다는 것을 확인할 수 있다. 


두 정규분포 모양의 모집단의 모평균의 차이에 대한 검정: Two sample t-test 


예를 들어, 서울사람과 부산사람의 몸무게의 평균의 차이가 있는지에 대해 검정하고자 한다. 이럴 때 흔히 쓸 수 있는 방법이 바로 two sample t-test 방법이다. 두 집단의 평균을 비교하여 차이가 유의한지를 보는 방법이다. 우선 언제 t-test 를 할 수 있는지 정리해보자.


t-test 는 연속형 확률 변수이고 평균이 mu 이고 표준편차가 sigma 인 분포를 갖는 모집단으로부터 n 개의 표본을 추출했을 때, n 이 크면 클 수록 중심극한정리에 의해 표본 평균이 평균이 mu 이고, 표준편차가 sigma/root(n) 인 정규분포에 근사된다는 사실을 이용한 방법이다. 따라서 t-test 의 가정은 아래 두 가지로 요약된다. 


t-test 의 가정


1. 연속형 변수 

2. 정규성

3. 임의추출 

4. 큰 샘플 사이즈 


변수가 연속형이어야하고 임의추출을 해야한다는것은 당연해 보이지만, t-test 에서 왜 정규성을 따라야할까? 처음 t-test 를 배울 때, 많이 질문하는 부분이 't-test 는 샘플이 정규분포를 따라야하나요?' 이다. 결론은 '정규분포를 따르면 좋다.' 라고 생각한다. 왜냐하면, 표본 크기가 크면, 어쨌든 표본 평균의 분포는 정규분포에 근사되기 때문에 t-test 가 불가능한 것은 아니기 때문이다. 하지만 문제는 평균이 대표성을 갖는가이다. 그리고 평균이 대표성을 갖기 위해서는 정규분포와 같은 모양이어야한다. 따라서 모평균을 비교하는 문제에 있어서, 모집단이 정규 분포와 비슷해야 타당한 것인데 모집단이 정규 분포임을 확인하는 방법은 샘플의 분포를 이용하는 길밖에 없다. 따라서 샘플의 정규성을 확인함으로써 평균이 대표성을 갖는지 확인해야한다. 만약 샘플이 정규성을 보이지 않는다면, t-test 의 타당성이 떨어지게 된다. 


예를 들어 모집단의 분포가 이항분포나 포아송 분포처럼 특정 상황에서 정규분포로 근사된다면, t-test 의 결과는 타당할 것이다. 하지만 정규성을 띄지 않는 분포에서 왔다면, t-test 보다는 비모수 검정이 더 적합하다. 왜냐하면 평균이 대표성을 띄지 않기 때문이다. 


Two sample t-test 의 종류



More about the basic assumptions of t-test: normality and sample size, statistical round, 2019


위 그림은 mu1-mu2 의 분포를 H0가 참일때와 H1 이 참일 때로 나누어서 보여주는 그림이다. 1종 오류와 2종 오류의 개념에 대해 잘 이해할 수 있는 그림이다. two sample t-test 에 개념에 대해 간단히 설명해보자. 두 집단의 모평균 (mu1, mu2) 에 차이가 있는지 검정하는 two-sample t-test 는 아래의 분포를 이용해 검정한다. 


$$ \tilde X - \tilde Y \sim N(\mu_1, \mu_2, \frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}) $$ 


H0가 참일 때 (mu1-mu2=0) 검정 통계량은 표본 분산을 이용해 아래와 같은 분포를 갖는다.


$$ \frac{\tilde X - \tilde Y}{\sqrt{\frac{s^2_1}{n_1}+\frac{s^2_2}{n_2}}} ... 식 1 $$


위 검정 통계량이 특정한 자유도를 갖는 t 분포를 갖는다는 사실을 이용해서 검정하는 것인데, t 분포의 자유도는 두 그룹의 분산에 대한 가정에 따라 다르기 계산된다.


1. 등분산 가정 (합동분산 추정량을 이용한 t-test)


합동분산추정량은 아래와 같이 계산된다. 


$$ s^2_p = \frac{n_1-1}{(n_1-1)+(n_2+1)}s^2_1 + \frac{n_2-1}{(n_1-1)+(n_2+1)}s^2_2 = \frac{(n_1-1)s^2_1 + (n_2-1)s^2_2}{n_1+n_2-2} $$ 


이 합동분산 추정량을 위 식1 의 s_1, s_2 대신에 넣으면 검정통계량을 계산할 수 있고, 이 검정통계량은 자유도가 n_1+n_2-2 인 t 분포를 따른다. 식 1을 합동 분산 추정량을 통해 쉽게 계산하기 위해 다음과 같이 정리할 수 있다. 


$$ \frac{\tilde X - \tilde Y}{\sqrt{s^2_p(\frac{1}{n_1}+\frac{1}{n_2})}} \sim t(n_1+n_2-2) $$


2. 등분산 가정을 하지 않을 때 (Welch's t-test)


등분산 가정을 하지 않으면 위 식1을 그대로 계산하면 된다. 다만, 자유도의 계산이 다소 복잡하다. 이 방법을 Welch's t-test 라 하는데 이 방법에 따라 자유도를 유도하는 것은 어려운 부분이므로, 단지 이렇게 계산한다는 것만 알고 넘어가도록 하자. Welch's t-test 에서 자유도는 아래와 같이 계산된다. 


$$ v=\frac{(s^2_1/n_1 + s^2_2/n_2)^2}{\frac{(s^2_1/n_1)^2}{n_1-1} + \frac{(s^2_2/n_2)^2}{n_2-1}} $$


즉, Welch's t-test 는 다음의 분포를 통해 검정한다. 


$$ \frac{\tilde X - \tilde Y}{\sqrt{\frac{s^2_1}{n_1}+\frac{s^2_2}{n_2}}} \sim t(v) $$



1, 2 두 방법 모두 표본수가 많다면 정규근사를 이용할 수 있다. 이 경우 t분포를 이용하지 않고 표본 수가 많을 때, 검정 통계량이 표준정규분포를 근사적으로 따른 다는 사실을 통해 검정한다. 

중심극한정리를 통한 표본 비율 분포 근사와 모비율의 검정


비율의 분포 


당신이 정책 담당자이고 전체 국민들 중 대상으로 새로운 교육 개혁에 대해 찬성하는 사람들의 비율을 알고 싶다고 해보자. 그리고 그 비율을 p 라고 할 때, p > 0.5 임을 보이고 싶다고 하자. 하지만 현실적으로 전체 국민에 대해 찬반 여부를 조사하는 것은 불가능에 가깝다. 따라서 전체 인구 집단에서 샘플링한 표본에서 구한 비율을 통해 위 가설 (p > 0.5) 을 보이고 싶다. 


예를 들어, 1000명을 조사해서 p=0.48 을 얻었다고 하자. 이 때, p > 0.5 라고 할 수 있을까? 이 표본에서는 p 가 0.5를 넘지 못했지만, 0.5 에 매우 근접했기 때문에 다른 1000 명을 조사하면 0.5가 넘을 수 있을 것을 기대할 수 있을 법하다. 한 가지 방법은, 1000 명을 여러번 뽑아서 p가 0.5를 넘는 비율을 보는 것일 것이다. 예를 들어, 1000 명을 100번 뽑았는데, 95번이 p>0.5 일 경우, p > 0.5 라고 결론내려볼 수 있을 것이다. 


하지만 1000 명을 100 번 뽑는 것도 마찬가지로 시간과 비용이 많이 든다. 따라서 현실적으로 비율의 분포를 이용한다. 표본을 통해 비율의 분포를 구할 수 있다면 더욱 쉬운 방법을 통해 가설을 검정할 수 있을 것이다. 표본 안에서 뽑은 값이 분포 안에서 어디에 위치해 있는 가를 보는 것이다. 가설이 분포 흔히 관찰할 수 있는 것이라면 그 가설이 맞는 것이고, 그렇지 않다면 가설이 틀린 것이다. 우선, 모집단의 일부를 샘플링하여 찬성 비율을 보는 것은 다음과 같이 모델링 해볼 수 있다. 모집단의 정책에 대한 찬반여부를 확률변수 X 라 할 때, 다음과 같은 베르누이 분포를 갖는다. 


$$ X \sim Be(P) $$ 


이 때, 1000의 크기를 갖는 표본 X1,...,X1000 이 존재하면, 찬성 비율은 다음과 같이 계산할 수 있다. 


$$ \tilde P = \frac{\sum_{i=1}^{i=1000}{X_i}}{n} $$ 


중심극한정리에 의해 p^ 은 정규근사가 가능하다. 왜냐하면 p^ 은 베르누이 분포의 표본 평균이기 때문이다. 


$$\tilde P \sim N(P, \frac{p(1-p)}{n}) $$ 


비율의 검정


중심극한정리를 통해 p^ 의 분포는 구했는데, 이것을 통해 p > 0.5 인지 아닌지를 어떻게 검정할 수 있을까? 사실 좋은 방법은 앞서 언급했듯, 1000 명을 여러번 뽑아서 p가 0.5를 넘는 비율을 보는 것일 것이다. 이 경우 분포를 모르더라도 결론을 내려볼 수 있다. 하지만 우리는 통계적, 근사적으로 이 문제를 해결하려고 한다. 


검정의 방법은 크게 두 가지로 나누어볼 수 있다. 


1. 가설검정을 이용한 방법

2. 신뢰구간을 이용한 방법 


가설검정을 이용한 방법


가설검정을 이용한 방법의 절차는 다음과 같다. 


1. 귀무가설 및 대립가설의 설정 


양측 검정, 단측 검정 둘 중 하나를 선택해 수행한다. 


$$ H_0 = 0.5 $$

$$ H_1 > 0.5 $$ 


2. 검정 통계량 계산 


아래분포를 이용해 검정통계량을 계산

$$ \tilde P \sim N(0.5, \frac{0.5*0.5}{1000}) = N(0.5, 0.00025) $$ 


> x = seq(0.01, 1, 0.001)

> y = dnorm(x=seq(0.01, 1, 0.001), mean=0.5, sd=sqrt(0.00025))

> plot(x, y)


검정통계량은 아래와 같다. 


$$ \frac{0.48-0.5}{\sqrt{0.00025}} = -1.26 $$ 


3. 유의수준 설정 


일반적으로 5 %의 유의수준을 설정한다. 


4. 기각역 (또는 p-value) 계산 


검정통계량과 유의수준을 통해 기각역을 결정한다.


기각역은 아래와 같은 절차로 계산한다. 


$$ P(\tilde p > R) = 0.05 $$ 


위 식에서 R 이 기각역인데, R = 0.526 이다. 


>qnorm(p=0.95, mean=0.5, sd=sqrt(0.00025)) 

0.526


유의수준을 통한 기각여부 결정은 아래와 같이한다. 유의수준 계산에는 검정통계량을 이용한다.


> 1-pnorm(-1.26)

0.89

 

$$ P(X > -1.26) = 0.89 $$ 


5. 기각역 포함 여부를 통해 의사 결정을 수행한다. 


p-value 가 0.05 보다 큰 0.89 이고, 기각역 (X > 0.526) 에 포함되지 않으므로, (사실 두 개는 같은 정보이다.) 귀무가설을 기각할 수 없다. 


신뢰구간을 이용하는 방법


신뢰구간을 이용한 검정 방법의 절차는 다음과 같다. 신뢰구간을 이용한 방법은 양측 검정을 위해 사용할 수 있다. 


귀무가설 및 대립가설의 설정 


$$ H_0 = 0.5 $$

$$ H_1 != 0.5 $$ 


신뢰수준을 정하고 p 의 구간추정을 통해 신뢰구간을 구한다. 95% 신뢰구간은 다음과 같이 구할 수 있다. 


$$ \tilde P \sim N(p, \frac{p(1-p)}{n}) $$ 

$$ P[\tilde p - 1.96* \sqrt{p(1-p)/n} < p < \tilde p+ 1.96 * \sqrt{p(1-p)/n}] = 0.95 $$

$$ (\tilde p - 1.96* \sqrt{p(1-p)/n}, \tilde p+ 1.96 * \sqrt{p(1-p)/n}) $$


위 식에서 n,p,sd 를 대입하면 다음과 같은 신뢰구간을 얻을 수 있다. p 대신의 p의 추정량 p^을 넣는다. 


(0.469, 0.531)


따라서 p^ 은 0.48 이고, 위 구간에 들어가기 때문에 귀무가설을 기각할 수 없다. 


모비율의 신뢰구간을 구할 때, 한가지 알아두면 좋은 테크닉은, p 에 표본 비율 p^을 대입하지 않고, 0.5를 대입하면 가장 보수적인 신뢰구간을 구할 수 있다는 것이다. 


$$ (\tilde p - 1.96* \sqrt{0.5(1-0.5)/n}, \tilde p+ 1.96 * \sqrt{0.5(1-0.5)/n}) $$


p를 0.5로 추정할 때, 신뢰구간의 길이가 가장 넓어지기 때문에 만약 표본 수가 충분한다면 신뢰구간의 타당성을 확보하기 좋은 방법이라고 할 수 있고, 표본수가 적어 표본비율의 변동이 큰 상황에서도 대략적인 신뢰구간을 구할 수 있다는 장점이 있다. 


Model Calibration


예측모형 (predicted model) 을 어떻게 평가할 수 있을까? 가장 직관적이면서 많이 쓰이는 평가 방법은 정확도 (accuracy), 즉, 예측한 것중 몇퍼센트나 맞았는가에 관한 지표일 것이다. 하지만 좋은 모델이란 정확해야할 뿐만아니라 잘 보정 (calibration) 되어야할 필요가 있다. 


calibration 을 평가하기 위해 사용되는 calibration plot 은 예측된 확률과, 실제 확률의 관계를 보여준다. 이를 통해 모델의 예측이 얼마나 "현실적인지" 를 측정하게 된다. 예를 들어, 이미지를 인풋으로 받아 개와 고양이를 분류하는 모델을 생각해보자. 어떤 이미지에 대해 0.8 의 확률로 개라고 반환했다면, "정말 이 이미지가 개일 확률이 0.8 인가?" 에 대한 답을 주는 것이 calibration plot 이다. 정확도란 (일반적으로 이진 분류에 관하여) 50 % 를 cutoff 로 사용하여, 예측을 A와 B 클래스로 나누어, 실제 값이랑 맞는지를 확인하는 것이지만, calibration 은 보다 면밀하게 모델의 결과값을 검증하는 과정이라고 볼 수 있다. 


Calibration plot 


만약 데이터의 실제 정답이 알려져 있다면, Calibration 을 평가하기 위해 Calibration plot 을 많이 그리게 되는데, 일반적인 방법은 다음과 같다. 


1. 모델의 예측값을 기준으로 [0,10%], (10,20%], (20,30%], … (90,100%] 에 맞게 데이터를 분할한다 (이를 binning 이라고도 한다).  

2. 각 카테고리에서 예측하고자 한 클래스의 비율 (event rate) 를 계산한다 (예제의 경우 개의 비율을 계산한다).

3. calibration plot을 그린다 : 각 카테고리에서의 중앙 값 (5 %, 15 %, 20 % ...) 를 x 로 놓고, event rate 를 y 로 놓고 그린 그림이다. 

4. calibration plot 의 선이 일직선 (45◦)임을 확인한다. 


Example


R을 통해 Calibration 을 실제로 해보자. 


diabetes.csv


실습 데이터는 Pima diabetes 데이터셋을 이용해보겠다. Pima diabetes 데이터셋은 사람들의 임상정보와 당뇨병 여부에 관한 정보를 갖고 있는 데이터셋이다. 이 때 당뇨병 여부를 예측하는 모형을 로지스틱 회귀분석 및 랜덤포레스트을 이용해 구축하고, 이 두 모델의 정확도 및 Calibration 을 평가해보자. 


데이터 로드 및 train/test split 

  • 이 데이터셋의 경우, missing value 가 많은 것이 특징이다. 
  • 아래 코드는 평균으로 missing 을 채워넣는 mean imputation 을 수행하고 train/test 를 50:50으로 나누는 코드이다. 
suppressPackageStartupMessages(library(tidyverse))
library(data.table) 
data <- readr::read_csv("../PimaDiabetes/diabetes.csv") 
data$Outcome <- factor(data$Outcome)

## Imputation
fix_missing <- function(x, missing_value) { 
  x[x == missing_value] <- NA 
  x 
} 
cols <- colnames(data)[1:8]
data[, cols] <- lapply(data[, cols], fix_missing, 0)

impute_mean <- function(x) {
    x[is.na(x)] <- mean(x, na.rm = TRUE)
    return(x)
}
data[, cols] <- lapply(data[, cols], impute_mean)
data %>% head

##  Train/Test Split
set.seed(123)
smp_size <- floor(0.5 * nrow(data))
train_ind <- sample(seq_len(nrow(data)), size = smp_size)
train <- data[trai n_ind, ]
test <- data[-train_ind, ]


로지스틱 회귀분석 모형 구축 및 test set 에 대한 예측

  • Pregnancies + Glucose + BloodPressure + Insulin + BMI + DiabetesPedigreeFunction + Age 를 통해 Outcome 을 예측하는 모형을 만든다. 
lrmodel <- glm(data = train, Outcome ~ Pregnancies + Glucose + BloodPressure + Insulin + BMI + DiabetesPedigreeFunction + Age, family = binomial("logit"))

x = predict(lrmodel, newdata = test)
p = (1 / (1+exp(-x)))
test <- test %>% mutate(lrmodel = p)

랜덤 포레스트 모형 구축 및 test set 에 대한 예측 

  • 랜덤 포레스트의 hyperparameter 인 mtry 와 ntree 는 적절한 값을 선택한다. 
library(randomForest)

rfmodel = randomForest(Outcome ~ Pregnancies + Glucose + BloodPressure + Insulin + BMI + DiabetesPedigreeFunction + Age
                      , data = train, mtry = floor(sqrt(7)), ntree = 500, importance = T)

p = predict(rfmodel, newdata = test, type = "prob")[, 2]
test <- test %>% mutate(rfmodel = p)

cutoff 정하기

  • 모형은 0~1사이의 확률을 의미하는 값을 내보내는데, 여기에 threshold 를 적용해서 0 또는 1로 변환한다. 
test <- test %>% mutate(
  rfclass = if_else(rfmodel >= 0.5, 1, 0),
  lrclass = if_else(lrmodel >= 0.5, 1, 0)
)
test$rfclass <- factor(test$rfclass)
test$lrclass <- factor(test$lrclass)


로지스틱 회귀분석 정확도 

  • 정확도는 최종 예측값 (0 또는 1) 을 기준으로, 예측한 값중 실제 정답으로 맞춘 비율을 의미하는 값이다. 
  • 로지스틱 회귀분석의 경우, 78.9 % 의 정확도를 보여준다. 

library(caret) confusionMatrix(test$lrclass, test$Outcome)

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 227  58
         1  23  76
                                          
               Accuracy : 0.7891          
                 95% CI : (0.7448, 0.8288)
    No Information Rate : 0.651           
    P-Value [Acc > NIR] : 2.494e-09       
                                          
                  Kappa : 0.5058          
                                          
 Mcnemar's Test P-Value : 0.0001582       
                                          
            Sensitivity : 0.9080          
            Specificity : 0.5672          
         Pos Pred Value : 0.7965          
         Neg Pred Value : 0.7677          
             Prevalence : 0.6510          
         Detection Rate : 0.5911          
   Detection Prevalence : 0.7422          
      Balanced Accuracy : 0.7376          
                                          
       'Positive' Class : 0  


로지스틱 회귀분석 Calibration plot

  • Calibration plot 을 그릴 수 있는 방법은 여러가지가 있지만, caret 패키지의 calibration 함수를 통해 쉽게 그려볼 수 있다. 

calibration 함수는 아래의 calibration plot 을 그릴 수 있는 정보를 dataframe으로 만들어 반환해준다. 
  • 모델의 예측값을 기준으로 [0,10%], (10,20%], (20,30%], … (90,100%] 에 맞게 데이터를 분할한다 (이를 binning 이라고도 한다).  
  • calibration plot을 그린다 : 각 카테고리에서의 중앙 값 (5 %, 15 %, 20 % ...) 를 x 로 놓고, event rate 를 y 로 놓고 그린 그림이다. 
library(caret)

cal_plot_data_lr = calibration(Outcome ~ lrmodel, 
  data = test, cuts = seq(0, 1, by=0.1), class = 1)$data 

ggplot() + xlab("Bin Midpoint") +
  geom_line(data = cal_plot_data_lr, aes(midpoint, Percent),
            color = "#F8766D") +
  geom_point(data = cal_plot_data_lr, aes(midpoint, Percent),
            color = "#F8766D", size = 3) +
  geom_line(aes(c(0, 100), c(0, 100)), linetype = 2, 
            color = 'grey50')

랜덤포레스트 정확도

  • 랜덤포레스트의 경우 로지스틱 회귀분석보다 조금 작은 0.77 % 의 정확도를 보인다. 
confusionMatrix(test$rfclass, test$Outcome)
Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 213  51
         1  37  83
                                          
               Accuracy : 0.7708          
                 95% CI : (0.7255, 0.8119)
    No Information Rate : 0.651           
    P-Value [Acc > NIR] : 2.421e-07       
                                          
                  Kappa : 0.4831          
                                          
 Mcnemar's Test P-Value : 0.1658          
                                          
            Sensitivity : 0.8520          
            Specificity : 0.6194          
         Pos Pred Value : 0.8068          
         Neg Pred Value : 0.6917          
             Prevalence : 0.6510          
         Detection Rate : 0.5547          
   Detection Prevalence : 0.6875          
      Balanced Accuracy : 0.7357          
                                          
       'Positive' Class : 0               
                                          


랜덤포레스트 Calibration plot

  • 랜덤포레스트에서도 같은 방법으로 calibration plot 을 그릴 수 있다. 
  • 정확도는 랜덤포레스트에서 약간 작았지만, Calibration 은 더 좋은 모습을 보인다.
  • 하지만 train/test 의 비율, hyperparameter 구성에 따라 calibration 이 달라지니, 다양한 세팅에서 검증해볼 필요가 있다. 
cal_plot_data_rf = calibration(Outcome ~ rfmodel, 
  data = test, class = 1)$data

ggplot() + xlab("Bin Midpoint") +
  geom_line(data = cal_plot_data_rf, aes(midpoint, Percent),
            color = "#F8766D") +
  geom_point(data = cal_plot_data_rf, aes(midpoint, Percent),
            color = "#F8766D", size = 3) +
  geom_line(aes(c(0, 100), c(0, 100)), linetype = 2, 
            color = 'grey50')


https://medium.com/optima-blog/model-calibration-4d710a76c54

http://appliedpredictivemodeling.com/blog?offset=1532965627474


Adaptive Histogram Equalization 이란 무엇인가? 


Adaptive Histogram Equalization (AHE) 이란 이미지 전처리 기법으로 이미지의 contrast 를 늘리는 방법이다. Histogram Equalization (HE) 방법을 조금 더 개선해서 만든 방법이라고 할 수 있는데, 그 차이점은 HE 는 하나의 histogram 을 통해 uniform 분포를 가정한 equalization 을 수행하는 반면, AHE 는 여러개의 histogram 을 만든다는 차이점이 있다. HE 의 문제점은 하나의 histogram 으로 pixel intensity 의 redistribution 을 수행하므로, pixel intensity 의 분포가 이미지의 전 범위에 걸쳐서 동일해야 적절히 수행될 수 있는 방법이다. 만약 이미지의 일부분이 다른 지역과 비교해 다른 분포를 갖고 있을 경우, HE 를 적용한 이미지는 왜곡된다. 




AHE 방법은 이미지를 grid 를 이용해서 여러 개로 쪼개고 이 쪼갠 이미지에 대해서 각각 HE 를 적용한다. 따라서 그 grid 안에서 이미지의 contrast 가 좋아지기 때문에 이미지의 local contrast 를 조정하기에 좋은 방법이라고 할 수 있다. 위 그림처럼 어떤 pixel 의 AHE 를 적용한 값은 각 box 안에서 HE 를 수행한 값이다. 경계에 위치한 pixel 에 대한 transformation 은 특별하게 다루어지는데, 그 반대편에 위치한 pixel 값들을 mirroring 함으로써 가상의 픽셀 값을 만든 후, AHE 를 적용하게 된다. 


AHE 의 단점은 어떤 grid 안에 있는 pixel intensity 들이 매우 좁은 지역에 몰려있는 경우에 발생한다. 이 경우 AHE 는 매우 좁은 지역에 위치해 있는 pixel intensity 들을 넓은 지역으로 퍼트리기 때문에 이 안에 다른 지역과 매우 작은 차이를 보이는 noise 가 있다고 하더라도 상당한 크기의 peak 를 만들어내게 된다. 예를 들어, 8x8 grid 안에서 하나의 pixel 이 1이고 나머지는 0인 경우, 1이 255로 변환이 되게 된다. 즉, 의미 없는 pixel 의 값이 크게 변화되는 문제가 생긴다. 이를 noise amplification 이라고 부르기도 한다. 또한 의미없는 지역의 pixel intensity 를 enhance 시키기도 한다. 즉, 이로 인해 원래 이미지에 왜곡이 생기게된다. 


Contrast Limited Adaptive Histogram Equalization 



Adaptive Histogram Equalization 과 늘 따라다니는 것이 Contrast limited adaptive histogram equalization (CLAHE) 방법이다. 이 방법은 AHE 의 변형으로 AHE 의 중대한 문제점인 noise amplification 을 해결하기 위해 contrast limit 을 활용한다. 만약 pixel intensity 가 좁은 구역에 몰려있다면, 아래와 같은 두 가지 현상이 생긴다. 


  1. 그 부근에서 slope of cdf 가 급격하게 증가한다.

  2. histogram의 높이가 매우 높다.


1, 2 는 결국 같다. slope of cdf 가 급격하게 증가하는 부분이 histogram 의 높이가 높은 지점이다. 따라서 CLAHE 방법은 위 그림과 같이 cdf 를 계산하기 전에 histogram 에 높이에 제한을 걸어 특정 높이 이상에 있는 pixel 값들을 redistribution 하는 과정을 먼저 거치게 된다. 이 때 histogram 의 높이를 제한하는 값을 clip limit 이라고 부른다. redistribution 이후에, 다시 clip limit 의 값을 초과할 수 있는데, 이를 그대로 놔두고 할 수 있지만, 제거하고 싶은 경우, 반복적으로 redistribution 을 실시한 후, 초과하는 pixel 값들이 무시할만 한 정도까지 반복할 수 있다. 결국 이러한 과정은 cdf function 의 기울기가 너무 높지 않게 제한한다. 최종적으로 가상의 cdf 함수를 구할 수 있고, 이를 통해 HE 를 적용하면 CLAHE 를 적용한 이미지를 얻을 수 있다. 


그러면 이 이미지는 좁은 범위에 pixel intensity 가 몰려 있는 경우 생기는 noise 에 robust 하도록 만든다. 왜냐하면 몰려 있는 픽셀의 값을 redistribution 을 하기 때문에 급격한 기울기가 없는 cdf를 통해 pixel intensity 를 변환하기 때문에 noise 의 pixel intensity 를 너무 큰 값으로 변환시키지 않도록 만든다. 그렇기 때문에 CLAHE 를 실시한 후의 local region 의 pixel intensity 가 0~255 사이에 위치해 있는 것이 아니라, 더 좁은 범위에 몰려 있게 만들 수 있다. 즉 CLAHE 를 사용했을 때 장점은, contrast 가 낮은 지역에 위치한 noise 에 robust 하도록 변환함으로써, 변환된 이미지가 실제 이미지와 비슷한 특성을 갖도록 만드는 것이다. 


참고

https://en.wikipedia.org/wiki/Adaptive_histogram_equalization