자기상관계수(Autocorrelation coefficients), 자기상관그림(Autocorrelation Plot) (a) 데이터가 시간에 의존하는 것 없이 무작위성 (randomness, no time dependence in the data) 을 띠는지 확인하는데 사용됩니다. 그리고 (b) 시계열분석에서 Box-Jenkins ARIMA (Autoregressive Integrated Moving Average) 모형 식별 단계(model identification stage)에서 AR 차수를 구하는데도 사용됩니다. 


이번 포스팅에서는 R을 사용해서 


(1) 자기상관계수(Autocorrelation coefficients)를 계산하고, 

(2) 자기상관계수의 95% 신뢰구간(95% confidence interval)을 계산하고, 

(3) 자기상관그림(Autocorrelation Plot)을 그리는 방법


을 소개하겠습니다. 


본론에 들어가기 전에 먼저, 상관계수(correlation coefficients)와 자기상관계수(autocorrelation coefficients)를 비교해보겠습니다. 


상관계수와 자기상관계수의 유사한 점은 두 연속형 변수 간의 관계를 분석하다는 점입니다. 기본 개념은 공분산을 표준편차로 나누어서 표준화해주며, (자기)상관계수가 -1 ~ 1 사이의 값을 가지게 됩니다. 


상관계수와 자기상관계수가 서로 다른 점은, 상관계수는 특정 동일 시점을 횡단면으로 해서 Y와 다른 X1, X2, ... 변수들 간의 관계를 분석합니다. 반면에 자기상관계수는 동일한 변수(Yt, Yt-1, Yt-2, ...)의 서로 다른 시간 차이 (time lag) 를 두고 관계를 분석하는 것입니다. 


기존에 cross-sectional 관점의 Y와 X 변수들 간의 상관관계 분석에 많이 익숙해져 있는 경우에, 시계열 분석을 공부할 때 보면 자기 자신의 시간 차이에 따른 자기상관관계 관점이 처음엔 헷갈리고 잘 이해가 안가기도 합니다. 관점이 바뀐것일 뿐 어려운 개념은 아니므로 정확하게 이해하고 가면 좋겠습니다. 




  (1) 자기상관계수(Autocorrelation coefficients) 구하고, 자기상관그림 그리기


자기상관계수는 아래의 공식과 같이  로 계산합니다. 



간단한 예로서, 1 ~ 50까지의 시간 t에 대해서 싸인곡선 형태의 주기적인 파동을 띠는 값에 정규확률분포 N(0, 1) 에서 추출한 난수를 더하여 생성한 Y 데이터셋에 대해서 R로 위의 공식을 이용해서 각 시차(time lag)별로 자기상관계수를 구해보겠습니다. 



# sine wave with noise

set.seed(123)

noise <- 0.5*rnorm(50, mean=0, sd=1)

t <- seq(1, 50, 1)

y <- sin(t) + noise


plot(y, type="b")





자기상관계수를 구하고 자기상관그림을 그리는 가장 간단한 방법은 R stats 패키지에 들어있는 acf() 함수를 사용하는 것입니다. 



z <- acf(y, type=c("correlation"), plot=TRUE)


# autocorrelation coefficients per time lag k

round(z$acf, 4)


         [,1]

 [1,]  1.0000

 [2,]  0.3753

 [3,] -0.3244

 [4,] -0.5600

 [5,] -0.4216

 [6,]  0.2267

 [7,]  0.5729

 [8,]  0.3551

 [9,] -0.0952

[10,] -0.5586

[11,] -0.4582

[12,]  0.0586

[13,]  0.3776

[14,]  0.4330

[15,]  0.0961

[16,] -0.4588

[17,] -0.4759


# autocorrelation plot



좀 복잡하기는 합니다만, 위의 자기상관계수를 구하는 공식을 이용해서 R로 직접 사용자 정의 함수를 만들고, for loop 반복문을 이용해서 시차(time lag) k 별로 자기상관계수를 구할 수도 있습니다. 자기공분산을 분산으로 나누어서 표준화해주었습니다. 시차 0(time lag 0) 일 경우에는 자기자신과의 상관관계이므로 자기상관계수는 '1'이 되며, 표준화를 해주었기 때문에 자기상관계수는 -1 <= autocorr(Y, Yt-k) <= 1  사이의 값을 가집니다. 



# User Defined Function of Autocorrelation

acf_func <- function(y, lag_k){

  # y: input vector

  # lag_k : Lag order of k

  

  N = length(y) # total number of observations

  y_bar = mean(y)

  

  # Variance

  var = sum((y - y_bar)^2) / N

  

  # Autocovariance

  auto_cov = sum((y[1:(N-lag_k)] - y_bar) * (y[(1+lag_k):(N)] - y_bar)) / N

  

  # Autocorrelation coefficient = Autocovariance / Variance

  r_k = auto_cov / var

  

  return(r_k)

}



# Compute Autocorrelation per lag (from 0 to 9 in this case)

acf <- data.frame()

for (k in 0:(length(y)-1)){

  acf_k <- round(acf_func(y, lag_k = k), 4)

  acf[k+1, 'lag'] = k

  acf[k+1, 'ACF'] = acf_k

}


> print(acf)

   lag     ACF

