Random walk model

아래와 같이 정의되는 X_t 를 random walk 이라고 한다. X_t 는 이전 time step 에서의 값 X_t-1 에 noise Z가 더해진 값이다. random sampling 과 다른점은 현재값이 이전값에 더해진다는것이다. 이는 랜덤하게 어떤 한 방향으로 걷는것과 비슷하다. 매번 시작점에서 한발짝 걷는 것이 아니라 한발짝 걸어서 도착한 곳에서 다시 한발짝을 간다. 

$$ X_t = X_{t-1} + Z_t $$

$$ Z_t \sim Normal(\mu, \sigma) $$

이러한 random walk 모델에서 X_t 는 이전 time step 에서의 값 X_t-1 과 매우 큰 연관성을 갖는다. 따라서 non-stationary time series 데이터이다. 

Random walk model simulation in R

R 로 random walk 모델을 만들어보자. 아래는 1000개의 random walk 데이터를 생성하는 예제이다. 시계열 그래프를 그려보면, 이 데이터는 non-stationary time series 데이터라는 것을 확인할 수 있다. 구간을 나눠서보면 트렌드를 보이기 때문이다. 

x <- NULL
x[1] <- 0
for(i in 2:1000){
  x[i] <- x[i-1]+rnorm(1)
}

random_walk <- ts(x)
plot(random_walk, main="A random walk", ylab="", xlab=" Days", col="black")

random walk 데이터에서 correlogram 을 그려보자. 인접한 time step 에서 auto correlation coefficient 가 큰 패턴을 보이기 때문에 non-stationary time series 라는 것을 다시 확인할 수 있다. 

acf_result <- acf(random_walk)

random walk 모델에서 noise Z는 stationary time series 라고 볼 수 있다. 

$$ Z_t \sim Normal(\mu, \sigma) $$

noise 가 stationary time series 라는 것을 데이터로 실제로 확인해보자. 

random_walk_diff <- diff(random_walk)
plot(random_walk_diff, main="A random walk diff", ylab="", xlab=" Days", col="black")

참고자료

https://people.duke.edu/~rnau/411rand.htm

Auto correlation coefficient 

앞선 포스팅에서 auto covariance coefficient 에 대해 설명하였다. auto covariance coefficient 은 time series 데이터에서의 각각의 time point 간 연관성을 의미하는데, stationary time series 에서는 k 라고하는 parameter 에 의해 달라진다. auto covariance coefficient 의 추정값 c_k 는 아래와 같이 계산된다. 

$$ c_k = \frac{\sum^{N-k}_{t=1}(x_t - \bar{x})(x_{t+k}-\bar{x})} {N-k} $$

이번에는 auto correlation coefficient 에 대해 정리해보려고 한다. auto correlation coefficient 도 auto covariance coefficient 와 마찬가지로 time series 데이터에서 time step 별 값의 연관성을 의미하는데 범위를 -1~1로 조정한 것으로 이해할 수 있다. 마치 공분산과 상관계수의 관계와 같다.

auto correlation coefficient 의 계산식은 아래와 같다. 

$$ -1 \leq \rho_k = \frac{\gamma_k}{\gamma_0} \leq 1 $$

이 때, rho_k 의 추정값 r_k 는 아래와 같다. c0 는 time step (=lag) 이 0일 때의 auto covariance coefficient 로 이는 분산과 같다.  

$$ \rho_k \sim r_k = \frac{c_k}{c_0} $$

$$ c_0 = \frac{\sum^{N}_{t=1}(x_t-\bar{x})^2}{N}  $$

Correlogram

rho_k 의 추정값을 k 에 따라 구한 뒤, 이를 시각화해서 표현한 것을 correlogram 이라고 한다. 이 때, rho_0 은 항상 1이다. (자기 자신과의 연관성이기 때문이다.) 

$$ r_0 = \frac{c_0}{c_0} , r_1 = \frac{c_1}{c_0} ... $$ 

R 에서 correlogram 은 acf 함수를 통해 쉽게 그려볼 수 있다.

예제 1)

아래 R 코드는 100개의 표준정규분포를 따르는 데이터를 만든 후, correlogram 을 그리는 코드이다. 파란선은 연관성이 유의한지에 대한 임계치를 의미한다. 유의한 데이터 포인트가 하나 밖에 없고, lag 에 따른 패턴이 보이지 않으므로, 전체적으로 시계열 데이터가 자기상관성이 없다고 결론 내릴 수 있다.  

purely_random_process <- ts(rnorm(100))
print(purely_random_process)
plot(purely_random_process)

auto_correlation_coef_by_lags <- acf(purely_random_process)
print(auto_correlation_coef_by_lags)

예제 2)

실제 데이터를 correlogram 을 그려보자. 다음은 모 어플리케이션의 월간 활성 이용자수 (MAU, monthly active user) 추이이다. 이 서비스는 점점 성장하는 추이를 보여주고 있다. 시계열 데이터의 관점에서는 시간에 따른 평균의 변화 (trend) 를 보이는 non-stationary time series 이다. 

