Summary : 다음을 R을 통하여 실습합니다. (1) Simple linear regression

Simple linear regression

단순 선형 회귀분석 - 1 : 시작하기

두 변수사이의 관계를 대략적으로 파악하기 위하여는 산점도 (scatter plot)를 이용하면 편리합니다. 그러나 두 변수 사이의 관계에 대하여 좀 더 정밀한 추론을 하려면, 추론의 대상을 명확히 하고 적절한 모형을 전제로 하여 개발된 통계적 방법을 필요로합니다.

두 변수 사이의 관계식을 파악하여 한 변수의 값으로부터 다른 변수의 값에 대한 예측을 필요로 하는 경우를 생각해봅시다. 이러한 경우에 관심있는 연속형 출력변수 (y)와 입력변수(X)의 함수 관계에 대한 통계적 분석을 하는 방법을 회귀분석 (regression analysis)라고 합니다.

\[ Y=f(X) + \epsilon \] 여기서 \(\epsilon\)은 평균 0, 분산 \(\sigma^2\)인 오차항을 나타내는 random variable입니다.

두 변수 사이의 함수 관계로서 가장 기본적인 형태인 일차함수관계, 즉 모형화하는 함수 \(f(X)\)가 선형관계일 때 선형회귀모형이라고 합니다. 수학적으로는 자료사이의 거리가 최소화되는 직선을 찾아낼 수 있다는 사실에 기반합니다. 이번 실습에서는 ’cars’라는 자료를 이용해서 분석방법을 알아보겠습니다.

여담) 단순선형회귀분석은 아버지와 아들의 키에 대한 조사에서 갈톤(Francis Galton: 1822-1911)에 의하여 처음으로 시도되었으며, 회귀분석이란 명칭도 그 조사의 결론으로부터 비롯된 것으로 알려져 있습니다.

data(cars)
str(cars)
## 'data.frame':    50 obs. of  2 variables:
##  $ speed: num  4 4 7 7 8 9 10 10 10 11 ...
##  $ dist : num  2 10 4 22 16 10 18 26 34 17 ...
plot(cars)

산포도에서 속도와 제동거리 사이에는 양의 상관관계가 있는 것으로 보입니다. 선형회귀분석을 해봅시다. R에서는 lm()이라는 함수를 이용해서 쉽게 계산할 수 있습니다. lm()이라는 함수는 ‘formula’ 라는 object를 인자로 받습니다. ’formula’라는 object는 R에서 수식을 읽기 쉽게 만든 객체입니다. (dependent variable ~ independent variable) 형태로 사용하며 intercept는 기본값으로 들어갑니다. 혹시 더 자세한 내용을 알고 싶다면 이 링크를 봐주세요~

fit <- lm(dist ~ speed, data=cars)
fit
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

intercept와 speed에 대한 계수를 확인할 수 있습니다. 실제로 우리가 계산하는 방식과 같을까요?

X <- cars$speed
Y <- cars$dist
S_xx <- sum((X-mean(X))^2)
beta_1 <- sum((X-mean(X))*(Y-mean(Y)))/S_xx 
beta_0 <- mean(Y) - beta_1*mean(X)
cat('lm :', coef(fit))
## lm : -17.57909 3.932409
cat('lse :', c(beta_0, beta_1))
## lse : -17.57909 3.932409

회귀분석의 결과를 abline()에 넘겨주어 그래프로 나타내봅시다. (이 결과 값들은 어떻게 계산되었을까요?, 이것으로 충분할까요?)

plot(cars)
abline(a=-5, b=3.5, col='red')
abline(fit, col='blue')
legend("topleft", legend = c('Random','Least squares'), col=c('red', 'blue'), lwd=1)

단순 선형 회귀분석 - 2 : 자세히 들여보기

summary(fit)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

summary()라는 함수를 이용해서 ‘lm’ class object의 내부를 좀 더 보도록 합시다! standard error, t-value, p-value를 확인할 수 있습니다.

residual <- (Y - beta_0 - beta_1*X)
N <- length(residual)

