4.1 Galton의 데이터 사용

  1. 먼저 UsingR 패키지를 설치하고 불러옵니다.
require(UsingR)
data(galton)
str(galton)
## 'data.frame':    928 obs. of  2 variables:
##  $ child : num  61.7 61.7 61.7 61.7 61.7 62.2 62.2 62.2 62.2 62.2 ...
##  $ parent: num  70.5 68.5 65.5 64.5 64 67.5 67.5 67.5 66.5 66.5 ...


  1. galton의 데이터는 928명의 부모의 키와 아이의 키에 대한 자료입니다. 이 자료에 포함되어 있는 부모의 키는 아빠의 키와 1.08 x 엄마의 키의 산술 평균 입니다. 이 자료들의 분포를 살펴보기 위해 화면을 둘로 나누고 히스토그램을 그려봅니다.
par(mfrow=c(1,2))
hist(galton$child, col="blue", breaks = 100)
hist(galton$parent, col="blue", breaks = 100)

par(mfrow=c(1,1)) # Default


  1. 자료의 상관관계는 다음과 같이 구합니다.
cor.test(galton$child, galton$parent)
## 
##  Pearson's product-moment correlation
## 
## data:  galton$child and galton$parent
## t = 15.711, df = 926, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4064067 0.5081153
## sample estimates:
##       cor 
## 0.4587624
  • 결과를 보면 상관계수 cor는 0.458이고 p-value는 0.05 미만입니다.

4.부모와 아이의 키를 표로 요약해봅시다.

xtabs(~child+parent, data=galton)
##       parent
## child  64 64.5 65.5 66.5 67.5 68.5 69.5 70.5 71.5 72.5 73
##   61.7  1    1    1    0    0    1    0    1    0    0  0
##   62.2  0    1    0    3    3    0    0    0    0    0  0
##   63.2  2    4    9    3    5    7    1    1    0    0  0
##   64.2  4    4    5    5   14   11   16    0    0    0  0
##   65.2  1    1    7    2   15   16    4    1    1    0  0
##   66.2  2    5   11   17   36   25   17    1    3    0  0
##   67.2  2    5   11   17   38   31   27    3    4    0  0
##   68.2  1    0    7   14   28   34   20   12    3    1  0
##   69.2  1    2    7   13   38   48   33   18    5    2  0
##   70.2  0    0    5    4   19   21   25   14   10    1  0
##   71.2  0    0    2    0   11   18   20    7    4    2  0
##   72.2  0    0    1    0    4    4   11    4    9    7  1
##   73.2  0    0    0    0    0    3    4    3    2    2  3
##   73.7  0    0    0    0    0    0    5    3    2    4  0


  1. 부모의 키와 자녀의 키 사이의 수학적 관계를 나타내는 공식은 회귀분석을 통하여 구할 수 있습니다. R에서 회귀분석은 선형모형(linear model)을 사용합니다. lm은 linear model의 약자로 사용법은 다음과 같습니다. lm(종속변수(결과) ~ 독립변수(원인), data=데이터) 즉, 부모의 키에 의해 자녀의 키가 결정되므로 아래와 같이 회귀공식을 구할 수 있습니다.
out = lm(child~parent, data=galton)
summary(out)
## 
## Call:
## lm(formula = child ~ parent, data = galton)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8050 -1.3661  0.0487  1.6339  5.9264 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 23.94153    2.81088   8.517   <2e-16 ***
## parent       0.64629    0.04114  15.711   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.239 on 926 degrees of freedom
## Multiple R-squared:  0.2105, Adjusted R-squared:  0.2096 
## F-statistic: 246.8 on 1 and 926 DF,  p-value: < 2.2e-16
  • 결과를 보면 y절편(Intercept)이 23.94이고 parent의 기울기는 0.65으로 나타납니다. 즉 y = 0.65x+23.94입니다.
plot(child~parent, data=galton)
abline(out,col='red')


4.2 단순 선형 회귀

  1. women 데이터는 30세부터 39세까지 미국 여성 15명의 키와 몸무게 데이터 입니다. 키는 인치, 몸무게는 파운드로 되어있습니다.
data(women)
str(women)
## 'data.frame':    15 obs. of  2 variables:
##  $ height: num  58 59 60 61 62 63 64 65 66 67 ...
##  $ weight: num  115 117 120 123 126 129 132 135 139 142 ...


  1. 선형모형을 사용하여 회귀분석을 해보면 다음과 같습니다.
