[회귀분석 3주차] 교란변수 및 t검정
오늘도 아니나 다를까 매우 내용이 많고 복잡하다. 수리통계학 내용도 다수 나온다. 이제는 몇 시간을 투자하여 깃헙을 정리하지는 않기로 했으므로, 나중에 이해를 도울 수 있도록 흐름과 맥락 위주로 서술한다.
교란변수(Confounding Variables)
MLR에서 각각의 회귀계수에 대한 영향은, 나머지 X값을 바꾸지 않고 고정한 채, 해당 회귀계수에 곱해지는 X값만 바꿔야 정확히 확인할 수 있다. 이처럼 나머지 X값을 고정(adjusted)해야하는 회귀계수를 partial regression coefficients라 한다.
그 예제로 차량의 무게와 마력에 따른 제로백 시간에 대한 SLR, MLR이 주어져있다. 차량이 무거우면 제로백이 오래 걸리는 것이 당연하니 positive한 적합선이 나와야할 것 같지만, 실제 둘 간의 관계는 negative였다. 그런데 마력까지 고려하여 다중 회귀를 하면 다시 positive가 나온다. 무엇이 문제일까?
앞선 SLR은 accelation과 weight에 모두 영향을 주는 마력(horsepower)이라는 변수를 모델에 고려하지 않았기 때문에 발생한 문제이다(마치 no-intercept 모델을 남용하면 오류가 발생하듯이). 이렇듯 계산에 고려되지 않았지만, 나머지 두 변수에 모두 영향을 주는 변수를 교란 변수(Confounding Variable)라 한다.
39페이지의 내용을 보면 알 수 있듯이, 제로백과 무게에 대한 전체적인 관계는 negative이지만, 마력에 따라 labeling(군집을 나눔)을 해보니 각각의 군집 내에서는 positive한 관계를 갖음을 볼 수 있다. 이게 옳은 해석이다. 이 교란변수의 영향을 무시하고 SLR을 하니 오류가 발생한 것이다.
이렇듯 세부적으로 보면 변수 간의 경향을 정확히 찾을 수 있으나, 전체 데이터를 나열함으로써 경향이 오히려 사라져버리는 현상을 심슨의 역설(Simpson’s Paradox)이라 한다.
R의 tree 데이터 모델링
data("trees")
라는 코드를 통해 R에 빌트인된 trees 데이터셋을 불러오자. 이 데이터셋은 각 나무들의 직경(Girth)과 높이(Height) 및 부피(Volume)에 대한 데이터를 저장하고 있다. 여기서 직경에 대한 부피의 산점도에 회귀선을 그어보자.
높이와 부피에 대한 산점도에서는 비교적 선형회귀가 잘 먹혔다. 물론 높이가 증가할수록 부피에 대한 변동(variability, 즉 퍼져있는 분산의 정도)이 증가하고 있긴 하다.
그런데 직경과 부피에 대한 산점도에서는 SLR의 선형 회귀보다, 이차곡선이 더 경향을 잘 설명하는 것을 볼 수 있다. 당연하다. 부피는 직경의 제곱에 비례할 테니까. 때로는 이러한 non-linearity 모델이 더 경향을 잘 반영할 수 있는 것이다.
두 경우 (즉, 선형회귀에서 점점 변동이 커지거나, 이차곡선이 경향을 더 잘 표현하는 경우)에는 x와 y에 로그를 취해, 더 나은 선형 회귀로 변환할 수 있다. 단 각 변수가 양수(positive)여야 로그를 취할 수 있을 것이다.
이 데이터셋에서도 교란변수 비슷한 상황을 볼 수 있는데, 다행히 심슨의 역설 상황까지는 아니다. 높이에 대한 부피의 산점도는 positive인데, R의 cut 함수를 통해 diameter를 5개의 구간으로 나누어 adjust 하더라도, 역시 각 직경의 구간 내에서 높이~부피가 positive함을 확인할 수 있다.
자 그렇다면
\[\text{Timber Volume}\approx \text{(constant)(Diameter)^2(Height)}\]라는 사실을 알기에, log transformation을 해주면
\[\log{\text{(Volume)}}=\beta_0+\beta_1\log\text{(Diameter)}+\beta_2\log\text{(Height)}+\epsilon\]라는 MLR 모델을 세울 수 있다. 이때 부피는 직경의 제곱에, 높이에 정비례하였으므로, \(\beta_1=2,\beta_2=1\)라는 것을 쉽게 추론할 수 있다. 그런데 이게 진짜 합당한 추리일까? 주어진 trees 데이터셋에서 이 값이 정말 합당한 값인지 구해보자. 먼저 코드로 두 intercept를 계산하면 다음과 같다.
data("trees")
trees$Diameter = trees$Girth
# 변수 Girth를 Diameter라고 rename하였다
lmtrees = lm(log(Volume) ~ log(Diameter) + log(Height), data=trees)
lmtrees$coef
## (Intercept) log(Diameter) log(Height)
## -6.631617 1.982650 1.117123
이 수치가 과연 합당한 값인지, LSE의 분산을 구함으로써 살펴보겠다.
\(\beta\)의 성질
\(\hat\beta\)는 consistent하다
먼저 LSE인 \(\hat\beta\)가 불편추정량(unbiased estimator)임을 보이자. r.v인 Y와 epsilon에 비해 X는 값이 관측된 fixed numbers(쉽게 말해, 상수 취급)이고 beta는 파라미터(모수이므로 상수)라는 것으로 보일 수 있다. 계산을 다 해보면
\[E[\hat\beta]=\beta\]라는 것을 쉽게 입증할 수 있다. 즉 LSE인 beta hat은 불편추정량, 다시 말해 consistent estimator(일치추정량)임을 알 수 있다. consistent하다는 것은 확률적으로 특정 값에 수렴한다는 뜻이다(수리통계학 내용).
X의 분산이 클수록 \(\hat\beta\)는 정확해진다
그렇다면 beta hat의 분산은 어떨까? beta hat은 회귀계수 추정량들을 모아놓은 벡터이므로, 그것의 분산은 행렬이 된다. 이러한 행렬을 Variance-Covariance Matrix(분산-공분산 행렬)이라 부른다. 주대각은 분산, 나머지는 공분산으로 이뤄진 행렬이다. 이때 이 행렬은 symmetric하다.
특히 SLR에서의 분산-공분산 행렬이 중요한데, 배부해주신 Algebra 참고자료와 회귀 노트의 내용을 종합하면
\[Var(\hat\beta)=\sigma^2(X^TX)^{-1}=\sigma^2\begin{align} \begin{bmatrix} \frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n(x_i-\bar{x})^2} & -\frac{\bar{x}^2}{\sum_{i=1}^n(x_i-\bar{x})^2} \\ -\frac{\bar{x}^2}{\sum_{i=1}^n(x_i-\bar{x})^2} & \frac{1}{\sum_{i=1}^n(x_i-\bar{x})^2} \end{bmatrix} \end{align}\]과 같이 정리되어 SLR에서의 회귀계수와 공분산 LSE를 바로 구할 수 있게 된다. 여기서 slope를 관장하는 beta hat 1의 표준편차를 구해보면
\[SD(\hat\beta_1)=\sqrt{Var(\hat\beta_1)}=\frac{\sigma}{\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2}}=\frac{\sigma}{\sqrt{n-1}SD_x}\]이다. beta hat 1의 표준편차가 작을수록, 모수인 beta 1에 가까워질 것이므로, 표본수 n이 커지거나 X의 표준편차가 커지면 \(\hat\beta_1\)은 더 좋은 추정량이 된다는 결론을 낼 수 있다.
이때 표본수는 그렇다 치고, 왜 X의 표준편차가 커지면 좋을까? 예측 변수 X의 값들이 모여있으면 올바른 직선을 만들기 어렵고, 그에따라 기울기도 산정하기 어렵다. 그러나 값들이 적절히 퍼져있으면 모델을 만들기에 유리하기 때문에, LSE인 beta hat 1이 모수 beta 1에 가까워지는 것이다!
\(\hat\beta_1\)의 상관성
두 가지를 기억해두면 되는데, 첫째로 y bar과 beta hat 1은 무관하다. 직선이 y축의 방향으로 평행이동시킨다고 해서 기울기가 바뀌지는 않듯이, 두 수치는 uncorrelated하다.
둘째로 beta hat 1과 beta hat 0, 즉 slope와 intercpet는 음의 상관관계를 가진다. 이는 식을 유도해보면 쉽게 확인할 수 있는데, beta hat 1과 beta hat 0간의 공분산을 구해보면 마이너스가 붙는 것을 볼 수 있다. 상식적으로, 기울기가 너무 가파르면 그에 맞춰 회귀선의 y절편은 낮아지기 마련이다.
표준 오차(Standard Errors)
이제 standard error에 대해 알아보자. 표준 오차(s.e.)란 표준 편차를 계산하는 수식에 알지 못하는 모수 sigma^2 대신 계산된 값인 MSE를 집어넣은 것이다. LSE의 표준 오차를 수식으로 표현하면
\[s.e.(\hat\beta_1)=\sqrt{\frac{MSE}{\sum_{i=1}^n(x_i-\bar{x})^2}}\] \[s.e.(\hat\beta_0)=\sqrt{MSE(\frac{1}{n}+\frac{\bar{x}^2}{\sum_{i=1}^n(x_i-\bar{x})^2})}\]이 표준 오차는 R 프로그램으로 계산된다.
lmtrees = lm(log(Volume)~log(Diameter)+log(Height), data=trees)
summary(lmtrees)$coef
와 같이 summary()$coef를 통하여 확인할 수 있다. 이떄 t value와 p-values가 같이 계산된다. 이제 수리통계학에서 배운 t-test를 할 수 있다. 이때 t와 자유도는
\[t=\frac{\hat\beta_j-\beta_{0j}}{s.e.(\hat\beta_j)}, df=n-p-1\]로 계산되며, 좌우측 혹은 양측 검정을 실시할 수 있다. 예를 들어 양측 검정은
# 양측 검정 시 p-values 계산
2*pt(abs(t), df=n-p-1, lower.tail=F)
# 단측 검정 시 2를 곱하지 않는다
로 p값을 계산할 수 있다. 이렇게 계산된 p values가 95% 신뢰구간 기준 0.05보다 작아야 귀무가설을 기각할 수 있다. 이 경우 귀무가설이 기각된다면 모수와 통계량은 서로 다른 수치인 것이고, 기각되지 않는다면 우리가 통계량을 모수에 알맞게 추론하여 모델을 세운 것이라 생각할 수 있다.
이외에도 R코드로 beta hat의 신뢰구간을 추정할 수 있다. beta hat과 표준 오차는 summmary() 함수를 통해, 임계 t값(그냥 t값과 다름에 유의)은 qt() 함수를 통해 구할 수 있다.
# 임계 t values 계산
# 두 코드는 값이 같게 나옴
qt(alpha/2, df=n-p-1, lower.tail=FALSE)
qt(1-alpha/2, df=n-p-1)
혹은 더 간단하게 confint() 함수를 통해 신뢰구간을 구할 수 있다.
# 신뢰구간 계산
confint(lmtrees)
# confint(lmtrees, level=0.95)
EX: Fire Damage Data
보험 회사에서 소방서로부터 떨어진 거리에 대한 화재 손실액에 대해 SLR을 구성하려고 한다. 계산해보면 기울기(단위: 마일 당 천 달러)가 약 4.9193이 나오는데, 이것이 기울기가 4라고 볼 수 있을까? 계산해보면 t값이 2.3498로, 단측 p-test를 돌려보면 0.018이 나와 0.05보다 작으므로, 손실액은 마일 당 4천 달러를 상회함을 알 수 있다.
전체 코드는 다음과 같다.
fire = data.frame(
dist = c(0.7, 1.1, 1.8, 2.1, 2.3, 2.6, 3.0, 3.1, 3.4, 3.8, 4.3, 4.6, 4.8, 5.5, 6.1),
damage = c(14.1, 17.3, 17.8, 24.0, 23.1, 19.6, 22.3, 27.5, 26.2, 26.1, 31.3, 31.3, 36.4, 36.0, 43.2)
)
lmfire = lm(damage ~ dist, data=fire)
summary(lmfire)$coef
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.277929 1.4202778 7.236562 6.585564e-06
## dist 4.919331 0.3927477 12.525421 1.247800e-08
t1 = (4.9193 - 4) / 0.3927
t1
## [1] 2.340973
pt(t1, df=13, lower.tail=F)
## [1] 0.01791179
# p값이 0.05보다 작으므로, beta1(slope) > 4이다!
출처: 회귀분석(STAT342) ㅊㅅㅂ 교수님