[회귀분석] Simple Linear Regression with R

오늘은 R을 이용하여

 

simple linear regression을 알아보겠습니다.

 

데이터는 wages.Rdata를 사용했습니다.

wages.Rdata

0.01MB

 

 

데이터 불러오기

setwd("경로")
load("wages.Rdata")
attach(wages)

setwd를 통하여 자신의 경로를 설정한 다음

 

load를 통해 경로 안에 있는 파일을 불러오면 됩니다.

 

attach를 통하여 데이터를 불러옴으로써

 

data.frame에서 column을

 

wages$logwage가 아닌 logwage라고 쓸 수 있게 됩니다.

 

 

데이터 구조 파악하기

str(wages)

>>> 'data.frame':	2178 obs. of  2 variables:
     $ education: num  16.8 15 10 12.7 15 ...
     $ logwage  : num  2.85 2.45 1.56 2.1 2.49 ...

data.frame 형식으로 되어있고

 

2178개의 observation이 2개의 변수로 되어있습니다.

 

하나의 변수는 education인 교육이고

 

또 다른 변수는 logwage로 월급과 같은 개념이라고 생각하시면 됩니다.

 

따라서 저희는 교육에 따라 임금이 어떻게 달라지는지를 알아볼 계획입니다.

 

summary(wages)

>>>    education        logwage     
     Min.   : 9.00   Min.   :0.240  
     1st Qu.:12.00   1st Qu.:1.950  
     Median :12.00   Median :2.241  
     Mean   :12.74   Mean   :2.240  
     3rd Qu.:14.00   3rd Qu.:2.524  
     Max.   :20.00   Max.   :4.265

summary를 통하여 요약 통계량 값을 알아봤습니다.

 

앞으로 독립변수(설명변수)로 사용할 education의 최솟값이

 

9인 것을 알 수 있습니다.

 

이것에 대한 얘기는 뒤에서 다루도록 하겠습니다.

 

 

회귀분석 line 그리기

wages.lm = lm(logwage ~ education)
plot(education, logwage, pch=23, bg='red', cex=2, cex.lab=3)
abline(wages.lm, lwd=4, col='black')

lm함수를 이용하여 종속변수를 logwage

 

독립변수를 education으로 설정한 

 

simple linear regression 모형을 만들었고

 

그것을 시각화한 것입니다.

 

어느정도 데이터의 값이 우상향 한다고 볼 수 있죠.

 

그렇다면 단순선형 회귀식인

에서 베타0과 베타1의 추정값을 직접 구해볼까요?

 

 

모수 베타0 베타1 구하기

beta.1.hat = cov(education, logwage) / var(education)
beta.0.hat = mean(logwage) - beta.1.hat * mean(education)
print(c(beta.0.hat, beta.1.hat))

>>> [1] 1.23919433 0.07859951

formula를 통하여 베타0와 베타1의 값을 추정해봤습니다.

 

각각 1.24, 0.08 정도가 되는것을 알 수 있습니다.

 

따라서 회귀식은 logwage = 1.24 + 0.08 * education이 됩니다.

 

공식을 통하여 구하지 않고 R에 내장된 coef를 통하여 

 

맞게 구했는지 알아볼까요?

 

print(coef(wages.lm))

>>> (Intercept)   education 
     1.23919433  0.07859951

제대로 구한것을 알 수 있습니다.

 

coef는 coefficient의 약자인 계수를 뜻하고

 

intercept는 절편을 뜻합니다.

 

이제 베타1의 표준오차와 test statistics를 구해볼까요?

 

 

베타1의 표준오차와 Test statistics 구하기

sigma.hat = sqrt(sum(resid(wages.lm)^2) / wages.lm$df.resid) 
SE.beta.1.hat = (sigma.hat * sqrt(1 / sum((education - mean(education))^2)))
Tstat = beta.1.hat / SE.beta.1.hat
data.frame(beta.1.hat, SE.beta.1.hat, Tstat)

>>>  beta.1.hat SE.beta.1.hat    Tstat
    1 0.07859951   0.004262471 18.43989

이렇게 베타1의 추정값, 베타1의 표준오차, test statistics 각각

 

0.08, 0.004, 18.44가 나옵니다.

 

그리고 이 모형의 결정계수인 R^2값을 구해볼까요?

 

 

결정계수 R^2 구하기

cor(logwage, education)^2 # == 결정계수

>>> [1] 0.1351453

결정계수는 0.135정도인 것을 확인할 수 있습니다.

 

그럼 지금까지 구한 값들을 summary를 통해 

 

한 번에 보여드리겠습니다.

 

 

