AI Analytics, by YD Lee
회귀분석 Regression, Part IV
Topics:
Part I:
iris 자료
단순 선형회귀
Part II
사영과 회귀분석
다중 선형회귀
Part III
가변수와 선형모형
모형의 다양성
Part IV
모형 선택
회귀 진단
1. 모형 선택
[모형선택]
위 그림의 여러 회귀 모형 중 가장 좋은 것은 어떤 것인가?
직선: 가장 단순한 모형
곡선: 중간 형태의 모형
불규칙 직선: 관측점을 완전 적합한 모형
좋은 모형이란 무엇일까?
관측점들을 잘 적합하는 모형 ?
단순한 모형 ?
적절한 타협이 필요하다면 어떤 기준으로 ?
가설검정에 의한 모형의 비교:
가설검정: 두 가지 모형, 귀무가설과 대립가설 중 우수한 모형 선택
전통적인 가설검정 방법에서,
귀무가설 \(~H_0~\) : 축약된 모형 reduced model
대립가설 \(~H_1~\) : 귀무가설을 포함한 넓은 모형 full model
즉, reduced model \(~\subset~\) full model
가설검정에서 단순성의 추구
가설검정에서는 기본적으로 단순한 모형 (\(~H_0~\)) 가 옳다고 가정함
관측된 자료가 이 가정(\(H_0\) 가 옳다는 가정) 에서 도저히 수용할 수 없을 때,
복잡한 모형 (\(~H_1~\)) 이 옳다고 결론을 냄
유의확률(p-value) 이란, \(H_0\) 가 옳다는 가정이 관측자료를 수용하는 정도를 말함
유의수준 (\(\alpha\)-value) 이란, \(H_0\) 가 옳다는 가정이 관측자료를 수용할 수 있는 정도의 하한
p-value 가 \(\alpha\)-value 보다 작아지면, \(~H_1~\) 이 옳다고 결론을 내림
참고: 일반적인 가설검정의 방법:
유의확률 p-value 와 유의수준 \(~\alpha\)-value 를 비교하여 결정함
p-value \(~\geq~\alpha\)-value : 귀무가설 \(~H_0~\) 이 옳다고 결론을 내림
p-value \(~\leq ~\alpha\)-value : 대립가설 \(~H_1~\) 이 옳다고 결론을 내림
등호의 경우는 발생하지 않는 무시할 만한 상황이므로 신경쓸 필요 없음
\(\alpha\)-value 는 보통 0.05 즉 5% 인 경우가 많이 사용되지만, 바꿀 수 있음
결국, p-value 가 낮을수록 대립가설 \(~H_1~\) 이 지지된다는 의미임
비유해서 말하자면, 가설검정은 \(~H_0~\) 이 시험을 치는 것과 동일
p-value: \(~H_0~\) 의 시험성적
\(\alpha\)-value: 커트라인
시험성적이 커트라인 보다 낮으면, 대립가설 \(~H_1~\) 이 합격
모형선택 방법으로써 가설검정의 한계:
포함관계가 없는 모형 사이에서의 모형 선택 불가 ( Bayesian 접근법에서는 가능함 )
단순한 모형에 대한 일방적 선호
적절한 모형선택 방법은 없을까? 있다.
그러면 기준은 ? 도대체 어떤 개념으로 ?
모형선택에 사용될 수 있는 기준들:
- \(R^2~\), \(~R_{adj}^2~\), MSE, F-value, p-value, \(~\log(\mbox{p-value})\)
수 많은 모형선택 방법들:
Mallow’s \(~C_p~\)
AIC, BIC
MAIC, AICb, AICc, AICi, …
EIC, DIC, GIC, TIC, FPE, HQ, …
모형선택과 관련하여 학습하여야 할 개념과 방법들:
수 많은 모형선택 방법들이 있음
먼저 Cp, AIC 와 BIC,
그리고 PRESS, CV, GCV 등에 대하여 살펴본다
일반적인 모형에 대한 논의는, likelihood 의 개념을 공부한 이후에 보기로 하고
먼저 회귀분석의 경우를 중심으로, 그 개념과 사용법을 살펴보기로 한다
AIC 와 Cp
AIC와 Cp:
[AIC & Cp]
정규모집단 모형에서의 기본적인 관계:
관측된 자료: \(~\mathbf{y} = (y_1 , \ldots, y_n )'\)
모형: \(~Y_j = \mu +\epsilon_j, ~~\epsilon_j ~\stackrel{iid}{\sim} ~N(0,\sigma^2), ~~j=1,\ldots,n\)
추정값: \(~\hat{\mathbf{y}}= \mathbf{1} \, \hat{\mu},~~\) 여기서 \(~\hat{\mu} = \bar{y}_n\)
관측되지 않은 자료: \(~Y_j^* = \mu +\epsilon_j^*,~~\epsilon_j^* ~\stackrel{iid}{\sim} ~N(0,\sigma^2), ~~j=1,\ldots,n\)
\(\hat{\mu}-\mu = \bar{y}_n -\mu = (1/n) \sum_j (y_j -\mu) = \bar{\epsilon}_n\)
\(A = E\left( n \, (\mu-\hat{\mu})^2 \right) = n \, Var( \bar{\epsilon}_n ) = \sigma^2\)
\(B= E \left( \sum_j (Y_j^* -\hat{\mu} )^2 \right) = n \, Var( \bar{\epsilon}_n ) +Var \left( \sum_j \epsilon_j^* \right) - 2\, Cov \left( \sum_j \epsilon_j^* , \bar{\epsilon}_n \right)\)
\(\qquad = \sigma^2 + n \, \sigma^2 - 0\)
\(C= E \left( \sum_j (Y_j -\hat{\mu} )^2 \right) = n \, Var( \bar{\epsilon}_n ) +Var \left( \sum_j \epsilon_j \right) - 2\, Cov \left( \sum_j \epsilon_j, \bar{\epsilon}_n \right)\)
\(\qquad = \sigma^2 + n \, \sigma^2 - 2 \, \sigma^2\)
\(A/\sigma^2~ \approx ~\sum_j (y_j -\hat{\mu})^2 / \sigma^2 -n +2\)
\(B/\sigma^2~ \approx ~\sum_j (y_j -\hat{\mu})^2 / \sigma^2 +2\)
회귀모형에서의 관계:
관측된 자료: \(~\mathbf{y} = (y_1 , \ldots, y_n )'\)
모형: \(~Y_j = \mu_j +\epsilon_j,~~\mu_j =\mathbf{x}_j ' \boldsymbol{\beta}, ~~\epsilon_j ~\stackrel{iid}{\sim} ~N(0,\sigma^2), ~~j=1,\ldots,n\)
추정값: \(~\hat{\mathbf{y}}= X \hat{\boldsymbol{\beta}},~~\) 여기서 \(~\hat{\boldsymbol{\beta}} = (X'X)^{-1} X'\mathbf{y}\)
관측되지 않은 자료: \(~Y_j^* = \mu_j +\epsilon_j^*,~~\mu_j =\mathbf{x}_j ' \boldsymbol{\beta},~~\epsilon_j^* ~\stackrel{iid}{\sim} ~N(0,\sigma^2), ~~j=1,\ldots,n\)
\(X \hat{\boldsymbol{\beta}} - X \boldsymbol{\beta} = P (\mathbf{Y}-E(\mathbf{Y})) = P \boldsymbol{\epsilon}, ~~\) 여기서 \(~P=X(X'X)^{-1} X'\)
\(\sum_j ( \mu_j - \hat{\mu}_j )^2 =\sum_j ( \mathbf{x}_j ' \boldsymbol{\beta} - \mathbf{x}_j ' \hat{\boldsymbol{\beta}} )^2 = \boldsymbol{\epsilon} ' P \boldsymbol{\epsilon}, ~~\) (why?)
\(A = E \left( \sum_j ( \mu_j - \hat{\mu}_j )^2 \right) = E(\boldsymbol{\epsilon} ' P \boldsymbol{\epsilon}) = k \, \sigma^2, ~~\) (why?)
\(B= E \left( \sum_j (Y_j^* -\hat{\mu}_j )^2 \right) = E(\boldsymbol{\epsilon} ' P \boldsymbol{\epsilon}) + E \left( \sum_j (\epsilon_j^*)^2 \right) - 2\, \sum_j Cov(\epsilon_j^* , (\hat{\mu}_j -\mu))\)
\(\qquad = k \, \sigma^2 + n \, \sigma^2 + 0\)
\(C= E \left( \sum_j (Y_j -\hat{\mu}_j )^2 \right) = E(\boldsymbol{\epsilon} ' P \boldsymbol{\epsilon}) + E \left( \sum_j \epsilon_j^2 \right) - 2\, E(\boldsymbol{\epsilon} ' P \boldsymbol{\epsilon})\)
\(\qquad = k \, \sigma^2 + n\, \sigma^2 - 2 k \, \sigma^2\)
\(A/\sigma^2~\approx~\sum_j (y_j -\hat{\mu}_j )^2 / \sigma^2 -n +2k\)
\(B/\sigma^2~\approx~\sum_j (y_j -\hat{\mu}_j )^2 / \sigma^2 +2k\)
\(A~\)와 \(~B~\)의 의미:
평균 제곱 모형추정오차: \(~A = E \left( \sum_j ( \mu_j - \hat{\mu}_j )^2 \right)\)
평균 제곱 예측오차: \(~B= E \left( \sum_j (Y_j^* -\hat{\mu}_j )^2 \right)\)
예측오차 (\(B\)) 와 모형추정오차 (\(A\)) 의 관심사항은 동일함
의미해석 측면에서, 예측오차 (\(B\)) 의 개념이 이해하기 쉽고,
CV 와 같은 다른 방법들과의 연관성을 찾기 쉽다는 장점이 있음
회귀모형에서, 오차분산 \(~\sigma^2~\)이 알려져 있는 경우:
\(C_p~=~SSE/ \sigma^2 -n +2k\)
\(AIC~=~ SSE/ \sigma^2 + 2k\)
\(k~\)는 독립변수로 사용된 변수의 개수, 추정할 모수의 개수
회귀모형에서, 오차분산 \(~\sigma^2~\)에 대한 추정이 필요한 경우:
\(C_p~=~SSE/ \hat{\sigma}^2 -n +2k\),
여기서 \(~\hat{\sigma}^2~\) 은 full model 에서 추정된 불편추정량
\(C_p~\) 는 전체 full model 의 sub model 중 우수한 모형을 선택하기 위하여 개발됨
가설에 대한 F-검정에서와 같이, full model, reduced model 관계가 발생함
\(AIC~=~ SSE/ \hat{\sigma}^2 + 2(k+1)\),
\(SSE/ \hat{\sigma}^2 = -2 \log(\mbox{likelihood})\) 이고,
여기서 \(~\hat{\sigma}^2~\) 은 최대우도추정량 (편의추정량)
\((k+1)~\) 이 되는 이유는 분산 모수가 추가되기 때문임
\(y\)-절편이 있는 모형이고, 독립변수의 개수가 \(~p~\) 라면,
\(k=(p+1)\) 이므로, \(~(k+1)~\)은 \(~(p+2)~\)를 의미함
회귀모형에서의 BIC
\(BIC~=~ SSE/ \hat{\sigma}^2 + \log(n)\, (k+1)\),
\(SSE/ \hat{\sigma}^2 = -2 \log(\mbox{likelihood})\) 이고,
\(n~\) 은 관측자료의 개수
BIC 는 AIC 에 비하여 단순한 모형을 최적모형으로 선정함 (why?)
BIC 는 자료에 대한 모형의 적합성을 지수화 함. 모형의 자료 설명력을 의미함
\(C_p~\) 와 F-value (\(~F~\)) 의 관계:
모형
M0 : \(~Y_j = \mu_{0j} +\epsilon_j, ~~~ \mu_{0j}= \mathbf{x}_{0j} ' \boldsymbol{\beta}_0\)
M1 : \(~Y_j = \mu_{1j} + \epsilon_j, ~~~ \mu_{1j}= \mathbf{x}_{0j} ' \boldsymbol{\beta}_0 + \mathbf{x}_{1j} ' \boldsymbol{\beta}_1\)
가설: \(~~H_0 : \boldsymbol{\beta}_2 = \mathbf{0}~(M0),~~\) vs. \(~~H_1 : \boldsymbol{\beta}_2 \neq \mathbf{0}~(M1)\)
rank, SSE & MSE
행렬 \(~X_0=(\mathbf{x}_{0j} ')_{n \times k_0 }~\), \(~X_1=(\mathbf{x}_{1j} ')_{n \times k_1 }~\) 는 full column matrix,
\(~rank(X_0)=k_0~\), \(~rank(X_1)=k_1~\)
\(SSE(M0)=\sum_j (y_j - \hat{\mu}_{0j})^2\)
\(SSE(M1)=\sum_j (y_j - \hat{\mu}_{1j})^2\)
\(MSR = (SSE(M0) - SSE(M1))/k_1\)
\(MSE(M1) = SSE(M1)/(n- k_0 - k_1)\)
\(F=MSR/MSE(M1)\)
\(\quad = [(n- k_0 - k_1)/k_1]\cdot (SSE(M0) - SSE(M1))/SSE(M1)\)
\(\quad = [(n- k_0 - k_1)/k_1]\cdot \left[ SSE(M0)/SSE(M1) - 1 \right]\)
\(C_p = SSE(M0)/MSE(M1) -n +2 k_0\)
\(\quad = (n-k_0 -k_1) SSE(M0)/SSE(M1) - n + 2 k_0\)
\(\quad = k_1 \cdot F - k_1 + k_0\)
\(M_0~\) 와 \(~M_1~\) 이 동일한 모형이 되면, \(~C_p = k_0 - k_1\)
AIC 와 Cp 의 비교:
AIC 와 Cp 모두 모형의 예측능력을 평가하는 지수임
AIC 의 장점:
full model 을 설정할 필요가 없어, 사용이 자유로움
우도 likelihood 이론을 적용할 수 있음
회귀분석 뿐만 아니라 다양한 통계모형에 적용이 가능함.
최근에는 Cp 를 사용하는 경우는 거의 없고, 주로 AIC 를 사용함
R 에서 Mallows Cp 계산:
R 에서 Cp 를 계산하는 함수가 들어 있는 패키지가 거의 없음
Cp 함수가 있다고 알려진 R 패키지
leaps: Thomas Lumley 가 2017년에 FORTRAN code 를 이식
olsrr: 2018년 말에 올라온 패키지
두 패키지 모두, 명목변수(factor 변수)를 제대로 처리하지 못함
leaps: factor 변수가 있는 경우 오류를 발생시킴
oslrr: factor 변수가 있는 경우 잘못된 Cp 값을 제시
Cp, PRESS, GCV :
현실 자료 분석에서는 AIC 를 사용하면 되므로
Cp 함수가 꼭 필요하지는 않지만, 교육 목적 상 필요하다.
그래서 그냥 내가 만들었다.
만드는 김에, 굳이 다른 패키지 찾아보기 뭐해서,
아래에서 사용할, PRESS, CV, GCV 함수도 함께 만들었고,
여러 모형들을 비교하려면,
주어진 모형에 대한 모든 sub model 들을
자동으로 생성하는 함수도 필요해서 함께 만들었다.
sub_models <- function(a_lm) {
colon <- function(a_colon){
tp = strsplit(a_colon,split=':')[[1]]
k = length(tp)
term = as.name(tp[1])
if( k > 1 ) for(j in 2:k ) term = call(':', term ,as.name(tp[j]))
term
}
integer_2_binary <- function(a_integer) {
twopower = -1 + 2^a_integer
t(matrix( as.integer(intToBits(0:twopower)), nrow=32 ))[,a_integer:1]
}
f_model = formula(a_lm$terms)
x_vars = attributes(a_lm$terms)$term.labels
offset = any(attributes(a_lm)$names == 'offset') ## offsset will be neglect here
one_zero = as.integer( any(names(coef(a_lm))=='(Intercept)') )
k = length(x_vars)
if(k == 0) {
x_model = 2*one_zero - 1
f_model[[3]] = x_model
sub_models = format(f_model)
}
else {
binary_sequences = integer_2_binary(k)
sub_models = NULL
for(i in 1:nrow(binary_sequences)) {
yes = binary_sequences[i,]
x_model = 2*one_zero - 1
for(j in 1:k ) if(!yes[j]) x_model = call('+', x_model, colon(x_vars[j]) )
f_model[[3]] = x_model
sub_models = c(sub_models, format(f_model) )
}
}
sub_models
}Cp <-function(R,F){
k1 <- F$rank-R$rank
k0 <- R$rank
if( is.na(anova(R,F)$F[2]) ) F<- 0 else F <- anova(R,F)$F[2]
k1*F -k1+k0 }
PRESS <- function(r){
e<- resid(r);
X<-model.matrix(r);
H <- X %*% solve(t(X) %*% X) %*% t(X);
mean(e*e/(1-diag(H))^2)
}
CV <- function(r){
X <- model.matrix(r)
rsx <- rep(0,nrow(X))
dfx <- eval(r$call[[3]])
d_var<- as.character(r$call[[2]][[2]])
for(i in 1:nrow(X)) {
rew <- eval(bquote(lm(.(r$call[[2]]),dfx[-i,])))
rsx[i]<-dfx[i,d_var] - sum(X[i,] * coef(rew) )
}
mean(rsx*rsx)
}
GCV <- function(r){
e<- resid(r);
X<-model.matrix(r);
H <- X %*% solve(t(X) %*% X) %*% t(X);
mean(e*e)/(1-mean(diag(H)))^2
}
ABC <- function(M, FULL_Model) {
rex <- eval(bquote( lm(.(M), .(FULL_Model$call[[3]]) )))
ONE <- eval(bquote( lm(.(FULL_Model$call[[2]][[2]])~1, .(FULL_Model$call[[3]]) )))
p0 <- anova(ONE, rex)$'Pr(>F)'[2]; if(is.na(p0)) p0 <- 1
p1 <- anova(rex, FULL_Model)$'Pr(>F)'[2] ; if(is.na(p1)) p1 <- 1
c( model_df = rex$rank,
Cp = Cp(rex, FULL_Model),
AIC = AIC(rex),
BIC = BIC(rex),
PRESS= PRESS(rex),
GCV = GCV(rex)
)
}iris 데이터에서의 AIC 와 Cp :
iris 자료에 대한 기본적인 모형에 대한 여러 가지 sub model 들을 구하고
각 sub model 들에 대하여 Cp 와 AIC 를 계산하였다
각 sub model 들에서 Cp 와 AIC 크기 순서가 일치함을 볼 수 있다
names(iris)<-c('sl','sw','pl','pw','sp')
res <- lm(sl~.,iris)
sub_models(res)## [1] "sl ~ 1 + sw + pl + pw + sp" "sl ~ 1 + sw + pl + pw"
## [3] "sl ~ 1 + sw + pl + sp" "sl ~ 1 + sw + pl"
## [5] "sl ~ 1 + sw + pw + sp" "sl ~ 1 + sw + pw"
## [7] "sl ~ 1 + sw + sp" "sl ~ 1 + sw"
## [9] "sl ~ 1 + pl + pw + sp" "sl ~ 1 + pl + pw"
## [11] "sl ~ 1 + pl + sp" "sl ~ 1 + pl"
## [13] "sl ~ 1 + pw + sp" "sl ~ 1 + pw"
## [15] "sl ~ 1 + sp" "sl ~ 1"abc_table <- do.call(rbind, lapply(sub_models(res), function(M) ABC(M,res) ))
abc_table <- as.data.frame(abc_table)
round(abc_table,3)## model_df Cp AIC BIC PRESS GCV
## 1 6 6.000 79.116 100.190 0.098 0.098
## 2 4 11.442 84.643 99.696 0.102 0.102
## 3 5 8.345 81.575 99.639 0.100 0.100
## 4 3 29.448 101.025 113.068 0.113 0.113
## 5 5 150.431 182.349 200.413 0.196 0.195
## 6 3 173.722 191.821 203.863 0.208 0.208
## 7 4 155.461 183.937 198.990 0.197 0.197
## 8 2 924.254 371.992 381.024 0.689 0.690
## 9 5 37.194 108.231 126.295 0.119 0.119
## 10 3 109.666 158.047 170.089 0.166 0.166
## 11 4 35.196 106.233 121.286 0.117 0.117
## 12 2 114.510 160.040 169.072 0.168 0.168
## 13 4 216.822 212.068 227.121 0.238 0.238
## 14 2 213.189 208.221 217.253 0.232 0.232
## 15 3 269.801 231.452 243.494 0.270 0.270
## 16 1 937.255 372.080 378.101 0.690 0.690library(ggplot2)
ggplot(abc_table) + geom_point(aes(x=Cp,y=AIC))
AIC 를 이용한 모형선택
R 에서의 AIC 와 BIC 계산:
R 에는 AIC 와 extractAIC, BIC 함수가 있음
AIC :
extractAIC : 값은 다르지만, AIC 와 같은 역할을 함.
BIC : 설명력이 높은 모형 선정. AIC 보다 단순한 모형을 선정
R 의 step 함수:
\(k~\) 개의 독립변수를 가진 데이터에서, 구성 가능한 회귀모형은 \(~2^k~\) 개임
독립변수의 개수가 20개가 되면, 1백만 개 이상의 회귀모형 구성이 가능함
많은 수의 회귀모형에 대하여 AIC 를 계산하여 최적의 모형을 찾는 함수가 필요함
R 에는 이를 위하여, step 이라는 함수와 Full Version 인 MASS::stepAIC 가 있음.
step 함수는 tree 탐색방법을 이용함 (모든 가능한 모형을 탐색하는 것은 아님)
step 함수에서 AIC 라고 제공하는 값은 extractAIC 임
왜 BIC 를 최적화하는 모형 탐색 함수는 없을까?
BIC 는 AIC 보다 단순한 모형을 최적모형으로 잡게 됨
따라서, step 을 이용하여 AIC 가 작은 모형을 선정하게 되면
선정된 모형의 sub model 중 BIC 가 작은 모형을 선정할 수
있으므로 굳이 별도의 함수가 필요하지 않을 수 있음.
step 함수의 인자 k 는 기정값이 2로 설정되어 있는데,
이를 k=log(n) 으로 대입하면, BIC 가 가장 작은 모형을 선정하게 됨.
아래에서는 step 으로 구한 모형의 모든 sub model 들의 BIC 를 구함
다음 예에서,
res1 의 경우 step 함수에 의한 탐색 과정이 길어, 그 과정을 출력하지 않고,
결과만 나타나도록 함. trace=0
처음 step 함수를 시작할 때, res1 이 res2 보다 더 많은 설명변수로 시작함
결과로 선정된 모형에서 res2 의 AIC 가 더 낮은 값을 가짐
step 함수가 모든 가능한 모형을 탐색하는 것이라면,
res1 의 경우가 res2 보다 더 낮은 AIC 를 갖는 모형을 선정하여야 하지만
결과는 그렇지 않음.
( res1<- step(lm(sl~ (sw+pl+pw+sp+log(sw)+log(pl)+log(pw))^2, iris),trace=0 ) )##
## Call:
## lm(formula = sl ~ sw + pl + pw + sp + log(sw) + log(pl) + log(pw) +
## sw:log(sw) + pl:pw + pl:log(pw) + pw:sp + sp:log(sw) + log(sw):log(pl) +
## log(pl):log(pw), data = iris)
##
## Coefficients:
## (Intercept) sw pl
## -35.823084 42.450255 2.147852
## pw spversicolor spvirginica
## 1.643967 -6.911207 -11.376370
## log(sw) log(pl) log(pw)
## -42.832818 5.503248 -3.160912
## sw:log(sw) pl:pw pl:log(pw)
## -12.839293 -1.493102 3.324021
## pw:spversicolor pw:spvirginica spversicolor:log(sw)
## -0.008256 2.173295 6.347488
## spvirginica:log(sw) log(sw):log(pl) log(pl):log(pw)
## 6.916618 -5.309039 -3.980517AIC(res1)## [1] 74.09104extractAIC(res1)## [1] 18.0000 -353.5905print('#### #### ####')## [1] "#### #### ####"res2<- step(lm(sl~ (sw+pl+pw+sp)^2, iris) )## Start: AIC=-348.13
## sl ~ (sw + pl + pw + sp)^2
##
## Df Sum of Sq RSS AIC
## - sw:sp 2 0.076800 12.135 -351.18
## - sw:pl 1 0.000414 12.059 -350.12
## - pl:pw 1 0.011712 12.070 -349.98
## - pl:sp 2 0.183280 12.242 -349.87
## - pw:sp 2 0.273788 12.332 -348.76
## <none> 12.059 -348.13
## - sw:pw 1 0.162166 12.221 -348.12
##
## Step: AIC=-351.18
## sl ~ sw + pl + pw + sp + sw:pl + sw:pw + pl:pw + pl:sp + pw:sp
##
## Df Sum of Sq RSS AIC
## - sw:pl 1 0.02651 12.162 -352.85
## - pl:pw 1 0.06908 12.204 -352.32
## - sw:pw 1 0.10389 12.239 -351.90
## <none> 12.135 -351.18
## - pw:sp 2 0.38586 12.521 -350.48
## - pl:sp 2 0.48227 12.618 -349.33
##
## Step: AIC=-352.85
## sl ~ sw + pl + pw + sp + sw:pw + pl:pw + pl:sp + pw:sp
##
## Df Sum of Sq RSS AIC
## - pl:pw 1 0.07691 12.239 -353.90
## <none> 12.162 -352.85
## - pw:sp 2 0.36206 12.524 -352.45
## - pl:sp 2 0.45639 12.618 -351.32
## - sw:pw 1 0.42641 12.588 -349.68
##
## Step: AIC=-353.9
## sl ~ sw + pl + pw + sp + sw:pw + pl:sp + pw:sp
##
## Df Sum of Sq RSS AIC
## <none> 12.239 -353.90
## - pw:sp 2 0.38440 12.623 -353.26
## - pl:sp 2 0.61094 12.850 -350.60
## - sw:pw 1 0.44186 12.681 -350.58print('**** **** ****')## [1] "**** **** ****"AIC(res2)## [1] 73.77855extractAIC(res2)## [1] 11.000 -353.903model_set <- sub_models(res2)
head(model_set, 40)## [1] "sl ~ 1 + sw + pl + pw + sp + sw:pw + pl:sp + pw:sp"
## [2] "sl ~ 1 + sw + pl + pw + sp + sw:pw + pl:sp"
## [3] "sl ~ 1 + sw + pl + pw + sp + sw:pw + pw:sp"
## [4] "sl ~ 1 + sw + pl + pw + sp + sw:pw"
## [5] "sl ~ 1 + sw + pl + pw + sp + pl:sp + pw:sp"
## [6] "sl ~ 1 + sw + pl + pw + sp + pl:sp"
## [7] "sl ~ 1 + sw + pl + pw + sp + pw:sp"
## [8] "sl ~ 1 + sw + pl + pw + sp"
## [9] "sl ~ 1 + sw + pl + pw + sw:pw + pl:sp + pw:sp"
## [10] "sl ~ 1 + sw + pl + pw + sw:pw + pl:sp"
## [11] "sl ~ 1 + sw + pl + pw + sw:pw + pw:sp"
## [12] "sl ~ 1 + sw + pl + pw + sw:pw"
## [13] "sl ~ 1 + sw + pl + pw + pl:sp + pw:sp"
## [14] "sl ~ 1 + sw + pl + pw + pl:sp"
## [15] "sl ~ 1 + sw + pl + pw + pw:sp"
## [16] "sl ~ 1 + sw + pl + pw"
## [17] "sl ~ 1 + sw + pl + sp + sw:pw + pl:sp + pw:sp"
## [18] "sl ~ 1 + sw + pl + sp + sw:pw + pl:sp"
## [19] "sl ~ 1 + sw + pl + sp + sw:pw + pw:sp"
## [20] "sl ~ 1 + sw + pl + sp + sw:pw"
## [21] "sl ~ 1 + sw + pl + sp + pl:sp + pw:sp"
## [22] "sl ~ 1 + sw + pl + sp + pl:sp"
## [23] "sl ~ 1 + sw + pl + sp + pw:sp"
## [24] "sl ~ 1 + sw + pl + sp"
## [25] "sl ~ 1 + sw + pl + sw:pw + pl:sp + pw:sp"
## [26] "sl ~ 1 + sw + pl + sw:pw + pl:sp"
## [27] "sl ~ 1 + sw + pl + sw:pw + pw:sp"
## [28] "sl ~ 1 + sw + pl + sw:pw"
## [29] "sl ~ 1 + sw + pl + pl:sp + pw:sp"
## [30] "sl ~ 1 + sw + pl + pl:sp"
## [31] "sl ~ 1 + sw + pl + pw:sp"
## [32] "sl ~ 1 + sw + pl"
## [33] "sl ~ 1 + sw + pw + sp + sw:pw + pl:sp + pw:sp"
## [34] "sl ~ 1 + sw + pw + sp + sw:pw + pl:sp"
## [35] "sl ~ 1 + sw + pw + sp + sw:pw + pw:sp"
## [36] "sl ~ 1 + sw + pw + sp + sw:pw"
## [37] "sl ~ 1 + sw + pw + sp + pl:sp + pw:sp"
## [38] "sl ~ 1 + sw + pw + sp + pl:sp"
## [39] "sl ~ 1 + sw + pw + sp + pw:sp"
## [40] "sl ~ 1 + sw + pw + sp"tail(model_set, 30)## [1] "sl ~ 1 + pw + sp + sw:pw + pw:sp" "sl ~ 1 + pw + sp + sw:pw"
## [3] "sl ~ 1 + pw + sp + pl:sp + pw:sp" "sl ~ 1 + pw + sp + pl:sp"
## [5] "sl ~ 1 + pw + sp + pw:sp" "sl ~ 1 + pw + sp"
## [7] "sl ~ 1 + pw + sw:pw + pl:sp + pw:sp" "sl ~ 1 + pw + sw:pw + pl:sp"
## [9] "sl ~ 1 + pw + sw:pw + pw:sp" "sl ~ 1 + pw + sw:pw"
## [11] "sl ~ 1 + pw + pl:sp + pw:sp" "sl ~ 1 + pw + pl:sp"
## [13] "sl ~ 1 + pw + pw:sp" "sl ~ 1 + pw"
## [15] "sl ~ 1 + sp + sw:pw + pl:sp + pw:sp" "sl ~ 1 + sp + sw:pw + pl:sp"
## [17] "sl ~ 1 + sp + sw:pw + pw:sp" "sl ~ 1 + sp + sw:pw"
## [19] "sl ~ 1 + sp + pl:sp + pw:sp" "sl ~ 1 + sp + pl:sp"
## [21] "sl ~ 1 + sp + pw:sp" "sl ~ 1 + sp"
## [23] "sl ~ 1 + sw:pw + pl:sp + pw:sp" "sl ~ 1 + sw:pw + pl:sp"
## [25] "sl ~ 1 + sw:pw + pw:sp" "sl ~ 1 + sw:pw"
## [27] "sl ~ 1 + pl:sp + pw:sp" "sl ~ 1 + pl:sp"
## [29] "sl ~ 1 + pw:sp" "sl ~ 1"abc_table <- do.call(rbind, lapply(model_set, function(M) ABC(M,res2) ))
abc_table <- as.data.frame(abc_table)
round(head(abc_table,40), 3)## model_df Cp AIC BIC PRESS GCV
## 1 11 11.000 73.779 109.906 0.095 0.095
## 2 9 11.366 74.417 104.524 0.095 0.095
## 3 9 13.939 77.085 107.192 0.097 0.097
## 4 7 15.387 78.584 102.669 0.098 0.098
## 5 10 14.018 77.099 110.216 0.097 0.097
## 6 8 14.419 77.613 104.709 0.098 0.097
## 7 8 15.740 78.942 106.038 0.098 0.098
## 8 6 15.965 79.116 100.190 0.098 0.098
## 9 9 24.698 87.755 117.862 0.105 0.104
## 10 7 24.649 87.490 111.575 0.105 0.104
## 11 7 22.085 85.077 109.162 0.102 0.102
## 12 5 22.890 85.569 103.633 0.103 0.102
## 13 8 23.828 86.833 113.929 0.104 0.103
## 14 6 23.923 86.676 107.750 0.104 0.103
## 15 6 21.018 83.960 105.034 0.101 0.101
## 16 4 22.060 84.643 99.696 0.102 0.102
## 17 10 12.912 75.942 109.059 0.097 0.096
## 18 8 11.138 74.260 101.356 0.095 0.095
## 19 8 13.381 76.560 103.656 0.097 0.097
## 20 6 13.753 76.946 98.020 0.097 0.097
## 21 10 14.018 77.099 110.216 0.097 0.097
## 22 7 16.003 79.193 103.278 0.098 0.098
## 23 8 15.740 78.942 106.038 0.098 0.098
## 24 5 18.610 81.575 99.639 0.100 0.100
## 25 8 23.817 86.823 113.919 0.105 0.103
## 26 6 22.750 85.585 106.660 0.103 0.102
## 27 6 21.381 84.302 105.377 0.102 0.102
## 28 4 21.339 83.982 99.035 0.101 0.101
## 29 8 23.828 86.833 113.929 0.104 0.103
## 30 5 33.393 94.942 113.005 0.110 0.109
## 31 6 21.018 83.960 105.034 0.101 0.101
## 32 3 41.450 101.025 113.068 0.113 0.113
## 33 11 11.000 73.779 109.906 0.095 0.095
## 34 9 11.366 74.417 104.524 0.095 0.095
## 35 8 172.997 186.633 213.729 0.203 0.201
## 36 6 171.988 184.087 205.162 0.201 0.197
## 37 10 14.018 77.099 110.216 0.097 0.097
## 38 8 14.419 77.613 104.709 0.098 0.097
## 39 7 171.155 184.710 208.795 0.198 0.198
## 40 5 170.529 182.349 200.413 0.196 0.195round(tail(abc_table,30),3)## model_df Cp AIC BIC PRESS GCV
## 99 7 202.114 199.115 223.200 0.221 0.218
## 100 5 207.977 199.428 217.492 0.222 0.219
## 101 9 49.490 109.788 139.894 0.121 0.121
## 102 7 49.086 108.730 132.816 0.120 0.120
## 103 6 237.031 212.659 233.733 0.238 0.239
## 104 4 241.652 212.068 227.121 0.238 0.238
## 105 8 53.017 112.287 139.383 0.123 0.122
## 106 6 53.465 111.813 132.887 0.123 0.122
## 107 5 215.073 202.456 220.520 0.224 0.223
## 108 3 214.295 199.811 211.853 0.221 0.219
## 109 7 76.857 129.701 153.786 0.139 0.137
## 110 5 76.563 128.290 146.354 0.137 0.136
## 111 4 241.823 212.135 227.188 0.237 0.238
## 112 2 238.045 208.221 217.253 0.232 0.232
## 113 9 39.082 100.929 131.036 0.114 0.114
## 114 7 46.447 106.576 130.661 0.118 0.118
## 115 6 201.432 197.699 218.773 0.220 0.216
## 116 4 211.458 199.772 214.825 0.220 0.219
## 117 9 49.490 109.788 139.894 0.121 0.121
## 118 6 47.131 106.767 127.842 0.119 0.118
## 119 6 237.031 212.659 233.733 0.238 0.239
## 120 3 298.436 231.452 243.494 0.270 0.270
## 121 7 52.157 111.199 135.284 0.122 0.122
## 122 5 76.080 127.955 146.018 0.137 0.136
## 123 4 214.243 200.949 216.003 0.222 0.221
## 124 2 212.710 197.985 207.016 0.217 0.216
## 125 7 76.857 129.701 153.786 0.139 0.137
## 126 4 76.058 127.321 142.374 0.136 0.135
## 127 4 241.823 212.135 227.188 0.237 0.238
## 128 1 1012.354 372.080 378.101 0.690 0.690AIC 와 BIC 측면에서 가장 좋은 모형?
AIC: sl ~ 1 + sw + pl + pw + sp + sw:pw + pl:sp + pw:sp
BIC: sl ~ 1 + sw + pl + sp + sw:pw
iris 자료에 대하여, 위 모형들이 최적 모형이라고 말할 수 있는 것은 아님
2차 교호작용까지를 포함한 모형들만 고려하였고,
step 으로 찾은 결과 내에서 최적결과를 찾은 것이고,
변수들을 함수적으로 변환하는 방법은 고려하지 않았음
aic <- abc_table$AIC
bic <- abc_table$BIC
aic[order(aic)[1:4]]## [1] 73.77855 73.77855 74.26028 74.26028model_set[order(aic)[1:4]]## [1] "sl ~ 1 + sw + pl + pw + sp + sw:pw + pl:sp + pw:sp"
## [2] "sl ~ 1 + sw + pw + sp + sw:pw + pl:sp + pw:sp"
## [3] "sl ~ 1 + sw + sp + sw:pw + pl:sp"
## [4] "sl ~ 1 + sw + pl + sp + sw:pw + pl:sp"best_aic_model <- eval(bquote(lm( .(model_set[order(aic)[1]]), iris)))
AIC(best_aic_model)## [1] 73.77855best_aic_model##
## Call:
## lm(formula = "sl ~ 1 + sw + pl + pw + sp + sw:pw + pl:sp + pw:sp",
## data = iris)
##
## Coefficients:
## (Intercept) sw pl pw
## 2.0896 0.7338 0.2325 1.0514
## spversicolor spvirginica sw:pw pl:spversicolor
## -1.0467 -2.6815 -0.2318 0.6599
## pl:spvirginica pw:spversicolor pw:spvirginica
## 0.7197 -1.1120 -0.4994bic[order(bic)[1:4]]## [1] 98.02025 99.03505 99.63873 99.69590model_set[order(bic)[1:4]]## [1] "sl ~ 1 + sw + pl + sp + sw:pw" "sl ~ 1 + sw + pl + sw:pw"
## [3] "sl ~ 1 + sw + pl + sp" "sl ~ 1 + sw + pl + pw"best_bic_model <- eval(bquote(lm( .(model_set[order(bic)[1]]), iris)))
BIC(best_bic_model)## [1] 98.02025best_bic_model##
## Call:
## lm(formula = "sl ~ 1 + sw + pl + sp + sw:pw", data = iris)
##
## Coefficients:
## (Intercept) sw pl spversicolor spvirginica
## 1.8064 0.6009 0.8401 -0.7275 -1.0381
## sw:pw
## -0.1041library( ggplot2 )
library( dplyr,warn.conflicts = FALSE, quiet=TRUE )
pp <- ggplot( filter(abc_table, AIC<150 ) ) + geom_point(aes(x=AIC, y=BIC))
pp +
geom_point(aes(x=AIC[order(aic)[1]], y=BIC[order(aic)[1]]), col='red',size=2)+
geom_point(aes(x=AIC[order(bic)[1]], y=BIC[order(bic)[1]]), col='green',size=2) 
AIC 와 BIC
AIC 는 예측능력이 높은 모형을 좋은 모형으로 선정
BIC 는 자료 설명력이 높은 모형을 좋은 모형으로 선정
BIC 는 AIC 에 비하여 단순한 모형을 좋은 모형으로 선정함
모형의 예측 능력 평가를 위한, PRESS 와 CV:
PRESS: Prediction Error Sum of Square
CV: Cross Validation
동일한 값에 대한 다른 이름
\(CV = (1/n) \sum_j (y_j - \hat{y}_{j, (-j)})^2\)
여기서, \(~\hat{y}_{j, (-j)}~\) 는 deleted fit (아래 쪽에서 설명 함)
\((y_j - \hat{y}_{j, (-j)}) = e_j /(1-h_{jj})~\) 이므로
\(PRESS= (1/n) \sum_j e_j^2 /(1-h_{jj})^2, ~~\)(식 11)
CV 의 식을 보면 여러번 회귀분석을 해야할 것처럼 보이지만,
PRESS 에서의 (식 11) 을 이용하면 한 번의 회귀분석으로 계산 가능함
모형의 적합기준이 아닌, 모형선택 기준으로 사용되어야 함
CV 의 개념은 보다 일반적인 경우로 확장 가능함
일반적인 CV
training 자료: 관측자료의 일부
testing 자료: training 에 사용되지 않은 관측자료의 일부
전체 관측자료에서 training 자료를 random 하게 뽑는 과정을 반복
여기서의 CV
traing 자료: n-1 개의 관측자료
test 자료: 1 개의 관측자료
전체 관측자료에서 돌아가며 하나씩을 test 자료로 사용
CV : Jackknife, Bootstrap 과 같은 재표본법을 이용한 모형 평가 방법
Bagging : 재표본법을 이용한 predictor 구성 방법
GCV: Generalized CV
PRESS (CV) 의 계산을 간단하게 근사한 것
\(GCV= (1/n) \sum_j e_j^2 /(1-\bar{h})^2\)
\(\bar{h} = (1/n) \sum_j h_{jj} = (1/n) trace(H)\)
\((1-\bar{h}) = (1/n) trace(I-H)\)
effective degree of freedom: 모형- \(~trace(H),~~\) 잔차-\(~trace(I-H)\)
- (식 11) 을 이용한 PRESS 계산과 CV 방법의 동일성:
res <- lm(sl~sw+pl,iris)
PRESS(res)## [1] 0.1130959CV(res)## [1] 0.1130959PRESS(res2)## [1] 0.09503959CV(res2)## [1] 0.09503959library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2abc_table %>% select(BIC, AIC, PRESS, GCV) %>% filter(AIC < 150 ) %>% ggpairs()
AIC 와 PRESS, GCV 는 그 순서가 일치하는데 반하여,
BIC 는 이와는 다른 특성을 보이고 있다
CV 와 deleted fit
최소자승법에 의하여 얻은 잔차의 성질:
관측된 자료 \(\mathbf{y}=\boldsymbol{\mu}+\boldsymbol{\epsilon}\) 을 \(\mathbf{y}=\hat{\boldsymbol{\mu}}+\mathbf{e}\) 이 되는 \(\hat{\boldsymbol{\mu}}\) 와 \(\mathbf{e}\) 로 분리
최소자승법은 \(\| \mathbf{e}\|^2\) 를 최소화 하는 방법
\(\mathbf{e} = (I-P) \boldsymbol{\epsilon}\), \(~~e_j = (1-h_{jj}) \epsilon_j +\cdots\),
\(\boldsymbol{\epsilon}~ \sim~ N_n (\mathbf{0}, \sigma^2 I)\), \(~~\mathbf{e}~ \sim~ N_n ( \mathbf{0} , \sigma^2 (I-P) )\),
\(Var(\epsilon_j ) = \sigma^2\), \(~~Var(e_j ) = \sigma^2 (1-h_{ii})\)
Deleted slope, Deleted fit, & Deleted variance
표현법: 음수첨자와 괄호 음수첨자
음수첨자: 행렬과 벡터에서 \(j\)-번째 행을 제거 한 것, 예) \(~X_{-j}\), \(\mathbf{y}_{-j}\)
괄호 음수첨자: \(~X_{-j}\), \(\mathbf{y}_{-j}\) 로부터 얻은 통계량, 예) \(~\hat{\boldsymbol{\beta}}_{(-j)}\)
Deleted slope:\(~~\hat{\boldsymbol{\beta}}_{(-j)} =(X_{-j}' X_{-j})^{-1} X_{-j} ' \mathbf{y}_{-j}\)
Deleted fit: \(~~\hat{\mathbf{y}}_{(-j)} =X \hat{\boldsymbol{\beta}}_{(-j)}\), \(~~\hat{\mathbf{y}}_{-j, (-j)} =X_{-j} \hat{\boldsymbol{\beta}}_{(-j)}\)
Deleted variance: \(~~ \hat{\tau}_j^2 = SSE_{(-j)}/(n-k-1), ~\) \(~(n-k-1)=(n-p-2)~\)
Deleted slope 계산 방법:
\(\hat{\boldsymbol{\beta}}- \hat{\boldsymbol{\beta}}_{(-j)}= (X'X)^{-1} \mathbf{x}_j \cdot e_j /(1-h_{jj}),~~\) (식7) \(~~\) (why?)
편의상, \(j=n\) 인 경우를 예로들어 살펴보자
\((X_{-n}' X_{-n})^{-1} \mathbf{x}_n = (X'X)^{-1} \mathbf{x}_n /(1-h_{nn}),~~\) (식8)
\(X'X' = X_{-n}' X_{-n} + \mathbf{x}_n \mathbf{x}_n '\) 이므로,
\((X_{-n}' X_{-n}) (X'X)^{-1} \mathbf{x}_n = (X'X - \mathbf{x}_n \mathbf{x}_n ' ) (X'X)^{-1} \mathbf{x}_n= (1-h_{nn}) \mathbf{x}_n\)
양변에 \((X_{-n}' X_{-n})^{-1}\) 을 곱하고, \((1-h_{nn})\) 으로 나누면 (식8)을 얻는다
\(\hat{\boldsymbol{\beta}}_{(-n)} = \hat{\boldsymbol{\beta}} - (X_{-n}' X_{-n})^{-1}\mathbf{x}_n (y_n -\hat{y}_n),~~\) (식9)
\((X_{-n}' X_{-n}) \hat{\boldsymbol{\beta}}_{(-n)} = X_{-n}' \mathbf{y}_{-n}\)
\((X_{-n}' X_{-n} + \mathbf{x}_n \mathbf{x}_n ' ) \hat{\boldsymbol{\beta}} = X_{-n}' \mathbf{y}_{-n} + \mathbf{x}_n y_n = (X_{-n}' X_{-n}) \hat{\boldsymbol{\beta}}_{(-n)}+ \mathbf{x}_n y_n\)
\(\hat{y}_n =\mathbf{x}_n ' \hat{\boldsymbol{\beta}}\) 이므로, 양변에 \((X_{-n}' X_{-n})^{-1}\) 을 곱하고 이항하여 정리하면, (식9)를 얻는다
(식8)을 (식9)에 대입하여 (식7) 을 얻는다
Deleted variance 계산 방법:
\(SSE_{(-j)} = SSE - e_j^2 /(1-h_{jj}), ~~\) (식10) \(~~\) (why?)
\(SSE_{(-j)}=\sum_{l \neq j} \left( y_l - \mathbf{x}_l ' \hat{\boldsymbol{\beta}}_{(-j)} \right)^2\)
\(y_l - \mathbf{x}_l ' \hat{\boldsymbol{\beta}}_{(-j)} = y_l - \mathbf{x}_l ' \hat{\boldsymbol{\beta}} + e_j h_{jl}/(1-h_{jj}) = e_l + e_j h_{jl}/(1-h_{jj})\)
\(y_j - \mathbf{x}_j ' \hat{\boldsymbol{\beta}}_{(-j)} = e_j /(1-h_{jj})\)
\(H \mathbf{y} = H \hat{\mathbf{y}}\) 이므로,\(\sum_l e_l h_{jl}=0\) 이고, \(H^2 =H\) 이므로, \(\sum_l h_{jl}^2 = h_{jj}\)
\(SSE_{(-j)} = \sum_l \left( e_l + e_j h_{jl}/(1-h_{jj}) \right)^2 - e_j^2 / ( 1-h_{jj} )^2\)
\(\qquad \quad ~~ = \sum_l e_l^2 + h_{jj} e_j^2 /(1-h_{jj})^2 - e_j^2 / ( 1-h_{jj} )^2\)
\(\qquad \quad ~~ = SSE - e_j^2 /(1-h_{jj})\)
2. 회귀진단
회귀분석 결과를 보고, 회귀모형에 대한 가정에 잘못이 없었는지 확인하는 과정
모든 통계분석 과정에서 꼭 거쳐야 할 단계
문제가 있는 경우 해결(치료)를 위한 전단계로써의 탐색 과정
최근에는 전형적인 해결 방법이 잘 알려져 있어, 진단의 중요성은 상대적으로 감소
이 강의에서는 회귀분석 방법의 이해와 관련된 중요 주제만 학습함
회귀진단의 8가지 검토사항
정규성 가정
등분산성
이상치
독립성
독립변수에 대한 선형성
다중공선성
지렛점
과대영향점
이전 강의자료
[가정과 진단]
[회귀진단의 주요 주제]
회귀진단
진단의 재료:
잔차
추정식/추정값
회귀계수
진단의 방법:
그림을 이용한 방법
검정을 이용한 방법
재표본법을 이용한 검정
잔차의 표준화:
standardized residuals: \(~~z_j = e_j /\hat{\sigma}\)
studentized residuals: \(~~t_j = e_j /(\hat{\sigma} \sqrt{1-h_{jj} })\)
externally standardized residuals: \(~~z_j^* = e_j /\hat{\tau}_j\)
externally studentized residuals: \(~~t_j^* = e_j /(\hat{\tau}_j \sqrt{1-h_{jj} })\)
Deleted variance: \(~~ \hat{\tau}_j^2 = SSE_{(-j)}/(n-k-1), ~\) \(~(n-k-1)=(n-p-2)~\)
PRESS residuals: \(~~r_j = y_j -\mathbf{x}_j' \hat{\boldsymbol{\beta}}_{(-j)}\)
PRESS =\(\sum_j r_j^2\)
\(r_j =e_j /(1-h_{jj})\), \(~~Var(r_j)= \sigma^2/(1-h_{jj})~~\) (why?)
정규성, 등분산성, 이상치, 독립성 문제
이전 강의자료
[Nonconstant Variance]
이전 강의자료
[Durbin-Watson]
함수의 형태에 대한 문제
plots
Resid vs. Predictor \(~~\mathbf{e}(\mathbf{y}|X) ~~vs. ~~\mathbf{x}_j\)
Partial regression plot (added variable plot) \(~~\mathbf{e}(\mathbf{y}|X_{j} ) ~~vs. ~~\mathbf{e}(\mathbf{x}_j|X_{j} )\)
Component-plus-residual plots (partial residual plot) \(~~\mathbf{e}(\mathbf{y}|X )+ \hat{\beta}_j \mathbf{x}_j ~~~vs. ~~\mathbf{x}_j\)
Augmented partial residual plots \(~~\mathbf{e}(\mathbf{y}|X, \mathbf{x}_j^2 ) + \hat{\beta}_j \mathbf{x}_j + \hat{\beta}_{jj} \mathbf{x}_j^2 ~~~vs. ~~\mathbf{x}_j\)
Nonparametric Regression
local regression
smoothing spline
GAM
MARS
Regression Tree
Neural Net
SVM
trees 자료, 함수 형태 진단을 위한 plot 예 :
residual vs. predictor plot
partial regression plot
partail residual plot
augmented partial residual plot
names(trees)<-c('g','h','v')
trees0 <- trees; trees0$g<-0
head(trees0)## g h v
## 1 0 70 10.3
## 2 0 65 10.3
## 3 0 63 10.2
## 4 0 72 16.4
## 5 0 81 18.8
## 6 0 83 19.7rex <- lm(v~g+h,trees)
rex2 <- lm(v~g+I(g^2)+h,trees)
vh <- resid(lm(v~h,trees))
gh <- resid(lm(g~h,trees))
vhat <- predict(rex,new=trees0)
vhat2 <- predict(rex2,new=trees0)
mydf <- data.frame(
g=trees$g,
resid=resid(rex),
p_res=trees$g-vhat, # partial residual
ap_res=trees$g-vhat2 # augmented partial residual
)
trz <- mydf %>%
reshape2::melt(id.var='g', measure.var=c('resid','p_res','ap_res')) %>%
rbind(data.frame(g=gh,variable='p_reg', value=vh))
levels(trz$variable)<-c( 'residual vs predictor',
'partial_residual' ,
'aug_partial_residual',
'partial_regression')
ggplot(trz)+ geom_point(aes(x=g, y=value, col=variable)) +
facet_wrap(~variable, ncol=2, scales='free') +xlab('')+
labs(title='Plots for Diagnostics')+theme(legend.position='')
다중공선성 muticollinearity
의미: 설명변수들이 선형독립이 아니거나, 거의 아닌 경우
원인
독립변수 선정의 부주의/오류
수치적인 오차
문제점
회귀계수 추정이 (거의) 불가능해지거나
추정된 회귀계수를 신뢰할 수 없게 됨
그러나 예측값에는 영향이 없음
해결책
Ridge 회귀 혹은 penalized regression 수행
VIF 를 점검하고 높은 경우를 제외
예측모형인 경우 굳이 신경쓰지 않아도 됨
예: 극단적인 경우로 다음 예를 고려하자
의도한 모형: (자식의 키) = a + c (아버지의 키) + d (어머니의 키) + 오차
변수 대입 오류: (자식의 키) = a + c (아버지의 키) + d (아버지의 키)
다중공선성 발생, 모수 추정 불가 (정규방정식의 근이 다양)
실제 적용된 모형: (자식의 키) = a + b (아버지의 키)+ 오차
c+d = b 이기만 하면, 합리적 추정이 된 것이고
다중공선성이 발생하더라도, 예측치에는 영향을 미치지 않음
즉, 모형 \(~\mathbf{Y} =X \boldsymbol{\beta} +\boldsymbol{\epsilon}~\) 에서 \(~X~\) 가 full column rank 가 아닌 경우,
\(\hat{\boldsymbol{\beta}}~\)는 여러 가능한 값을 가질 수 있으나, \(~\hat{\mathbf{y}}~\)는 영향을 받지 않음
인공신경망과 같이 회귀계수에는 관심이 없는 예측 모형에서는,
독립변수로 동일한 변수를 사용하는 과포화 모형이 되더라도 문제가 되지 않음
VIF 의 계산
독립변수들 사이에 다중회귀를 해서 독립변수들 사이의 종속성을 지수화함
독립변수 중 하나, 예를 들어, \(~\mathbf{x}_1~\) 을 다른 독립변수들을 이용하여 다중회귀를 함
그 결과로 얻은 결정계수를 \(~R_1^2~\) 이라 하면, \(~(1-R_1^2)^{-1}~\) 을 \(~VIF_1~\) 이라 함.
독립변수 각각에 대하여 VIF 를 계산하여, 그 값 중에 높은 값이 있다면 다중공선성이 있다고 말함
VIF 값을 크다고 하는 기준이 얼마인지에 대하여는 정해진 바 없지만, 보통 10 이상을 기준으로 사용함
이전 강의자료
[VIF]
지렛점, 과대영향점, 이상점 (이상치)
개념의 구별
지렛점 leverage point : 독립변수들만으로 그 여부 확인할 수 있음
과대영향점 influential point : 독립변수와 종속변수가 모두 있어야 그 여부 확인 가능
이상점 outlier : 잔차로 판단 가능함
의미
지렛점: 개별 관측값이 회귀식 전체에 큰 영향력을 미칠 수 있는 관측점
과대영향점: 개별 관측값이 회귀식 전체에 영향을 미쳐 그 결과를 크게 바꾸는 경우
이상점: 잔차값이 매우 큰 경우
이전 강의자료
[Leverage]
지렛점, 과대영향점
과대영향점은 보통 지렛점 혹은 leverage 가 큰 경우에 나타난다
지렛점이라고 해서 모두 과대영향점이 되는 것은 아니다
지렛점은 leverage 를 보고 판별할 수 있다
leverage : 사영행렬의 대각원소들 \(~h_{jj}~\)
leverage 들의 합은 \(~r~\) 이다 ( why? )
leverage 들은 비음수이다 ( why? )
정해진 기준은 없지만 보통 \(~h_{jj} > 2 (r/n)~\) 이면 지렛점으로 판단한다.
이전 강의자료
[Leverage Point & Influentail Point]
과대영향점의 판별
과대영향점의 판단은, 그 관측점이 있는 경우와 없는 경우를 시험하여
회귀결과에서의 차이가 큰 지 아닌지를 보아 판별한다
판별척도
개별관측값이 예측값 \(~\hat{\mathbf{y}}~\) 에 미치는 영향에 대한 관점
COOk’s distance:
DFFITS
개별관측값이 예측값 \(~\hat{\boldsymbol{\beta}}~\) 에 미치는 영향에 대한 관점
DFBETAS
COVRATIO
Cook’s distance: \(~D_j ~\)
\(j~\) 번째 관측값을 제거한 이후 각 값들에 대한 예측벡터와 \(~\hat{\mathbf{y}}_{(-j)}~\) 와
전체 자료를 다 이용한 예측벡터 \(~\hat{\mathbf{y}}~\) 사이의 거리를 \(~k \, \hat{\sigma}^2~\) 으로 나눈 값
즉, \(~D_j = Q_j /(k \, \hat{\sigma}^2)),~~\) 여기서 \(~Q_j = (\hat{\mathbf{y}}- \hat{\mathbf{y}}_{(-j)})'(\hat{\mathbf{y}}- \hat{\mathbf{y}}_{(-j)})\)
\(Q_j = (\hat{\mathbf{y}}- \hat{\mathbf{y}}_{(-j)})'(\hat{\mathbf{y}}- \hat{\mathbf{y}}_{(-j)}) = \sum_l^n (\hat{y}_l -\hat{y}_{l,(-j)})^2\)
\((\hat{y}_l -\hat{y}_{l,(-j)}) =(y_l -\hat{y}_{l,(-j)})- (y_l -\hat{y}_{l}),~~\) 이고
\((y_l -\hat{y}_{l,(-j)})= y_l - \mathbf{x}_l ' \hat{\boldsymbol{\beta}}_{(-j)}= e_l + e_j h_{jl}/(1-h_{jj}) ~~\) 이므로,
\(Q_j = e_j^2 h_{jj}/(1-h_{jj})^2, ~~\) (why?)
즉, \(~D_j = (1/k) \, t_j^2 \, h_{jj}/(1-h_{jj}),~~\) 여기서 \(~t_j~\) 는 studentized residual
명확한 기준이 있는 것은 아니나, \(~D_j >1~\) 이면,
\(j~\) 번째 관측점이 influential point 로 판정하기도 한다
influential point 를 판정하는 절대적 기준은 아니고 참고적으로 사용한다
DFFITS
\(DFFITS_j = (\hat{y}_j - \hat{y}_{j,(-j)} )/(\hat{\tau}_j \sqrt{h_{jj}})\)
\((\hat{y}_j -\hat{y}_{j,(-j)}) = e_j h_{jj}/(1-h_{jj}) ~~\) 이므로,
\(DFFITS_j = t_j^* \, (h_{jj} /(1-h_{jj}))^{1/2}, ~~\) 여기서 \(~t_j^*~\) 는 ext. studnt. residual
\(|DFFITS_j | > 2 \sqrt{k/n}~\) 이면 influential point 로 보기도 하고,
\(|DFFITS_j | > 2~\) 인 경우를 확실한 influential point 로 보지만, 참고적인 판단 기준임
DFBETAS
\(DFBETAS_{k,j} = (\hat{\beta}_k - \hat{\beta}_{k,(-j)} )/(\hat{\tau}_j \sqrt{c_{kk}}), ~~\) 여기서 \(~\{c_{kl}\}= (X'X)^{-1}\)
\(| DFBETAS_{k,j} | > 2 /\sqrt{n}~~\) 이면, influential point 일 가능성이 있다고 봄
COVRATIO
\(Cov(\hat{\boldsymbol{\beta}}) = \sigma^2 (X'X)^{-1}\)
COVRATIO 는, 관측점 \(~j~\) 를 제거하기 전후의 행렬식 \(~| \sigma^2 (X'X)^{-1} |~\) 의 변화 비율
즉, \(~\mbox{COVRATIO} = | \hat{\tau}_j^2 (X_{-j}'X_{-j})^{-1} |/ \hat{\sigma}^2 (X'X)^{-1}|~\)
\(|\mbox{COVRATIO}-1 | > 3 \sqrt{k/n}~~\) 이면, influential point 일 가능성이 있다고 봄
과대영향점의 처리
그 관측값이 정확하게 얻어지고 기록된 것인지 검토할 필요가 있음
과대영향점 자료라고 하여, 그 자료를 제거할 이유는 없음
과대영향점 자료의 경우, 바른 모형에 대한 더욱 많은 정보를 암시할 수 있음
과대영향점 자료에 영향을 덜 받는 통계적 방법의 적용을 고려할 필요가 있음
예) Bagging predictor (배깅방법에 대하여는 이후에 별도로 언급함)
과대영향점 자료와 bagging 의 효과
[bagging & influential point]
- trees 자료, 회귀진단 함수의 적용
names(trees)<-c('g','h','v')
rex <-lm(v~g,trees)
cbind(
cooks_D = cooks.distance(rex),
COvratio = covratio(rex),
dfbetas(rex)
)## cooks_D COvratio (Intercept) g
## 1 1.098362e-01 1.0752524 0.4491311407 -0.402632292
## 2 4.924417e-02 1.1309468 0.2931502713 -0.260665055
## 3 2.223564e-02 1.1579623 0.1939601381 -0.171433663
## 4 4.160457e-05 1.1384345 0.0073239827 -0.005960195
## 5 3.971050e-03 1.1234178 0.0697847174 -0.055880467
## 6 6.044249e-03 1.1152450 0.0849099405 -0.067386189
## 7 1.528743e-02 1.0830615 -0.1313710071 0.102194280
## 8 5.099320e-04 1.1268851 -0.0237564672 0.018480296
## 9 1.602747e-02 1.0776660 0.1320984607 -0.101613721
## 10 1.583849e-05 1.1249202 0.0040237058 -0.003057700
## 11 2.080153e-02 1.0561809 0.1448584520 -0.108631917
## 12 4.922094e-05 1.1217008 0.0067760165 -0.005008391
## 13 4.656599e-04 1.1202953 0.0208482702 -0.015409686
## 14 1.278988e-03 1.1131523 -0.0318138501 0.022304086
## 15 2.524848e-02 1.0172266 -0.1300100598 0.084691630
## 16 3.718828e-02 0.9468992 -0.0933107180 0.031261544
## 17 2.809168e-02 0.9853522 0.0802962791 -0.026901365
## 18 8.762270e-03 1.0686563 -0.0276494752 -0.002194424
## 19 4.450995e-02 0.9183785 -0.0256864567 -0.044487915
## 20 6.408063e-02 0.8430924 -0.0196272183 -0.066274006
## 21 2.754799e-04 1.1094629 -0.0002281443 0.005457506
## 22 1.137428e-02 1.0648496 0.0104904246 -0.044132322
## 23 5.014444e-05 1.1143703 0.0015313385 -0.003697407
## 24 6.088904e-02 0.9882962 0.1700264620 -0.236377261
## 25 1.847520e-02 1.1033207 0.0997734793 -0.134044426
## 26 6.459356e-02 1.0696731 -0.2304097889 0.287654563
## 27 5.008387e-02 1.1052970 -0.2071149431 0.255972776
## 28 7.597544e-02 1.0950295 -0.2687867964 0.326475494
## 29 2.844384e-02 1.1650876 0.1637461789 -0.198127619
## 30 3.976321e-02 1.1499744 0.1942385377 -0.235022394
## 31 8.880581e-01 0.8244374 -1.2037539861 1.370062034standardized 잔차와 studentized 잔차의 비교:
studentized 잔차의 경우,
상대적으로 큰 값은 더 커지고 작은 값은 더 작아지는 경향이 있다.
ggplot()+geom_point(aes(x=rstandard(rex),y=rstudent(rex)))+
geom_abline(slope=1,intercept=0,col='red')
leverage_influence_plot <- function(a_df){
ggplot(a_df) +
# geom_point(aes(x=a_df[,1], y=a_df[,2], size=3),col='white')+
geom_point(aes(x=a_df[,1], y=a_df[,2], col=high, size=value ) ) +
# scale_color_manual(values = c('FALSE' = 'skyblue4', 'TRUE' = 'orange')) +
geom_point(aes(x=a_df[,1], y=a_df[,2], size=value+1.5) , shape=1 )+
geom_smooth(aes(x=a_df[,1], y=a_df[,2]), method='lm', se=FALSE, col='orange') +
labs(x=names(a_df)[1], y=names(a_df)[2] ) +
facet_wrap(~variable,ncol=2, scales='free') # +
# guides(color=FALSE, size=FALSE)
}
lv = hatvalues(rex); h_lv = ( lv > 0.1)
df = dffits(rex); h_df = ( abs(df) > 0.4 )
trz <- data.frame(trees, lv, h_lv, df , h_df ) %>%
reshape2::melt(id.var=c('g','v'), measure.var=c('lv', 'df')) %>%
cbind(high=c(h_lv, h_df) )
levels(trz$variable) =c('leverage', 'influence(dffits)' )
# trz$variable = relevel(trz$variable,ref='influence(dffits)')
trz <- trz %>% arrange(high) # FALSE first
leverage_influence_plot(trz)## `geom_smooth()` using formula 'y ~ x'
- setosa 자료에서의 지렛점(leverage)과 과대영향점(dffit)
setosa <- iris[1:50,1:4]
names(setosa) <- c('sl','sw','pl','pw')
rout <- lm(sl ~sw, setosa)
infres<- influence.measures(rout)
infres$infmat %>% head()## dfb.1_ dfb.sw dffit cov.r cook.d hat
## 1 -0.002173923 0.005092115 0.02702353 1.063531 0.0003726311 0.02073628
## 2 0.145834973 -0.133858780 0.17802279 1.063380 0.0159601031 0.04601750
## 3 -0.064141657 0.054674293 -0.10529357 1.054625 0.0056142730 0.02738325
## 4 -0.107390546 0.096008164 -0.14588482 1.054864 0.0107346814 0.03528008
## 5 0.026090323 -0.034481211 -0.08275397 1.056426 0.0034765395 0.02420180
## 6 -0.048123423 0.053025730 0.06774178 1.095868 0.0023390989 0.05164186lv <-infres$infmat[,'hat']
cd <- infres$infmat[,'cook.d']
# sort(lv)
# sort(cd)
h_lv = (lv > 0.08 )
h_cd = (cd > 0.2 )
trz <- data.frame(sl=setosa$sl, sw=setosa$sw, lv, h_lv, cd , h_cd ) %>%
reshape2::melt(id.var=c('sl','sw'), measure.var=c('lv', 'cd')) %>%
cbind(high=c(h_lv, h_cd) )
levels(trz$variable) =c('leverage', 'influence(Cooks D)' )
leverage_influence_plot(trz)## `geom_smooth()` using formula 'y ~ x'
실습 과제
setosa 자료에서, 종속변수를 pl 이라 할 때,
step 함수로 좋은 모형을 찾고, 그 모형의 모든 하위 모형들에 대하여
AIC 와 BIC 값을 조사하시오
mtcars 자료에서, mpg 변수를 종속변수로 하는 가장 좋은 모형을 찾으시오
먼저 어떤 모형을 좋은 모형이라고 할 지에 대하여 설정하시오
setosa 자료에서, 종속변수를 pl 이라 할 때,
회귀진단에서의 검토하여야 할 8 가지 사항을 검토하시오
setosa 자료에서, 종속변수를 pl 이라 할 때,
회귀진단 방법으로 지렛점이 있는지, 과대영향점이 있는지 검토하시오
#mpg
mpg %>% group_by(class) %>% summarise(n=n(),hwy=mean(hwy))## # A tibble: 7 × 3
## class n hwy
## <chr> <int> <dbl>
## 1 2seater 5 24.8
## 2 compact 47 28.3
## 3 midsize 41 27.3
## 4 minivan 11 22.4
## 5 pickup 33 16.9
## 6 subcompact 35 28.1
## 7 suv 62 18.1