위 데이터에서 correlogram 을 그리면 아래와 같이 나타난다. lag 에 따른 auto correlation coef 의 패턴이 보이며 (점점 감소), 인접한 데이터 포인트에서는 유의한 상관성을 보이고 있는 것을 확인할 수 있다. 

Time series data (시계열 데이터) 

어떤 종류의 데이터이든 상관 없으며, 그저 시간에 따라 수집된 데이터를 시계열 데이터 (timeseries data) 라고 한다. 

예를 들어, 일별 코로나 확진자수는 1일이라고 하는 time step 으로 수집된 시계열 데이터의 한 종류이다. 

Week stationary time series 

week stationary time series 란 다음의 조건을 만족한다. 

1) 시간에 따른 평균 (mean) 에 변화가 없다. 

2) 시간에 따른 분산 (variance) 의 변화가 없다. 

3) 주기적인 등락 (flucation) 이 없다.

이러한 조건을 만족하긴 위해서는 time series 의 한 섹션 (A 섹션) 고른 후, 다른 섹션 (B 섹션) 을 골랐을 때, A, B 섹션이 비슷하면 된다. 

Stochastic process

random variable 의 collection - X1,X2,X3 .. 가 있다고 하자. 이들이 각각 다른 모수를 가진 분포를 따를 때, 이를 stochastic process 라고 한다. stochastic process 의 반대개념은 deterministic process 이다. deterministic process 는 모든 step (t) 에 대해서 예측 가능하다. 예를 들어, 어떤 함수에 대한 미분함수는 특정 X 에서의 Y 값을 정확하게 알 수 있다. 이와 반대로 stochastic process 는 매 step 이 random 이기 때문에 어떤 확률 분포에서 왔다는 것만을 알 수 있을 뿐, 값을 정확하게 예측할 수 없다.

$$ X_t \sim distribution(\mu_t, \sigma_t) $$

예를 들어, 다음과 같은 시계열 데이터가 있다고 해보자. 

$$ X_1 = 30, X_2 = 29, X_3 = 57 ... $$ 

시계열 데이터를 바라보는 한 가지 관점은 stochastic process 의 실현 (realization) 으로 보는것이다. 매 timestep 별로 어떤 확률 변수가 정해지고 우리는 그 확률변수에서 나온 하나의 샘플값을 관찰하는 것이다. 

Autocovariance function

stationary time series 라는 가정을 하자. 두 가지 timestep s,t 에서의 covariance 를 정의해볼 수 있다. (확률 변수이기 때문)

$$ \gamma(s,t) = Cov(X_s, X_t) $$

$$ \gamma(s,s) = Var(X_s)   $$  

또한 아래처럼 covariance function 을 정의할 수 있는데, 이 함수는 stationary time series 라는 가정 하에 t 에 따라서는 값이 바뀌지 않으며, k가 결정하는 함수가 된다. (아래 식에서 c 는 추정값이다.) 이러한 time step (k) 에 따른 공분산의 식을  autocovariance function  이라고 한다. 

$$ \gamma_k = \gamma(t, t+k) \sim c_k $$ 

즉, stationary time series 에서는 Cov(X1, X2) 나 Cov(X10,X11) 이나 기댓값은 같다고 할 수 있다. 그 이유는 데이터에서 두 가지 섹션을 선택했을 때, 그 모습이 똑같다고 기대하는것이 stationary time series 이기 때문이다.  

또한 gamma(t_k, t) 는 autocovariance coefficient 라고 하며, stochastic process 에서의 실제 autocovariance 값이다. 데이터를 통해 구한 c_k 를 통해 autocovariance coefficient 를 추정한다. 

Autocovariance coefficient

그러면 Autocovariance coefficient 의 추정값은 어떻게 구할까? timestep 을 k 라고 할 때, 추정값은 아래와 같다. 

$$ c_k = \frac{\sum^{N-k}_{t=1}(x_t - \bar{x})(x_{t+k}-\bar{x})} {N} $$

이는 time series 데이터에서 k time step 만큼 차이나는 점들의 묶음 (x_t, x) 을 확률 변수의 관찰값으로 놓고 추정한 covariance 와 같다.

참고로 확률 변수 X,Y 의 covariance 의 추정값은 아래와 같이 구할 수 있다. 

$$ s_{xy} =  \frac{\sum^{N}_{t=1}(x_t - \bar{x})(y_t-\bar{y})} {N-1} $$

R 에서는 acf 함수를 통해 auto covariance coefficient 추정값을 계산할 수 있다. 

purely_random_process <- ts(rnorm(100))
print(purely_random_process)
plot(purely_random_process)

auto_covariance_coef_by_lags <- acf(purely_random_process, type = "covariance")
print(auto_covariance_coef_by_lags)

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

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

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

전략

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

다중 비교 (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} = 118.8 $$

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

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

p-value 가 0.085로 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개월 간 관찰아혀 사망한 환자의 수를 관찰한다. 이 때, 처리집단이 대조집단에 비해 더 사망자가 적다는 것을 비율의 차이를 검정하는 방법을 통해 확인해보자.  

$$ 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분포를 이용하지 않고 표본 수가 많을 때, 검정 통계량이 표준정규분포를 근사적으로 따른 다는 사실을 통해 검정한다.