sigma_squared <- sum(residual^2) / (N-2)
beta_1_se <- sqrt( sigma_squared / S_xx )
beta_0_se <- sqrt( sigma_squared*(1/N + mean(X)^2 / S_xx) )
cat('Srd. Error', c(beta_0_se, beta_1_se))
## Srd. Error 6.75844 0.4155128
cat('residual standard error', sqrt(sigma_squared))
## residual standard error 15.37959

일치하는 것을 확인할 수 있습니다. 그리고 t-test에 대한 것도 아래와 같이 계산할 수 있습니다. 여기서 사용되는 t-test에 대한 식은

\[ \frac{\hat{\beta}}{ \hat{\rm s.e.}(\hat{\beta})} \] 와 같습니다. (왜 그럴까요?)

t_statistics <- c(beta_0/beta_0_se, beta_1/beta_1_se)
2*(1-pt(abs(t_statistics), df=(N-2))) # p-value
## [1] 1.231882e-02 1.489919e-12

R-squared 값 계산하기 (simple linear regression에서 partition of sum of squares)

SST <- sum((Y-mean(Y))^2)
SSR <- sum((beta_0 + beta_1*X-mean(Y))^2)
SSE <- sum(residual^2)

R_squared <- SSR/SST
cat('partition of sum of squares')
## partition of sum of squares
cat('SSR+SSE:', (SSR+SSE), ', SST:',SST) # check !!
## SSR+SSE: 32538.98 , SST: 32538.98
cat('R_squared:', R_squared)
## R_squared: 0.6510794

추가: summary() 읽기

결정계수 (Multiple R-squared) 는 회귀직선이 얼마나 자료를 잘 설명하는지에 대한 척도입니다. 그러나 결정계수는 설명 변수가 늘어나면 그 값이 커지는 성질이 있으므로 이를 조정한 조정결정게수 (Adjusted R-squared)가 더 많이 사용됩니다. F-statistic은 모형에서 기울기가 유의한지에 대한 척도입니다. 앞의 p-value값과 비교해봅시다. (나오는 수치를 완벽히 이해할 수 있도록 합시다~~)

좀 더 음미 해봅시다! attributes() 함수를 사용하면 ‘lm’ class object가 무엇을 결과 값으로 가지고 있는지 확인할 수 있습니다.

coef(fit)
## (Intercept)       speed 
##  -17.579095    3.932409
vcov(fit)
##             (Intercept)      speed
## (Intercept)   45.676514 -2.6588234
## speed         -2.658823  0.1726509
cars$dist[1:5]
## [1]  2 10  4 22 16
fitted(fit)[1:5]
##         1         2         3         4         5 
## -1.849460 -1.849460  9.947766  9.947766 13.880175
residuals(fit)[1:5]
##         1         2         3         4         5 
##  3.849460 11.849460 -5.947766 12.052234  2.119825
confint(fit)
##                  2.5 %    97.5 %
## (Intercept) -31.167850 -3.990340
## speed         3.096964  4.767853
deviance(fit) # residual sum of sqaures
## [1] 11353.52

Visualization

oldpar <- par(mfrow=c(2,2))
plot(fit)

표를 보는 방법 !

  1. residuals vs fitted : 선형 회귀분석에서 잔차는 평균이 0 이고 분산이 일정하며 fitted value에 독립이므로 잔차가 random 으로 뿌려져 있어야 정상입니다.

  2. normal Q-Q : 잔차 (residual)를 이론값 (정규분포)에 가까운지 확인하는 plot입니다. 기울기가 1인 직선이 되어야 이상적입니다.

  3. scale-location : 이후에 자세히 하도록 합시다!

  4. residuals vs leverage : 이상치 (outlier)의 유무를 검사하는데 유용합니다. (이후에 자세히)

ggplot 코드는 다음과 같이 사용할 수 있습니다!

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.3.2
plt <- ggplot(cars, aes(speed, dist)) +
  geom_point() +
  stat_smooth(method='lm')

plt

What’s next?

이제 선형회귀분석은 끝 ?? 무엇을 더 할 수 있을까요 ??

Disclaimer: 서울대학교 통계학과 원중호 선생님의 통계자료분석 과정과 통계학 주교재인 일반통계학을 참고하였습니다.