Autoregressive Processes

Autoregressive process 란 history 가 현재 값에 직접적인 영향을 주는 time series 를 말한다. 식으로 표현하면 아래와 같다. 

현재시점 (t) 의 값은 과거 시점들 (t-1 ~ t-p) 의 값을 가중치를 두고 합한 값에 error term (Z) 을 더한 값이다. 

$$ X_t = Z_t + \phi_1(X_{t-1})  + \phi_2(X_{t-2}) ... + \phi_p(X_{t-p}) $$

Example

p=2 이고, 가중치가 0.7, 0.2 인 autoregressive process 는 아래와 같다.

$$ X_t = Z_t + 0.7X_{t-1} + 0.2X_{t-2} $$

이를 R 코드로 구현하면 아래와 같다. 선차트를 통해 보면 현재 값이 과거 값과 높은 상관성이 있다는것을 확인할 수 있다. correlogram 을 통해 가까운 시간에 측정된 값이 현재값과 더 높은 상관성이 있다는것을 확인할 수 있다. (가중치가 0.7, 0.2 이므로)

set.seed(2017)
X.ts <- arima.sim(list(ar=c(0.7,0.2)), n=1000)
par(mfrow=c(2,1))
plot(X.ts, main="AR(2) Time series, phi1=0.7, phi2=0.2")


X.acf <- acf(X.ts, main="Autocorrelation of AR(2) Time series")
X.acf

Moving average process 와의 관계 

Autoregressive process 는 moving average process 의 무한 수열로 나타낼 수 있다. 

차수 (p) 가 1인 AR 을 생각해보자.

아래와 같이 식을 쓸 수 있다. (Z는 평균이 0, 분산이 sigma^2 을 따른다고 가정하고, phi 를 theta 로 치환하자.) 

$$ X_t = Z_t + \phi X_{t-1} = Z_t + \phi Z_{t-1} + \phi^2 X_{t-2} ... = Z_t + \theta_1 Z_{t-1} + \theta_2 X_{t-2} ... $$

AR (1) 의 통계량

위와 같이 식을 써서 moving average 처럼 표현하면 AR process 에서의 X(t) 의 기댓값과 분산을 쉽게 구해볼 수 있다. 

$$ E(X_t) = 0 $$

$$ V(X_t) = \sigma^2 \sum^{\infty}_{i=0}\theta_i^2  $$ 

time series 가 stationarity 를 만족하기 위해서는 분산이 t 에따라 바뀌지 않고, 일정해야한다. 그렇기 때문에 phi 의 절댓값이 1보다 작은것은 stationarity 에 대한 필요 조건이라고 할 수 있다. 

AR(1) 의 auto covariance function 

MA process 의 acf 와 유사한 형태라는것을 확인할 수 있다. (link)

$$ \gamma(k) = \sigma^2 \sum^{\infty}_{i=0} \theta_i \theta_{i+k} $$ 

AR(1) 의 auto correlation coefficient 

$$ \rho(k) = \frac{\sum^{\infty}_{i=0} \theta_i \theta_{i+k}}{\sum^{\infty}_{i=0} \theta_i \theta_i} $$ 

theta 를 phi 로 치환하여 AR(1) 의 auto covariance 와 auto covariance coefficient 를 구해보자. (무한등비수열 공식 사용하여 정리) 

$$ \gamma(k) = \sigma^2 \frac{\phi^k}{1-\phi^2} $$

$$ \rho(k) = \phi^k $$

AR Process 의 Stationarity 를 확인하는 방법

Example: AR(1) process

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

위 AR(1) process 에서 Z 를 제외한 나머지 텀들을 한쪽으로 옮겨서 아래와 같은 식을 만들 수 있다.

$$ \phi(B) = 1-\phi B $$ 

이 때, 우변을 0으로 만드는 B 의 해를 찾는다. 해는 B = 1/phi 이다. B 의 해가 단위원 (unit circle) 바깥에 있는 것이  stationarity 를 만족하기 위한 조건이 된다. 따라서 AR(1) 모델에서는 phi 의 절댓값이 1 미만이어야 stationarity 를 만족한다. 

Example: AR(2) process 

$$ X_t = \frac{1}{3} X_{t-1} + \frac{1}{2} X_{t-2} + Z_t $$

$$ \phi(B) = 1-\frac{1}{3} B-\frac{1}{2} B^2 $$

위 식에서 B 의 해를 찾으면 (-2+sqrt(76))/6, (-2-sqrt(76))/6 이 된다. 위 두 값이 모두 단위원 바깥에 있기 때문에 위 AR(2) process 는 stationarity 를 만족한다. 

Week stationarity

이전 포스팅에서 개념을 직관적으로 소개하였다. 이번엔 좀 더 포멀한 정의와 함께 예시를 통해 설명해보려고한다.

아래 조건을 만족할 때 weekly stationary 라고 부른다. 

1) 시간에 따른 평균이 같다. 

$$ \mu(t) = \mu  $$ 

2) Auto covariance function 이 time spacing 에만 의존한다. (t2=t1+tau 라고 생각하면 이해가 쉽다.)

