다중회귀모형에서 설명변수들의 상대적 중요성

다중회귀모형에서 여러 개의 설명변수가 있는 경우 각 변수들의 상대적인 중요성을 평가하기 위한 방법은 모든 가능한 서브모형을 만들고 설명변수 하나를 추가하였을 때 \(R^2\) 값이 평균적으로 어느 정도 증가하는지 구해 봄으로써 알 수 있다. 다음 함수는 Robert I. Kabacoff의 R in action(2nd edition)책에 소개되어 있는 함수로 회귀모형에서 상대적인 중요성을 보여준다.

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)
}

사용 예는 다음과 같다. 자동차의 연비에 관한 데이터인 mtcars 데이터로 자동차의 연비(mile per gallon, mpg)를 결정하는 요인들을 회귀분석 해보면 다음과 같다.

fit=lm(mpg~wt+hp+cyl+am,data=mtcars)
summary(fit)

Call:
lm(formula = mpg ~ wt + hp + cyl + am, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4765 -1.8471 -0.5544  1.2758  5.6608 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 36.14654    3.10478  11.642 4.94e-12 ***
wt          -2.60648    0.91984  -2.834   0.0086 ** 
hp          -0.02495    0.01365  -1.828   0.0786 .  
cyl         -0.74516    0.58279  -1.279   0.2119    
am           1.47805    1.44115   1.026   0.3142    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.509 on 27 degrees of freedom
Multiple R-squared:  0.849, Adjusted R-squared:  0.8267 
F-statistic: 37.96 on 4 and 27 DF,  p-value: 1.025e-10

다중회귀분석 결과 차의 중량(wt)이 가장 유의하고 마력(hp), 엔진수(cyl), 자동변속여부(am) 순으로 유의한 것으로 나타난다. 이 모형의 \(R^2\) 값은 0.849, adjusted \(R^2\) 값은 0.8267이다. 이 모형에서 설명변수들의 상대적 중요성은 다음과 같이 구할 수 있다.

relweights(fit,col="blue")    

     Weights
am  15.53312
hp  26.75247
cyl 27.00191
wt  30.71250

이 모형의 설명변수 중 차중량(wt)의 상대적 중요성이 30.7%로 나타났다.

설명변수 중 범주형/문자형 변수가 있는 경우

다중회귀모형에서 설명변수가 범주형/문자형 변수가 있어도 회귀분석에는 아무런 문제가 없다. 저자가 만든 moonBook패키지에 포함되어 있는 radial 데이타에는 동맥경화 정도를 나타내는 NTAV(표준화된 동맥경화반 부피)와 함께 여러 임상적인 특징들이 포함되어 있다. NTAV를 반응변수로 하고 나이와 성별, 당뇨병 유무,고혈압유무를 설명변수로 다중회귀분석을 해보면 다음과 같이 할 수 있다.

require(moonBook)
Loading required package: moonBook
fit1=lm(NTAV~age+sex+DM+HBP, data=radial)
summary(fit1)

Call:
lm(formula = NTAV ~ age + sex + DM + HBP, data = radial)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.054 -10.957  -2.535   7.263  85.306 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.2295    14.1442   1.218   0.2258    
age           0.5371     0.2144   2.505   0.0137 *  
sexM         21.9737     4.1654   5.275  6.7e-07 ***
DM            0.2187     4.2774   0.051   0.9593    
HBP          10.2796     4.3041   2.388   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.45 on 110 degrees of freedom
Multiple R-squared:  0.2376,    Adjusted R-squared:  0.2099 
F-statistic: 8.571 on 4 and 110 DF,  p-value: 4.653e-06

다중회귀분석 결과 나이와 성별, 고혈압이 유의한 설명변수로 나타났고 당뇨병 유무는 통계적으로 유의하지 않았다. 이 모형에서 상대적인 중요성을 보기 위해 relweight()함수를 써보면 다음과 같은 에러가 난다.

relweights(fit1)
Error in cor(fit$model): 'x' must be numeric

숫자가 아닌 설명변수는 relweights() 함수로 처리가 안된다. 문자형인 sex를 factor()gkatnfh 범주형변수로 바꾸어 주면 어떻게 될까?

fit2=lm(NTAV~age+factor(sex)+DM+HBP, data=radial)
summary(fit2)

Call:
lm(formula = NTAV ~ age + factor(sex) + DM + HBP, data = radial)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.054 -10.957  -2.535   7.263  85.306 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   17.2295    14.1442   1.218   0.2258    
age            0.5371     0.2144   2.505   0.0137 *  
factor(sex)M  21.9737     4.1654   5.275  6.7e-07 ***
DM             0.2187     4.2774   0.051   0.9593    
HBP           10.2796     4.3041   2.388   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.45 on 110 degrees of freedom
Multiple R-squared:  0.2376,    Adjusted R-squared:  0.2099 
F-statistic: 8.571 on 4 and 110 DF,  p-value: 4.653e-06
relweights(fit2)
Error in cor(fit$model): 'x' must be numeric

여전히 에러가 난다. 범주형 변수인 factor(sex)를 다시 숫자로 바꾸어주면 다음과 같다.

fit3=lm(NTAV~age+as.numeric(factor(sex))+DM+HBP, data=radial)
summary(fit3)

Call:
lm(formula = NTAV ~ age + as.numeric(factor(sex)) + DM + HBP, 
    data = radial)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.054 -10.957  -2.535   7.263  85.306 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -4.7442    16.1717  -0.293   0.7698    
age                       0.5371     0.2144   2.505   0.0137 *  
as.numeric(factor(sex))  21.9737     4.1654   5.275  6.7e-07 ***
DM                        0.2187     4.2774   0.051   0.9593    
HBP                      10.2796     4.3041   2.388   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.45 on 110 degrees of freedom
Multiple R-squared:  0.2376,    Adjusted R-squared:  0.2099 
F-statistic: 8.571 on 4 and 110 DF,  p-value: 4.653e-06
relweights(fit3)

                           Weights
DM                       0.4525551
age                     15.2713467
HBP                     15.4709620
as.numeric(factor(sex)) 68.8051361

사실 radial 데이터에는 문자형변수인 sex이외에 숫자형변수인 male 이라는 변수가 있다. 이 변수는 남자는 1,여자는 0으로 입력되어 있다. sex 대신 male을 써서 다중회귀를 해보면 다음과 같다.

fit4=lm(NTAV~age+male+DM+HBP, data=radial)
summary(fit4)

Call:
lm(formula = NTAV ~ age + male + DM + HBP, data = radial)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.054 -10.957  -2.535   7.263  85.306 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.2295    14.1442   1.218   0.2258    
age           0.5371     0.2144   2.505   0.0137 *  
male         21.9737     4.1654   5.275  6.7e-07 ***
DM            0.2187     4.2774   0.051   0.9593    
HBP          10.2796     4.3041   2.388   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.45 on 110 degrees of freedom
Multiple R-squared:  0.2376,    Adjusted R-squared:  0.2099 
F-statistic: 8.571 on 4 and 110 DF,  p-value: 4.653e-06
relweights(fit4)

        Weights
DM    0.4525551
age  15.2713467
HBP  15.4709620
male 68.8051361

데이터프레임의 문자형 변수를 범주형변수로 바꾸는 방법

데이터프레임에 있는 문자형변수를 모두 숫자형 변수로 바꾸려면 어떻게 해야 할까? 저자가 쓴 “의학논문 작성을 위한 R 통계와 그래프”(한나래,2015) 책에 소개한 방법은 다음과 같다. acs데이터에는 문자형변수가 많이 있다. acs데이터의 구조는 다음과 같다.

str(acs)
'data.frame':   857 obs. of  17 variables:
 $ age             : int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex             : chr  "Male" "Female" "Female" "Female" ...
 $ cardiogenicShock: chr  "No" "No" "Yes" "No" ...
 $ entry           : chr  "Femoral" "Femoral" "Femoral" "Femoral" ...
 $ Dx              : chr  "STEMI" "STEMI" "STEMI" "STEMI" ...
 $ EF              : num  18 18.4 20 21.8 21.8 22 24.7 26.6 28.5 31.1 ...
 $ height          : num  168 148 NA 165 162 153 167 160 152 168 ...
 $ weight          : num  72 48 NA 50 64 59 78 50 67 60 ...
 $ BMI             : num  25.5 21.9 NA 18.4 24.4 ...
 $ obesity         : chr  "Yes" "No" "No" "No" ...
 $ TC              : num  215 NA NA 121 195 184 161 136 239 169 ...
 $ LDLC            : int  154 NA NA 73 151 112 91 88 161 88 ...
 $ HDLC            : int  35 NA NA 20 36 38 34 33 34 54 ...
 $ TG              : int  155 166 NA 89 63 137 196 30 118 141 ...
 $ DM              : chr  "Yes" "No" "No" "No" ...
 $ HBP             : chr  "No" "Yes" "Yes" "No" ...
 $ smoking         : chr  "Smoker" "Never" "Never" "Never" ...

먼저 데이타프레임 중 문자형 변수를 찾아 select에 저장한다.

df<-acs
select<-sapply(df,function(x) {is.character(x)})
select
             age              sex cardiogenicShock            entry 
           FALSE             TRUE             TRUE             TRUE 
              Dx               EF           height           weight 
            TRUE            FALSE            FALSE            FALSE 
             BMI          obesity               TC             LDLC 
           FALSE             TRUE            FALSE            FALSE 
            HDLC               TG               DM              HBP 
           FALSE            FALSE             TRUE             TRUE 
         smoking 
            TRUE

select에 저장된 문자형변수를 범주형변수로 바꾼다.

df[,select]<-lapply(df[,select],factor)
str(df)
'data.frame':   857 obs. of  17 variables:
 $ age             : int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex             : Factor w/ 2 levels "Female","Male": 2 1 1 1 2 1 2 2 1 2 ...
 $ cardiogenicShock: Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
 $ entry           : Factor w/ 2 levels "Femoral","Radial": 1 1 1 1 2 2 2 1 2 1 ...
 $ Dx              : Factor w/ 3 levels "NSTEMI","STEMI",..: 2 2 2 2 1 3 3 2 3 2 ...
 $ EF              : num  18 18.4 20 21.8 21.8 22 24.7 26.6 28.5 31.1 ...
 $ height          : num  168 148 NA 165 162 153 167 160 152 168 ...
 $ weight          : num  72 48 NA 50 64 59 78 50 67 60 ...
 $ BMI             : num  25.5 21.9 NA 18.4 24.4 ...
 $ obesity         : Factor w/ 2 levels "No","Yes": 2 1 1 1 1 2 2 1 2 1 ...
 $ TC              : num  215 NA NA 121 195 184 161 136 239 169 ...
 $ LDLC            : int  154 NA NA 73 151 112 91 88 161 88 ...
 $ HDLC            : int  35 NA NA 20 36 38 34 33 34 54 ...
 $ TG              : int  155 166 NA 89 63 137 196 30 118 141 ...
 $ DM              : Factor w/ 2 levels "No","Yes": 2 1 1 1 2 2 2 2 2 2 ...
 $ HBP             : Factor w/ 2 levels "No","Yes": 1 2 2 1 2 2 2 2 2 1 ...
 $ smoking         : Factor w/ 3 levels "Ex-smoker","Never",..: 3 2 2 2 3 2 1 1 2 3 ...

바뀐 데이터프레임을 보면 모두 범주형변수로 바뀐 것을 알 수 있다. 이상의 과정을 함수로 만들면 다음과 같이 만들 수 있다.

chr2factor=function(df){
    select=sapply(df,function(x) {is.character(x)})
    df[,select]<-lapply(df[,select],factor)
    df
}

함수가 잘 동작하는지 시험해보자.

df2<-chr2factor(acs)
str(df2)
'data.frame':   857 obs. of  17 variables:
 $ age             : int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex             : Factor w/ 2 levels "Female","Male": 2 1 1 1 2 1 2 2 1 2 ...
 $ cardiogenicShock: Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
 $ entry           : Factor w/ 2 levels "Femoral","Radial": 1 1 1 1 2 2 2 1 2 1 ...
 $ Dx              : Factor w/ 3 levels "NSTEMI","STEMI",..: 2 2 2 2 1 3 3 2 3 2 ...
 $ EF              : num  18 18.4 20 21.8 21.8 22 24.7 26.6 28.5 31.1 ...
 $ height          : num  168 148 NA 165 162 153 167 160 152 168 ...
 $ weight          : num  72 48 NA 50 64 59 78 50 67 60 ...
 $ BMI             : num  25.5 21.9 NA 18.4 24.4 ...
 $ obesity         : Factor w/ 2 levels "No","Yes": 2 1 1 1 1 2 2 1 2 1 ...
 $ TC              : num  215 NA NA 121 195 184 161 136 239 169 ...
 $ LDLC            : int  154 NA NA 73 151 112 91 88 161 88 ...
 $ HDLC            : int  35 NA NA 20 36 38 34 33 34 54 ...
 $ TG              : int  155 166 NA 89 63 137 196 30 118 141 ...
 $ DM              : Factor w/ 2 levels "No","Yes": 2 1 1 1 2 2 2 2 2 2 ...
 $ HBP             : Factor w/ 2 levels "No","Yes": 1 2 2 1 2 2 2 2 2 1 ...
 $ smoking         : Factor w/ 3 levels "Ex-smoker","Never",..: 3 2 2 2 3 2 1 1 2 3 ...

하지만 이 방법은 치명적인 단점이 하나 있다. 이 함수는 문자형 변수가 하나만 있을 때는 제대로 동작하지 않는다. 예를 들어 acs데이터 중 첫번째 열인 age와 두번째 열인 sex 만으로 이루어진 acs2를 만들고 테스트해보자.

acs2=acs[1:2]
str(acs2)
'data.frame':   857 obs. of  2 variables:
 $ age: int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex: chr  "Male" "Female" "Female" "Female" ...
df3<-chr2factor(acs2)
Warning in `[<-.data.frame`(`*tmp*`, , select, value =
list(structure(1L, .Label = "Male", class = "factor"), : provided 857
variables to replace 1 variables
str(df3)
'data.frame':   857 obs. of  2 variables:
 $ age: int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex: Factor w/ 1 level "Male": 1 1 1 1 1 1 1 1 1 1 ...

뭔가 경고가 생기고 제대로 변환이 이루어지지 않았다.이 점을 보완한 함수는 다음과 같다.

chr2factor=function(df){
    select=sapply(df,function(x) {is.character(x)})
    if(sum(select)==1) df[[which(select)]]<-factor(df[[which(select)]])
    else if(sum(select)>1) df[,select]<-lapply(df[,select],factor)
    df
}

이번에는 제대로 동작하는지 테스트해보자.

str(acs2)
'data.frame':   857 obs. of  2 variables:
 $ age: int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex: chr  "Male" "Female" "Female" "Female" ...
df4<-chr2factor(acs2)
str(df4)
'data.frame':   857 obs. of  2 variables:
 $ age: int  62 78 76 89 56 73 58 62 59 71 ...
 $ sex: Factor w/ 2 levels "Female","Male": 2 1 1 1 2 1 2 2 1 2 ...

함수가 제대로 동작하는 것을 확인했으니 다시 범주형변수를 숫자형변수로 만드는 함수를 만들어보면 다음과 같다.

factor2number=function(df){
    select=sapply(df,function(x) {is.factor(x)})
    if(sum(select)==1) df[[which(select)]]<-as.numeric(df[[which(select)]])
    else if(sum(select)>1) df[,select]<-lapply(df[,select],as.numeric)
    df
}

relweight()함수의 개선

처음 만들었던 relweights()함수를 개선하여 회귀모형에 문자형 변수, 범주형 변수가 포함되어 있으면 숫자형변수로 바꾸어 주고 상대적 중요성을 구해주는 함수를 만들어 보면 다음과 같다.

relweights <- function(fit,plot=FALSE,...){
    df<-fit$model
    df1<-chr2factor(df)
    df2<-factor2number(df1)
    
    R <- cor(df2)
    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]
    if(plot) {
        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)),
                 ...)
    }  
    attr(import,"R-square")<-round(rsquare, digits=3)
    return(import)
}

제대로 동작하는지 확인해보면 다음과 같다.

fit=lm(NTAV~age+sex+DM+HBP, data=radial)
summary(fit)

Call:
lm(formula = NTAV ~ age + sex + DM + HBP, data = radial)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.054 -10.957  -2.535   7.263  85.306 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.2295    14.1442   1.218   0.2258    
age           0.5371     0.2144   2.505   0.0137 *  
sexM         21.9737     4.1654   5.275  6.7e-07 ***
DM            0.2187     4.2774   0.051   0.9593    
HBP          10.2796     4.3041   2.388   0.0186 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.45 on 110 degrees of freedom
Multiple R-squared:  0.2376,    Adjusted R-squared:  0.2099 
F-statistic: 8.571 on 4 and 110 DF,  p-value: 4.653e-06
result=relweights(fit)
result
       Weights
DM   0.4525551
age 15.2713467
HBP 15.4709620
sex 68.8051361

ggplot을 이용한 relweights의 시각화

위에서 구한 ggplot2를 이용하여 가로막대그래프를 그려주는 함수는 다음과 같이 만들 수 있다.

require(ggplot2)
Loading required package: 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
   }

위의 함수를 이용하여 plot을 그려보면 다음과 같다.

plotRelWeights(fit1)