다중회귀분석 multiple regression analysis (통계 R 초급 - 12)

다중회귀분석을 실행하기 위하여 R에서 제공하는 데이터를 아래와 같이 불러들인다.

> st77 = data.frame(state.x77)
> str(st77)
'data.frame':	50 obs. of  8 variables:
 $ Population: num  3615 365 2212 2110 21198 ...
 $ Income    : num  3624 6315 4530 3378 5114 ...
 $ Illiteracy: num  2.1 1.5 1.8 1.9 1.1 0.7 1.1 0.9 1.3 2 ...
 $ Life.Exp  : num  69 69.3 70.5 70.7 71.7 ...
 $ Murder    : num  15.1 11.3 7.8 10.1 10.3 6.8 3.1 6.2 10.7 13.9 ...
 $ HS.Grad   : num  41.3 66.7 58.1 39.9 62.6 63.9 56 54.6 52.6 40.6 ...
 $ Frost     : num  20 152 15 65 20 166 139 103 11 60 ...
 $ Area      : num  50708 566432 113417 51945 156361 ...

1977년도 미국의 50개 주에 대한 여러 데이터를 모아 놓은 것이다. 

Population은 인구, Income은 1인당 수입, Illiteracy는 문맹도, Life.Exp는 평균수명, HS.Grad는 고등교육졸업률, Frost는 추운정도, Area는 면적을 각각 나타낸다. 

데이터에 인구밀도를 계산해서 넣으려면 아래와 같이 하면 된다.

> st77$Density = st77$Population * 1000 / st77$Area
> summary(st77)
   Population        Income       Illiteracy       Life.Exp    
 Min.   :  365   Min.   :3098   Min.   :0.500   Min.   :67.96  
 1st Qu.: 1080   1st Qu.:3993   1st Qu.:0.625   1st Qu.:70.12  
 Median : 2838   Median :4519   Median :0.950   Median :70.67  
 Mean   : 4246   Mean   :4436   Mean   :1.170   Mean   :70.88  
 3rd Qu.: 4968   3rd Qu.:4814   3rd Qu.:1.575   3rd Qu.:71.89  
 Max.   :21198   Max.   :6315   Max.   :2.800   Max.   :73.60  
     Murder          HS.Grad          Frost             Area       
 Min.   : 1.400   Min.   :37.80   Min.   :  0.00   Min.   :  1049  
 1st Qu.: 4.350   1st Qu.:48.05   1st Qu.: 66.25   1st Qu.: 36985  
 Median : 6.850   Median :53.25   Median :114.50   Median : 54277  
 Mean   : 7.378   Mean   :53.11   Mean   :104.46   Mean   : 70736  
 3rd Qu.:10.675   3rd Qu.:59.15   3rd Qu.:139.75   3rd Qu.: 81163  
 Max.   :15.100   Max.   :67.30   Max.   :188.00   Max.   :566432  
    Density        
 Min.   :  0.6444  
 1st Qu.: 25.3352  
 Median : 73.0154  
 Mean   :149.2245  
 3rd Qu.:144.2828  
 Max.   :975.0033

그러면 Density라는 변수는 Population과 Area로부터 계산된 것이므로 redundant한 변수가 된다. 심할경우 회귀분석을 하는동안 역행렬을 구하지못하는 경우가 발생할때도 있다. 이럴때에는 서로 linear dependent한 변수를 제거해 주어야만 된다. 만약 임의로 이렇게 변수를 만들어 넣지 않았는데에도 역행렬을 구하지 못한다고 나온다면 linear dependent 변수를 제거하는 작업은 꽤 골치아플수도 있다. 

correlation matrix를 아래와 같이 구할수도 있고 모든 변수를 다 scatter plot을 그려서 한눈에 어떤 변수끼리 관계가 있는지 알아보는 것도 가능하다.

> cor(st77)
            Population     Income   Illiteracy    Life.Exp     Murder
