12.2 4가지 유형의 제곱합 (Four Types of Sum of Square)
SAS회사에서 ANOVA시에 계산하는 Sum of Square(SS)를 4가지 유형으로 구분하여 제시하였다. 대부분의 교과서에는 SAS의 방식으로 설명하고 있지 않지만, 실무에서 피해갈 수 없는 문제이므로 여기에서 설명한다. 그런데, SAS가 계산하는 방식으로 SS를 구하려면 먼저 g2 type inverse matrix를 이해해야 하므로 이를 먼저 설명한다.
12.2.1 g2 type inverse matrix
g2 inverse란 행렬 generalized inverse의 type 2를 의미한다. 보통 역행렬은 정방행렬에 대해서 정의되어 있고, determinant가 0이면 역행렬이 존재하지 않는다.
역행렬이 존재하지 않는 행렬을 singular matrix (비정칙 행렬)라 하고, 존재하면 non-singular matrix (정칙 행렬)라고 한다.
A라는 정방행렬이 있을 때 A의 역행렬은 다음의 조건을 만족하는 행렬로서, A-1라 표기하고 존재한다면 이것이 유일하다. 이런 행렬 A를 non-singular matrix라고 한다.
AA-1 = A-1A = I
여기서 I는 A와 dimension이 같은 identity matrix (좌상 우하 대각선의 원소만 1이고 나머지는 0인 행렬)를 의미한다. I의 dimension이 \(n \times n\)일 때 \(I_n\)으로 표기하기도 한다.
그런데, A가 정방행렬이 아니거나, 정방행렬이어도 역행렬이 존재하지 않는 경우(singular matrix)가 종종 있다. 이 경우에 좀 더 일반화된 역행렬을 생각할 수 있는데, 이것은 다음의 조건을 만족하는 행렬이다.
AA-A = A
이때 이런 역행렬을 보통의 역행렬과 구분하기 위해 위 첨자로 -1을 쓰지 않고, - 만 쓰고 이를 일반화 역행렬(generalized inverse)라 부른다.
다음과 같이 보통의 역행렬이 존재하지 않는 경우(singular)에도 일반화 역행렬은 존재한다.
A = matrix(1:9, nrow=3)solve(A) # not solved due to singularityG2SWEEP(A) # generalized inverse [,1] [,2] [,3]
[1,] -1.6666667 1.3333333 0
[2,] 0.6666667 -0.3333333 0
[3,] 0.0000000 0.0000000 0
attr(,"rank")
[1] 2일반화 역행렬은 정방행렬이 아니어도 구할 수 있지만, 단점은 이것이 unique하지 않다는 것이다. 처리하는 사람에 따라서 결과가 달라지고 (일관되지 않다면) 연구자들 간의 의사소통에 문제가 생긴다.
따라서, 누가 처리하든 동일한 결과를 얻기 위해 계산 방법을 지정해둔 것들이 있는데, 이것 중의 하나가 g2 inverse이다. R에서는 sasLM::G2SWEEP 함수로 g2 type의 inverse를 구할 수 있다. 알고리듬이 포함된 source는 다음과 같다.
G2SWEEP = function(A, Augmented=FALSE, eps=1e-8)
{
idx = abs(diag(A)) > eps
p = sum(idx)
if (p == 0) {A[,] = 0 ; return(A)}
B = A[idx, idx, drop=F]
p0 = ifelse(Augmented, p - 1, p)
r = 0
for (k in 1:p0) {
d = B[k, k]
if (abs(d) < eps) { B[k, ] = 0 ; B[, k] = 0 ; next }
B[k, ] = B[k, ]/d
r = r + 1
for (i in 1:p) {
if (i != k) {
c0 = B[i, k]
B[i, ] = B[i, ] - c0*B[k, ]
B[i, k] = -c0/d
}
}
B[k, k] = 1/d
}
A[!idx, !idx] = 0
A[idx, idx] = B
attr(A, "rank") = r
return(A)
}선형 회귀 분석에서도 model matrix (design matrix, 설계행렬)의 역행렬이 필요한데, 이것이 종종 singular하다. SAS, SPSS 진영에서는 이 경우에도 부족하지만 (unique하지는 않지만, 유일한 해는 아니지만) 어떤 일관된 해를 제시하고자 하여 g2 inverse를 사용한다.
하지만, 지금의 R 진영 (W. N. Venables 와 같은 분이 대표적인데, 사실 R이 만들어지기 전부터 있었던 분들이다.)에서는 이것의 사용을 반대했거나, 적어도 찬성하지는 않았다.
대가들 중에 Shayle R. Searle (1928-2013) 과 같은 분이 대표적인 찬성자이며, Norman R Draper같은 분은 온건한 반대론자이다.
Venables 는 ’Exegeses on Linear Models’이라는 글(1998년 10월 S-PLUS User’s Conference, 2000년 5월 13일자 파일)을 발표했는데, 지금도 인터넷에서 찾을 수 있다. 내용은 SAS 의 Type 3 SS에 대한 신랄한 비판이다.
Type 3 SS를 구하려면 g2 inverse가 핵심이고 필요하기 때문에, 이에 대한 반대는 g2 inverse의 사용에 대한 반대라고도 볼 수 있다.
Draper는 singular matrix를 만나게 되면, 이것을 non-singular하게 만들라고 권하고 있는데(Applied regression analysis 3e 1998. p444), 이것은 data와 model이 맞지 않는 것이니, data를 더 모으거나, model을 변경하라고 권하고 있다. (ibid., p126)
그래서, SAS, SPSS 진영에서는 Type 3 SS를 사용하고, S-PLUS, R 진영은 Type 3 SS를 사용하지 않게 (또는 권장하지 않게) 되었고, R에서는 SAS와 동일한 결과를 얻기 어렵게 되었다.
R에는 car라는 package의 Anova 함수에서 인자를 type=3이라고 주면 된다. 다만 lm을 쓰기 전에 contrast option을 다음과 같이 바꾸어 주어야 한다.
options(contrasts=c("contr.sum", "contr.poly"))하지만, 이렇게 해도 종종 SAS와 다른 결과를 얻게 되는데, 다음은 간단한 예이다.
weight = c(8,13,9,12,7,11,6,12,12,14,9,7,14,16,10,14,11,13)
treatment = c("ta","ta","ta","ta","ta","ta","tb","tb","tb","tb","tc","tc","tc",
"tc","tc","tc","tc","tc")
variety = c("va","va","va","vc","vd","vd","va","va","vb","vb","vb","vb","vc",
"vc","vd","vd","vd","vd")
d1 = data.frame(weight, treatment, variety)
options(contrasts = c("contr.sum", "contr.poly"))
Anova(lm(weight ~ treatment*variety, d1), type=3, singular.ok=TRUE) # NOT OKNote: model has aliased coefficients
sums of squares computed by model comparisonAnova Table (Type III tests)
Response: weight
Sum Sq Df F values Pr(>F)
treatment 0.000 0
variety 0.000 0
treatment:variety 34.714 2 3.0995 0.08965 .
Residuals 56.000 10
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1하지만, 다음과 같은 R script로 SAS와 동일한 결과를 얻을 수 있다.
require(sasLM)
aov3(weight ~ treatment*variety, d1)Response : weight
Df Sum Sq Mean Sq F value Pr(>F)
MODEL 7 82.000 11.7143 2.0918 0.13995
treatment 2 12.471 6.2353 1.1134 0.36595
variety 3 34.872 11.6240 2.0757 0.16719
treatment:variety 2 34.714 17.3571 3.0995 0.08965 .
RESIDUALS 10 56.000 5.6000
CORRECTED TOTAL 17 138.000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1위에서 결과가 다르게 나온 이유는 다음과 같이 design matrix의 cross product(이 행렬이 \(\beta\) 계산에 필요하다. 앞 chapter 참고)가 singular 하기 때문이다.
X = model.matrix(~ treatment*variety, d1)
solve(crossprod(X))역행렬이 존재하지 않는 이유는 여러 가지가 있을 수 있는데, 대표적인 것인 aliased coefficient의 존재이다. Aliased coefficient란 다른 변수의 조합으로 설명되는 x 변수가 함께 들어 있다는 뜻이다. 즉 perfect correlation을 갖게 되면 역행렬이 존재하지 않는다.
R에서는 다음과 같이 aliased coefficient를 찾을 수 있다.
alias(weight ~ treatment*variety, d1)Model :
weight ~ treatment * variety
Complete :
(Intercept) treatment1 treatment2 variety1 variety2 variety3
treatment1:variety2 -1/6 -1/3 2/3 -1/2 -1/2 1/2
treatment2:variety2 5/12 -2/3 4/3 -3/4 -3/4 1/4
treatment1:variety3 -5/6 4/3 -2/3 -7/6 3/2 1/2
treatment2:variety3 -5/12 2/3 -1/3 -7/12 3/4 -1/4
treatment1:variety1 treatment2:variety1
treatment1:variety2 1 0
treatment2:variety2 1 -1
treatment1:variety3 2/3 8/3
treatment2:variety3 1/3 4/3 이 외에도, nest design, split-plot design, hierarchical model 등 좀 복잡한 모형의 경우에 car::Anova의 Type 3 SS와 SAS는 다르게 나온다. 대표적인 것이 생물학적 동등성 시험의 자료(bioequivalence data)이다.
Type 1 ~ 4 SS라는 용어는 SAS회사에서 만들어낸 용어로서, 일반적인 통계교과서에는 잘 나오지 않는다.
Type 1 ~ 4 SS에 대하여 가장 먼저 참고할 문서는 SAS User’s Guide에 있는 “The Four Types of Estimable Functions” chapter이다. Chapter 번호는 SAS version에 따라 다르다.
내용은 상당히 난해하여, 이것만으로 R에서 구현하기는 매우 힘들다. 대표적인 원인이 용어의 뜻을 명확히 알기 어려운데, ‘contained’, ‘containing’, ‘associated’ factor/column이 그 예들이다.
Type 1 - 3 SS의 구체적인 계산 실례는 ’SAS Techinical Report R-101. Tests of Hypotheses in Fixed-Effects Linear Models’를 참고한다.
12.2.2 Type 1 SS
R에서 기본적으로 출력해주는 것은 Type 1 SS이다.
Type 1 SS 는 sequential sum of square라고도 하는데, 이는 불균형 자료(unbalanced data)의 경우에 설명 변수를 넣어 주는 순서에 따라 결과가 달라진다다. 균형 자료(balanced data)인 경우에는 설명 변수의 나열 순서와 무관하게 결과는 동일한데, 이 때문에 대부분의 기초교과서에서는 balanced data만 다룬다.
다음은 R에 내장된 CO2 자료를 첫 행을 제거하여 unbalanced data로 만든 뒤에 설명 변수의 순서에 따라 type 1 SS가 달라지는 것을 보인 것이다.
d2 = CO2[-1,]
anova(lm(uptake ~ Type + Treatment, d2))Analysis of Variance Table
Response: uptake
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 3553.5 3553.5 58.050 4.490e-11 ***
Treatment 1 1129.0 1129.0 18.443 4.879e-05 ***
Residuals 80 4897.2 61.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm(uptake ~ Treatment + Type, d2))Analysis of Variance Table
Response: uptake
Df Sum Sq Mean Sq F value Pr(>F)
Treatment 1 1080.5 1080.5 17.651 6.849e-05 ***
Type 1 3602.0 3602.0 58.843 3.558e-11 ***
Residuals 80 4897.2 61.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SAS에서 Type 1 SS를 구하기 위한 contrast matrix를 만드는 방법은 여러 가지가 있는데, 대표적인 것은 forward Doolittle method이다.
sasLM package의 e1 함수가 그것을 쓰고 있다.
ce1 = e1(uptake ~ Type + Treatment, d2) ; zapsmall(ce1) (Intercept) TypeQuebec TypeMississippi Treatmentnonchilled Treatmentchilled
L1 1 0.4939759 0.5060241 0.4939759 0.5060241
L2 0 1.0000000 -1.0000000 -0.0121951 0.0121951
L3 0 0.0000000 0.0000000 0.0000000 0.0000000
L4 0 0.0000000 0.0000000 1.0000000 -1.0000000
L5 0 0.0000000 0.0000000 0.0000000 0.0000000zapsmall은 0에 가까운 값은 0으로 표기하도록 해주어 사람이 좀 더 읽기 쉽게 한다.
이제 REG 결과와 위에서 구한 e1 contrast를 이용하여 SS를 구하면 다음과 같다.
rx = REG(uptake ~ Type + Treatment, d2, summarize=FALSE)
cSS(ce1[2:3, ], rx) # SS for 'Type' varaible Df Sum Sq Mean Sq F value Pr(>F)
1 1 3553.544 3553.544 58.05039 4.489786e-11cSS(ce1[4:5, ], rx) # SS for 'Treatment' variable Df Sum Sq Mean Sq F value Pr(>F)
1 1 1128.998 1128.998 18.44321 4.878752e-05ce1의 첫째 줄은 intercept를 위한 것인데, 일반적으로 ANOVA 결과표에 나오지 않는다. 구하면 다음과 같다.
cSS(ce1[1, , drop=F], rx) # SS for Intercept Df Sum Sq Mean Sq F value Pr(>F)
1 1 62077.66 62077.66 1014.095 0cSS 함수는 contrast를 row vector로 받는다. 그런데, 하나의 행을 고르게 되면 자동으로 column vector 로 바꾸어 버리기 때문에, 이를 방지하기 위해 drop=F 라는 옵션을 사용하였다.
일반적으로 분산 분석에서는 intercept가 있는 모형이며 corrected SS를 보여준다. 따라서, corrected sum of square total (corrected SST)의 자유도는 n - 1이다. Intercept가 없는 모형을 지정해 줄 수도 있는데, 이때는 uncorrected SST가 출력되며, 자유도는 n이다. 이를 수식으로 나타내면 다음과 같다.
corrected sum of squares total = cSST = \(\sum (y_i - \bar{y})^2\) = \((y - \bar{y})^T (y - \bar{y})\) = crossprod(\(y - \bar{y}\))
uncorrected sum of squares total = uSST = \(\sum y_i^2\) = \(y^T y\) = crossprod(y)
12.2.3 Type 2 SS
Type 2 SS를 위한 contrast는 e2 함수로 구할 수 있다.
Type 2 SS는 동일한 차수의 설명 항(term)들에 대해서는 순서와 상관없이 같게 나온다. 즉, 동일한 차수끼리는 서로 보정해주며, 하위 차수의 SS를 먼저 계산한 뒤, 상위 차수를 계산한다.
ce2 = e2(uptake ~ Type + Treatment, d2) ; ce2 (Intercept) TypeQuebec TypeMississippi Treatmentnonchilled Treatmentchilled
L1 1 0 0 0 0
L2 0 1 -1 0 0
L3 0 0 0 0 0
L4 0 0 0 1 -1
L5 0 0 0 0 0cSS(ce2[2:3,], rx) # SS for 'Type' Df Sum Sq Mean Sq F value Pr(>F)
1 1 3602.033 3602.033 58.8425 3.558154e-11cSS(ce2[4:5,], rx) # SS for 'Treatment' Df Sum Sq Mean Sq F value Pr(>F)
1 1 1128.998 1128.998 18.44321 4.878752e-05Regression에 의한 \(\beta\) 계수는 type 1 ~ 4 SS와는 무관하므로 refit할 필요는 없다.
Type 2 SS는 비록 g2 inverse를 쓰기는 하지만 행렬 연산으로 잘 정의되어 있기 때문에 (물론 contained, containing, associated column/factor의 뜻을 먼저 알아야 한다.) software가 달라도 일정한 결과를 보여준다. 즉, car::Anova 함수의 type=2 option은 비교적 만족스러운 결과를 보여준다.
R 진영에서도 Type 2 SS에 대한 비판은 보이지 않으며, Type 3 SS보다 Type 2 SS를 더 추천한다는 의견들이 있다.
12.2.4 Type 3 SS
Type 3 SS는 차수와도 상관없이 포함된 모든 설명 항(term)들에 대하여 대등한 자격으로 서로를 보정해 준 것이다. SAS에서 우선적으로 이것을 사용하라고 추천하여, 여러 통계학자들(특히 R 진영)로부터 비판을 받았다.
Type 3 SS는 Type 1이나 Type 2에 비하여 계산하기 쉽지 않은 면이 있다. Unbalanced data에 대하여 SS의 합들이 total SS와도 맞지 않는다.
차수가 다른 것을 동등한 입장에서 다루어야 할지 말지는, 또는 설명 변수를 넣는 순서에 따라 다르게 나오는 것이 좋을지 그렇지 않을지 등은, 그 분야 전문가의 domain knowledge를 바탕으로 판단해야 하기 때문에 blind하게 Type 3 SS를 사용하는 것은 비판을 받는 것이 당연하다.
만약 A, B, C의 설명 변수가 있고, 연구자는 A, B의 효과를 배제하고도 C가 효과가 있는지 보고자 한다면, Type 1 SS를 사용하되, A와 B를 C보다 앞에 두고 ANOVA분석을 해야한다. 하지만, 그렇지 않은 경우를 흔히 보게 된다.
Type 1 SS는 누가 봐도 이견이 없고, 계산이 매우 명확(나열한 설명 변수의 순서에 따라 SS를 제거해 나간다. 그래서 최종적으로 남은 것이 residual SS이다.)하기 때문에 장점이 있다. 사실 Type 2 SS와 Type 3 SS의 하단에 나오는 Residual SS는 Type 1 SS로 계산한 것이다.
sasLM package의 e3 함수로 contrast (estimable function)를 볼 수 있다.
다음은 GLM 함수로 다음 모형의 Type 1, 2, 3 SS가 모두 다름을 볼 수 있다.
아래 e3 contrast (etimable function)는 보기 좋도록 transpose 하였다.
rx2 = REG(uptake ~ Type*Treatment, d2, summarize=F)
ce3 = e3(uptake ~ Type*Treatment, d2) ; t(zapsmall(ce3)) L1 L2 L3 L4 L5 L6 L7 L8 L9
(Intercept) 1 0.0 0 0.0 0 0 0 0 0
TypeQuebec 0 1.0 0 0.0 0 0 0 0 0
TypeMississippi 0 -1.0 0 0.0 0 0 0 0 0
Treatmentnonchilled 0 0.0 0 1.0 0 0 0 0 0
Treatmentchilled 0 0.0 0 -1.0 0 0 0 0 0
TypeQuebec:Treatmentchilled 0 0.5 0 -0.5 0 1 0 0 0
TypeQuebec:Treatmentnonchilled 0 0.5 0 0.5 0 -1 0 0 0
TypeMississippi:Treatmentchilled 0 -0.5 0 -0.5 0 -1 0 0 0
TypeMississippi:Treatmentnonchilled 0 -0.5 0 0.5 0 1 0 0 0cSS(ce3[2:3,], rx2) # SS for Type Df Sum Sq Mean Sq F value Pr(>F)
1 1 3582.646 3582.646 59.77216 2.905065e-11cSS(ce3[4:5,], rx2) # SS for Treatment Df Sum Sq Mean Sq F value Pr(>F)
1 1 1118.29 1118.29 18.65733 4.505622e-05cSS(ce3[6:9,], rx2) # SS for Type:Treatment Df Sum Sq Mean Sq F value Pr(>F)
1 1 162.0548 162.0548 2.703691 0.1040902GLM(uptake ~ Type*Treatment, d2)$ANOVA
Response : uptake
Df Sum Sq Mean Sq F value Pr(>F)
MODEL 3 4844.6 1614.87 26.942 4.208e-12 ***
RESIDUALS 79 4735.1 59.94
CORRECTED TOTAL 82 9579.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$Fitness
Root MSE uptake Mean Coef Var R-square Adj R-sq
7.741987 27.34819 28.30895 0.5057134 0.4869431
$`Type I`
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 3553.5 3553.5 59.2866 3.344e-11 ***
Treatment 1 1129.0 1129.0 18.8360 4.177e-05 ***
Type:Treatment 1 162.1 162.1 2.7037 0.1041
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type II`
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 3602.0 3602.0 60.0956 2.646e-11 ***
Treatment 1 1129.0 1129.0 18.8360 4.177e-05 ***
Type:Treatment 1 162.1 162.1 2.7037 0.1041
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$`Type III`
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 3582.6 3582.6 59.7722 2.905e-11 ***
Treatment 1 1118.3 1118.3 18.6573 4.506e-05 ***
Type:Treatment 1 162.1 162.1 2.7037 0.1041
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1