[생명과학을 위한 통계학] (5) Confidence Interval, Effect Size, Power Calculation

[생명과학을 위한 통계학] 시리즈 글들은 Data Analysis for the Life Sciences by Rafael A Irizarry & Michael I Love와 오픈 강의 플랫폼들 (Coursera, edX)에서 제공하는 강의들을 바탕으로 작성되었습니다. 더 자세한 내용을 확인하고 싶으신 경우 글 하단의 링크를 확인해주세요.

---

이전 게시물에서도 예시로 들었던 일반식 먹은 쥐와 고지방식을 먹은 쥐의 몸무게를 비교하는 실험을 생각해보자.

앞으로 일반식을 먹은 쥐들을 그룹 N, 고지방식을 먹은 쥐들을 그룹 H라고 부르겠다. 만약 실험자가 "그룹 H의 평균 몸무게가 그룹 N의 평균 몸무게보다 크다 (p-value <0.05)"라고 결론을 내린다면 이는 올바른 statistical summary일까?

특정 데이터를 분석함에 있어서 p-value만을 제공하는 것은 데이터에 대한 제대로 된 statistical summary라고 할 수 없다. 통계학적으로 유의미하다는 것이 항상 과학적으로 유의미하다는 것을 나타내지는 않기 때문이다. 

이게 무슨 말인가 하면... 쥐는 평균 몸무게가 49g 정도 된다. 이러한 상황에서 "그룹 H의 평균 몸무게가 그룹 N의 평균 몸무게보다 0.1g 크고, 이것이 통계학적으로 유의미한 p-value를 가진다는 사실"은 별다른 과학적 의미를 가지지 않는다는 것이다.

따라서 statistical summary에 있어서 p-value 못지않게 중요한 것이 바로 confidence interval과 effect size이다.

---

1. Confidence Interval (신뢰구간)

95% confidence interval이란 같은 population으로부터 추출된 sample의 경우, 95%의 확률로 population의 평균이 sample의 confidence interval 안에 들어간다는 의미이다. 

만일 sample의 confidence interval 안에 population의 평균이 들어오지 않는다면 sample이 population으로부터 추출되었을 확률이 5% 미만이라는 것을 의미한다. 이는 p-value < 0.05와 의미하는 바가 같다. 

과연 sample의 confidence interval이 population의 평균을 포함하는 확률이 얼마나 되는지 chowPopulation 데이터를 가지고 시뮬레이션을 돌려보자.

# Library loading
library(rafalib)

# Data processing
dat <- read.csv("mice_pheno.csv")
chowPopulation <- dat[dat$Sex == "F" & dat$Diet == "chow", 3]

# Parameter setting
B <- 250
N <- 30
Q <- qnorm(1 - 0.05/2)

# Plotting Graph
mypar()
plot(mean(chowPopulation) + c(-3,3), c(1,1), type = "n", xlab = "weight", ylab = "interval", ylim = c(1,B))
abline(v = mean(chowPopulation))

# Plotting Confidence Intervals
for (i in 1:B){
        chow <- sample(chowPopulation, N)
        se <- sd(chow)/sqrt(N)
        interval <- c(mean(chow) - Q*se, mean(chow) + Q*se)
        covered <- mean(chowPopulation) <= interval[2] & mean(chowPopulation) >= interval[1]
        color <- ifelse(covered, 1, 2)
        lines(interval, c(i, i), col = color)
}

Confidence interval 안에 chowPopulaiton weight mean이 들어올 경우 녹색, 들어오지 못할 경우 붉은색으로 표기.

위 그래프를 보면 각 선들이 한번 sampling을 바탕으로 한 95% confidence interval을 나타낸다. chowPopulation의 몸무게 평균이 포함되는 confidence interval의 경우 녹색으로 표시했고, 그렇지 않은 경우 붉은색으로 표시했다.

예상한 대로 대략 5%의 경우 sample confidence interval에 chowPopulation의 평균이 포함되지 않는 것을 볼 수 있다.

Confidence interval은 예측하고자 하는 값이 나올 수 있는 구간을 제시해주기 때문에 단순히 p-value만 제시하는 것보다 유용하다.

Confidence interval과 effect size에 대해서는 아래 링크를 참조하면 더 자세히 설명되어 있으니 참고하길 바란다.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5133225/

Alternatives to P value: confidence interval and effect size

---

2. Effect Size

Effect size란 두 집단 사이의 차이의 크기를 비교 가능하게 수치화 하는 것을 의미한다. p-value는 sample size가 커질 경우 필연적으로 작아지는데, effect size는 sample size에 의한 p-value 왜곡을 제거해주고, 단순히 Yes or No의 분석이 아니라 How much의 수치적인 분석을 가능하게 해 준다.

$$ Effect\; Size = \frac{[Mean \;of \;experimental \;group] - [Mean \;of \;control \;group]}{Standard \;Deviation}$$

역시 effect size에 대해 잘 설명해놓은 링크를 첨부한다.

https://www.leeds.ac.uk/educol/documents/00002182.htm

