- 참고
- R 응용 다변량분석, 나종화 저
- SAS를 이용한 다변량 통계 분석, 김재희 저
1. 다변량 정규분포
- 일변량 정규분포를 다차원으로 확장한 다변량 정규분포(Multivariate Normal Distribution)는 다변량 분석 절차에서 중요한 역할을 한다.
- 실제로 관측된 데이터가 다변량 정규분포를 따르지 않더라도, 관측
개수(개체 수)가 충분히 크면
표본평균벡터는 중심극한정리(Central Limit Theorem)에 의해 근사적으로 다변량 정규분포를 따른다.- 중심극한정리 : 모집단 분포에 상관없이 표본의 크기가 충분히 크면 표본평균은 근사적으로 정규분포를 따른다.
- 다변량 데이터 분석에서 데이터의 분포는 주로 다변량 정규분포를 가정하며, 이렇게 가정한 경우 수학적으로 다루기 쉬운 상황이 되어 통계적 추론에서 많이 이점이 있다.
- 다변량 데이터의 복잡성을 다루는 연구자들은 종종 다변량 데이터가 다변량 정규분포를 따르는 지에 대한 검증을 요구한다.
1-1. 다변량 정규분포의 확률밀도함수
- 일변량 확률변수 \(X\)가 평균이 \(\mu\)이고 분산이 \(\sigma^2\)인 정규분포를 따를 때, \(X \sim N(\mu, \sigma^2)\)로 표기하며 확률변수 \(X\)는 다음과 같은 확률밀도함수(Probability Density Function)를 가진다. \[ \begin{align*} f(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{ -\frac{(x-\mu)^2}{2\sigma^2} }, \;\; -\infty < x < \infty. \end{align*} \]
- 확률벡터 \(\mathbf{X}=(X_1, \ldots,
X_p)^T\)가 평균이 \(\boldsymbol{\mu}\)이고 공분산행렬이 \(\boldsymbol{\Sigma}\)인 다변량 정규분포를
따를 때, 확률밀도함수는 다음과 같으며 \(\mathbf{X}\sim N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma})\)로 표기한다. \[
\begin{align*}
f(\mathbf{x}) = \frac{1}{(2\pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}}e^{
-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T
\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})}
\end{align*}
\]
- 만약 \(p=2\)인 경우, 즉, \(\boldsymbol{\mu}=(\mu_1,\; \mu_2),
\;\boldsymbol{\Sigma}= \begin{pmatrix}\sigma^2_1 &
\rho\sigma_1\sigma_2\\\rho\sigma_1\sigma_2 &
\sigma^2_2\end{pmatrix}\)일 때, 이변량 정규분포(Bivariate Normal
Distribution)의 확률밀도함수는 다음과 같다. \[
\begin{align*}
&f(x_1, x_2) =
\frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}\exp{\left[-\frac{1}{2(1-\rho^2)}\left\lbrace
\left(\frac{x_1-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x_1-\mu_1}{\sigma_1}\right)\left(\frac{x_2-\mu_2}{\sigma_2}\right)+\left(\frac{x_2-\mu_2}{\sigma_2}\right)^2 \right\rbrace \right]},\\
\;\;\; &-\infty < x_1 < \infty, \;\; -\infty < x_2 <
\infty.
\end{align*}
\]
- 여기서, \(-\infty < \mu_1, \;\mu_2<\infty,\; \sigma_1^2 >0,\; \sigma_2^2>0,\; -1<\rho<1\)이다.
- 예를 들어, \(\boldsymbol{\mu}=(0, 0), \;\boldsymbol{\Sigma} = \begin{pmatrix} 1 & \rho\\\rho & 1 \end{pmatrix}\)인 경우의 확률밀도함수를 그리는 방법은 다음과 같다.
- 만약 \(p=2\)인 경우, 즉, \(\boldsymbol{\mu}=(\mu_1,\; \mu_2),
\;\boldsymbol{\Sigma}= \begin{pmatrix}\sigma^2_1 &
\rho\sigma_1\sigma_2\\\rho\sigma_1\sigma_2 &
\sigma^2_2\end{pmatrix}\)일 때, 이변량 정규분포(Bivariate Normal
Distribution)의 확률밀도함수는 다음과 같다. \[
\begin{align*}
&f(x_1, x_2) =
\frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}\exp{\left[-\frac{1}{2(1-\rho^2)}\left\lbrace
\left(\frac{x_1-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x_1-\mu_1}{\sigma_1}\right)\left(\frac{x_2-\mu_2}{\sigma_2}\right)+\left(\frac{x_2-\mu_2}{\sigma_2}\right)^2 \right\rbrace \right]},\\
\;\;\; &-\infty < x_1 < \infty, \;\; -\infty < x_2 <
\infty.
\end{align*}
\]
# 이변량 정규분포의 확률밀도 함수 그리기(rho = 0)
mu1 <- 0 # x1의 평균NAmu2 <- 0 # x2의 평균NAs11 <- 1 # x1의 분산rho <- 0 # Rho
s22 <- 1 # x2의 분산
x1 <- seq(-3, 3, length = 50) # x1값값
x2 <- seq(-3, 3, length = 50) # x2값값
# 이변량 정규분포의 확률밀도함수gaussian_func <- function(x1, x2){
term1 <- 1/(2*pi*sqrt(s11*s22*(1-rho^2)))
term2 <- -1/(2*(1-rho^2))
term3 <- (x1 - mu1)^2/s11
term4 <- (x2 - mu2)^2/s22
term5 <- 2*rho*((x1 - mu1)*(x2 - mu2))/(sqrt(s11)*sqrt(s22))
term1*exp(term2*(term3 + term4 - term5))
}
# 확률밀도함수 그림림
persp(x1, x2,
outer(x1, x2, gaussian_func), # 확률밀도함수 값 계산산
zlab = "",
main = expression(paste("Bivariate Normal Density for ", rho==0)),
theta = 30, phi = 10)
# 등고선 그림그림
contour(x1, x2,
outer(x1, x2, gaussian_func), # 확률밀도함수 값 계산산
xlab = "x1", ylab = "x2",
main = expression(paste("Contour of Bivariate Normal Density for ", rho==0)))
# 이변량 정규분포의 확률밀도 함수 그리기(rho = 0.75)
mu1 <- 0 # x1의 평균NAmu2 <- 0 # x2의 평균NAs11 <- 1 # x1의 분산rho <- 0.75 # Rho
s22 <- 1 # x2의 분산
x1 <- seq(-3, 3, length = 50) # x1값값
x2 <- seq(-3, 3, length = 50) # x2값값
# 이변량 정규분포의 확률밀도함수gaussian_func <- function(x1, x2){
term1 <- 1/(2*pi*sqrt(s11*s22*(1-rho^2)))
term2 <- -1/(2*(1-rho^2))
term3 <- (x1 - mu1)^2/s11
term4 <- (x2 - mu2)^2/s22
term5 <- 2*rho*((x1 - mu1)*(x2 - mu2))/(sqrt(s11)*sqrt(s22))
term1*exp(term2*(term3 + term4 - term5))
}
# 확률밀도함수 그림림
persp(x1, x2,
outer(x1, x2, gaussian_func), # 확률밀도함수 값 계산
zlab = "",
main = expression(paste("Bivariate Normal Density for ", rho==0.75)),
theta = 30, phi = 10)
# 등고선 그림그림
contour(x1, x2,
outer(x1, x2, gaussian_func), # 확률밀도함수 값 계산산
xlab = "x1", ylab = "x2",
main = expression(paste("Contour of Bivariate Normal Density for ", rho==0.75)))
1-2. 다변량 정규분포의 성질
1-2-1. 선형결합식의 분포
- \(\mathbf{X}\sim N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma})\)일 때,
- 상수벡터 \(\mathbf{a}=(a_1, \ldots, a_p)^T\)에 대해 선형결합식 \(\mathbf{a}^T\mathbf{X} = a_1X_1+ \ldots, a_pX_p\)는 평균이 \(\mathbf{a}^T\boldsymbol{\mu}\)이고 분산이 \(\mathbf{a}^T\boldsymbol{\Sigma}\mathbf{a}\)인 일변량 정규분포를 따른다. 즉, \(\mathbf{a}^T\mathbf{X} \sim N(\mathbf{a}^T\boldsymbol{\mu}, \mathbf{a}^T\boldsymbol{\Sigma}\mathbf{a})\).
- \(q\times p\) 상수행렬 \(\mathbf{A}\), \(p\times 1\) 상수벡터 \(\mathbf{d}\)에 대해 \(\mathbf{AX}\sim N_q(\mathbf{A}\boldsymbol{\mu}, \mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^T)\)이고, \(\mathbf{X+d}\sim N_p(\boldsymbol{\mu}+\mathbf{d}, \boldsymbol{\Sigma})\)이다.
1-2-2. 표준화 변수
- \(\mathbf{X}\sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)일 때, 표준화벡터 \(\mathbf{Z}=\boldsymbol{\Sigma}^{-1/2}(\mathbf{X}-\boldsymbol{\mu})\)는 \(N_p(\mathbf{0}, \mathbf{I})\)를 따른다.
1-2-3. 다변량 정규분포의 주변 및 조건부확률밀도함수
- \(p\)-변량 정규분포를 따르는
확률벡터 \(\mathbf{X}\)를 \(q\times 1\)벡터 \(\mathbf{X}_{1}\)과 \((p\times q)\times 1\)벡터 \(\mathbf{X}_2\)로 분할할 때,
- \(\mathbf{X}_{p\times 1}=\begin{pmatrix}\mathbf{X}_1 \\ \mathbf{X}_2 \end{pmatrix}\)는 평균벡터 \(\boldsymbol{\mu}_{p\times 1}= \begin{pmatrix}\boldsymbol{\mu}_1 \\ \boldsymbol{\mu}_2 \end{pmatrix}\)와 공분산행렬 \(\boldsymbol{\Sigma}_{p\times p}=\begin{pmatrix}\boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12}\\ \boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22} \end{pmatrix}\)을 가지며 \(\mathbf{X}_1\)은 \(q\)-변량 정규분포를 따른다. 즉, \(\mathbf{X}_1 \sim N_q(\boldsymbol{\mu_1}, \boldsymbol{\Sigma}_{11})\).
- \(\mathbf{X}_{2}=\mathbf{x}_{2}\)가 주어졌을 때 \(\mathbf{X}_{1}\)의 조건부 분포 \(\mathbf{X}_{1}|\mathbf{X}_{2}=\mathbf{x}_{2}\)는 \(N(\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}(\mathbf{x}_2-\boldsymbol{\mu}_2)\),\(\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21})\)이다.
1-2-4. 다변량 정규분포의 독립성
- 확률벡터 \(\mathbf{X}\)의
부분벡터인 \(\mathbf{X}_{1}\)과 \(\mathbf{X}_2\)에 대해, 공분산행렬 \(\boldsymbol{\Sigma}_{12}=\mathbf{0}\)이면
\(\mathbf{X}_{1}\)과 \(\mathbf{X}_2\)는 서로 독립이다.
\(\Leftrightarrow\) \(\mathbf{X}_{1}\)과 \(\mathbf{X}_2\)이 독립이면, 공분산행렬 \(\boldsymbol{\Sigma}_{12}=\mathbf{0}\)이다.
1-2-5. 정규 확률벡터 합의 분포
- \(q\times 1\) 크기의 확률벡터 \(\mathbf{X}_{1} \sim N_q(\boldsymbol{\mu}_{1},
\boldsymbol{\Sigma}_{11})\)과 \(\mathbf{X}_2\sim N_q(\boldsymbol{\mu}_{2},
\boldsymbol{\Sigma}_{22})\)가 서로 독립일 때,
- 두 벡터의 합 : \(\mathbf{X}_1 +\mathbf{X}_2 \sim N_q(\boldsymbol{\mu}_{1}+\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_{11} + \boldsymbol{\Sigma}_{22})\)
- 두 벡터의 차 : \(\mathbf{X}_1 -\mathbf{X}_2 \sim N_q(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_{11} + \boldsymbol{\Sigma}_{22})\)
1-2-6. 표본평균벡터와 표본공분산행렬의 분포
- \(\mathbf{X}_1,
\ldots,\mathbf{X}_i=(X_{i1}, \ldots, X_{ip})^T, \ldots,
\mathbf{X}_n\)을 \(N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma})\)로부터 추출한 확률표본이라고 할 때, 다음이
성립한다.
- 표본평균벡터 \(\bar{\mathbf{X}}=(\bar{X}_1, \ldots,
\bar{X}_p)^T\)는 \(N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma}/n)\)을 따른다.
- 여기서 \(\bar{X}_k=\frac{1}{n}\sum_{j=1}^n X_{jk}\)이다.
- 표본공분산행렬 \(\boldsymbol{S}\)에
대해 \((n-1)\boldsymbol{S}\)는 \(W_p(n-1, \boldsymbol{\Sigma})\)를 따른다.
- 여기서 \(W_p(n-1, \boldsymbol{\Sigma})\)는 자유도가 \((n-1)\)인 위샤트 분포이다.
- \(\boldsymbol{Z}_1, \ldots, \boldsymbol{Z}_n\)을 \(N_p(\boldsymbol{0}, \boldsymbol{\Sigma})\)로부터 추출한 확률표본일 때, 확률행렬 \(\sum_{i=1}^n \boldsymbol{Z}_i \boldsymbol{Z}_i^T\)은 자유도가 \(n\)인 위샤트 분포를 따르며, 기호로는 \(W_p(n, \boldsymbol{\Sigma})\)로 나타낸다.
- 위샤트 분포는 카이제곱분포의 다변량 형태로 이해할 수 있다.
- \(\bar{\mathbf{X}}\)와 \(\boldsymbol{S}\)는 서로 독립이다.
- \((\bar{\mathbf{X}}-\boldsymbol{\mu})^T\left(\frac{\boldsymbol{\Sigma}}{n}\right)^{-1}(\bar{\mathbf{X}}-\boldsymbol{\mu})\)는 자유도가 \(p\)인 카이제곱분포를 따른다.
- \(n-p\to\infty\)이면 \((\bar{\mathbf{X}}-\boldsymbol{\mu})^T\left(\frac{\boldsymbol{S}}{n}\right)^{-1}(\bar{\mathbf{X}}-\boldsymbol{\mu})\)는 근사적으로 자유도가 \(p\)인 카이제곱분포를 따른다.
- 표본평균벡터 \(\bar{\mathbf{X}}=(\bar{X}_1, \ldots,
\bar{X}_p)^T\)는 \(N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma}/n)\)을 따른다.
1-3. 다변량 데이터의 정규성 평가
1-3-1. 카이제곱그림
- \(\mathbf{X}\sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)이고 \(|\boldsymbol{\Sigma}|>0\)일 때, \((\mathbf{X}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{X}-\boldsymbol{\mu}) \sim \chi^2(p)\)이다.
- \(\mathbf{X}_1, \ldots,
\mathbf{X}_n\)을 \(N_p(\boldsymbol{\mu},
\boldsymbol{\Sigma})\)로부터 추출한 확률표본이라고 할 때, \(\mathbf{X}_i=(X_{i1}, \ldots, X_{ip})^T\)의
\(\bar{\mathbf{X}}=(\bar{X}_1, \ldots,
\bar{X}_p)^T\)로부터의 마할라노비스 제곱거리(Mahalanobis Squared
Distance) 또는 일반화 제곱거리(Generalized Squared Distance) \(d_i^2\)은
근사적으로 \(\chi^2(p)\)를 따른다.- 즉, \(d_i^2 =
(\mathbf{X}_i-\bar{\mathbf{X}})^T\boldsymbol{S}^{-1}(\mathbf{X}_i-\bar{\mathbf{X}})
\overset{\underset{\mathrm{.}}{}}{\sim} \chi^2(p)\).
- 여기서 \(\boldsymbol{S}\)는 표본공산분행렬이다.
- 즉, \(d_i^2 =
(\mathbf{X}_i-\bar{\mathbf{X}})^T\boldsymbol{S}^{-1}(\mathbf{X}_i-\bar{\mathbf{X}})
\overset{\underset{\mathrm{.}}{}}{\sim} \chi^2(p)\).
- 마할라노비스 제곱거리 \(d_i^2\)는 독립이거나 정확한 카이제곱분포를 따르지는 않으나, 관측값이 다변량 정규분포를 따르는 지를 판단하는 카이제곱그림(Chi-square Plot)을 제공한다.
- 카이제곱그림은 \(d_i^2\)를 순서대로 늘어놓고 누적확률이 \(\frac{i-0.5}{n}\)에 해당하는 자유도가 \(p\)인 카이제곱분포로부터의 분위수를 구해 (분위수, 거리)의 이차원 공간에 점을 찍고 기울기가 “1”인 직선과 비교하여 다변량 정규성을 검토하는 그림이다.
- 카이제곱그림은 다음의 두 단계로 그려진다.
- 마할라노비스 제곱거리 \(d_i^2\)를
구하고 크기순으로 정렬한다.
- 즉, \(d_{(1)}^2\le d_{(2)}^2\le \ldots \le d_{(n)}^2\).
- \(\left\lbrace \chi^2_p
\left(\frac{i-0.5}{n} \right), d^2_{(i)}\right\rbrace\) 점을
그린다.
- 즉, \(d_{(i)}^2\)와 이론적으로 정규분포를 따를 때의 사분위수를 구해 비교하여 다변량 정규성을 검토한다.
- 마할라노비스 제곱거리 \(d_i^2\)를
구하고 크기순으로 정렬한다.
- 관측된 다변량 데이터 \(\mathbf{x}_1,
\ldots, \mathbf{x}_n\)가 다변량 정규분포를 만족한다면
마할라노비스 제곱거리 \(d_i^2\)는
근사적으로 카이제곱분포를 따르기 때문에 카이제곱그림은 기울기가 “1”인
직선이 된다.
- 즉, 카이제곱그림이 직선에 가까우면 관측 데이터가 다변량 정규분포를 잘 따른다고 판단한다.
- 또한, 직선으로부터 크게 벗어나는 점을 이상치 데이터로 판단한다.
# 예제 11
pacman::p_load("MVA") # For Data
There are binary versions available but the source versions
are later:
binary source needs_compilation
HSAUR2 1.1-18 1.1-19 FALSE
MVA 1.0-7 1.0-8 FALSE'data.frame': 41 obs. of 7 variables:
$ SO2 : int 46 11 24 47 11 31 110 23 65 26 ...
$ temp : num 47.6 56.8 61.5 55 47.1 55.2 50.6 54 49.7 51.5 ...
$ manu : int 44 46 368 625 391 35 3344 462 1007 266 ...
$ popul : int 116 244 497 905 463 71 3369 453 751 540 ...
$ wind : num 8.8 8.9 9.1 9.6 12.4 6.5 10.4 7.1 10.9 8.6 ...
$ precip : num 33.36 7.77 48.34 41.31 36.11 ...
$ predays: int 135 58 115 111 166 148 122 132 155 134 ...USairpollution
SO2 temp manu popul wind precip predays
Albany 46 47.6 44 116 8.8 33.36 135
Albuquerque 11 56.8 46 244 8.9 7.77 58
Atlanta 24 61.5 368 497 9.1 48.34 115
Baltimore 47 55.0 625 905 9.6 41.31 111
Buffalo 11 47.1 391 463 12.4 36.11 166
Charleston 31 55.2 35 71 6.5 40.75 148
Chicago 110 50.6 3344 3369 10.4 34.44 122
Cincinnati 23 54.0 462 453 7.1 39.04 132
Cleveland 65 49.7 1007 751 10.9 34.99 155
Columbus 26 51.5 266 540 8.6 37.01 134
Dallas 9 66.2 641 844 10.9 35.94 78
Denver 17 51.9 454 515 9.0 12.95 86
Des Moines 17 49.0 104 201 11.2 30.85 103
Detroit 35 49.9 1064 1513 10.1 30.96 129
Hartford 56 49.1 412 158 9.0 43.37 127
Houston 10 68.9 721 1233 10.8 48.19 103
Indianapolis 28 52.3 361 746 9.7 38.74 121
Jacksonville 14 68.4 136 529 8.8 54.47 116
Kansas City 14 54.5 381 507 10.0 37.00 99
Little Rock 13 61.0 91 132 8.2 48.52 100
Louisville 30 55.6 291 593 8.3 43.11 123
Memphis 10 61.6 337 624 9.2 49.10 105
Miami 10 75.5 207 335 9.0 59.80 128
Milwaukee 16 45.7 569 717 11.8 29.07 123
Minneapolis 29 43.5 699 744 10.6 25.94 137
Nashville 18 59.4 275 448 7.9 46.00 119
New Orleans 9 68.3 204 361 8.4 56.77 113
Norfolk 31 59.3 96 308 10.6 44.68 116
Omaha 14 51.5 181 347 10.9 30.18 98
Philadelphia 69 54.6 1692 1950 9.6 39.93 115
Phoenix 10 70.3 213 582 6.0 7.05 36
Pittsburgh 61 50.4 347 520 9.4 36.22 147
Providence 94 50.0 343 179 10.6 42.75 125
Richmond 26 57.8 197 299 7.6 42.59 115
Salt Lake City 28 51.0 137 176 8.7 15.17 89
San Francisco 12 56.7 453 716 8.7 20.66 67
Seattle 29 51.1 379 531 9.4 38.79 164
St. Louis 56 55.9 775 622 9.5 35.89 105
Washington 29 57.3 434 757 9.3 38.89 111
Wichita 8 56.6 125 277 12.7 30.58 82
Wilmington 36 54.0 80 80 9.0 40.25 114# 마할라노비스 제곱거리 계산계산
S <- cov(USairpollution) # 표본공분산행렬
D2 <- mahalanobis(USairpollution, # 데이터 행렬NAcolMeans(USairpollution), # 평균벡터 S) # 표본공분산행렬
D2
Albany Albuquerque Atlanta Baltimore
4.081128 8.063093 1.448871 3.421276
Buffalo Charleston Chicago Cincinnati
12.880983 7.888784 26.891450 7.265585
Cleveland Columbus Dallas Denver
11.489013 3.145760 6.368290 5.400965
Des Moines Detroit Hartford Houston
4.310154 7.222633 7.760199 8.684843
Indianapolis Jacksonville Kansas City Little Rock
4.270100 5.286583 3.060650 6.207429
Louisville Memphis Miami Milwaukee
3.060861 3.573384 14.277037 5.593824
Minneapolis Nashville New Orleans Norfolk
4.945830 2.564959 5.394046 3.754730
Omaha Philadelphia Phoenix Pittsburgh
3.052545 6.722708 20.258912 7.955753
Providence Richmond Salt Lake City San Francisco
18.176040 2.492672 4.431380 4.539210
Seattle St. Louis Washington Wichita
6.535345 4.767204 1.538300 9.060638
Wilmington
2.156830 # 카이제곱그림
qqplot(qchisq(ppoints(41, # 데이터 개수
a = 1/2), # (1:n-a)/(n+1-2a)
df = 7), # 카이제곱분포 자유도는 변수 개수
D2)
abline(a = 0, b = 1, col = "grey") # a+bx 직선을 그리는 함수
Caution! 함수
mahalanobis(데이터 행렬, 평균벡터, 표본공분산행렬)를 통해
관측된 데이터의 마할라노비스 거리를 계산할 수 있다. 계산된 마할라노비스
거리가 카이제곱분포를 따르는 지 파악하기 위해 카이제곱그림을 그렸으며,
이때 함수 qqplot(x, y)을 이용하였다. 함수
qqplot(x, y)는 두 데이터 셋 “x”와 “y”가 같은 분포로부터
왔는지를 확인할 때 사용하는 함수로서 일반적으로 관측된 데이터가 이론적인
분포를 따르는지 파악할 때 쓰인다. 함수 qqplot()에 입력된
벡터 qchisq(ppoints(41, a = 1/2), df = 7)의 함수
ppoints(41, a = 1/2)는 \(\left(\frac{1:41-0.5}{41} \right)\)을
계산하며, 함수 qchisq(a, df)는 누적확률이 “a”에 해당하는
자유도가 “df”인 카이제곱분포로부터의 분위수를 구해준다. 여기서 총 7개의
변수가 있기 때문에 자유도는 7이 되었다.
Result! 큰 상위 3개 정도의 데이터가 다변량 정규분포로부터
벗어난 점으로 보인다.
'data.frame': 32 obs. of 11 variables:
$ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
$ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
$ disp: num 160 160 108 258 360 ...
$ hp : num 110 110 93 110 175 105 245 62 95 123 ...
$ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
$ wt : num 2.62 2.88 2.32 3.21 3.44 ...
$ qsec: num 16.5 17 18.6 19.4 17 ...
$ vs : num 0 0 1 1 0 1 0 1 1 1 ...
$ am : num 1 1 1 0 0 0 0 0 0 0 ...
$ gear: num 4 4 4 3 3 3 3 4 4 4 ...
$ carb: num 4 4 1 1 2 1 4 2 2 4 ...mtcars
mpg cyl disp hp drat wt qsec vs am gear
Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4
Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4
Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4
Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3
Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3
Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3
Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3
Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4
Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4
Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4
Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4
Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3
Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3
Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3
Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3
Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3
Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3
Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4
Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4
Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4
Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3
Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3
AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3
Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3
Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3
Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4
Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5
Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5
Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5
Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5
Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5
Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4
carb
Mazda RX4 4
Mazda RX4 Wag 4
Datsun 710 1
Hornet 4 Drive 1
Hornet Sportabout 2
Valiant 1
Duster 360 4
Merc 240D 2
Merc 230 2
Merc 280 4
Merc 280C 4
Merc 450SE 3
Merc 450SL 3
Merc 450SLC 3
Cadillac Fleetwood 4
Lincoln Continental 4
Chrysler Imperial 4
Fiat 128 1
Honda Civic 2
Toyota Corolla 1
Toyota Corona 1
Dodge Challenger 2
AMC Javelin 2
Camaro Z28 4
Pontiac Firebird 2
Fiat X1-9 1
Porsche 914-2 2
Lotus Europa 2
Ford Pantera L 4
Ferrari Dino 6
Maserati Bora 8
Volvo 142E 2mtcars.d <- mtcars[,c("mpg", "disp", "hp", "drat", "wt", "qsec")]
# 마할라노비스 제곱거리 계산계산
S <- cov(mtcars.d) # 표본공분산행렬
D2 <- mahalanobis(mtcars.d, # 데이터 행렬NAcolMeans(mtcars.d), # 평균벡터 S) # 표본공분산행렬
D2
Mazda RX4 Mazda RX4 Wag Datsun 710
4.078929 3.198245 2.050829
Hornet 4 Drive Hornet Sportabout Valiant
3.498982 4.481779 6.263927
Duster 360 Merc 240D Merc 230
4.073871 3.511888 14.717512
Merc 280 Merc 280C Merc 450SE
3.694698 4.253093 3.752894
Merc 450SL Merc 450SLC Cadillac Fleetwood
1.816520 2.007100 5.887196
Lincoln Continental Chrysler Imperial Fiat 128
6.223139 10.653411 7.244391
Honda Civic Toyota Corolla Toyota Corona
8.984821 8.069408 5.234610
Dodge Challenger AMC Javelin Camaro Z28
5.782990 3.878396 3.904593
Pontiac Firebird Fiat X1-9 Porsche 914-2
6.255946 1.952287 4.889832
Lotus Europa Ford Pantera L Ferrari Dino
8.017199 11.448015 7.294298
Maserati Bora Volvo 142E
16.345103 2.534100 # 카이제곱그림
qqplot(qchisq(ppoints(32, # 데이터 개수
a = 1/2), # (1:n-a)/(n+1-2a)
df = 6), # 카이제곱분포 자유도는 변수 개수
D2)
abline(a = 0, b = 1, col = "grey") # a+bx 직선을 그리는 함수
Result! 전반적으로 점들이 직선에 가깝기 때문에 다변량
정규분포를 잘 따르는 것으로 보인다.
'data.frame': 150 obs. of 5 variables:
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...# 마할라노비스 제곱거리 계산계산
S <- cov(iris[,1:4]) # 표본공분산행렬
D2 <- mahalanobis(iris[,1:4], # 데이터 행렬NAcolMeans(iris[,1:4]), # 평균벡터 S) # 표본공분산행렬
D2
[1] 2.1344679 2.8491187 2.0813387 2.4523816 2.4621545
[6] 3.8834177 2.8621081 1.8333003 3.3840731 2.3752179
[11] 3.2831069 2.7747975 2.6132975 3.6034324 8.7375184
[16] 9.7127899 5.7605877 2.3213894 4.4996899 3.4388658
[21] 2.6360071 2.9292496 3.6134114 2.2371731 5.3023607
[26] 2.4453103 1.7658286 2.1971806 2.5027712 2.4643980
[31] 1.9849638 4.5911380 8.3583413 7.2213139 1.9820679
[36] 3.4173031 5.3372175 3.4513350 3.1549793 1.8926197
[41] 2.5485013 11.4240288 3.3144697 3.7085855 4.4840560
[46] 2.9786602 4.4333077 2.3594309 3.0076970 1.9285641
[51] 4.4528311 0.6273996 3.0186280 3.6125278 2.0404590
[56] 3.3195858 1.2763441 4.3093305 2.7426578 3.1665229
[61] 7.6832930 0.4323176 7.5614023 1.5482870 1.0372548
[66] 2.7865631 3.4732716 3.3289796 7.3923779 2.1837666
[71] 3.2615033 1.3601100 2.3885294 4.5073307 1.6095193
[76] 2.3362612 4.1982905 1.4219664 0.3194730 2.1112223
[81] 2.8147732 3.1338323 0.8579871 2.5178289 6.2164814
[86] 2.7910313 1.7971852 6.5544318 1.6677135 2.0056841
[91] 4.9630348 1.1821941 1.3713881 4.5698858 1.4932609
[96] 2.6869276 1.1870542 0.6575354 4.8207427 0.6341019
[101] 8.9395988 2.9682833 2.4451456 3.5551146 2.3541837
[106] 6.0617111 10.1378044 7.5880086 3.3301655 7.3453376
[111] 2.2611037 1.2342978 2.4316524 4.7588845 11.4105735
[116] 5.8815414 1.9743869 12.8130732 7.1768159 4.3974981
[121] 4.0399248 5.3595987 8.7973122 1.7967966 2.8868280
[126] 5.8835660 1.4473698 1.1348285 1.6688603 7.0710485
[131] 5.6935238 13.1010925 2.5441171 2.4087367 12.8803310
[136] 9.6569355 8.2028004 3.2575803 1.4687445 3.4924068
[141] 5.9713581 12.4413843 2.9682833 3.1069367 7.4592052
[146] 9.0639497 4.0366487 1.7670035 7.6824724 3.4787688# 카이제곱그림
qqplot(qchisq(ppoints(150, # 데이터 개수
a = 1/2), # (1:n-a)/(n+1-2a)
df = 4), # 카이제곱분포 자유도는 변수 개수
D2)
abline(a = 0, b = 1, col = "grey") # a+bx 직선을 그리는 함수
Result! 전반적으로 점들이 직선에 가깝기 때문에 다변량
정규분포를 잘 따르는 것으로 보인다.
1-3-2. Shapiro Wilk 정규성 검정
- 위의 카이제곱그림은 시각적으로 정규성 가정을 확인하기 때문에 정확한 판단이 어려울 때가 있다.
- 그래서, 가설 검정으로 정규성을 검토하기 위해 Shapiro-Wilk 정규성 검정을 수행할 수 있다.
- Shapiro-Wilk 검정은 주어진 데이터가 있을 때, 그 데이터들이
정규분포로부터 추출된 표본인지 검정하는 방법이다.
- Shapiro-Wilk 검정의 귀무가설(\(H_0\))은 “주어진 데이터의 모집단은 다변량 정규분포를 따른다.”이고 대립가설(\(H_1\))이 “주어진 데이터의 모집단은 다변량 정규분포를 따르지 않는다.”이다.
- 다변량 데이터에 대한 정규성 검정은 Package
mvnormtest에 내장되어 있는 함수mshapiro.test()를 통해 Shapiro-Wilk 검정을 수행할 수 있다.
pacman::p_load("mvnormtest")
package 'mvnormtest' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\User\AppData\Local\Temp\RtmpExJq88\downloaded_packages
Shapiro-Wilk normality test
data: Z
W = 0.59549, p-value = 2.025e-09Result! Shapiro-Wilk 검정통계량 값이 0.59549이며, \(p\)-값이 2.025e-09로 0에 가깝다. 이에
근거하여, 유의수준 5%에서 \(p\)-값이
\(\alpha = 0.05\)보다 작으므로
귀무가설을 기각한다. 즉, “USairpollution” 데이터는 다변량 정규분포를
따르지 않는다.
# 예제 22
data(mtcars)
mtcars.d <- mtcars[,c("mpg", "disp", "hp", "drat", "wt", "qsec")]
mshapiro.test(t(mtcars.d)) # 한 행이 하나의 변수값NA
Shapiro-Wilk normality test
data: Z
W = 0.80886, p-value = 6.008e-05Result! Shapiro-Wilk 검정통계량 값이 0.80886이며, \(p\)-값이 6.008e-05로 0에 가깝다. 이에
근거하여, 유의수준 5%에서 \(p\)-값이
\(\alpha = 0.05\)보다 작으므로
귀무가설을 기각한다. 즉, “mtcars.d” 데이터는 다변량 정규분포를 따르지
않는다.
Shapiro-Wilk normality test
data: Z
W = 0.97935, p-value = 0.02342Result! Shapiro-Wilk 검정통계량 값이 0.97935이며, \(p\)-값이 0.02342이다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha = 0.05\)보다 작으므로 귀무가설을
기각한다. 즉, “iris” 데이터는 다변량 정규분포를 따르지 않는다.
1-4. 정규화 변환: 박스-콕스 변환
- 정규성이 성립한다고 할 수 없을 때, 대체로 두 가지 방법이 있다.
- 정규성 검정 결과를 무시하고 무작정 정규성을 가정하는 것
- 표본의 크기가 충분히 클 때는 큰 문제가 되지 않을 수도 있으나 결론에 오류 가능성을 배제할 수 없다.
- 변수변환을 통해 정규성을 만족하는 변수로 만드는 것
- 대표적인 방법이 박스-콕스 변환(Box-Cox Transformation)이다.
- 정규성 검정 결과를 무시하고 무작정 정규성을 가정하는 것
일변량의 경우
- 박스-콕스 변환은 분산안정화 및 정규화를 위한 변환을 수행한다.
- 이 변환은 원자료(\(x\))에 대해
다음의 변환을 수행한다. \[
\begin{align*}
x^* = \begin{cases}
\frac{x^\lambda-1}{\lambda}, \;\; &\lambda \ne 0\\
\log(x), \;\;\; &\lambda = 0
\end{cases}
\end{align*}
\]
- 위의 변환에서 \(\lambda\)값은 원자료(\(x\))로부터 추정되며, 변환된 자료(\(x^*\))는 정규분포에 가까운 형태를 취하게 된다.
- R에서 박스-콕스 변환은 다음 함수를 이용할 수 있다.
- Package
MASS에 내장되어 있는 함수boxcox() - Package
car에 내장되어 있는 함수bcPower() - Package
car에 내장되어 있는 함수powerTransform()
- Package
# 박스-콕스 변환변환
p <- boxcox(lm(x~1),
plotit = FALSE,
lambda = seq(-10, 10, by = 0.0001)) # 후보 Lambda 값값
lambda <- p$x[which.max(p$y)] # 로그 우도함수가 최대가 되는 Lambda값a값
lambda
[1] 0.0296x_star <- if(lambda==0){
log(x)
}else{
(x^lambda - 1) / lambda
}
x_star
[1] 3.185921 3.185921 3.276024 3.206574 3.059200 3.023653 2.767803
[8] 3.350498 3.276024 3.087988 3.005448 2.916352 2.974435 2.833899
[15] 2.424878 2.424878 2.797663 3.663509 3.592950 3.713706 3.211679
[22] 2.855089 2.833899 2.689452 3.087988 3.474145 3.420376 3.592950
[29] 2.875885 3.116057 2.819546 3.206574
Result! 원자료 “x”에 대한 히스토그램과 정규확률그림을
보면, “x”가 정규분포를 따르고 있지 않다는 것을 알 수 있다. 왜냐하면
히스토그램은 왼쪽으로 치우쳐져 있으며, Q-Q plot은 직선에 가깝지 않기
때문이다. 하지만, 박스-콕스 변환을 적용한 “x_star”는 히스토그램도
좌우대칭 종모양이며, Q-Q plot도 직선에 가깝기 때문에 정규분포를 따르고
있음을 알 수 있다.
pacman::p_load("car")
# 원자료set.seed(100)
x <- rexp(1000) # 지수분포로부터 난수발생NApar(mfrow = c(1,2))
hist(x) # 히스토그램
qqnorm(x) # 정규확률그림림
# 박스-콕스 변환변환
p <- powerTransform(x) # 람다 추정p
Estimated transformation parameter
x
0.2568089 
Result! 원자료 “x”에 대한 히스토그램과 정규확률그림을
보면, “x”가 정규분포를 따르고 있지 않다는 것을 알 수 있다. 왜냐하면
히스토그램은 왼쪽으로 치우쳐져 있으며, Q-Q plot은 직선에 가깝지 않기
때문이다. 하지만, 박스-콕스 변환을 적용한 “x_star”는 히스토그램도
좌우대칭 종모양이며, Q-Q plot도 직선에 가깝기 때문에 정규분포를 따르고
있음을 알 수 있다.
다변량의 경우
- 몇몇 다변량 데이터 분석의 이론은 모집단에 대해 다변량 정규분포의 가정을 요구한다.
- 그러나, 실제로 데이터 분석을 수행할 때는 이 가정에 대한 검토를
생략하는 경우가 매우 빈번하다.
- 그 이유는 다변량 데이터에 대한 정규화 변환의 절차가 다소 복잡한 면이 있으며, 표본의 크기가 충분히 큰 경우에는 비정규성이 완화되는 측면이 있기 때문이다.
- 일변량의 경우에서와 마찬가지로 박스-콕스 변환은 다변량의 경우에도 유사하게 적용된다.
- 다만, 각 변수별로 변환의 차수를
동시에 결정하는 과정에 반복식을 통한 계산 과정이 필요하다.- 이 과정을 간단히 소개하면 다음과 같다.
- 변환된 변수들이 다변량 정규분포를 따른다고 가정하면, 원 변수들의 분포는 변환식에 의해 구할 수 있다.
- 이 분포에서 변환모수 \(\boldsymbol{\lambda}=(\lambda_1, \ldots, \lambda_p)^T\)에 대한 프로파일 로그-가능도함수를 정의하고, 이를 최대로 하는 \(\boldsymbol{\lambda}\)를 추정한다.
- 이 과정을 간단히 소개하면 다음과 같다.
2. 모집단 평균벡터에 대한 추론
- 다변량 정규모집단으로부터 얻은 확률표본에 대해 모집단의 평균벡터에 대한 추론을 알아보고자 한다.
- 모평균벡터의 검정에는 일변량 \(t\)-검정의 확장된 개념으로 Hotelling의 \(T^2\) 검정을 사용한다.
- 일변량 \(t\)-검정과 달리 Hotelling의 \(T^2\) 검정에서는 변수들 간의 공분산을 이용하여 검정을 수행한다.
- Hotelling의 \(T^2\) 검정은 Package
ICSNP에 내장되어 있는 함수HotellingsT2()를 통해 수행할 수 있다.- 함수
HotellingsT2()는 두 집단에 대한 검정일 경우, 공분산행렬은 같다고 가정한다.
- 함수
2-1. Hotelling의 \(T^2\) 검정
- 다변량 정규모집단의 가정하에서 모평균벡터에 대한 검정통계량으로
Hotelling의 \(T^2\) 통계량이 사용된다.
- 이 통계량은 일변량의 \(t\)-통계량을
다변량으로 확장한 것이다.
- 일변량의 \(t\)-통계량은 소표본(관측 개수 \(n<30\))이고, 모분산을 모를 때 모평균 추론에서 사용되는 통계량이다.
- 이 통계량은 일변량의 \(t\)-통계량을
다변량으로 확장한 것이다.
2-1-1. 일표본 \(t\)-검정과 다변량인 경우로의 확장
- 일변량 확률표본 \(X_1, \ldots,
X_n\)이 서로 독립이고 \(N(\mu_0,
\sigma^2)\)를 따를 때, \(t\)-통계량은 \(t=\frac{\bar{X}-\mu_0}{s/\sqrt{n}}\sim
t(n-1)\)이다.
- 여기서 \(s^2=\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})\)은 \(\sigma^2\)의 불편추정량이다.
- 일변량 확률표본에 대한 \(t\)-통계량의 제곱 형태 \(t^2=n(\bar{X}-\mu_0)(s^2)^{-1}(\bar{X}-\mu_0)\)는 다변량의 경우로 확장할 수 있다.
- 다변량 확률표본 \(\mathbf{X}_{1}, \ldots,
\mathbf{X}_n\)이 서로 독립이고 \(p\)-변량 정규분포 \(N_p(\boldsymbol{\mu}_0, \Sigma)\)를
따른다고 하자.
- 여기서 \(\mathbf{X}_i=(X_{i1}, \ldots, X_{ip})^T\)는 \(i\)번째 개체에 대한 확률벡터이며, 관측 개수는 \(n\)이고 변수 개수는 \(p\)개이다.
- 만약 모공분산행렬 \(\boldsymbol{\Sigma}\)를 모른다면, \(\boldsymbol{\Sigma}\)의 추정통계량으로 표본공분산행렬 \(\hat{\boldsymbol{\Sigma}}=\boldsymbol{S}=\frac{1}{n-1}\sum_{i=1}^n (\mathbf{X}_i-\bar{\mathbf{X}})(\mathbf{X}_i-\bar{\mathbf{X}})^T\)을 얻는다.
- 모평균벡터에 대한 추론을 위한 통계량은 Hotelling의 \(T^2\) 통계량으로 \(T^2=n(\bar{\mathbf{X}}-\boldsymbol{\mu}_0)^T\boldsymbol{S}^{-1}(\bar{\mathbf{X}}-\boldsymbol{\mu}_0)\)이다.
- Hotelling의 \(T^2\) 통계량은 \(T^2\sim \frac{(n-1)p}{n-p}F_{p,n-p}\), 즉, 자유도가 \(p,\;n-p\)인 \(F\)분포를 따른다.
2-1-2. \(\boldsymbol{\Sigma}\)를 모를 때 Hotelling \(T^2\) 검정
- 일변량의 경우, 한 개의 정규분포 \(N(\mu_0,
\sigma^2)\)를 따르는 모집단으로부터 얻은 확률표본 \(X_1, \ldots, X_n\)에 대해 특정한
모평균(\(H_0: \mu = \mu_0\))에 대한
가설검정을 하기 위한 검정통계량 \(t=\frac{\bar{X}-\mu_0}{s/\sqrt{n}}\)는
\(H_0\)하에서 \(t_{n-1}\)분포를 따른다.
- 이때 \(n<30\)으로 소표본이고 모분산 \(\sigma^2\)을 모른다고 가정한다.
- 다변량의 경우, 서로 독립인 \(n\)개의 \(p\times
1\) 확률벡터 \(\mathbf{X}_1, \ldots,
\mathbf{X}_n\)를 모평균벡터가 \(\boldsymbol{\mu}\)이고 공분산행렬이 \(\boldsymbol{\Sigma}\)인 \(p\)-변량 정규분포로부터 얻었을 때, \(\boldsymbol{\mu}=\boldsymbol{\mu}_0\)에
대한 검정법은 다음과 같다.
- 통계적 가설 : \(H_0 : \boldsymbol{\mu}=\boldsymbol{\mu}_0\) vs \(H_1 :\boldsymbol{\mu}\ne\boldsymbol{\mu}_0\).
- 검정법 : 유의수준 \(\alpha\)에서
\(T^2=n(\bar{\mathbf{X}}-\boldsymbol{\mu}_0)^T\boldsymbol{S}^{-1}(\bar{\mathbf{X}}-\boldsymbol{\mu}_0)
\ge \frac{(n-1)p}{n-p}F_{p,n-p}(\alpha)\)이면 \(H_0\)를 기각한다.
- \(F_{p,n-p}(\alpha)\)는 자유도가 \(p,\; n-p\)인 \(F\)분포를 따르는 확률변수 \(F\)에 대해 \(P(F\ge F_{p, n-p}(\alpha))=\alpha\)를 만족하는 값이다.
pacman::p_load("ICSNP")
There is a binary version available but the source version
is later:
binary source needs_compilation
survey 4.0 4.1-1 FALSE
package 'mitools' successfully unpacked and MD5 sums checked
package 'ICS' successfully unpacked and MD5 sums checked
package 'ICSNP' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\User\AppData\Local\Temp\RtmpExJq88\downloaded_packagesdata(pulmonary)
pulmonary
FVC FEV CC
1 -0.11 -0.12 -4.3
2 0.02 0.08 4.4
3 -0.02 0.03 7.5
4 0.07 0.19 -0.3
5 -0.16 -0.36 -5.8
6 -0.42 -0.49 14.5
7 -0.32 -0.48 -1.9
8 -0.35 -0.30 17.3
9 -0.10 -0.04 2.5
10 0.01 -0.02 -5.6
11 -0.10 -0.17 2.2
12 -0.26 -0.30 5.5# 일표본 Hotelling T^2 test
HotellingsT2(X = pulmonary, # p개의 변수로 이루어진 데이터 행렬NA= c(0, 0, 0), # True mean vector
test = "f") # "f" : F-distribution , "chi" : Chi-sqaured approximation
Hotelling's one sample T2-test
data: pulmonary
T.2 = 3.8231, df1 = 3, df2 = 9, p-value = 0.05123
alternative hypothesis: true location is not equal to c(0,0,0)Result! \(H_0 :
\boldsymbol{\mu}=\begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\) vs
\(H_1 :\boldsymbol{\mu}\ne\begin{pmatrix} 0 \\
0 \\ 0 \end{pmatrix}\)일 때, \(\frac{n-p}{(n-1)p}T^2\) 값은 3.8231이고
\(p\)-값은 0.05123이다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 크기 때문에 귀무가설을
기각하지 못한다. 즉, 평균벡터가 \(\begin{pmatrix} 0 \\ 0 \\ 0
\end{pmatrix}\)이라고 할 수 있다.
# 일표본 Hotelling T^2 test
HotellingsT2(X = pulmonary, # p개의 변수로 이루어진 데이터 행렬NA= c(0, 0, 2), # True mean vector
test = "f") # "f" : F-distribution , "chi" : Chi-sqaured approximation
Hotelling's one sample T2-test
data: pulmonary
T.2 = 6.6204, df1 = 3, df2 = 9, p-value = 0.01178
alternative hypothesis: true location is not equal to c(0,0,2)Result! \(H_0 :
\boldsymbol{\mu}=\begin{pmatrix} 0 \\ 0 \\ 2 \end{pmatrix}\) vs
\(H_1 :\boldsymbol{\mu}\ne\begin{pmatrix} 0 \\
0 \\ 2 \end{pmatrix}\)일 때, \(\frac{n-p}{(n-1)p}T^2\) 값은 6.6204이고
\(p\)-값은 0.01178이다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 평균벡터가 \(\begin{pmatrix} 0
\\ 0 \\ 2 \end{pmatrix}\)라고 할 수 없다.
2-1-3. 두 모집단에 대한 Hotelling \(T^2\) 검정
\(\boldsymbol{\Sigma}_1=\boldsymbol{\Sigma}_2=\boldsymbol{\Sigma}\)인 경우
- 두 개 모집단으로부터 서로 독립이고 \(p\)-변량 정규분포를 따르는 \(n_1\)개의 \(p\times 1\) 확률벡터로 구성된 \(\mathbf{X}_{11}, \ldots,
\mathbf{X}_{1n_1}\) 확률표본과 \(n_2\)개의 \(p\times 1\) 확률벡터로 구성된 \(\mathbf{X}_{21}, \ldots,
\mathbf{X}_{2n_2}\) 확률표본을 얻었다고 가정하자.
- 즉, \(\mathbf{X}_{1j} \overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)\), \(\mathbf{X}_{2j} \overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)\)일 때 두 집단의 모평균벡터가 같은지 검정하고자 한다.
- 이때 공분산행렬은 \(\boldsymbol{\Sigma}_1=\boldsymbol{\Sigma}_2=\boldsymbol{\Sigma}\)이며, \(\boldsymbol{\Sigma}\)는 알려지지 않았다고 가정한다.
- 그러면, \(\boldsymbol{\mu}_1=\boldsymbol{\mu}_2\)에
대한 검정법은 다음과 같다.
- 통계적 가설 : \(H_0 : \boldsymbol{\mu}_1=\boldsymbol{\mu}_2\) vs \(H_1 :\boldsymbol{\mu}_1\ne\boldsymbol{\mu}_2\).
- 검정법 : 유의수준 \(\alpha\)에서
다음을 만족하면 \(H_0\)를 기각한다.
\[
\begin{align*}
T^2&=(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)^T\left[\left(\frac{1}{n_1}
+
\frac{1}{n_2} \right)\boldsymbol{S}_{pl}\right]^{-1}(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)\\
&=
\frac{n_1n_2}{n_1+n_2}(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)^T\boldsymbol{S}_{pl}^{-1}(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)\ge\frac{(n_1+n_2-2)p}{n_1+n_2-p-1}F_{p,n_1+n_2-p-1}(\alpha)
\end{align*}
\]
- 여기서 \(\boldsymbol{S}_{pl}=\frac{(n_1-1)\boldsymbol{S}_1+(n_2-1)\boldsymbol{S}_2}{n_1+n_2-2}\)은 \(p\times p\) 합동공분산행렬(Pooled Covariance Matrix)이고, \(\boldsymbol{S}_1=\hat{\boldsymbol{\Sigma}}_1\), \(\boldsymbol{S}_2=\hat{\boldsymbol{\Sigma}}_2\)이다.
# 두 집단에 대한 Hotelling T^2 testset.seed(100)
x1 <- matrix(rnorm(50, 0, 1), nrow = 10)
x2 <- matrix(rnorm(100, 2, 1), nrow = 20)
x1
[,1] [,2] [,3] [,4] [,5]
[1,] -0.50219235 0.08988614 -0.4380900 -0.09111356 -0.10162924
[2,] 0.13153117 0.09627446 0.7640606 1.75737562 1.40320349
[3,] -0.07891709 -0.20163395 0.2619613 -0.13792961 -1.77677563
[4,] 0.88678481 0.73984050 0.7734046 -0.11119350 0.62286739
[5,] 0.11697127 0.12337950 -0.8143791 -0.69001432 -0.52228335
[6,] 0.31863009 -0.02931671 -0.4384506 -0.22179423 1.32223096
[7,] -0.58179068 -0.38885425 -0.7202216 0.18290768 -0.36344033
[8,] 0.71453271 0.51085626 0.2309445 0.41732329 1.31906574
[9,] -0.82525943 -0.91381419 -1.1577295 1.06540233 0.04377907
[10,] -0.35986213 2.31029682 0.2470760 0.97020202 -1.87865588x2
[,1] [,2] [,3] [,4] [,5]
[1,] 1.5529378 2.44890327 1.4428777 2.1725935 1.5822056
[2,] 0.2614021 0.93564433 3.4283014 2.2546013 1.1496192
[3,] 2.1788648 0.83758068 1.1070426 1.3854662 2.6890462
[4,] 3.8974657 3.64852175 0.8424288 0.5707849 1.5398038
[5,] -0.2719255 -0.06209602 1.4697035 1.6690246 3.3481844
[6,] 2.9804641 2.01274972 4.4456828 2.1283861 2.4430714
[7,] 0.6011744 0.91247165 1.1675042 3.0181200 1.8490738
[8,] 3.8248724 2.27053949 2.4135198 1.7444263 2.4555489
[9,] 3.3812987 3.00845187 0.8213169 1.6974590 1.9598453
[10,] 1.1611481 -0.07440475 0.8259652 3.6151907 2.4561210
[11,] 1.7380042 2.89682227 1.6670766 1.2262866 1.5915750
[12,] 1.9311560 1.95000423 3.3631137 2.4240024 -0.1364939
[13,] 1.6211164 0.65465069 1.5308527 1.4160530 2.1568219
[14,] 4.5819589 0.06878847 2.8428756 2.4150357 2.6600489
[15,] 2.1298341 2.70958158 0.5420063 0.4547383 1.0181656
[16,] 1.2869750 1.84209497 1.5996941 1.4812505 0.8863563
[17,] 2.6379942 2.21636787 1.2235827 1.7202084 1.5626523
[18,] 2.2016916 2.81736208 1.6307035 3.0074574 1.4838888
[19,] 1.9300831 3.72717575 3.2401015 1.5304300 2.4189960
[20,] 1.9075101 1.89622971 1.8925662 2.2978970 2.1341554dt <- rbind(x1, x2)
g <- factor(rep(c(1,2), c(10, 20))) # 집단을 식별하는 변수
HotellingsT2(dt ~ g, # X ~ g -> X : p개의 변수로 이루어진 데이터 행렬, g : 데이터 행렬의 관측값에 대응되는 두 개의 그룹을 나타내는 범주형 벡터NA= rep(0, 5))
Hotelling's two sample T2-test
data: dt by g
T.2 = 21.076, df1 = 5, df2 = 24, p-value = 4.542e-08
alternative hypothesis: true location difference is not equal to c(0,0,0,0,0)Result! \(H_0 :
\boldsymbol{\mu}_1-\boldsymbol{\mu}_2=\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0
\\ 0 \end{pmatrix}\) vs \(H_1
:\boldsymbol{\mu}_1-\boldsymbol{\mu}_2 \ne\begin{pmatrix} 0 \\ 0 \\ 0 \\
0 \\ 0 \end{pmatrix}\)일 때, \(\frac{n_1+n_2-p-1}{(n_1+n_2-2)p}T^2\) 값은
21.076이고 \(p\)-값은 4.542e-08이다.
이에 근거하여, 유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 두 평균벡터의 차이가 \(\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0
\end{pmatrix}\)라고 할 수 없다.
pacman::p_load("mvtnorm")
set.seed(100)
x1 <- rmvnorm(20, mean = c(0, 0, 1, 1), sigma = diag(1:4))
x2 <- rmvnorm(30, mean = c(2, 2, 2, 2), sigma = diag(1:4))
x1
[,1] [,2] [,3] [,4]
[1,] -0.50219235 0.18601316 0.863311591 2.7735696
[2,] 0.11697127 0.45061099 -0.007691025 2.4290654
[3,] -0.82525943 -0.50892191 1.155687368 1.1925489
[4,] -0.20163395 1.04629247 1.213699564 0.9413666
[5,] -0.38885425 0.72245985 -0.582772598 5.6205936
[6,] -0.43808998 1.08054489 1.453730266 2.5468092
[7,] -0.81437912 -0.62006274 -0.247460318 1.4618891
[8,] -1.15772946 0.34941822 0.842186682 4.5147512
[9,] -0.13792961 -0.15725135 -0.195139861 0.5564115
[10,] 0.18290768 0.59018425 2.845330961 2.9404040
[11,] -0.10162924 1.98442940 -2.077465669 2.2457348
[12,] -0.52228335 1.86991695 0.370502888 3.6381315
[13,] 0.04377907 -2.65682063 0.225665587 -2.4771959
[14,] 0.17886485 2.68342173 -2.935090373 2.9609283
[15,] -1.39882562 2.58075933 3.392479581 -0.6777038
[16,] -0.26199577 -0.09736016 0.343754430 6.1639179
[17,] 0.12983414 -1.00836960 2.105038444 1.4033832
[18,] -0.06991695 -0.13080044 1.777523277 -1.1287113
[19,] -1.16241932 2.33136181 -2.571655076 1.0254994
[20,] -1.08752835 0.38260062 2.746689881 -3.1488095x2
[,1] [,2] [,3] [,4]
[1,] 2.8968223 1.92929531 -0.33021336 -1.8624231
[2,] 2.7095816 1.77668856 2.37476015 3.6347242
[3,] 3.7271758 1.85324665 1.03503589 4.8566029
[4,] 1.1070426 0.36294705 1.08149960 6.8913655
[5,] 1.1675042 2.58480538 -0.04153909 -0.3480695
[6,] 1.6670766 3.92773389 1.18741297 3.6857513
[7,] 0.5420063 1.43388194 0.65520581 1.2614070
[8,] 3.2401015 1.84806565 2.29894072 2.5092025
[9,] 1.3854662 -0.02121537 1.42673373 2.2567721
[10,] 3.0181200 1.63856422 1.47598360 5.2303814
[11,] 1.2262866 2.59962995 0.98857416 2.8300714
[12,] 0.4547383 1.26637741 1.51538681 4.0149148
[13,] 1.5304300 2.42129003 1.27635881 0.2992384
[14,] 2.6890462 1.34918430 4.33512384 2.8861428
[15,] 1.8490738 2.64424337 1.93045005 2.9122421
[16,] 1.5915750 -1.02145859 2.27162353 3.3200978
[17,] 1.0181656 0.42506997 1.24249160 0.9677775
[18,] 2.4189960 2.18972444 3.79212951 5.3070065
[19,] 1.9820532 1.96577134 2.43344035 1.3257509
[20,] 1.8866463 1.86015844 2.45741179 2.2779674
[21,] 1.7577305 2.08348298 1.69295612 3.5893605
[22,] 2.0067378 1.10934203 1.56267487 0.6191557
[23,] 2.2025421 3.19696411 3.09478439 2.4028270
[24,] 1.9089294 2.40939238 1.90528291 -2.0836997
[25,] 2.3583692 1.47306282 4.19677535 6.3372006
[26,] 0.7602772 2.83420765 2.21480771 0.9525844
[27,] 2.6202280 3.00157657 1.83857572 1.4096066
[28,] 0.9141848 1.11637807 1.59642083 1.4983663
[29,] 2.9538953 1.62385807 5.28271423 1.1400183
[30,] 3.5755470 2.22901944 0.11994042 3.1538746dt <- rbind(x1, x2)
g <- factor(rep(c(1,2), c(20, 30))) # 집단을 식별하는 변수
HotellingsT2(dt ~ g, # X ~ g -> X : p개의 변수로 이루어진 데이터 행렬, g : 데이터 행렬의 관측값에 대응되는 두 개의 그룹을 나타내는 범주형 벡터NA= c(-2, -2, -1, -1))
Hotelling's two sample T2-test
data: dt by g
T.2 = 2.1642, df1 = 4, df2 = 45, p-value = 0.08837
alternative hypothesis: true location difference is not equal to c(-2,-2,-1,-1)Result! \(H_0 :
\boldsymbol{\mu}_1-\boldsymbol{\mu}_2=\begin{pmatrix} -2 \\ -2 \\ -1 \\
-1 \end{pmatrix}\) vs \(H_1
:\boldsymbol{\mu}_1-\boldsymbol{\mu}_2 \ne\begin{pmatrix} -2 \\ -2 \\ -1
\\ -1 \end{pmatrix}\)일 때, \(\frac{n_1+n_2-p-1}{(n_1+n_2-2)p}T^2\) 값은
2.1642이고 \(p\)-값은 0.08837이다. 이에
근거하여, 유의수준 5%에서 \(p\)-값이
\(\alpha=0.05\)보다 크기 때문에
귀무가설을 기각하지 못한다. 즉, 두 평균벡터의 차이가 \(\begin{pmatrix} -2 \\ -2 \\ -1 \\ -1
\end{pmatrix}\)라고 할 수 있다.
\(\boldsymbol{\Sigma}_1\ne\boldsymbol{\Sigma}_2\)인 경우
- 두 개 모집단으로부터 서로 독립이고 \(p\)-변량 정규분포를 따르는 \(n_1\)개의 \(p\times 1\) 확률벡터로 구성된 \(\mathbf{X}_{11}, \ldots,
\mathbf{X}_{1n_1}\) 확률표본과 \(n_2\)개의 \(p\times 1\) 확률벡터로 구성된 \(\mathbf{X}_{21}, \ldots,
\mathbf{X}_{2n_2}\) 확률표본을 얻은 경우, 즉, \(\mathbf{X}_{1j}
\overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_1,
\boldsymbol{\Sigma}_1)\), \(\mathbf{X}_{2j}
\overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_2,
\boldsymbol{\Sigma}_2)\)를 따르는 경우, 두 집단의 모평균벡터가
같은지 검정하고자 한다.
- 여기서 공분산행렬은 \(\boldsymbol{\Sigma}_1\ne\boldsymbol{\Sigma}_2\)이며, \(\boldsymbol{\Sigma}_1\)과 \(\boldsymbol{\Sigma}_2\)는 알려지지 않았다고 가정한다.
- 그러면, \(\boldsymbol{\mu}_1=\boldsymbol{\mu}_2\)에
대한 검정법은 다음과 같다.
- 통계적 가설 : \(H_0 : \boldsymbol{\mu}_1=\boldsymbol{\mu}_2\) vs \(H_1 :\boldsymbol{\mu}_1\ne\boldsymbol{\mu}_2\).
- 검정법 : 유의수준 \(\alpha\)에서
근사적인 검정법으로 \(n_1-p,\;
n_2-p\)가 충분히 클 때, \(T^2=(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)^T\left(\frac{\boldsymbol{S}_1}{n_1}
+ \frac{\boldsymbol{S}_2}{n_2}
\right)^{-1}(\bar{\mathbf{X}}_1-\bar{\mathbf{X}}_2)\ge
\chi^2_p(\alpha)\)이면 \(H_0\)를
기각한다.
- 여기서 \(\boldsymbol{S}_1=\hat{\Sigma}_1\)이고 \(\boldsymbol{S}_2=\hat{\Sigma}_2\)이다.
2-1-4. 짝지어진 두 집단에 대한 검정
- 관측값 간에 짝지어진 관계로 인해 두 개의 표본이 서로 독립이 아닐 때
두 집단에 대한 \(T^2\) 검정은 적절하지 않다.- 예를 들어, 한 학생에게 특별한 교육방법 시행 전과 후 2회 3과목 시험이 치러진 경우와 한 환자에게 약을 투여한 뒤 1시간 후, 3시간 후에 효과와 혈압을 측정한 경우를 들 수 있다.
- 이러한 경우에는 짝지어진 관계로 인해 각
측정치 간의 차이를 구해한 개 모집단 문제로 전환하여 보는 것이 타당하다. - 각 개체에 대해 \(p\)-변량
정규분포를 따르는 \(p\times 1\)
확률벡터를 쌍(Pair)으로 측정하여 \((\mathbf{X}_1, \mathbf{Y}_1), \ldots,
(\mathbf{X}_n, \mathbf{Y}_n)\)로 확률표본을 얻은 경우 두 집단의
모평균벡터가 같은지 검정하고자 한다.
- 통계적 가설 : \(H_0 : \mathbf{d}=\boldsymbol{\mu}_1-\boldsymbol{\mu}_2=0\) vs \(H_1 :\mathbf{d}=\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\ne 0\).
- 검정법 : 유의수준 \(\alpha\)에서
\(T^2=
\bar{\mathbf{d}}^T\left(\frac{\boldsymbol{S}_d}{n}\right)^{-1}\bar{\mathbf{d}}=n\bar{\mathbf{d}}^T\boldsymbol{S}_d^{-1}\bar{\mathbf{d}}\ge\frac{(n-1)p}{n-p}F_{p,
n-p}(\alpha)\)이면 \(H_0\)를
기각한다.
- 여기서 \(\bar{\mathbf{d}}=\frac{1}{n}\sum_{i=1}^n \mathbf{d}_i\), \(\boldsymbol{S}_d=\frac{1}{n-1}\sum_{i=1}^n(\mathbf{d}_i-\bar{\mathbf{d}})(\mathbf{d}_i-\bar{\mathbf{d}})^T\)이다.
# 짝지어진 두 집단에 대한 검정NAx <- cbind(x1 = c(1, 5, 3, 6, 9, 3, 4, 3, 2), # 특별한 교육방법 시행 전 2과목 시험 등급NA= c(1, 3, 2, 5, 9, 6, 2, 1, 5))
y <- cbind(y1 = c(1, 2, 4, 4, 8, 1, 5, 2, 4), # 특별한 교육방법 시행 후 2과목 시험 등급NA= c(2, 4, 5, 3, 6, 3, 6, 3, 2))
d <- x-y
d
x1 x2
[1,] 0 -1
[2,] 3 -1
[3,] -1 -3
[4,] 2 2
[5,] 1 3
[6,] 2 3
[7,] -1 -4
[8,] 1 -2
[9,] -2 3# 일표본 Hotelling T^2 test
HotellingsT2(X = d, # p개의 변수로 이루어진 데이터 행렬NA= c(0, 0), # True mean vector
test = "f") # "f" : F-distribution , "chi" : Chi-sqaured approximation
Hotelling's one sample T2-test
data: d
T.2 = 0.46483, df1 = 2, df2 = 7, p-value = 0.6463
alternative hypothesis: true location is not equal to c(0,0)Result! \(H_0 :
\mathbf{d}=\boldsymbol{\mu}_1-\boldsymbol{\mu}_2=\begin{pmatrix} 0 \\ 0
\end{pmatrix}\) vs \(H_1
:\mathbf{d}=\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\ne \begin{pmatrix} 0
\\ 0 \end{pmatrix}\)일 때, \(\frac{n-p}{(n-1)p}T^2\) 값은 0.46483이고
\(p\)-값은 0.6463이다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 크기 때문에 귀무가설을
기각하지 못한다. 즉, 두 평균벡터의 차이가 \(\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)라고
할 수 있으며, 이는 특별한 교육방법 시행 전과 후에 시험 등급은 변화가
없다는 것을 의미한다.
3. 다변량 분산분석
- 다변량 분산분석(Multivariate Analysis of Variance, MANOVA)는
세 개 이상의 다변량 표본평균이 같은 지를 검정하는 방법이다.- 위에서 소개한 Hotelling \(T^2\) 검정은 한 집단 또는 두 집단에 대한 모평균벡터에 대한 추론을 위해 사용된다.
- 다변량 분산분석은 두 개 또는 그 이상의 변수가 있을 때 사용하는 방법으로, 통상적으로 개별 변수에 대한 유의성 검정을 수행한 후 적용한다.
- 다변량 분산분석은 일변량 분산분석(ANOVA)의 일반화된 형태이나, 일변량 분산분석(One-way ANOVA)과는 달리 평균의 차이에 대한 유의성 검정에서 변수 간의 공분산을 사용한다.
3-1. 일원배치 다변량 분산분석
- \(g(\ge 3)\)개 모집단에 대해 모평균벡터에 대한 비교를 위해 일원배치 다변량 분산분석(One-way MANOVA)를 수행하고자 한다.
공분산행렬이 같은다변량 정규분포를 따르는 \(g\)개 모집단에서 각각 \(n_1, \ldots, n_g\)개의 \(p \times 1\) 관측벡터가 얻어졌으며 이들은 서로 독립이라 가정한다. \[ \begin{align*} \text{모집단 1 : } &\mathbf{X}_{11}, \ldots, \mathbf{X}_{1n_1} \overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_1, \Sigma)\\ \text{모집단 2 : } &\mathbf{X}_{21}, \ldots, \mathbf{X}_{2n_2} \overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_2, \Sigma)\\ &\vdots\\ \text{모집단 g : } &\mathbf{X}_{g1}, \ldots, \mathbf{X}_{gn_g}\overset{\underset{\mathrm{iid}}{}}{\sim}N_p(\boldsymbol{\mu}_g, \Sigma)\\ \end{align*} \]- 데이터의 변동을 처리에 대한 변동과 오차에 의한 변동으로 분해하면 다음과 같다. \[ \begin{align*} X_{ij}-\bar{X} = (\bar{\mathbf{X}}_i -\bar{\mathbf{X}}) + (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i) \end{align*} \]
- 마찬가지로 총제곱합행렬을 처리제곱합행렬과 오차제곱합행렬로 분해하면
다음과 같다. \[
\begin{align*}
&\sum_{i=1}^{g}\sum_{j=1}^{n_i}
(\mathbf{X}_{ij}-\bar{\mathbf{X}})(\mathbf{X}_{ij}-\bar{\mathbf{X}})^T\\
&= \sum_{i=1}^g n_i(\bar{\mathbf{X}}_i
-\bar{\mathbf{X}})(\bar{\mathbf{X}}_i -\bar{\mathbf{X}})^T+
\sum_{i=1}^{g}\sum_{j=1}^{n_i} (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i)
(\mathbf{X}_{ij}-\bar{\mathbf{X}}_i)^T\\
& = \mathbf{B}+\mathbf{E}
\end{align*}
\]
- 여기서 처리간 제곱합(Between Sum of Squares)은 \(\mathbf{B}=\sum_{i=1}^g n_i(\bar{\mathbf{X}}_i -\bar{\mathbf{X}})(\bar{\mathbf{X}}_i -\bar{\mathbf{X}})^T\)이며 처리평균과 전체 평균간의 변동으로 처리제곱합행렬을 나타낸다.
- 또한, 처리내 제곱합(Within Sum of Squares)은 \(\mathbf{E}=\sum_{i=1}^{g}\sum_{j=1}^{n_i} (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i) (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i)^T\)이며 오차제곱합행렬을 나타낸다.
| 일원배치 다변량 분산분석표 | ||
|---|---|---|
| 요인 | 제곱합과 교차곱 행렬 | 자유도 |
| 처리 | \(\mathbf{B}=\sum_{i=1}^g n_i(\bar{\mathbf{X}}_i -\bar{\mathbf{X}})(\bar{\mathbf{X}}_i -\bar{\mathbf{X}})^T\) | \(g-1\) |
| 오차 | \(\mathbf{E}=\sum_{i=1}^{g}\sum_{j=1}^{n_i} (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i) (\mathbf{X}_{ij}-\bar{\mathbf{X}}_i)^T\) | \(\sum_{i=1}^g n_i -g\) |
| 총 | \(\mathbf{B}+\mathbf{E}\) | \(\sum_{i=1}^g n_i-1\) |
- 일원배치 다변량 분산분석에서 관심 있는 귀무가설(\(H_0\))은 “집단(처리)간의 모평균벡터에
차이가 없다. 즉, 집단(처리)에 따른 효과가 없다.”이고, 대립가설(\(H_1\))은 “적어도 하나 이상의 집단에서
모평균벡터에 차이가 있다. 즉, 집단에 따른 효과가 있다.”이다.
- 통계적 가설 : \(H_0 : \boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_g\) vs \(H_1 : \text{Not } H_0\)
- 검정법 : 유의수준 \(\alpha\)에서
Wilks Lambda 통계량 \(\Lambda\)를
이용하여 \(\Lambda =
\frac{|\mathbf{E}|}{|\mathbf{B}+\mathbf{E}|}\le \Lambda_{p, g-1,
N-g}(\alpha)\)이면 \(H_0\)를
기각한다.
- 여기서 \(N=\sum_{i=1}^g n_i\)이며, \(\Lambda_{p, g-1, N-g}\)는 Wilks Lambda의 분포이다.
- Wilks Lambda의 분포 : \(H_0\)하에서
\(N\)이 충분히 크면 근사적으로 \(\Lambda^*=-\frac{(N-1-p-g)}{2}\ln{\Lambda}\)는
자유도가 \(p(g-1)\)인 카이제곱분포를
따른다.
- 즉, 근사적인 검정법은 \(\Lambda^*\ge \chi^2_{p(g-1)}(\alpha)\)이면 \(H_0\)를 기각한다.
# 일원배치 다변량 분산분석 with Wilks statisticscs
data(iris)
fit <- manova(cbind(Sepal.Length, Sepal.Width) ~ Species, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= iris)
summary(fit, test = "Wilks")
Df Wilks approx F num Df den Df Pr(>F)
Species 2 0.16654 105.88 4 292 < 2.2e-16 ***
Residuals 147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “Species”의 자유도는 집단 내의 범주 개수(3)
- 1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(150) - 집단 내의
범주 개수(3)으로 147이 된다. 집단 간의 평균벡터에 차이가 있는 지
검정하기 위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Wilks Lambda
검정통계량 값은 0.16654이고 근사적인 \(F\)값은 105.88이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 붓꽃 종류(“Species”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
# 일원배치 다변량 분산분석 with Wilks statisticscs
pacman::p_load("car")
data(Baumann)
fit <- manova(cbind(pretest.1, post.test.1, pretest.2, post.test.2) ~ group, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= Baumann)
summary(fit, test = "Wilks")
Df Wilks approx F num Df den Df Pr(>F)
group 2 0.50701 6.066 8 120 1.518e-06 ***
Residuals 63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Caution! Package car에 내장되어 있는 데이터
Baumann은 Baumann와 Jones에 의해 실시된 실험연구로부터
수집된 데이터이다. 총 66개의 행과 6개의 변수로 이루어져 있다. 변수
“group”은 교수법의 종류, “pretest”은 교수법 실시전 시험 점수,
“post.test”는 교수법 실시후 시험점수를 나타낸다.
Result! 집단 “group”의 자유도는 집단 내의 범주 개수(3) -
1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(66) - 집단 내의 범주
개수(3)으로 63이 된다. 집단 간의 평균벡터에 차이가 있는 지 검정하기
위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Wilks Lambda
검정통계량 값은 0.50701이고 근사적인 \(F\)값은 6.066이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 교수법(“group”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
3-1-1. 그외의 통계량
- 일원배치 다변량 분산분석의 가설에 대해 널리 사용되는 검정통계량은
다음과 같다.
- Wilks Lambda 검정 통계량
- Pillai 검정통계량
- Lawley-Hotelling 검정통계량
- Roy의 검정통계량
- 위의 통계량들은 행렬 \(\mathbf{E}^{-1}\mathbf{B}\)의 고유값 \(\lambda_1>\ldots>\lambda_s\)을
이용하여 검정통계량을 계산한다.
- 여기서, \(s=\text{min}(g-1, p)\)이다.
Pillai 검정
- Pillai 검정통계량은 \(V^{(S)}=\sum_{i=1}^{s}\frac{\lambda_{i}}{1+\lambda_i}\)로 정의되며 \(V^{(S)}\)가 클수록 \(H_0\)를 기각한다.
- 유의수준 \(\alpha\)에서 검정법은
\(V^{(S)}\ge V_{m, N}(\alpha)\)이면
\(H_0\)를 기각한다.
- 여기서 \(m=\frac{1}{2}(|g-1-p|-1)\), \(N=\frac{1}{2}(\sum_{i=1}^g n_i-g-p-1)\)이다.
# 일원배치 다변량 분산분석 with Pillaiai
data(iris)
fit <- manova(cbind(Sepal.Length, Sepal.Width) ~ Species, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= iris)
summary(fit, test = "Pillai")
Df Pillai approx F num Df den Df Pr(>F)
Species 2 0.94531 65.878 4 294 < 2.2e-16 ***
Residuals 147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “Species”의 자유도는 집단 내의 범주 개수(3)
- 1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(150) - 집단 내의
범주 개수(3)으로 147이 된다. 집단 간의 평균벡터에 차이가 있는 지
검정하기 위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Pillai
검정통계량 값은 0.94531이고 근사적인 \(F\)값은 65.878이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 붓꽃 종류(“Species”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
# 일원배치 다변량 분산분석 with Pillaiai
pacman::p_load("car")
data(Baumann)
fit <- manova(cbind(pretest.1, post.test.1, pretest.2, post.test.2) ~ group, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= Baumann)
summary(fit, test = "Pillai")
Df Pillai approx F num Df den Df Pr(>F)
group 2 0.56225 5.9637 8 122 1.892e-06 ***
Residuals 63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “group”의 자유도는 집단 내의 범주 개수(3) -
1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(66) - 집단 내의 범주
개수(3)으로 63이 된다. 집단 간의 평균벡터에 차이가 있는 지 검정하기
위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Pillai
검정통계량 값은 0.56225이고 근사적인 \(F\)값은 5.9637이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 교수법(“group”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
Lawley-Hotelling 검정
- Lawley-Hotelling 검정통계량은 \(U^{(S)}=\sum_{i=1}^s \lambda_i\)로 정의되며 \(U^{(S)}\)가 클수록 \(H_0\)를 기각한다.
- 유의수준 \(\alpha\)에서 검정법은 \(U^{(S)}\ge U_{g-1, \sum_{i=1}^g n_i-g}(\alpha)\)이면 \(H_0\)를 기각한다.
# 일원배치 다변량 분산분석 with Lawley-Hotellingng
data(iris)
fit <- manova(cbind(Sepal.Length, Sepal.Width) ~ Species, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= iris)
summary(fit, test = "Hotelling-Lawley")
Df Hotelling-Lawley approx F num Df den Df Pr(>F)
Species 2 4.3328 157.06 4 290 < 2.2e-16 ***
Residuals 147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “Species”의 자유도는 집단 내의 범주 개수(3)
- 1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(150) - 집단 내의
범주 개수(3)으로 147이 된다. 집단 간의 평균벡터에 차이가 있는 지
검정하기 위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때,
Lawley-Hotelling 검정통계량 값은 4.3328이고 근사적인 \(F\)값은 157.06이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 붓꽃 종류(“Species”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
# 일원배치 다변량 분산분석 with Lawley-Hotellingng
pacman::p_load("car")
data(Baumann)
fit <- manova(cbind(pretest.1, post.test.1, pretest.2, post.test.2) ~ group, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= Baumann)
summary(fit, test = "Hotelling-Lawley")
Df Hotelling-Lawley approx F num Df den Df Pr(>F)
group 2 0.83573 6.1635 8 118 1.237e-06 ***
Residuals 63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “group”의 자유도는 집단 내의 범주 개수(3) -
1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(66) - 집단 내의 범주
개수(3)으로 63이 된다. 집단 간의 평균벡터에 차이가 있는 지 검정하기
위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때,
Lawley-Hotelling 검정통계량 값은 0.83573이고 근사적인 \(F\)값은 6.1635이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 교수법(“group”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
Roy 최대근 검정
- Roy 검정통계량은 \(\theta= \frac{\lambda_1}{1+\lambda_1}\)로 정의되며 \(\theta\)가 클수록 \(H_0\)를 기각한다.
- 유의수준 \(\alpha\)에서 검정법은
\(\theta\ge R_{m, N}(\alpha)\)이면
\(H_0\)를 기각한다.
- 여기서 \(m=\frac{1}{2}(|g-1-p|-1)\), \(N=\frac{1}{2}(\sum_{i=1}^g n_i-g-p-1)\)이다.
# 일원배치 다변량 분산분석 with Royoy
data(iris)
fit <- manova(cbind(Sepal.Length, Sepal.Width) ~ Species, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= iris)
summary(fit, test = "Roy")
Df Roy approx F num Df den Df Pr(>F)
Species 2 4.1718 306.63 2 147 < 2.2e-16 ***
Residuals 147
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “Species”의 자유도는 집단 내의 범주 개수(3)
- 1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(150) - 집단 내의
범주 개수(3)으로 147이 된다. 집단 간의 평균벡터에 차이가 있는 지
검정하기 위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Roy
검정통계량 값은 4.1718이고 근사적인 \(F\)값은 306.63이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 붓꽃 종류(“Species”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
# 일원배치 다변량 분산분석 with Royoy
pacman::p_load("car")
data(Baumann)
fit <- manova(cbind(pretest.1, post.test.1, pretest.2, post.test.2) ~ group, # p개의 변수로 이루어진 데이터 행렬 ~ 각 관측에 해당되는 그룹변수NA= Baumann)
summary(fit, test = "Roy")
Df Roy approx F num Df den Df Pr(>F)
group 2 0.6128 9.3453 4 61 5.872e-06 ***
Residuals 63
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Result! 집단 “group”의 자유도는 집단 내의 범주 개수(3) -
1로 2가 되고, 잔차에 대한 자유도는 전체 관측 개수(66) - 집단 내의 범주
개수(3)으로 63이 된다. 집단 간의 평균벡터에 차이가 있는 지 검정하기
위해, 가설이 \(H_0 :
\boldsymbol{\mu}_1=\ldots=\boldsymbol{\mu}_3\) vs \(H_1 : \text{Not } H_0\)일 때, Roy
검정통계량 값은 0.6128이고 근사적인 \(F\)값은 9.3453이며, \(p\)-값은 0에 가깝다. 이에 근거하여,
유의수준 5%에서 \(p\)-값이 \(\alpha=0.05\)보다 작기 때문에 귀무가설을
기각한다. 즉, 교수법(“group”)의 적어도 한 그룹 이상은 평균벡터가
다르다.
3-2. 공분산의 동질성 검정
- 다변량 분산분석(MANOVA)과 선형판별분석 등에서는 분산-공분산행렬의 동질성을 가정한다.
- 이 가정에 대해서는 일반적으로 Box의 M-검정이 사용된다.
- 그러나, 이 검정은 정규성이 위배될 경우 매우 민감하여, 많은 경우 기각을 하게 된다.
- Box의 M-검정은 Package
biotools에 내장되어 있는 함수boxM()을 통해 수행할 수 있다.- 분산-공분산행렬의 동질성 검정에서 귀무가설(\(H_0\))은 “모든 집단의 공분산행렬은 동일하다.”이고 대립가설(\(H_1\))이 “적어도 하나 이상의 집단의 분산-공분산행렬에 차이가 있다.”이다.
# 분산-공분산행렬 검정pacman::p_load("biotools")
There is a binary version available but the source version
is later:
binary source needs_compilation
biotools 4.1 4.2 FALSEboxM(data = iris[, 1:4], # X_1, ..., X_p p개의 변수에 대한 값으로 이루어진 데이터 행렬 NA= iris$Species) # 범주형 그룹의 벡터NA
Box's M-test for Homogeneity of Covariance Matrices
data: iris[, 1:4]
Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16Result! “iris” 데이터의 붓꽃 종류(“Species”) 간에
분산-공분산행렬 동질성 검정 결과, 근사적인 카이제곱 검정통계량 값은
140.94이고 \(p\)-값은 0에 가깝다. 이에
근거하여, 유의수준 5%에서 \(p\)-값이
\(\alpha = 0.05\)보다 작으므로
귀무가설을 기각한다. 즉, 붓꽃 종류(“Species”)간에 분산-공분산행렬은
적어도 하나 이상의 집단에서 차이가 있다.
# 분산-공분산행렬 검정pacman::p_load("car")
data(Baumann)
boxM(data = Baumann[, 2:5], # X_1, ..., X_p p개의 변수에 대한 값으로 이루어진 데이터 행렬 NA= Baumann$group) # 범주형 그룹의 벡터NA
Box's M-test for Homogeneity of Covariance Matrices
data: Baumann[, 2:5]
Chi-Sq (approx.) = 20.585, df = 20, p-value = 0.4219Result! “Baumann” 데이터의 교수법(“group”) 간에
분산-공분산행렬 동질성 검정 결과, 근사적인 카이제곱 검정통계량 값은
20.585이고 \(p\)-값은 0.4219이다. 이에
근거하여, 유의수준 5%에서 \(p\)-값이
\(\alpha = 0.05\)보다 크기 때문에
귀무가설을 기각하지 못한다. 즉, 교수법(“group”)간에 분산-공분산행렬은
모두 같다.