Summary를 통하여 한번에 보기

summary(wages.lm)

이렇게 지금까지 직접 구한 추정 값들을

 

summary함수를 통해 한번에 구할 수 있습니다.

 

 

Test 하기

H_0 : 베타1 = 0

H_A : 베타1 =/ 0

 

위의 귀무가설과 대립 가설을 검정하는 방법은 간단합니다.

 

유의수준 알파가 0.05라고 했을 때 위의 표에서

 

유의확률인 p-value가 <2e-16으로 되어있으므로

 

소수점 뒤의 16자리에서 처음 2가 나오는 값보다

 

작다는 것을 통해 알파인 0.05보다 작으므로 귀무가설을 기각하게 됩니다.

 

만약

 

H_0 : 베타1 = 0

H_A : 베타1 > 0

 

위와 같이 단측검정를 하게 되면 위의 <2e-16은 양측 검정일 때의 p-value이므로

 

2로 나눠준 <e-16이 유의확률인 p-value값이 됩니다.

 

따라서 위와 같은 가설검정에서도 귀무가설을 기각하게 됩니다.

 

즉, 베타1이 0보다 크므로 우상향 한다는 결과를 얻을 수 있게 됩니다.

 

사실 엄밀히 말하면 위의 데이터는 결정계수 값이 너무 낮아

 

유의미하다고 보기 힘들고

 

observation이 많아 표준오차값이 작게나와

 

어지간하면 귀무가설을 기각시킬 확률이 높지만

 

공부를 하기 위해 사용한 데이터임을 참고해주시기 바랍니다.

 

 

모수의 신뢰구간 구하기

L = beta.1.hat - qt(0.975, wages.lm$df.resid) * SE.beta.1.hat
U = beta.1.hat + qt(0.975, wages.lm$df.resid) * SE.beta.1.hat
data.frame(L, U)

>>>            L          U
    1 0.07024057 0.08695845

베타1의 신뢰구간이 다음과 같이 나옵니다.

 

R의 내장함수를 이용하여 빠르게 구해볼까요?

 

confint(wages.lm)

>>>                  2.5 %     97.5 %
    (Intercept) 1.13138690 1.34700175
    education   0.07024057 0.08695845

confident interval을 줄인 confint함수를 이용하여

 

95%의 신뢰구간을

 

베타0와 베타1에 대하여 각각 구해줬습니다.

 

마지막으로 회귀분석을 통해

 

education의 값에 따른 logwage를 예측해 볼까요?

 

 

Prediction

edu16 <- data.frame(education = 16)
predict(wages.lm, edu16, interval='confidence', level=0.95)

>>>        fit      lwr      upr
    1 2.496786 2.464661 2.528912

위에서 구한 것은 education의 값이 16일 때

 

logwage의 평균에 대한 구간입니다.

 

따라서 education의 값이 16일 때

 

logwage는 평균적으로 2.46 ~ 2.53 사이에 있다고 할 수 있습니다.

 

그렇다면 만약 education이 16일 때 

 

특정 어떤 한 사람의 logwage가 몇이냐고 물어보면 어떻게 해야 될까요?

 

이때는 특정 한 명을 추정해야 하므로 에러텀이 더해져서 구간의 폭이 늘어나게 됩니다.

 

predict(wages.lm, edu16, interval='predict', level=0.95)

>>>        fit      lwr      upr
    1 2.496786 1.704295 3.289278

이처럼 1.70 ~ 3.29로 폭이 늘어나는 것을 확인할 수 있죠.

 

마지막으로 centering에 대해 말씀드리겠습니다.

 

 

Centering

우리는 위에서 베타0이 아닌 베타1 위주로 살펴봤습니다.

 

Simple Linear Regression에서 원래도 우리가 궁금한 모수는

 

베타1이기도 하지만 위에서 말씀드린 대로

 

education의 최솟값이 9인 것을 확인할 수 있죠.

 

대한민국은 의무교육이 있어 교육이 0이 될 수가 없습니다.

 

그렇다면 education의 값이 0일 때의 절편 값인 베타0이 사실상 의미가 없다고 볼 수 있죠.

 

이럴 때 베타0에 대해서 알아보려면 사용하는 방법이 centering입니다.

 

edu.center <- education - mean(education)
cen.lm <- lm(logwage ~ edu.center)
summary(cen.lm)

이렇게 centering을 해도 베타1에 대해서는 변한 것이 없고

 

베타0에만 영향이 간 것을 확인할 수 있습니다.

 

만약 앞으로 베타0의 추정값에 대해 궁금하시다면

 

centering을 해줘야 하는지도 고려해봐야 합니다.