$$ \gamma(t_1, t_2) = \gamma(t_2-t_1) = \gamma(\tau) $$

: 이는 시간에 따른 분산이 같다는 조건을 포함하는 조건이다. 

Examples

실제 시계열 데이터의 예시를 통해 stationarity 에 대해 더 이해해보려고한다. 

1) White noise 는 stationarity 를 만족한다.

White noise model

$$ X_t \sim N(0, \sigma) $$

White noise model은 시간에 따른 평균이 같다. 

$$ \mu = 0 $$

White noise model은 Auto covariance function 이 time spacing 에만 의존한다.

$$ \gamma(t_1, t_2) = 0, \ when \ t_1 = t_2 $$

$$ \gamma(t_1, t_2) = \sigma^2, \ when \ t_1 \neq t_2 $$

2) Random walk 는 stationarity 를 만족하지 않는다.

Random walk model

아래 식에서 X(t) 를 random walk model 이라고 한다. 

$$ Z_t \sim iid(\mu, \sigma^2) $$

$$ X_t = X_{t-1} + Z_t = \sum^{t}_{i=1}Z_i $$

따라서 Random walk model 의 시간에 따른 평균은 t*mu 이고, 분산은 t*sigma^2 이다. 만약 Z의 평균이 0이라고 가정하더라도 분산이 시간에 따라 점점 커진다는 것을 알 수 있다. 따라서 Random walk model 은 stationarity 를 만족하지 않는다. 

3) Moving average model 는 stationarity 를 만족한다. 

moving average model 

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} ... + \theta_q Z_{t-q} $$

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

moving average 의 parameter q 와 가중치 theta 를 고정해놓고 계산을 하면, 평균과 분산은 t 와는 관계 없이 고정된다는 것을 알 수 있다. 따라서 moving average model 은 stationarity 를 만족한다. 

추가적으로 Moving average model 의 auto covariance function 을 구해보자. 

moving average model 은 stationarity 를 만족하기 때문에 auto covariance function 은 time spacing 에만 의존한다. 또한 이전 포스팅에서 time spacing 이 최대 q 인 경우에만 자기상관성이 존재한다는 것을 correlogram 을 통해 확인할 수 있었다. moving average model 의 노이즈의 평균이 0일 때를 가정하고 covaraicne 를 구해보자.

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

$$ Cov(X_t, X_{t+k}) = E(X_t X_{t+k}) - E(X_t)E(X_{t+k}) =  E(X_t X_{t+k}) $$

위 기댓값을 정리하면 아래와 같은 식이 된다.

$$ Cov(X_t, X_{t+k}) = E(X_t X_{t+k}) = \sigma^2 \sum_{i=0}^{q-k} \beta_i \beta_{i+k}, \ when \ k \leq q $$

$$ Cov(X_t, X_{t+k}), \ when \ k \gt q $$ 

참고) 위 식의 정리에는 아래 기댓값과 분산의 성질을 이용하면 된다.

$$ V(X) = E(X^2) - E(X)^2 $$

$$ E(XY) = E(X)E(Y), \ when \ X, \  Y \ is \ independent $$ 

Moving average process

특정 시점 t에서의 주가를 X_t 라고하자. 또한 특정 시점 t 에서의 회사의 공지 Z_t (noise) 가 주가에 영향을 미친다고 하자. 그런데 과거 시점 (t-1, t-2...) 에 회사의 공지도 주가에 영향을 미친다. 이런 경우에 X_t 를 다음과 같이 모델링할 수 있다.

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} ... + \theta_q Z_{t-q} $$

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

이 때, q는 어떤 시점의 noise 까지 현재값에 영향을 미치는지를 의미하며, MA(2) 는 이와 같이 정의된다. 

MA(2) Process

$$ X_t = Z_t + \theta_1 Z_{t-1} + \theta_2 Z_{t-2} $$

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

MA(2) process 의 예를 들면 아래와 같다. 

$$ X_t = Z_t + 0.7 Z_{t-1} + 0.2 Z_{t-2} $$

MA(2) simulation

위에 예시로 든 MA(2) process 를 R 을 통해 simulation 해보자.

# noise 생성
noise <-rnorm(10000)

ma_2 = NULL

# ma(2) 생성을 위한 loop
for (i in 3:10000) {
  ma_2[i] = noise[i] + 0.7*noise[i-1]+0.2*noise[i-2]
}

# shift 
moving_average_process <- ma_2[3:10000]
moving_average_process <- ts(moving_average_process)

par(mfrow=c(2,1))

plot(moving_average_process, main = "A moving average process of order 2", ylab = "")
acf(moving_average_process, main = "Correlogram of ma (2)", ylab = "")

correlogram 을 보면 time step 이 0,1,2 인 경우에만 상관성이 있는 것을 확인할 수 있다. 우선, time step 이 0 인 경우는 항상 auto correlation coef 1이다. 또한 현재값에는 최대 2 time step 전의 noise 까지 반영이 되기 때문에, 최대 2 time step 의 값과 상관성이 있다는 것을 확인할 수 있따. 

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 $$

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