fit = lm(weight~height, data=women)
summary(fit)
## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14
  • 회귀분석 결과, y절편(Intercept)는 -87.52, 기울기는 3.45 입니다. p-value는 매우 낮습니다. 즉 몸무게와 키는 다음의 관계가 성립합니다. weight = -87.52 + 3.45 x height
  • 결과의 마지막 세 줄을 보면, Residual standard error는 이 모형을 사용하여 키로부터 몸무게를 예측했을 때 평균 1.525 파운드의 오차가 생긴다는 뜻입니다. **Multiple R-*squared가** 0.991 이라는 것은 이 모형은 몸무게 분산의 99.1%를 설명해주는다는 뜻 입니다.
  1. 이를 그림으로 그리면 다음과 같습니다. 실제수치와 예측치가 가장 많이 차이를 보이는 경우는 키가 아주 클때와 작을때 입니다. 이 그림을 보면 회귀선이 한 번 구부러지면 훨씬 정확한 에측이 가능할 수 있을 것 같습니다.
plot(weight~height, data=women)
abline(fit, col='blue')


4.3 다항 회귀

  1. 보다 정확한 예측을 위해 이차방정식을 도입하여 회귀분석을 해보겠습니다.
fit2 = lm(weight~height+I(height^2), data=women)
summary(fit2)
## 
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50941 -0.29611 -0.00941  0.28615  0.59706 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
## height       -7.34832    0.77769  -9.449 6.58e-07 ***
## I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared:  0.9995, Adjusted R-squared:  0.9994 
## F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16
  • lm 함수를 호출할 때 I() 함수를 사용하면 괄호 안의 수식을 R의 보통 수식처럼 계산하라는 명령어 입니다.
  • 이차 방정식의 회귀 결과를 보면 height와 height^2의 회귀게수가 모두 유의미 합니다. 또한 Multiple R-squared 값이 0.999로 증가하였습니다.
  • weight = 261.9 - 7.35 x height + 0.0831 x height^2
  1. 이 회귀모형의 예측치를 실제 값과 그려보면 다음과 같습니다.
plot(weight~height,data=women)
lines(women$height, fitted(fit2), col="red")


  1. car 패키지의 scatterplot() 함수는 두 개의 변수 간의 관계를 보는 데 매우 유용합니다. 다음의 그래프를 그려봅시다.
require(car)
scatterplot(weight~height, data=women)

- scatterplot() 함수는 scatterplot과 함께 가로, 세로축의 변수의 boxplot을 그려주고 가장 잘 맞는 직선회귀선을 그려줍니다.

scatterplot(weight~height, data=women, pch=19, # pch=점을 속이 차있는 circle
            spread=FALSE, smoother.args=list(lty=2),
            main="Women Age 30~39",
            xlab="Height(inches)", ylab="Weight(lbs.)")


4.4 다중 회귀

  1. 다중 회귀는 단순회귀분석과는 달리 예측변수가 두 개 이상인 경우에 사용합니다. 다중회귀분석에서 예측변수의 수가 많을 경우 모형선택의 과정을 거치는데, 회귀모형에 포함시킬 예측변수들의 선택 기준은 다음과 같습니다.
  • 반응변수(종속변수)와 상관관계가 높아야 합니다.
  • 선택된 예측변수(독립변수)들 사이에는 상관관계가 낮아야 합니다.
  • 소수 정예의 원칙을 따릅니다.
  1. 두번째 원칙에 있는 예측변수들 간의 상관관계가 높다는 것은 중복된 정보를 모형에 넣어준다는 뜻으로 변수의 낭비일 뿐만 아니라 계산상으로도 문제를 일으켜 영향력있는 예측변수를 의미없게 만들 수도 있습니다(다중공선성, multicolinearity).

  2. 모형선택에는 가능한 모든 예측변수들의 조합으로부터 모형을 만들어 비교하는 all possible regressions, 가장 유의한 변수부터 덜 유의한 변수 순으로 하나씩 추가시키는 forward selection, 모든 변수를 넣고 가장 기여도가 낮은 변수부터 하나씩 제거해가는 backward elimination등이 사용됩니다.

  3. MASS 패키지의 birthwt 데이터를 사용해보겠습니다.

require(MASS)
str(birthwt)
## 'data.frame':    189 obs. of  10 variables:
##  $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : int  2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
  • 이 자료를 쓸 때 유의할 점은 race 변수는 엄마가 백인인 경우 1, 흑인인 경우 2, 그 외의 경우 3으로 코딩되어있는 범주형 변수라는 점입니다. 그냥 race를 쓰게되면 연속형 변수로 인식하므로 factor(race)로 사용하여야 합니다.
# age(엄마의 나이), lwt(마지막 생리기간의 엄마의 몸무게), race(인종), smoke(임신 중 흡연상태), ptl(과거 조산 횟수), ht(고혈압의 기왕력), ui(uterine irritability 존재 여부)
out=lm(bwt~age+lwt+factor(race)+smoke+ptl+ht+ui,data=birthwt)
anova(out)
## Analysis of Variance Table
## 
## Response: bwt
##               Df   Sum Sq Mean Sq F value    Pr(>F)    
## age            1   815483  815483  1.9380 0.1656022    
## lwt            1  2967339 2967339  7.0519 0.0086284 ** 
## factor(race)   2  4750632 2375316  5.6450 0.0041901 ** 
## smoke          1  6291918 6291918 14.9529 0.0001538 ***
## ptl            1   732501  732501  1.7408 0.1887130    
## ht             1  2852764 2852764  6.7796 0.0099900 ** 
## ui             1  5817995 5817995 13.8266 0.0002676 ***
## Residuals    180 75741025  420783                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


