library(MVA)
library(mvtnorm)
library(gapminder)
library(GGally)
library(pheatmap)
library(mmrm)
library(nlme)
library(tidyverse)
library(here)
library(knitr)
library(purrr)
library(rmarkdown)
library(kableExtra)
library(flextable)4 공분산 행렬의 추정
4.1 공분산 행렬의 정의
공분산행렬은 다변량 분석에서 여러 변수들 간의 관계를 살펴볼 수 있는 가장 기본적이고 중요한 요소이다. 공분산행렬의 구조를 이해하고 적절히 추정하는 것은 다변량 데이터 분석에서 매우 중요하다. 본 장에서는 공분산 행렬의 다양한 구조적 형태와 이를 추정하는 방법들을 간단히 소개한다.
먼저 \(p\)-차원 확률벡터 \(\pmb X\) 가 다변량 정규분포를 따른다고 가정하자.
\[ \pmb X = \begin{bmatrix} X_1 \\ X_2 \\ \vdots \\ X_p \end{bmatrix} \sim N_p({\pmb \mu}, {\pmb \Sigma}) \]
공분산 행렬 \(\Sigma\) 는 다음과 같이 정의되며
\[ \pmb \Sigma = \mathrm{Var}(\pmb{X}) = \pmb{E}\big[(\pmb{X}- \pmb \mu)(\pmb{X}-\pmb \mu)^\top\big], \quad \pmb{X} \in \RR^p \] 다음과 같은 성질을 가지고있다.
- 대각원소는 각 변수의 분산, 비대각원소는 변수 간의 공분산을 나타낸다.
- 대칭 행렬(symmeric matrix): \(\pmb \Sigma^\top = \pmb \Sigma\)
- 양반정치 행렬(semi-positive matrix):
\[ \pmb{a}^{\top} {\pmb \Sigma} \pmb{a} \ge 0 \quad \text{ for all } \pmb{a} \in \RR^p\]
- 양반정치를 확인하는 방법은 공분산 행렬의 고유값이 모두 0 보다 같거나 커야한다.
이제 다변량 정규분포를 따르는 예제 데이터를 생성하고, 표본 공분산 행렬을 계산해 보자. 평균는 모두 0 이고 다음과 같은 공분산행렬을 가지는 6차원 다변량 정규분포를 고려한다. 다변량 정규분포에서 100개의 표본을 임의로 추출한 다음 표본 공분산 행렬을 구해보자.
\[ \Sigma = \begin{bmatrix} 1.0 & 0.4 & 0.2 & 0.5 & 0.1 & -0.2 \\ 0.4 & 1.0 & 0.3 & 0.4 & -0.3 & 0.01 \\ 0.2 & 0.3 & 1.0 & 0.3 & 0.2 & -0.1 \\ 0.5 & 0.4 & 0.3 & 1.0 & 0.4 & -0.2 \\ 0.1 & -0.3 & 0.2 & 0.4 & 1.0 & 0.2 \\ -0.2 & 0.01 & -0.1 & -0.2 & 0.2 & 1.0 \end{bmatrix} \]
# 예제 데이터 (다변량 정규)
# 공분산 행렬
Sigma_true <- matrix(c(
1.0, 0.4, 0.2, 0.5, 0.1, -0.2,
0.4, 1.0, 0.3, 0.4, -0.3, 0.01,
0.2, 0.3, 1.0, 0.3, 0.2, -0.1,
0.5, 0.4, 0.3, 1.0, 0.4, -0.2,
0.1, -0.3, 0.2, 0.4, 1.0, 0.2,
-0.2, 0.01, -0.1, -0.2, 0.2, 1.0
), 6, 6, byrow = TRUE)
Sigma_true [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0 0.40 0.2 0.5 0.1 -0.20
[2,] 0.4 1.00 0.3 0.4 -0.3 0.01
[3,] 0.2 0.30 1.0 0.3 0.2 -0.10
[4,] 0.5 0.40 0.3 1.0 0.4 -0.20
[5,] 0.1 -0.30 0.2 0.4 1.0 0.20
[6,] -0.2 0.01 -0.1 -0.2 0.2 1.00# 고유값 확인- 공분산행렬의 양반정치 점검
eigen(Sigma_true)$values [1] 2.1507035 1.3633588 1.0108471 0.8176237 0.4860156 0.1714512# 다변량 정규분포에서 표본 생성
set.seed(123) # 재현 가능성
X <- rmvnorm(n = 100, mean = c(0, 0, 0, 0, 0, 0), sigma = Sigma_true)
colnames(X) <- c("x1", "x2", "x3", "x4", "x5", "x6")
# 표본의 일부 보기
head(X) x1 x2 x3 x4 x5 x6
[1,] -0.60697291 -0.016012988 1.3739779 -0.09555296 0.5483513 1.6489336
[2,] 0.06807651 -1.514672553 -0.7643606 -0.39888978 1.2909460 0.4933002
[3,] 0.98823956 0.285412440 -0.1300729 2.00665652 0.5701723 -2.0653381
[4,] 0.46989680 -0.333468735 -1.1591075 -0.42188577 -1.0883766 -0.8543093
[5,] -0.99065257 -1.222874588 0.3448429 -0.68099634 -0.4556374 1.0022332
[6,] 0.58628469 -0.005235656 1.0113511 1.05145508 1.2414359 0.5797083# 표본 공분산 행렬
S_hat <- cov(X)
S_hat x1 x2 x3 x4 x5 x6
x1 0.88057500 0.33490828 0.1872037 0.3207651 -0.05887987 -0.32665921
x2 0.33490828 0.88314376 0.3453332 0.2363243 -0.37813166 -0.06704301
x3 0.18720366 0.34533319 1.0006850 0.1587787 0.10083233 -0.01456920
x4 0.32076508 0.23632432 0.1587787 0.8100167 0.28207315 -0.44437452
x5 -0.05887987 -0.37813166 0.1008323 0.2820731 0.82579538 0.04660305
x6 -0.32665921 -0.06704301 -0.0145692 -0.4443745 0.04660305 1.06484897# 추정량과 실제 공분산 행렬 비교: 각 분산과 공분산의 상대적인 오차(%)
round(100 * (S_hat - Sigma_true) / (Sigma_true), 2) x1 x2 x3 x4 x5 x6
x1 -11.94 -16.27 -6.40 -35.85 -158.88 63.33
x2 -16.27 -11.69 15.11 -40.92 26.04 -770.43
x3 -6.40 15.11 0.07 -47.07 -49.58 -85.43
x4 -35.85 -40.92 -47.07 -19.00 -29.48 122.19
x5 -158.88 26.04 -49.58 -29.48 -17.42 -76.70
x6 63.33 -770.43 -85.43 122.19 -76.70 6.484.2 공분산 행렬의 형태
- 일반형 공분산
p-차원 공분산 행렬의 일반적인 구조(general 또는 unstructured) 는 다음과 같으며 대칭 행렬이기 떄문에 모수(parameter)의 개수는 \(p(p-1)/2\) 개이다.
\[ \pmb \Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \cdots & \sigma_{1p} \\ \sigma_{12} & \sigma_{22} & \cdots & \sigma_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{1p} & \sigma_{2p} & \cdots & \sigma_{pp} \end{bmatrix} \tag{4.1}\]
- 독립과 등분산
다변량 정규분포에서 모든 변수 \(X_1, X_2, \dots, X_p\) 가 독립인 경우 공분산이 0인 것과 동일하기 떄문에 다음과 같은 대각행렬의 형태를 가지게 된다.
\[ \pmb \Sigma = \begin{bmatrix} \sigma_1^2 & 0 & \cdots & 0 \\ 0 & \sigma_2^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma_p^2 \end{bmatrix} \]
이 경우 다음과 같이 각 확률변수가 독립적으로 분산이 다른 일변량 정규분포를 따른다고 할 수 있다.
\[ X_i \sim_{indep} N(\mu_i, \sigma_i^2), \quad i=1,2,\dots,p \] 더 나아가 모든 변수의 분산이 동일한 경우(등분산)에는 다음과 같은 구형(spherical) 형태의 공분산 행렬을 가진다.
\[ \pmb \Sigma = \begin{bmatrix} \sigma^2 & 0 & \cdots & 0 \\ 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^2 \end{bmatrix} = \sigma^2 {\pmb I}_p \]
- 균등 상관 구조
특별한 구조를 고려하는 경우, 공분산 행렬의 형태로 가장 자주 사용되는 구조가 균등 상관 구조(Compound Symmetry, CS) 이다. 이 형태는 분산이 모두 동일하고 공분산도 모두 동일한 형태이며 분산과 공분산은 다르게 설정된다. 즉 모든 변수의 분산이 \(\sigma^2\) 이고 모든 변수 쌍들의 공분산이 \(\rho\sigma^2\) 인 경우이다. 따라서 \(\rho\) 는 상관 계수이다. 따라서 모든 변수의 상관계수도 동일한 형태를 가지는 구조이다.
\[ cor(X_i, X_j) = \rho \quad \text{ for all } i,j \]
\[ \pmb \Sigma_{\text{CS}} = \begin{bmatrix} \sigma^2 & \rho\sigma^2 & \cdots & \rho\sigma^2 \\ \rho\sigma^2 & \sigma^2 & \cdots & \rho\sigma^2 \\ \vdots & \vdots & \ddots & \vdots \\ \rho\sigma^2 & \rho\sigma^2 & \cdots & \sigma^2 \end{bmatrix} = \sigma^2 \big [ (1-\rho) \pmb I_p + \rho \pmb J_p \big ] \tag{4.2}\]
위의 식에서 \(\pmb J_p\) 는 모든 원소가 1인 \(p\times p\) 행렬이다.
4.2.1 AR(1) 구조
공분산 행렬의 또 다른 구조적 형태로는 AR(1) 구조가 있다. 이 형태는 시계열 자료에서 자주 사용되는 구조로서, 인접한 변수들 간의 상관관계가 멀어질수록 지수적으로 감소하는 특징을 가진다.
\[ cor(X_i, X_j) = \rho^{|i-j|} \]
AR(1) 구조의 공분산 행렬은 다음과 같은 형태를 가진다.
\[ \pmb \Sigma_{AR} = \begin{bmatrix} \sigma^2 & \rho\sigma^2 & \rho^2\sigma^2 & \cdots & \rho^{p-1}\sigma^2 \\ \rho\sigma^2 & \sigma^2 & \rho\sigma^2 & \cdots & \rho^{p-2}\sigma^2 \\ \rho^2\sigma^2 & \rho\sigma^2 & \sigma^2 & \cdots & \rho^{p-3}\sigma^2 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{p-1}\sigma^2 & \rho^{p-2}\sigma^2 & \rho^{p-3}\sigma^2 & \cdots & \sigma^2 \end{bmatrix} \]
따라서 AR(1) 구조를 가진 다변량 확률 벡터의 상관 계수 행렬은 다음과 같이 나타낼 수 있다.
\[ \pmb R = \begin{bmatrix} 1 & \rho & \rho^2 & \cdots & \rho^{p-1} \\ \rho & 1 & \rho & \cdots & \rho^{p-2} \\ \rho^2 & \rho & 1 & \cdots & \rho^{p-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \rho^{p-1} & \rho^{p-2} & \rho^{p-3} & \cdots & 1 \end{bmatrix} \]
4.2.2 블록 대각 구조
공분산 행렬의 또 다른 구조적 형태로는 블록 대각(Block Diagonal) 구조가 있다. 이 형태는 변수들이 여러 개의 그룹으로 나누어져 있고, 각 그룹 내에서는 변수들 간에 상관관계가 존재하지만, 그룹 간에는 상관관계가 없는 경우에 적합하다. 블록 대각 구조의 공분산 행렬은 다음과 같은 형태를 가진다.
\[ \pmb \Sigma = \begin{bmatrix} \pmb \Sigma_{1} & \pmb 0 & \pmb 0\\ \pmb 0 & \pmb \Sigma_{2} & \pmb 0 \\ \pmb 0 & \pmb 0 & \pmb \Sigma_{3} \end{bmatrix} \]
4.3 공분산의 추정
4.3.1 표본 공분산 행렬
만약 \(n\) 개의 표본 \(\pmb X_1, \pmb X_2, \dots, \pmb X_n\) 이 관측되었다고 하자.
만약 분포의 가정에서 공분산이 제약이 없는 일반적인 형태 식 4.1 라고 한다면 다음과 같은 표본 공분산 행렬을 이용하여 추정한다. 이 추정량은 불편추정량(unbiased estimator)이고 동시에 \(n\) 이 충분히 크면 최대가능도 추정량과 동일하다고 볼 수 있다.
\[ \hat{\Sigma} = \frac{1}{n-1}\sum_{i=1}^n (\pmb{X}_i-\bar{\pmb{X}})(\pmb{X}_i-\bar{\pmb{X}})^\top \tag{4.3}\]
특별한 구조를 가진 공분산 행렬은 최대가능도 추정을 이용하여 추정할 수 있다.
4.3.2 예제: 반복측정자료
이 장에서는 의학통계에 자주 사용되는 반복측정자료 형태의 다변량 확률벡터에 대하여 균등 상관 구조 형태의 공분산 행렬을 추정하는 예제에 대해서 살펴 보자.
확률변수 \(X_i\) 는 \(i\) 시점에서 순서대로 5번 관측한 자료이며 한 명의 개체가 각 시점에서 반응값을 측정한다고 하자. 따라서 반복으로 측정한 확률 변수 \(X_1, X_2, \dots, X_5\) 는 독립이 아니다.
이제 확률 벡터 \(\pmb X = (X_1, X_2, \dots, X_5)^t\) 가 다변량 정규분포를 따른다고 가정하자. 5-차원의 다변량 정규 확률 벡터를 고려하고 각 변수의 평균은 시간을 나타내는 시점 (\(i\))에 비례하게 다음과 같이 정의한다.
\[ \mu_i = E(X_i) = 0.5 + 0.1 (i - 1) \]
공분산 행렬은 식 4.2 의 균등 상관 구조를 가지며 분산은 모두 1 이고 상관계수는 0.6으로 가정한다. 표본의 개수는 100 개이다.
먼저 주어진 평균벡터와 공분산 행렬에 대하여 분포를 정의하고 가상의 자료를 임의로 추출하는 R 코드를 고려한다.
set.seed(121)
n <- 100 # 표본 수
p <- 5 # 확률 벡터의 차원
# 각 개체의 id 와 시간 변수 생성
id <- factor(rep(1:n, each = p))
time <- rep(1:p, times = n)
mu <- 0.5 + 0.1 * (rep(1:p, times = 1) - 1) # 평균 벡터
# p- 차원 CS 공분산을 만드는 함수
make_Sigma_CS <- function(p, sigma2 = 1, rho = 0.3) {
if (rho <= -1/(p - 1) || rho >= 1) {
stop("rho must be in (-1/(p-1), 1) to ensure positive definiteness.")
}
J <- matrix(1, p, p) # 모든 원소가 1인 행렬
Sigma <- sigma2 * ((1 - rho) * diag(p) + rho * J)
return(Sigma)
}
# 균등 상관 구조의 공분산 행렬
Sigma_true <- make_Sigma_CS(p, sigma2 = 1, rho = 0.6)
Sigma_true [,1] [,2] [,3] [,4] [,5]
[1,] 1.0 0.6 0.6 0.6 0.6
[2,] 0.6 1.0 0.6 0.6 0.6
[3,] 0.6 0.6 1.0 0.6 0.6
[4,] 0.6 0.6 0.6 1.0 0.6
[5,] 0.6 0.6 0.6 0.6 1.0이제 100 개의 표본을 다변량 정규분포에서 임의로 추출하고 wide 형식의 자료로 변환해 보자.
Xmat <- rmvnorm(n = n, mean = mu, sigma = Sigma_true) # 다변량 표본 추출
colnames(Xmat) <- paste0("X", 1:p)
df_wide <- as.data.frame(Xmat) # wide 형식의 자료
head(df_wide) X1 X2 X3 X4 X5
1 0.1412804 0.47132670 0.5836000 0.55150916 0.25141129
2 1.9259813 1.19236353 1.0678925 0.50890233 1.88213477
3 1.4160815 0.08403608 0.4753524 0.47093252 1.27335173
4 1.1549582 -0.30705522 1.4909395 0.08919076 1.23971146
5 1.2631167 1.31918169 1.2023846 1.90247109 1.83556681
6 -0.6336252 0.47159709 -0.6432958 -0.37359215 -0.06331619추출된 표본을 이용하여 식 4.3 에 주어진 표본 공분산 행렬을 계산해 보자. 표본 공분산 행렬은 제약조건이 없는 형태로 나타나기 때문에 모든 분산과 공분산이 각각 다르게 추정된다.
cov(df_wide) X1 X2 X3 X4 X5
X1 0.7180133 0.4216849 0.3886342 0.3553769 0.4661159
X2 0.4216849 0.8704348 0.4896497 0.4478777 0.4441199
X3 0.3886342 0.4896497 0.8766630 0.4369850 0.5013652
X4 0.3553769 0.4478777 0.4369850 0.7887407 0.4043893
X5 0.4661159 0.4441199 0.5013652 0.4043893 0.9236021이제 균등 상관 구조를 가지는 공분산을 추정할 수 있는 방법을 알아보자. 먼저 자료는 측정 시간을 변수로 하는 긴 형식의 자료로 생성한다.
# pivot_longer 함수를 이용하여 긴 형식의 자료를 생성
df_long <- df_wide %>%
mutate(id = factor(1:n)) %>%
pivot_longer(cols = starts_with("X"),
names_to = "time",
values_to = "y",
names_prefix = "X") %>%
mutate(time = as.integer(time))
head(df_long, 10)# A tibble: 10 × 3
id time y
<fct> <int> <dbl>
1 1 1 0.141
2 1 2 0.471
3 1 3 0.584
4 1 4 0.552
5 1 5 0.251
6 2 1 1.93
7 2 2 1.19
8 2 3 1.07
9 2 4 0.509
10 2 5 1.88 공분산 추정을 위하여 nlme 패키지에 있는 다변량 확률변수의 회귀식을 추정하는 함수 nlme 를 사용하여고 하며, 균등 상관 구조를 가지는 공분산을 가정한다.
다음 추정된 다변량 회귀모형의 결과를 보면 식 4.2 의 균등 상관 구조에서 분산 \(\sigma\) 의 추정값은 \(0.8362 = (0.9144)^2\), 상관계수 \(\rho\) 의 추정값은 \(0.52078\) 로 나타난다. 참고로 아래 결과에서 rho 는 공통 상관계수의 추정값,Residual standard error 는 표준편차 추정값이다.
fit_gls_cs <- gls(y ~ time, data = df_long,
correlation = corCompSymm(form = ~ 1 | id))
summary(fit_gls_cs) Generalized least squares fit by REML
Model: y ~ time
Data: df_long
AIC BIC logLik
1163.245 1180.088 -577.6226
Correlation Structure: Compound symmetry
Formula: ~1 | id
Parameter estimate(s):
Rho
0.5207752
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.3750164 0.09360729 4.006274 1e-04
time 0.0817480 0.02001772 4.083779 1e-04
Correlation:
(Intr)
time -0.642
Standardized residuals:
Min Q1 Med Q3 Max
-2.76158379 -0.71458351 0.04453439 0.76834584 3.26193595
Residual standard error: 0.9144187
Degrees of freedom: 500 total; 498 residual# rho 는 공통 상관계수, Residual standard error 는 표준편차 추정값위에서 주어진 결과를 가지고 공분산 행렬의 추정값을 다음과 같이 구할 수 있다.
getVarCov(fit_gls_cs) Marginal variance covariance matrix
[,1] [,2] [,3] [,4] [,5]
[1,] 0.83616 0.43545 0.43545 0.43545 0.43545
[2,] 0.43545 0.83616 0.43545 0.43545 0.43545
[3,] 0.43545 0.43545 0.83616 0.43545 0.43545
[4,] 0.43545 0.43545 0.43545 0.83616 0.43545
[5,] 0.43545 0.43545 0.43545 0.43545 0.83616
Standard Deviations: 0.91442 0.91442 0.91442 0.91442 0.91442