Population  1.00000000  0.2082276  0.107622373 -0.06805195  0.3436428
Income      0.20822756  1.0000000 -0.437075186  0.34025534 -0.2300776
Illiteracy  0.10762237 -0.4370752  1.000000000 -0.58847793  0.7029752
Life.Exp   -0.06805195  0.3402553 -0.588477926  1.00000000 -0.7808458
Murder      0.34364275 -0.2300776  0.702975199 -0.78084575  1.0000000
HS.Grad    -0.09848975  0.6199323 -0.657188609  0.58221620 -0.4879710
Frost      -0.33215245  0.2262822 -0.671946968  0.26206801 -0.5388834
Area        0.02254384  0.3633154  0.077261132 -0.10733194  0.2283902
Density     0.24622789  0.3299683  0.009274348  0.09106176 -0.1850352
               HS.Grad        Frost        Area      Density
Population -0.09848975 -0.332152454  0.02254384  0.246227888
Income      0.61993232  0.226282179  0.36331544  0.329968277
Illiteracy -0.65718861 -0.671946968  0.07726113  0.009274348
Life.Exp    0.58221620  0.262068011 -0.10733194  0.091061763
Murder     -0.48797102 -0.538883437  0.22839021 -0.185035233
HS.Grad     1.00000000  0.366779702  0.33354187 -0.088367214
Frost       0.36677970  1.000000000  0.05922910  0.002276734
Area        0.33354187  0.059229102  1.00000000 -0.341388515
Density    -0.08836721  0.002276734 -0.34138851  1.000000000
> pairs(st77)

 

이제 dependent variable로 Life.Exp를 사용하고 나머지 모든 변수를 independent variable로 사용하여 선형회귀분석을 하여보자.

> model1 = lm(Life.Exp ~ Population + Income + Illiteracy + Murder + 
+                        HS.Grad + Frost + Area + Density, data = st77)
> summary(model1)

Call:
lm(formula = Life.Exp ~ Population + Income + Illiteracy + Murder + 
    HS.Grad + Frost + Area + Density, data = st77)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.47514 -0.45887 -0.06352  0.59362  1.21823 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.995e+01  1.843e+00  37.956  < 2e-16 ***
Population   6.480e-05  3.001e-05   2.159   0.0367 *  
Income       2.701e-04  3.087e-04   0.875   0.3867    
Illiteracy   3.029e-01  4.024e-01   0.753   0.4559    
Murder      -3.286e-01  4.941e-02  -6.652 5.12e-08 ***
HS.Grad      4.291e-02  2.332e-02   1.840   0.0730 .  
Frost       -4.580e-03  3.189e-03  -1.436   0.1585    
Area        -1.558e-06  1.914e-06  -0.814   0.4205    
Density     -1.105e-03  7.312e-04  -1.511   0.1385    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7337 on 41 degrees of freedom
Multiple R-squared:  0.7501,	Adjusted R-squared:  0.7013 
F-statistic: 15.38 on 8 and 41 DF,  p-value: 3.787e-10

R-squared가 0.7정도인데 독립변수로 쓰인 것들 중에 평균수명을 예측하는데 중요하지 않은 변수들이 보인다. 그중 하나가 Area로 ANOVA table에서 F테스트의 probability 0.96이나 된다. 반면 Income, Illiteracy, Murder, HS.Grad는 아주 significant한 변수로 간주할수있다. 

 

이것을 고려하여 Area를 모델에서 삭제하여 다시 만들수도 있다.

model2 = update(model1, .~. - Area)
summary(model2) 
summary.aov(model2)
anova(model1, model2)

두 모델 즉 Area를 predictor로 사용한 모델과 사용하지 않는 모델 사이의 ANOVA 결과를 볼수도 있는데 p-value가 약 0.42로 두 모델의 차이가 significant하지 않음을 보여준다. 

 

이렇게 하나씩 하나씩 비교해가면서 변수를 더하고 빼는 것은 참으로 번거로운 일이다. 

그런데 이것을 자동으로 해주는 것이 있는데 stepwise regression 기능이다. 

modelf = step(model1, direction = "backward") 
summary(modelf)
summary.aov(modelf)

Adjusted R-squared가 0.70이었지만 지금 0.71로 약간 증가하였다.

사용된 predictor는 Population, Murder, HS.Grad, Frost이다. 이들 변수의 coefficient는 모두 significant하다.