# coefplot 패키지를 이용하면 계수 플롯을 만들 수 있습니다.
# 신뢰구간이 0을 포함하지 않으면 통계적으로 유의하다고 볼 수 있습니다.
library(coefplot)
coefplot(out)


  1. 분산분석표에서 age와 ptl에 유의하지 않으므로 두 변수를 제거한 모형을 만들고 두 변수의 중요도를 평가하기 위해 F-test를 실시합니다. F-test는 두 모형의 설명력을 비교하여 첫번째 모형에서 제거된 변수들의 기여도를 평가한다고 생각하면 됩니다. anova(작은모형, 큰모형)
out2=lm(bwt~lwt+factor(race)+smoke+ht+ui,data=birthwt)
anova(out2,out)
## Analysis of Variance Table
## 
## Model 1: bwt ~ lwt + factor(race) + smoke + ht + ui
## Model 2: bwt ~ age + lwt + factor(race) + smoke + ptl + ht + ui
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1    182 75937505                           
## 2    180 75741025  2    196480 0.2335  0.792
  • F-test 결과 p-value가 0.792로 두 변수를 모두 제거할 수 있습니다.
anova(out2)
## Analysis of Variance Table
## 
## Response: bwt
##               Df   Sum Sq Mean Sq F value    Pr(>F)    
## lwt            1  3448639 3448639  8.2654 0.0045226 ** 
## factor(race)   2  5076610 2538305  6.0836 0.0027701 ** 
## smoke          1  6281818 6281818 15.0557 0.0001458 ***
## ht             1  2871867 2871867  6.8830 0.0094402 ** 
## ui             1  6353218 6353218 15.2268 0.0001341 ***
## Residuals    182 75937505  417239                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • 두 변수를 제거한 모든 변수가 유의합니다.
  • 최종모형을 살펴보면 R^2 = 0.2404로 모형이 반응변수의 24.04%를 설명합니다. 설명력이 너무 작다고 볼수도 있지만 통제 불가능한 변인들이 많은 연구에서 이 정도의 R^2는 결코 작다고 할 수 없습니다.

4.5 상호작용이 있는 다중 회귀

  1. 다중회귀에서 예측인자 사이에 상호작용이 있는 경우에 대해 살펴봅시다. 자동차의 연비에 관한 mtcars 데이터에서 자동차의 공차 중량과 마력이 연비에 미치는 영향에 대해 회귀분석을 해보겠습니다. 이 때 공창 중량과 마력을 예측인자에 넣고 공차 중량과 마력의 상호작용을 같이 넣어보겠습니다.
data(mtcars)
fit=lm(mpg~hp*wt, data=mtcars)
summary(fit)
## 
## Call:
## lm(formula = mpg ~ hp * wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0632 -1.6491 -0.7362  1.4211  4.5513 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.80842    3.60516  13.816 5.01e-14 ***
## hp          -0.12010    0.02470  -4.863 4.04e-05 ***
## wt          -8.21662    1.26971  -6.471 5.20e-07 ***
## hp:wt        0.02785    0.00742   3.753 0.000811 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared:  0.8848, Adjusted R-squared:  0.8724 
## F-statistic: 71.66 on 3 and 28 DF,  p-value: 2.981e-13
  • 여기에서 hp*wt로 표현한 것은 all possible interaction을 뜻합니다. \[A*B = A + B + A:B\]
  • Coefficients 열을 보면 hp, wt 뿐만 아니라 hp와 wt의 상호작용(interaction)도 유의한 것을 알 수 있습니다. 즉, 두 개의 예측변수들 사이에 유의한 상호작용이 있다는 말은 한 예측인자와 반응인자 사이의 관계가 다른 예측 인자의 값에 따라 달라진다는 뜻 입니다.
  • mpg = 49.8 - 0.12 x hp - 8.22 x wt + 0.027 x hp x wt
plot(mpg~hp, data=mtcars, main="Interaction of hp:wt")
curve(31.41-0.06*x,add=TRUE)
curve(23.37-0.03*x,add=TRUE,lty=2,col=2)
curve(15.33-0.003*x,add=TRUE,lty=3,col=3)
legend("topright",c("2.2","3.2","4.2"),title="wt",lty=1:3,col=1:3)

- 상호작용을 이해하기 위해서 wt에 여러값을 대입한 식을 보도록 합시다. wt의 평균이 3.2이고 표준편차가 1이므로 2.2와 4.2값을 각각 공식에 넣으면 mpg와 hp 사이의 관계식에서 회귀계수가 (0.06, 0.03, 0.003)으로 변화하는 것을 볼 수 있습니다.