1    0  1.0000

2    1  0.3753

3    2 -0.3244

4    3 -0.5600

5    4 -0.4216

6    5  0.2267

7    6  0.5729

8    7  0.3551

9    8 -0.0952

10   9 -0.5586

11  10 -0.4582

12  11  0.0586

13  12  0.3776

14  13  0.4330

15  14  0.0961

16  15 -0.4588

17  16 -0.4759

18  17 -0.0594

19  18  0.2982

20  19  0.3968

21  20  0.1227

22  21 -0.2544

23  22 -0.3902

24  23 -0.1981

25  24  0.1891

26  25  0.3701

27  26  0.1360

28  27 -0.1339

29  28 -0.2167

30  29 -0.1641

31  30  0.0842

32  31  0.2594

33  32  0.1837

34  33  0.0139

35  34 -0.1695

36  35 -0.1782

37  36 -0.0363

38  37  0.1157

39  38  0.1754

40  39  0.0201

41  40 -0.1428

42  41 -0.1014

43  42  0.0000

44  43  0.0750

45  44  0.0641

46  45 -0.0031

47  46 -0.0354

48  47 -0.0405

49  48 -0.0177

50  49 -0.0054

 



이번 예의 경우 시간의 흐름에 따라 싸인 파동 (sine wave) 형태를 보이는 데이터이므로, 당연하게도 시차를 두고 주기적으로 큰 양(+)의 자기상관, 큰 음(-)의 자기 상관을 보이고 있습니다. 




  (2) 자기상관계수의 95% 신뢰구간(95% confidence interval) 계산


자기상관계수의 95% 신뢰구간은 아래의 공식으로 간단하게 구할 수 있습니다. 


자기상관계수의 95% 신뢰구간 = 


이때 N은 샘플 개수, z 는 표준정규분포의 누적분포함수, 는 유의수준 5% 입니다. 


이미 익숙하게 알고 있다시피 유의수준 =0.05 에서의 z = 1.96 입니다. 이번 예에서의 관측치 개수 N = 50 입니다. 따라서 자기상관계수의 95% 신뢰구간은  이 됩니다. 



> # z quantile of 95% confidence level

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)

[1] 1.959964

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)

[1] -1.959964

> # sample size

> N <- length(y)

> print(N)

[1] 50

> # 95% confidence interval of autocorrelation

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] 0.2771808

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] -0.2771808




위의 acf() 함수로 그린 자기상관그림(autocorrelation plot)에서 ACF 축 +0.277 과 -0.277 위치에 수평으로 그어진 점선이 바로 95% 신뢰구간의 상, 하한 선입니다. 


위의 자기상관그림에서 보면 95% 신뢰구간의 상, 하한 점선을 주기적으로 벗어나므로 이 Y 데이터셋은 무작위적이지 않으며 (not random), 시계열적인 자기상관(autocorrelation)이 존재한다고 판단할 수 있습니다. 


만약, 무작위적인 데이터셋에 대해서 자기상관계수를 계산하고 자기상관그림을 그린다면, 시차(time lag)에 따른 자기상관계수가 모두 95% 신뢰구간 상, 하한 점선의 안쪽에 위치할 것입니다. 


많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 안녕하세요 2020.05.12 15:07  댓글주소  수정/삭제  댓글쓰기

    안녕하세요. 블로그에서 정말 큰 도움받고 있다는 감사의 말씀부터 전해드리고 싶습니다.

    다름이 아니라, 저희가 혼란 변수에 대해 분석을 해서 테이블을 작성해야하는데

    약의 종류에 따라서 신생아의 생존에 차이가 있는지를 분석하려고 합니다.

    이때 성별, 체중을 혼란 변수의 후보군으로 설정하였습니다.

    그러면 성별에 따라서 약의 종류가 다르게 설정되었는지와 성별에 따라서 생존에 차이가 있는지//(아마 Chi-square test를 수행하면 될까요?)

    체중에 따라서 약의 종류가 다르게 설정되었는지와 체중에 따라서 생존에 차이가 있는지
    //(이것은 T-test를 수행하면 될까요?)

    분석을 어떻게 시도해야할 지 정말 모르겠습니다..

    • R Friend R_Friend 2020.05.12 23:02 신고  댓글주소  수정/삭제

      안녕하세요.

      범주형 자료 분석(categorical data analysis) 중에서 '코크란-맨틀-핸챌 검정 (CMH test)'에 대해서 한번 찾아보시기 바랍니다.
      2 x 2 x K 분할표 (즉, 신생아 생존여부(survive, death) x 약(A, B) x 성별(남, 여) 분할표)에서 Z(성별)가 주어졌을 때 X(약)와 Y(생존여부)의 조건부 독립이라는 귀무가설을 검정할 때 사용합니다.

      혹은 로지스틱 회귀모형을 적합해보시는 것도 방법일 것 같네요. (y 는 생존여부 1, 0, X는 성별, 체중, 약 종류) 모형 적합 후에 모형 전체가 유의미할 경우 각 변수별 유의미 여부를 검정하고, 각 변수의 y에 대한 영향력도 파악할 수 있겠네요. (해석하는 방법 참고 => https://rfriend.tistory.com/401)