[R] Multiple linear regression (다중회귀분석)

[R] Multiple linear regression (다중회귀분석) Start

BioinformaticsAndMe

Multiple linear regression (다중회귀분석)

---

: 다중회귀분석은 예측변수(독립변수;설명변수)가 2개 이상인 회귀분석

: 다중회귀분석에서 예측변수 개수가 많을 경우, 적절한 회귀모형 선택이 필요함

: 회귀모형에 포함되는 예측변수의 선정 기준

ㄱ) 반응변수(종속변수)와 높은 상관관계

ㄴ) 선택된 예측변수들은 서로 낮은 상관관계를 보임 (다중공선성 문제를 회피)

ㄷ) 예측변수의 개수는 적을수록 유리

Feature selection (변수 선택법)

---

1) All possible regressions

- 변수들의 가능한 모든 조합들로부터 최적의 모형을 찾아냄

- 유의한 변수가 누락되지 않는 안전한 방법

- 변수가 많을수록 탐색 시간이 급증함

2) Forward stepwise selection (Forward selection)

- 기여도가 높은 유의한 변수부터 하나씩 추가하는 기법

- 빠른 계산이 장점

- 이미 선택된 변수는 다시 제거되지 않음

3) Backward stepwise selection (Backward elimination)

- 모든 변수를 포함한 상태에서 불필요한 변수를 제거해나가는 방법

- 중요한 변수가 제외될 가능성이 매우 적음

- 이미 제외된 변수는 다시 선택되지 않음

         

1. 실습데이터 로딩

---

# MASS 패키지의 birthwt 데이터
library(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 ...

2. 다중회귀분석 수행

---

■반응변수 - bwt(출생시 체중)

■예측변수 - age(엄마나이), lwt(엄마몸무게), race(인종), smoke(임신중흡연상태), ptl(과거조산횟수), ht(고혈압기왕력), ui(uterine irritability 존재여부) 

# race는 연속형 변수로 인식하므로, factor(race)로 사용해야함
model_1=lm(bwt ~ age+lwt+factor(race)+smoke+ptl+ht+ui, data=birthwt)
summary(model_1)
Call:
lm(formula = bwt ~ age + lwt + factor(race) + smoke + ptl + ht + 
    ui, data = birthwt)

Residuals:
    Min      1Q  Median      3Q     Max 
-1838.7  -454.5    57.6   465.1  1711.0 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2934.132    311.450   9.421  < 2e-16 ***
age             -4.093      9.440  -0.434 0.665091    
lwt              4.301      1.722   2.497 0.013416 *  
factor(race)2 -488.196    149.604  -3.263 0.001318 ** 
factor(race)3 -353.334    114.319  -3.091 0.002314 ** 
smoke         -351.314    106.180  -3.309 0.001132 ** 
ptl            -47.423    101.663  -0.466 0.641443    
ht            -586.836    200.841  -2.922 0.003925 ** 
ui            -514.937    138.483  -3.718 0.000268 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 648.7 on 180 degrees of freedom
Multiple R-squared:  0.2424,	Adjusted R-squared:  0.2087 
F-statistic: 7.197 on 8 and 180 DF,  p-value: 2.908e-08

3. 유의하지 않은 변수 제거

---

: 위 다중회귀분석 결과에서 age와 ptl 변수가 유의하지 않은 것으로 확인됨

: 아래의 3가지 방법 중 하나로, 유의하지 않은 변수를 제거할 수 있음

1) 유의하지 않은 변수를 누락하고 회귀모형을 새롭게 정의
model_2 = lm(bwt~ lwt+factor(race)+smoke+ht+ui, data=birthwt)

2) update 함수를 사용하여, 기존 회귀모형에서 유의하지 않은 변수 제거
model_2 = update(model_1, .~. -age-ptl)

3) step 함수를 사용하여, 기존 회귀모형에서 유의하지 않은 변수를 제거해나감
model_2 = step(model_1, direction = "backward")

summary(model_2)
Call:
lm(formula = bwt ~ lwt + factor(race) + smoke + ht + ui, data = birthwt)

Residuals:
     Min       1Q   Median       3Q      Max 
-1842.14  -433.19    67.09   459.21  1631.03 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2837.264    243.676  11.644  < 2e-16 ***
lwt              4.242      1.675   2.532 0.012198 *  
factor(race)2 -475.058    145.603  -3.263 0.001318 ** 
factor(race)3 -348.150    112.361  -3.099 0.002254 ** 
smoke         -356.321    103.444  -3.445 0.000710 ***
ht            -585.193    199.644  -2.931 0.003810 ** 
ui            -525.524    134.675  -3.902 0.000134 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 645.9 on 182 degrees of freedom
Multiple R-squared:  0.2404,	Adjusted R-squared:  0.2154 
F-statistic:   9.6 on 6 and 182 DF,  p-value: 3.601e-09

4. 회귀모델 비교 및 평가

---

: F-test 분산분석으로 두 회귀모형의 설명력을 비교하여, 첫 번째 회귀모형에서 제거된 변수들의 기여도를 평가함

# anova( 제거 전 모형, 제거 후 모형 )
anova(model_1, model_2)
Analysis of Variance Table

