R활용 회귀분석_기본가정
백승훈
2021 11 9
1 회귀분석을 위한 기본가정의 위배
1.1 오차항의 비정규분포
최소제곱 추정량은 오차항이 비정규분포를 갖는 경우에도 BLUE가 된다. (Best Linear Unbiased Estimator) 그러나 오차항이 비정규분포인 경우 회귀모수에 대한 유의성 검정이나 신뢰구간 추정이 유효하지 못하다.
1.1.1 잔차의 정규성 탐지
1.1.1.1 잔차의 히스토그램에 의한 탐지
오차항이 정규분포한다면 잔차항의 히스토그램은 다음의 특징을 가진다.
- 좌우대칭
- 종모양
- 하나의 봉우리
\[\divideontimes \; 오차항: 모집단, \; 잔차항 : 표본\]
1.1.1.1.1 오차항의 정규성 히스토그램으로 확인하기
mtcars data 중 mpg(연비)와 wt(무게)를 회귀분석하여 무게와 연비간의 반응함수를 추정하고 이 결과를 바탕으로 오차항이 정규분포하는지 여부를 잔차의 히스토그램을 통해 살펴보자.
summary(mtcars)## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000model <- lm(mtcars$mpg ~ mtcars$wt)
hist(model$residuals)
히스토그램을 그려본 결과, 정규분포를 따르는 것으로 보이지는 않는다.
1.1.1.1.2 Q-Q plot으로 정규성 확인
점들이 직선에 모여있을수록 대상이 정규분포를 따른다.
qqnorm(model$residuals)
qqline(model$residuals)
QQ plot을 그려본 결과로도 정규성을 따르지 않는 것으로 판단된다.
1.1.1.2 자크-베라(Jarque-Bera)의 검정통계량
히스토그램은 시각적으로 오차항의 분포형태를 쉽게 알 수 있으나 정확하게 판단할 수 없으므로 오차항의 정규분포 검정이 필요하다.
왜도와 초과첨도의 추정량
\[ skew = \frac{\mu^3}{\sigma^3} \\ kurt = \frac{\mu^4}{\sigma^4}−3 \]
자크−베라의 검정통계량은 점근적으로 \(\chi^2\)−분포를따른다.
\[ LM = n(\frac{skew^2}{6}+\frac{kurt^2}{24}) \]
- 검정통계량 LM > 임계값 : 귀무가설기각 (오차항은 정규분포를 따르지 않는다)
- 검정통계량 LM < 임계값 : 귀무가설채택 (오차항은 정규분포를 따른다)
- 자크-베라 검정을 위한 R 함수 : jarque.bera.test (tseries 패키지)
- jaeque.bera.test(x)에서 인수x는 수치형 데이터 혹은 시계열 자료형태
#자크베라검정
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoojarque.bera.test(model$residuals) ##
## Jarque Bera Test
##
## data: model$residuals
## X-squared = 2.399, df = 2, p-value = 0.3013qchisq(0.05,df=2,lower.tail=F) # 5% 유의수준의 임계값## [1] 5.991465하지만 자크-베라검정 결과는 p값 > 0.05이므로 귀무가설을 유지하며 정규성을 따른다고 얘기하고 있다.
그 밖의 함수들
- shapiro-wilk 정규성검정, 기본 내장함수, shapiro.test()
- anderson-darling 정규성검정, nortest패키지 내장, ad.test()
- shapiro-wilk 정규성검정, 기본 내장함수, shapiro.test()
#샤피로검정
shapiro.test(model$residuals) ##
## Shapiro-Wilk normality test
##
## data: model$residuals
## W = 0.94508, p-value = 0.1044샤피로 검정 결과도 p > 0.05 이므로 정규성을 따른다고 결론짓고 있다.
#앤더슨검정
library(nortest)
ad.test(model$residuals) ##
## Anderson-Darling normality test
##
## data: model$residuals
## A = 0.46842, p-value = 0.233앤더슨 검정 결과도 p > 0.05 이므로 정규성을 따른다고 결론짓고 있다.
1.1.2 오차항이 비정규분포를 갖는 경우의 추정방법
1.1.2.1 최우추정량
오차항이 비정규분포이면서 분포를 알 때 효율적인 추정량.
1.1.2.1.1 최우추정량과 최소제곱추정량 비교
mtcars data 중 mpg(연비)와 wt(무게)를 회귀분석하여 무게와 연비간의 반응함수를 추정하고, 오차항이 다변량 t-분포를 따른다는 가정 하에 최우추정량을 이용하여 다음 모형을 추정해보고, 동일한 회귀모형을 최소제곱추정량을 이용하여 분석한 결과와 비교해보자.
y <- mtcars$mpg #종속변수(mpg)를 식별자 y에 할당
x <- mtcars$wt #독립변수(wt)를 식별자 x에 할당
ls.print(com.ls.res<-lsfit(x,y)) #최소제곱추정결과 출력## Residual Standard Error=3.0459
## R-Square=0.7528
## F-statistic (df=1, 30)=91.3753
## p-value=0
##
## Estimate Std.Err t-value Pr(>|t|)
## Intercept 37.2851 1.8776 19.8576 0
## X -5.3445 0.5591 -9.5590 0str(mtcars) # 데이터 갯수 확인## '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 ...com.fn<-function(b0,b1){ #대수우도함수 입력(자유도31)
-sum(dt(y-b0-b1*x,df=31,log=T))
}
library(bbmle)## 필요한 패키지를 로딩중입니다: stats4summary(mle2(com.fn,start=list(b0=-12,b1=0.03))) #최우추정 결과출력## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = com.fn, start = list(b0 = -12, b1 = 0.03))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## b0 36.51960 0.81018 45.076 < 2.2e-16 ***
## b1 -5.18931 0.24253 -21.396 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 275.1463위의 결과를 정리하면 다음과 같다.
| 구분 | 최소제곱추정량 | 최대우도추정량 |
|---|---|---|
| \(\beta_0\) | 37.2851 | 36.51960 |
| \(\beta_1\) | -5.3445 | -5.18931 |
1.1.2.2 로버스트 추정량
오차항이 두터운 꼬리분포를 가져 무한한 분산을 갖게 되면 최소제곱추정량이나 최우추정량 대신 로버스트 추정량을 사용한다. 오차항이 두터운 꼬리분포를 갖는 이유는 극단적인 이상치가 빈번히 관측되기 때문이다. 최소제곱추정량의 문제는 이상치에 높은 가중치를 주기 때문에 이상치의 영향을 많이 받는다는 것이다. 이런 경우에 최소제곱추정량은 더이상 불편추정량이 아니다. 이런 문제를 해결하기 위해 이상치에 작은 가중치를 부여하는 로버스트 추정량이 필요하다. 이를 위한 구체적인 방법은 최소절대편차추정량, 절사최소제곱추정량, M-추정량등이 있고 최근에는 코엔커와 바셋의 분위수회귀추정량이 있다.
- 로버스트성(robustness): 만약 어떤 추정량의 표본분포가 가정사항이 맞지 않는 것에 대한 영향을 심각하게 받지 않으면 그 추정량은 로버스트 하다고 한다.
\[\textbf{M-추정량} \\ e=y−X\beta\]
최소제곱합이 잔차의 제곱합을 최소화하는 방법이라면 M-추정량은 목적함수의 합을 최소화하는 추정량을 구하는 것으로 정의된다.
\[\sum \rho(e) = \sum \rho (y−X \hat{\beta_M})\]
1.1.2.2.1 로버스트 추정량 산출
mtcars data 중 mpg(연비)와 wt(무게)를 회귀분석하여 무게와 연비간의 회귀모형을 로버스트 M-추정량을 이용하여 추정해보자(목적함수는 Huber의 함수를 사용하도록 한다.)
library(MASS)
M.res<-rlm(mpg ~ wt ,data=mtcars,method="M",psi=psi.huber) #M-추정 실행
summary(M.res) #M-추정 결과 출력##
## Call: rlm(formula = mpg ~ wt, data = mtcars, psi = psi.huber, method = "M")
## Residuals:
## Min 1Q Median 3Q Max
## -4.29431 -2.06044 0.08458 1.73940 7.21287
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 36.7379 1.8730 19.6143
## wt -5.2503 0.5577 -9.4138
##
## Residual standard error: 3.056 on 30 degrees of freedompt(19.6143,df=30,lower.tail=F) #절편계수 추정치에 대한 P-값## [1] 5.813877e-19pt(-9.4138,df=30) #기울기계수 추정치에 대한 P-값## [1] 9.162802e-11p값이 모두 0.05보다 작으므로 유의미하다고 판단된다.
1.1.2.3 분위수회귀추정량
회귀모형의 종속변수에 대한 기대치나 중앙값을 추정하는 기존 방법과는 다르게 다양한 분위수 값을 추정하여 종속변수분포의 특성을 알 수 있게 해주는 강점을 가진다.
1.1.2.3.1 분위수회귀추정량 산출
mtcars data 중 mpg(연비)와 wt(무게)를 회귀분석하여 무게와 연비간의 회귀모형을 분위수회귀추정량을 이용하여 추정해보자(분석에 사용될 분위수는 0.1, 0.25, 0.5, 0.75, 0.9로 한다.)
library(quantreg)## 필요한 패키지를 로딩중입니다: SparseM##
## 다음의 패키지를 부착합니다: 'SparseM'## The following object is masked from 'package:base':
##
## backsolvevt <- c(0.1,0.25,0.5,0.75,0.9) #분위수 벡터 생성
quant.res <- rq (mpg ~ wt, data = mtcars, tau = vt) #분위수회귀분석 실행
quant.res #결과 반환## Call:
## rq(formula = mpg ~ wt, tau = vt, data = mtcars)
##
## Coefficients:
## tau= 0.10 tau= 0.25 tau= 0.50 tau= 0.75 tau= 0.90
## (Intercept) 37.509794 36.913333 34.232237 38.137152 44.391047
## wt -6.494845 -6.083333 -4.539474 -5.113782 -6.266786
##
## Degrees of freedom: 32 total; 30 residual1.1.2.3.2 분위수 회귀분석 시각화
mtcars의 mpg와 wt간 산포도를 그려보고, 각 분위수에 대한 회귀직선을 그래프에 표시해보자.
plot(x = mtcars$wt, y = mtcars$mpg, xlab = "Weight", ylab = "miles per gallon", main = "Scatterplot and Quantile Regression Fit")
for( i in 1:length(vt)){
vt <- c(0.1,0.25,0.5,0.75,0.9)
color <- c("blue","yellow","pink","grey","red")
abline(rq(mpg ~ wt, data = mtcars, tau = vt[i]), col = color[i])
} 
1.2 다중공선성
독립 변수들 간의 상관관계가 심각하게 존재하느냐의 문제이다.
1.2.1 다중공선성 탐지
1.2.1.1 상관계수, 결정계수, t-검정에 의한 판단
독립변수 간의 상관계수 추정값이 큰 경우 다중공선성의 존재 가능성을 의심할 수 있다.(상관계수와 행렬산점도를 그려서 파악가능) 또한 결정계수가 매우 커서 0.7이상 임에도 불구하고 t-검정을 수행했을 때 부분 회귀계수가 통계적 유의성이 없는 경우 다중공선성이 존재할 확률이 매우 높다.
1.2.1.1.1 다중공선성 파악
mtcars의 mpg와 disp, hp, drat, wt, qsec간의 회귀분석 모형을 최소제곱추정량을 이용하여 추정해보고 다중공선성 존재 여부를 살펴보자.
mpg_reg <- lm(mpg ~ disp+hp+drat+wt+qsec, data = mtcars)
summary(mpg_reg) ##
## Call:
## lm(formula = mpg ~ disp + hp + drat + wt + qsec, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5404 -1.6701 -0.4264 1.1320 5.4996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.53357 10.96423 1.508 0.14362
## disp 0.00872 0.01119 0.779 0.44281
## hp -0.02060 0.01528 -1.348 0.18936
## drat 2.01578 1.30946 1.539 0.13579
## wt -4.38546 1.24343 -3.527 0.00158 **
## qsec 0.64015 0.45934 1.394 0.17523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.558 on 26 degrees of freedom
## Multiple R-squared: 0.8489, Adjusted R-squared: 0.8199
## F-statistic: 29.22 on 5 and 26 DF, p-value: 6.892e-10cor(mtcars)## mpg cyl disp hp drat wt
## mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 0.68117191 -0.8676594
## cyl -0.8521620 1.0000000 0.9020329 0.8324475 -0.69993811 0.7824958
## disp -0.8475514 0.9020329 1.0000000 0.7909486 -0.71021393 0.8879799
## hp -0.7761684 0.8324475 0.7909486 1.0000000 -0.44875912 0.6587479
## drat 0.6811719 -0.6999381 -0.7102139 -0.4487591 1.00000000 -0.7124406
## wt -0.8676594 0.7824958 0.8879799 0.6587479 -0.71244065 1.0000000
## qsec 0.4186840 -0.5912421 -0.4336979 -0.7082234 0.09120476 -0.1747159
## vs 0.6640389 -0.8108118 -0.7104159 -0.7230967 0.44027846 -0.5549157
## am 0.5998324 -0.5226070 -0.5912270 -0.2432043 0.71271113 -0.6924953
## gear 0.4802848 -0.4926866 -0.5555692 -0.1257043 0.69961013 -0.5832870
## carb -0.5509251 0.5269883 0.3949769 0.7498125 -0.09078980 0.4276059
## qsec vs am gear carb
## mpg 0.41868403 0.6640389 0.59983243 0.4802848 -0.55092507
## cyl -0.59124207 -0.8108118 -0.52260705 -0.4926866 0.52698829
## disp -0.43369788 -0.7104159 -0.59122704 -0.5555692 0.39497686
## hp -0.70822339 -0.7230967 -0.24320426 -0.1257043 0.74981247
## drat 0.09120476 0.4402785 0.71271113 0.6996101 -0.09078980
## wt -0.17471588 -0.5549157 -0.69249526 -0.5832870 0.42760594
## qsec 1.00000000 0.7445354 -0.22986086 -0.2126822 -0.65624923
## vs 0.74453544 1.0000000 0.16834512 0.2060233 -0.56960714
## am -0.22986086 0.1683451 1.00000000 0.7940588 0.05753435
## gear -0.21268223 0.2060233 0.79405876 1.0000000 0.27407284
## carb -0.65624923 -0.5696071 0.05753435 0.2740728 1.000000001.2.1.2 분산팽창인자(VIF)
독립변수들의 상관행렬에 의하면 분산팽창인자는 계수의 완전한(다중공선성이 존재하지 않는 경우)분산에 대산 계수의 실제 분산 비율로 정의된다. 즉 분산팽창인자는
\[VIF(\hat{\beta_k})= \frac{1}{1−R^2_k}\]
이다. 여기서 \(R^2_k\)은 k번째 독립변수를 종속변수로, 나머지 모든 독립변수들을 독립변수로 하는 회귀모형의 결정계수이다.
\(R^2_k\)의 값이 1에 가까우면 분산팽창 인자의 값은 커지게 되고 다중공선성이 존재할 가능성이 매우 높은 것으로 판단할 수 있다. (일반적으로 VIF가 10 이상이면 다중공선성이 존재한다고 판단한다.)
mtcars의 분산팽창인자를 계산해보라.
library(car)## 필요한 패키지를 로딩중입니다: carDatampg_reg <- lm(mpg ~ disp+hp+drat+wt+qsec, data = mtcars)
vif(mpg_reg) ## disp hp drat wt qsec
## 9.110869 5.201833 2.322343 7.012686 3.191939고려했던 모든 독립변수들의 VIF값이 10보다 작으므로 다중공선성 문제는 없다고 판단할 수 있다.
1.2.2 다중공선성의 해결방안
다중공선성이 존재하는 경우 가장 쉽게 문제를 해결하는 방법은 상관계수가 높은 독립변수 중 하나 혹은 일부를 회귀모형에서 제거하는 것이다. 그러면 나머지 독립변수의 유의성을 높여줄 수 있으나 필요한 변수의 누락으로 인한 모형설정의 오류 때문에 새로운 문제가 발생할 수 있다. 즉, 독립변수의 제거로 인해 최소제곱추정량이 더 이상 최우수 선형불편 추정량(BLUE)이 되지 못할 수 있다. 따라서 경제이론에 근거하여 회귀모형에서 중요한 역할을 하는 변수라면 단순히 다중공선성의 문제로 변수를 제거하는 방법은 지양하는 것이 바람직하다. 이런 경우 다음 세 가지 방법을 고려해볼 수 있다.
- 관측값을 늘려 표본크기를 증가시킨다.
- 원자료에 대해 차분 혹은 로그변환을 하거나 명목변수는 실질변수를 사용한다.
- 사전정보를 이용하여 변수를 상관관계가 높은 다른 변수로 대체한다.
1.2.3 능선추정량
다중공선성의 문제를 해결하기 위한 통계적 기법으로 호얼과 케나드가 제안한 능선회귀방법이 있다.
\(\beta\)의 능선추정량은 다음과 같이 정의된다.
\[ \hat{\beta_R} = (X′X+ \lambda I)^{−1}X′y\]
mtcars의 데이터를 이용하여 mpg와 disp, hp, drat, wt, qsec간의 능선추정량을 이용하여 추정해보고, 최소제곱추정 결과와 비교해보라.
library(MASS)
select(lm.ridge(mpg ~ disp+hp+drat+wt+qsec,data = mtcars, lambda=seq(0.01,3,0.0001))) #람다의 추정량 확인## modified HKB estimator is 0.8427224
## modified L-W estimator is 0.657137
## smallest value of GCV at 2.2622ridge.res <- lm.ridge(mpg ~ disp+hp+drat+wt+qsec,data = mtcars, lambda=2.2622) #능선추정(GCV 기준 값 사용)
lm(mpg ~ disp+hp+drat+wt+qsec,data = mtcars) #최소제곱추정 결과 출력##
## Call:
## lm(formula = mpg ~ disp + hp + drat + wt + qsec, data = mtcars)
##
## Coefficients:
## (Intercept) disp hp drat wt qsec
## 16.53357 0.00872 -0.02060 2.01577 -4.38546 0.64015ridge.res## disp hp drat wt qsec
## 19.449163151 -0.002281332 -0.020611520 1.867272544 -3.062039294 0.4104965821.3 이분산
오차항 입실론의 분산이 모든 엑스에 관계없이 동일하다는 기본가정인 등분산이 위배되었을 때 이분산이 존재한다고 한다. 자기상관이 시계열 자료에서 흔히 발생하는 것과는 달리 이분산은 횡단면 자료에서 자주 발생하게된다.
\[Var(\epsilon^2) = E(\epsilon^2_i) = \sigma^2_i \]
이분산이 존재해도 최소제곱(OLS)추정량은 여전히 불편 추정량이고 일치 추정량이긴 하지만 최소분산을 갖지 못하여 소표본이나 대표본에서 더 이상 효율적 추정량이 아니다. 따라서 독립변수의 통계적 설명력이 떨어지고 추정결과를 이용하여 종속변수에 대한 예측을 할 경우 오류가 발생될 확률이 높아지게 된다. 또한 계수 추정값의 표준오차가 편의되어 t-검정통계량 및 신뢰구간등이 정확하게 계산되지 못한다.
1.3.1 이분산 탐지
1.3.1.1 \(\epsilon_i\)와 \(\hat{y_i}\)의 산포도를 이용
이분산의 존재여부를 판단할 수 있는가장 쉬 운방법은 \(\epsilon_i\)와 \(\hat{y_i}\)의 산포도를 그려보는 것이다. 등분산이 충족되면 잔차들 \(\epsilon_i\)는 일정한 범위 내에 동일하게 분포하게 되나 이분산이 존재하는 경우 잔차들의 퍼짐은 증가 혹은 감소 추세를 나타내게 된다. 물론 잔차들의 퍼짐이 증가하다 감소하거나 혹은 감소하다 증가할 때도 이분산이 존재하게 된다. 따라서 잔차의 변동을 나타내주는 \(\epsilon_i\)와 \(\hat{y_i}\)의 산포도를 그려보고 이 두 변수 사이에 함수관계가 성립된다고 보이면 이분산이 존재할 가능성이 높다고 판단할 수도 있다.
1.3.1.2 \(\epsilon_i\)와 \(\hat{y_i}\)의 산포도 시각화
mtcars의 mpg와 wt간 회귀 모형을 다음과 같이 설정한 후 분석하여 결과를 도출해보고, 모형에 이분산이 존재하는지 여부를 \(\epsilon_i\) 와 \(\hat{y_i}\)의 산포도를 통해 살펴보라. 산포도를 그려보고, 각 분위수에 대한 회귀직선을 그래프에 표시해보자.
model1 <- lm(mpg~wt,data=mtcars) #회귀분석 수행
summary(model1) #결과출력##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10residuals<-resid(model1) #잔차
yhat <- predict(model1,interval="none") #종속변수의 예측값
plot(x=yhat,y=residuals) #산포도 생성
그려본 결과, \(\epsilon_i\)와 \(\hat{y_i}\)간 상관관계가 존재한다고 보기 어려우므로 이분산이 아니라고 보인다.
1.3.1.3 Breausch-Pagan-Godfrey 검정통계량
\(\epsilon_i\)와 \(\hat{y_i}\)의 산포도를 그려보면 이분산의 존재 여부를 쉽게 판단할 수 있으나 판단의 통계적 정확성이 보장되지 못하는 한계가 있다. 따라서 이분산 탐지를 위한 여러 가지 검정통계량 중에서 브로쉬-파간-가드프레이 검정통계량을 설명하고자 한다. 회귀모형 \(y=X \beta+\epsilon\)을 추정하여 얻은 잔차의 제곱 \(\epsilon^2_i\)을 종속변수로, \(X\)를 독립변수(절편항 포함)로 한 선형회귀모형
\[\epsilon^2_i = b_0+b_1X_{1i}+ \cdots +b_K X_{Ki}+ν_i\]
을 다시 추정한다. 이 때 얻은 결정계서 \(R^2\)와 \(n\)의 곱에 의해 브로쉬-파간-가드프레이의 카이제곱 검정통계량(\(BPG\))을 계산할 수 있다.
\[BPG=nR^2\]
BPG 검정통계랑이 \(\chi^2_K\)분포의 임계값을 초과하면 등분산을 가정하는 귀무가설은 기각된다.
- 귀무가설 : 등분산이다
- 대립가설 : 이분산이다
1.3.1.4 \(BPG\) 검정통계량 산출
mtcars의 mpg와 wt간 회귀 분석 결과를 이용하여 BPG 검정통계량을 계산하고 그 결과를 해석해보라.
library(lmtest) #lmtest 패키지 실행## 필요한 패키지를 로딩중입니다: zoo##
## 다음의 패키지를 부착합니다: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(model1) #BPG 검정통계량 계산##
## studentized Breusch-Pagan test
##
## data: model1
## BP = 0.040438, df = 1, p-value = 0.8406\(p\)값 > 0.05 이므로 등분산이라고 판단된다.
1.3.2 일반최소제곱(GLS)추정량
이분산이 존재하는 경우 OLS추정량은 불편추정량이지만 가장 효율적인 선형추정량이 아니므로 BLUE가 아니다. 하지만 오차항의 표준편차를 이용한 가중치를 구하여 변형한 회귀모형에 최소제곱법을 적용함으로써 BLUE를 구할 수 있는데 이를 일반최소제곱(GLS)추정량이라고 한다.
\[\hat{\beta_G}=(X^{*'}X^*)^{−1}X^{*'}y^*= (X^{'} \Omega X)^{−1}X^{'} \Omega^{−1}y\]
이는 다음 속성을 만족하게 되어 BLUE이다.
\[E(\hat{\beta_G})= \beta \\ E(\hat{\beta_G}−\beta)(\hat{\beta_G}−\beta)^{'}=(X^{*'}X^*)^{−1}=(X^{'} \Omega X)^{−1}\]
mtcars의 mpg와 wt간 회귀모형을 GLS추정량을 이용하여 추정해보자. (GLS추정량을 이용하여 분석을 수행하는 과정에서 이분산 구조를 결정하는 함수는 고정된 가중치를 부여하는 varFixed()를 사용한다.)
library(nlme)
model_gls <- gls(model=mpg~wt,data=mtcars,weights=varFixed(~wt)) #GLS분석실행
model_gls## Generalized least squares fit by REML
## Model: mpg ~ wt
## Data: mtcars
## Log-restricted-likelihood: -80.13914
##
## Coefficients:
## (Intercept) wt
## 39.25858 -5.95787
##
## Variance function:
## Structure: fixed weights
## Formula: ~wt
## Degrees of freedom: 32 total; 30 residual
## Residual standard error: 1.7798931.4 자기상관
오차항이 서로 상관되지 않는다는(오차항이 서로독립) 가정이 위배되어
\[Cov(\epsilon_i,\epsilon_j) = E(\epsilon_i,\epsilon_j) \ne 0 \quad (i \ne j)\]
인 경우 자기상관(autocorrelation)이 서로 존재한다고 한다. 이런 자기상관은 시계열 자료가 가지는 관성 혹은 느린 조정과정, 거미줄 현상, 모형설정의 오류, 자기회귀 조건부 이분산(ARCH)효과 등에 의해서 발생하게 된다.
예를 들어, 종합주가지수(KOSPI) 수익률은 이자율 \(I_i\), 경재성장률 \(G_i\), 미국의 다우존스지수 수익률 \(D_i\)에 의해 결정된다고 가정하면 종합지수 수익률 회귀모형은
\[R_i= \beta_0 + \beta_1 I_i + \beta_2 G_i + \beta_3 D_i + \epsilon_i\]
에 의해 설정된다.여기서 오차항 \(\epsilon_i\)는 자기상관이 없다. 하지만 다우존스지수 수익률 변수를 누락하여 모형
\[R_i= \beta_0 + \beta_1 I_i + \beta_2 G_i + \nu_i \]
을 설정하게 되면 이 모형의 오차항 \(\nu_i = \beta_3 D_i + \epsilon_i\)은 \(D_i\)변수와 함수 관계를 갖게 되면서 \(D_i\)변수의 자기상관으로 인해 \(\nu_i\)에 자기상관이 발생하게 된다. 이와 같이 변수의 누락은 자기상관을 유발하게 된다. 자기상관이 존재하는 경우에도 이분산의 경우처럼 OLS추정량은 여전히 불편추정량이기는 하지만 최소분산을 갖지 못하여 소표본이나 대표본에서 더 이상 효율적 추정량이 아니다. 따라서 자기상관이 존재할 경우 독립변수의 통계적 설명력이 떨어지고 추정 결과를 이용하여 종속변수에 대한 예측을 할 경우 오류가 발생될 확률이 높아지게 된다. 또한 계수 추정값의 표준오차가 편의되어 t-검정통계량 및 신뢰구간 등이 정확하게 계산되지 못한다.
1.4.1 자기상관의 탐지
1.4.1.1 \(\epsilon_i\)와 \(i\)의 산포도를 이용
자기상관이 존재하는 경우에는 잔차 상호 간에 어떠한 함수관계가 존재하게 된다. 산점도가 구형인 경우 잔차들 간에 함수관계가 성립되지 않아 자기상관이 없다고 볼 수 있으나, 산점도에서 잔차들 간에 선형 혹은 비선형의 함수관계가 성립되는 경우에는 자기상관이 존재한다고 볼 수 있다.
1990년부터 2020년까지 GDP 증가율과 실업률에 대한 대한민국의 자료를 이용하여 다음과 같은 회귀모형을 설정하고
\[Unemp_i = \beta_0 + \beta_1 GDP_i + \epsilon_i\]
회귀분석을 수행해보자. 그리고 잔차의 산포도를 통해 위의 모형에 자기상관이 존재하는지 검토해보자.
library(readxl)
unemployment <- read_excel("unemployment.xlsx", col_types = c("numeric", "numeric", "numeric"))
model_unemploy <-lm(unemploy~GDP_GR,data=unemployment) #회귀분석 수행
summary(model_unemploy) #회귀분석 결과 출력##
## Call:
## lm(formula = unemploy ~ GDP_GR, data = unemployment)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3328 -0.4319 -0.0882 0.1485 3.7399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.07729 0.24772 16.459 2.99e-16 ***
## GDP_GR -0.18159 0.05822 -3.119 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9103 on 29 degrees of freedom
## Multiple R-squared: 0.2512, Adjusted R-squared: 0.2254
## F-statistic: 9.729 on 1 and 29 DF, p-value: 0.004077plot(resid(model_unemploy)) #잔차의 산포도 출력
잔차의 산포도 상으로 다소 자기 상관이 있어 보이나 판단하기 어려운 상황이다.
1.4.1.2 더빈-왓슨의 d 통계량
잔차의 산포도를 이용하는 경우보다 명확하게 자기상관의 존재 여부를 판단하기 위해 사용되는 일반적인 검정통계량으로 더빈-왓슨의 \(d\) 통계량이 있다.
\[\frac{\sum(e_i−e_{i−1})^2}{\sum e^2_i}\]
\(d\)가 0이나 4에 가까울수록 잔차들은 매우 높은 자기상관을 갖게 됨을 알 수 있다.
- 자기상관검정 가설
- 귀무가설:자기상관이 없다
- 대립가설:자기상관이 있다
- 검정통계량: \(\frac{\sum(e_i−e_{i−1})^2}{\sum e^2_i}\)
- 검정결과 판정
- 검정통계량이 \((4-d_L)\)보다 크거나 \(d_L\)보다 작으면 \(H_0\) 기각
- 검정통계량이 \(d_U\)보다 크고\((4-d_U)\)보다 작으면 \(H_0\) 채택
- 검정통계량이 \(d_L\)과 \(d_U\)사이 혹은 \((4-d_U)\)와 \((4-d_L)\)사이에 놓이면 미결정
위의 실업률 회귀분석 결과를 이용하여 더빈-왓슨의 \(d\) 검정통계량을 구해복 자기상관 검정을 수행해보자.
library(car)
dwt(model_unemploy) #더빈-왓슨 검정## lag Autocorrelation D-W Statistic p-value
## 1 0.3897787 1.21837 0.02
## Alternative hypothesis: rho != 0p 값 < 0.05이므로 자기상관이 있다고 보인다.
1.4.1.3 더빈의 \(h\) 통계량
독립변수가 오차항과 독립이라는 가정이 위배되어 더빈-왓슨의 \(d\) 통계량이 정확하지 못할 때, 자기상관의 존재에도 불구하고 자기상관이 없다는 귀무가설이 채택되는 문제가 발생할 수 있다. 이런 경우 표준정규분포를 이루는 더빈의 \(h\) 통계량을 사용한다.
\[ h = \hat{\rho} \sqrt{\frac{n}{1−n Var(\hat{\beta_1})}}\]
여기서 \(Var(\hat{\beta_1})\)은 \(\hat{\beta_1}\)의 OLS 추정값의 추정된 분산이고 \(n\)은 표본크기이다.
위의 실업률 회귀분석 결과를 이용하여 시차종속변수를 독립변수로 포함하는 새로운 회귀모형이 다음과 같을 때, 5% 유의수준에서 모형에 자기상관이 존재하는지 더빈의 \(h\)통계량을 이용하여 검정해보자.
\[Unemp_i = \beta_0 + \beta_1 Unemp_{t−1} + \beta_2 GDP_i + \epsilon_i\]
attach(unemployment)
unemploy1 <- c(NA,unemployment$unemploy[1:(length(unemployment$unemploy)-1)]) #시차변수생성
model2 <- lm(unemploy ~ unemploy1 + GDP_GR, data=unemployment) #회귀분석 실행
n <- length(resid(model2)) #관측값의 수
rho.hat <- 1-0.5*dwt(model2)[[2]][1] #상관계수 추정치 계산
var.b1.hat <- summary(model2)$coefficients[2,2]^2 #베타1의 분산 추정치 계산
h.stat <- rho.hat*sqrt(n/(1-n*var.b1.hat)) #더빈의 h통계량 계산
detach(unemployment)
h.stat;qnorm(0.05,0,1,lower.tail=F) #h통계량과 임계값.## [1] NaN## [1] 1.644854