It's the effect size, stupid: what effect size is and why it is important

 

---

3. Power Calculation

population의 평균 간에 실제로 차이가 존재하더라도 가설 검정을 통해 항상 귀무가설을 기각할 수 있는 것이 아니다. 아래의 시뮬레이션을 보자.

# Library loading
library(dplyr)

# Data processing
dat <- read.csv("mice_pheno.csv")
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% select(Bodyweight) %>% unlist
hfPopulation <- filter(dat, Sex == "F" & Diet == "hf") %>% select(Bodyweight) %>% unlist
mu_hf <- mean(hfPopulation)
mu_control <- mean(controlPopulation)

# Calculating actual difference
print((mu_hf - mu_control)/mu_control * 100)

# Simulation
set.seed(1)
N <- 5
hf <- sample(hfPopulation, N)
control <- sample(controlPopulation, N)
t.test(hf, control)$p.value

실제로는 controlPopulation과 hfPopulation 사이에 9.94%의 평균 차이가 존재하는 것을 알 수 있다. 하지만 5개씩 sampling 하여 t-test를 진행하면 0.58의 p-value가 나오기 때문에 이를 바탕으로 두 집단 간에 차이가 존재하지 않는다는 귀무가설을 기각할 수 없다.

이 결과에서 볼 수 있듯, 귀무가설을 기각하지 못했다는 것이 귀무가설이 꼭 사실임을 의미하지는 않는다. 이 simluation에서 귀무가설이 기각되지 않은 이유는 우리가 사용한 가설 검정이 충분한 power를 가지지 못했기 때문이다.

 

통계에는 두 종류의 error가 존재한다.

Type 1 error (false positive) : 귀무가설을 기각해서는 안되는데 기각할 확률 =  alpha (유의 수준)

Type 2 error (false negative) : 귀무가설을 기각해야 하는데 기각하지 않을 확률

Power : 귀무가설을 기각해야 할 때 기각할 확률 = 1 - Type 2 error

 

아래 시뮬레이션으로 위 실험의 power를 계산해보면 0.214가 나온다. 즉 sample size가 12이고, 유의 수준을 0.05로 놓았을 때, 귀무가설을 기각해야 할 때 기각할 확률이 대락 22% 밖에 되지 않는다는 것을 알 수 있다.

# Parameter Setting
N <- 12
alpha <- 0.05
B <- 2000

# Reject function
reject <- function(N, alpha = 0.05){
        hf <- sample(hfPopulation, N)
        control <- sample(controlPopulation, N)
        pval <- t.test(hf, control)$p.value
        pval < alpha
}

# Power calculation
rejections <- replicate(B, reject(N))
mean(rejections)

 

Power는 sample size에 의해 영향을 받는데, 시뮬레이션을 해보면 sample size가 커질수록 power 역시 같이 커지는 것을 확인할 수 있다.

 

Ns <- seq(5, 50, 5)
power <- sapply(Ns, function(N){
        rejections <- replicate(B, reject(N))
        mean(rejections)
})
plot(Ns, power, type = 'b')

Sample size (Ns) vs Power

 

alpha (유의 수준)이 커질 경우, 즉 type 1 error가 증가할수록 power 역시 커지는 것을 시뮬레이션을 통해 확인할 수 있다.

N <- 30
alphas <- c(0.1, 0.05, 0.01, 0.001, 0.0001)
power <- sapply(alphas, function(alpha){
        rejections <- replicate(B, reject(N, alpha = alpha))
        mean(rejections)
})
plot(alphas, power, xlab = "alpha", type = "b", log = "x")

alpha vs power

 

모든 경우에 일률적으로 적합한 power나 유의 수준은 존재하지 않는다. power와 유의 수준, sample size간의 상관관계를 잘 이해하고 각 실험에 맞게 적합한 sample size, 유의수준을 설정하는 것이 중요하다.

 

아래 그래프는 sample size, 유의수준, power 간의 상관관계를 하나의 그래프에 그린 것이다.

alpha vs power. Different sample sizes are represented as different line colors.

 

---

[생명과학을 위한 통계학] 시리즈 글들은 Data Analysis for the Life Sciences by Rafael A Irizarry & Michael I Love와 오픈 강의 플랫폼들 (Coursera, edX)에서 제공하는 강의들을 바탕으로 작성되었습니다. 더 자세한 내용을 확인하고 싶으신 경우 글 하단의 링크를 확인해주세요.

---

[참고서적]
 - [PDF] Data Analysis for the Life Sciences by Rafael A Irizarry & Michael I Love

[오픈 강의]
1. [Coursera] R Programming from Johns Hopkins University
2. [edX] Statistics and R from Harvard University
3. [edX] Introduction to Linear Models and Matrix Algebra from Harvard University
4. [edX] Statistical Inference and Modeling for High-throughput Experiments from Harvard University
5. [edX] High-Dimensional Data Analysis from Harvard University
6. [edX] Introduction to Bioconductor from Harvard University
7. [edX] Case Studies in Functional Genomics from Harvard University
8. [edX] Advanced Bioconductor from Harvard University

---