Model 1: bwt ~ age + lwt + factor(race) + smoke + ptl + ht + ui
Model 2: bwt ~ lwt + factor(race) + smoke + ht + ui
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    180 75741025                           
2    182 75937505 -2   -196480 0.2335  0.792

 : F-test 결과 p-value가 0.792로 매우 크므로, 앞서 제거된 두 변수는 회귀모형에 대한 기여도가 적음을 알 수 있음

# summary 함수로 최종 회귀모형을 평가
summary(model_2)
Call:
lm(formula = bwt ~ lwt + factor(race) + smoke + ht + ui, data = birthwt)

Residuals:
     Min       1Q   Median       3Q      Max 
-1842.14  -433.19    67.09   459.21  1631.03 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2837.264    243.676  11.644  < 2e-16 ***
lwt              4.242      1.675   2.532 0.012198 *  
factor(race)2 -475.058    145.603  -3.263 0.001318 ** 
factor(race)3 -348.150    112.361  -3.099 0.002254 ** 
smoke         -356.321    103.444  -3.445 0.000710 ***
ht            -585.193    199.644  -2.931 0.003810 ** 
ui            -525.524    134.675  -3.902 0.000134 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 645.9 on 182 degrees of freedom
Multiple R-squared:  0.2404,	Adjusted R-squared:  0.2154 
F-statistic:   9.6 on 6 and 182 DF,  p-value: 3.601e-09

 : R² 값이 0.2404 → 최종 회귀모형이 예측변수들의 24.04% 설명함

5. 다중공선성(Multicollinearity) 확인

---

: 다중공선성은 분산팽창지수(Variation Inflation Factor;VIF)라는 통계량을 사용하여 계산 가능

: 한 예측변수에 대해 VIF의 제곱근은 다중공선성의 정도를 나타냄

: VIF가 10을 넘으면(= GVIF^(1/(2*Df)) 값이 2를 넘으면), 다중공선성 문제를 보임

# VIF값은 car 패키지의 vif() 함수로 계산
install.packages('car')
library(car)
vif(model_2)
                 GVIF Df GVIF^(1/(2*Df))
lwt          1.182676  1        1.087509
factor(race) 1.257687  2        1.058993
smoke        1.154755  1        1.074595
ht           1.073548  1        1.036122
ui           1.036844  1        1.018255

: vif 함수로 완성된 회귀모형을 검사했을 때, 다중공선성 문제는 없는 것으로 확인됨

6. 변수의 상대적 중요성

---

: 현재까지 사용한 변수들 중에서, 결과 예측에 가장 중요한 변수를 서열화할 수 있음

# 아래 코드는 Dr. Johnson 논문에 소개된 스크립트로, 변수의 상대적 중요도를 출력함
relweights <- function(fit,...){
         R <- cor(fit$model)
         nvar <- ncol(R)
         rxx <- R[2:nvar, 2:nvar]
         rxy <- R[2:nvar, 1]
         svd <- eigen(rxx)
         evec <- svd$vectors
         ev <- svd$values
         delta <- diag(sqrt(ev))
         lambda <- evec %*% delta %*% t(evec)
         lambdasq <- lambda ^ 2
         beta <- solve(lambda) %*% rxy
         rsquare <- colSums(beta ^ 2)
         rawwgt <- lambdasq %*% beta ^ 2
         import <- (rawwgt / rsquare) * 100
         import <- as.data.frame(import)
         row.names(import) <- names(fit$model[2:nvar])
         names(import) <- "Weights"
         import <- import[order(import),1, drop=FALSE]
         dotchart(import$Weights, labels=row.names(import),
            xlab="% of R-Square", pch=19,
            main="Relative Importance of Predictor Variables",
            sub=paste("Total R-Square=", round(rsquare, digits=3)),
            ...)
         return(import)
}

# 변수의 상대적 중요도를 시각화
model_3 = lm(bwt~ lwt+race+smoke+ht+ui, data=birthwt)
result = relweights(model_3, col="blue")
result

# ggplot2을 사용하여, 변수의 상대적 중요도를 시각화
library(ggplot2)
plotRelWeights=function(fit){
  data<-relweights(fit)
  data$Predictors<-rownames(data)
  p<-ggplot(data=data,aes(x=reorder(Predictors,Weights),y=Weights,fill=Predictors))+ 
    geom_bar(stat="identity",width=0.5)+
    ggtitle("Relative Importance of Predictor Variables")+
    ylab(paste0("% of R-square \n(Total R-Square=",attr(data,"R-square"),")"))+
    geom_text(aes(y=Weights-0.1,label=paste(round(Weights,1),"%")),hjust=1)+
    guides(fill=FALSE)+
    coord_flip()
  p
}

plotRelWeights(model_3)

#Reference

1) http://rstudio-pubs-static.s3.amazonaws.com/189354_277dfb3a83a34a2abaae855b90fcf269.html

2) http://brokerstir.com/multiple-linear-regression-intuition/

3) https://quantifyinghealth.com/stepwise-selection/

4) https://slidesplayer.org/slide/11708770/

5) https://rpubs.com/cardiomoon/190997

6) http://leoslife.com/archives/3874

7) Johnson(2000, Mutilvariate BehavioralResearch, 35,1-19)

[R] Multiple linear regression (다중회귀분석) End

BioinformaticsAndMe