메뉴 건너뛰기

웹에서 하는 R 통계

다중회귀모형에서 상호작용의 시각화

단순회귀모형

먼저 단순회귀모형에서 회귀계수를 추출하고 회귀선을 시각화하는 방법은 다음과 같다. 예를 들어 자동차의 연비에 관한 mtcars 데이터에서 연비(mile per gallon, mpg)를 반응변수(종속변수)로 자동차의 공차중량(weight, wt)를 영향변수(독립변수)를 하는 단순회귀모형을 만들어보면 다음과 같다.

fit=lm(mpg~wt,data=mtcars)
summary(fit)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

회귀분석 결과에서 mpg와 wt의 관계는 기울기가 -5.3445이고 y절편이 37.2851인 직선이라는 것을 알 수 있다. 즉

mpg=-5.3445 * wt + 37.2851

이다. 회귀계수는 다음의 명령으로 알 수 있다.

coef(fit)
(Intercept)          wt 
  37.285126   -5.344472 
str(coef(fit))
 Named num [1:2] 37.29 -5.34
 - attr(*, "names")= chr [1:2] "(Intercept)" "wt"

즉, y절편은 coef(fit)[1],기울기는 coef(fit)[2]이다. 이를 plot함수로 그려보면 다음과 같다.

intercept=coef(fit)[1]
slope=coef(fit)[2]
plot(mpg~wt,data=mtcars)
abline(a=intercept,b=slope,col="red")  #(1)

위의 코드에서 (1)은 abline(fit, col=“red”)로 입력해도 같은 결과를 보이나 여기서는 y절편과 기울기를 이용한 직선을 보여주기 위해 위의 코드를 사용하였다. 같은 plot을 ggplot으로 그려보면 다음과 같다.

require(ggplot2)
Loading required package: ggplot2
p<-ggplot(data=mtcars,aes(x=wt,y=mpg))+geom_point()
p+geom_abline(intercept=intercept,slope=slope,col="red") #(2)

물론 (2)대신에 ggplot의 stat_smooth를 사용하면 일일이 y절편과 기울기를 지정하지 않아도 회귀선을 그릴 수 있다. (95% 신뢰구간을 표시하지 않으려면 인수로 se=FALSE를 주면 된다.)

p+stat_smooth(method=lm)

상호작용이 있는 다중회귀모형

하나의 연속형 변수와 하나의 범주형 변수를 갖는 경우

다음으로 하나의 연속형 변수와 하나의 범주형 변수및 상호작용을 갖는 다중회귀모형에서 회귀계수를 추출하고 회귀선을 시각화해보자. mtcars 데이타에서 am은 수동변속기는 0, 자동변속기는 1로 입력되어 있다. 새로 am1이라는 변수를 만들고 수동은 “manual”, 자동은 “automatic”으로 입력해보자.

mtcars$am1=ifelse(mtcars$am==1,"automatic","manual")

mpg를 반응변수로, am1과 wt를 영향변수로 갖는 다중회귀모형을 만들어보자.

fit1=lm(mpg~wt*am1,data=mtcars)
summary(fit1)

Call:
lm(formula = mpg ~ wt * am1, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6004 -1.5446 -0.5325  0.9012  6.0909 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    46.294      3.010  15.379 3.49e-15 ***
wt             -9.084      1.212  -7.493 3.68e-08 ***
am1manual     -14.878      4.264  -3.489  0.00162 ** 
wt:am1manual    5.298      1.445   3.667  0.00102 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.591 on 28 degrees of freedom
Multiple R-squared:  0.833, Adjusted R-squared:  0.8151 
F-statistic: 46.57 on 3 and 28 DF,  p-value: 5.209e-11

am1은 “automatic”, “manual”의 두개의 고유값을 가진다. 우리가 따로 reference를 정해주지 않으면 알파벳순으로 빠른 “automatic”을 reference로 취급한다. 회귀 결과 요약을 보면 먼저 자동변속기의 경우(am1이 “automatic”인 경우)의 wt와 mpg의 회귀식은 다음과 같다.

mpg=-9.084 * wt + 46.294

수동변속기의 경우(am1이 “manual”인 경우) 회귀식은 다음과 같다.

mpg=(-9.084+5.298) * wt+(46.294-14.878)

이 두개의 회귀식을 직선으로 plot해 보면 다음과 같다.

plot(mpg~wt,data=mtcars,col=factor(mtcars$am1))
abline(a=fit1$coef[1],b=fit1$coef[2])
abline(a=fit1$coef[1]+fit1$coef[3],b=fit1$coef[2]+fit1$coef[4],col="red")

ggplot을 이용하여 plot을 하고 stat_smooth()를 사용하면 회귀직선을 알아서 그려준다.

p<-ggplot(data=mtcars,aes(x=wt,y=mpg,colour=am1))+geom_point()+geom_smooth(method="lm")
p

ggplot의 회귀선이 맞는지 보기 위해 우리가 계산한 회귀직선을 덧그려보면 다음과 같다..

p<-p+geom_abline(intercept=fit1$coef[1],slope=fit1$coef[2],colour="orange",linetype=2)
p<-p+geom_abline(intercept=fit1$coef[1]+fit1$coef[3],slope=fit1$coef[2]+fit1$coef[4],
                 colour="green",linetype=2)
p

당연한 이야기이지만 회귀직선이 정확히 일치함을 알수 있다.

두개의 연속형 변수의 상호작용이 있는 경우

이번에는 두개의 연속형 변수의 상호작용을 갖는 다중회귀모형의 경우를 생각해보자. mtcars데이터를 이용하여 mpg를 반응변수로, 공차중량(wt)과 마력수(hp) 및 상호작용을 설명변수로 하는 회구모형을 만들면 다음과 같이 할 수 있다.

fit2=lm(mpg~wt*hp,data=mtcars)
summary(fit2)

Call:
lm(formula = mpg ~ wt * hp, 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 ***
wt          -8.21662    1.26971  -6.471 5.20e-07 ***
hp          -0.12010    0.02470  -4.863 4.04e-05 ***
wt:hp        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와 mpg 간의 회귀공식은 wt의 값에 따라 달라진다. wt 값의 사분위수를 알아보면 다음과 같다

summary(mtcars$wt)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.513   2.581   3.325   3.217   3.610   5.424 

즉 mtcars데이타에 있는 자동차들 중 가장 가벼운 차는 1.513, 가장 무거운 차는 5.424임을 알 수 있다. wt와 mpg사이의 회귀공식은 hp에 따라 달라진다. 즉 회귀계수를 소수점 두자리에서 반올림하면 다음과 같은 회귀식을 얻을 수 있다.

mpg=49.81 - 0.12 * hp- 8.22 * wt + 0.03 * hp * wt 

wt의 값이 평균인 3.2일때 회귀공식은 다음과 같이 구할 수 있다.

mpg=49.81-0.12 * hp - 8.22 * 3.2 + 0.03 * hp * 3.2 

mpg=23.51 - 0.03 * hp 

wt의 평균과 3.2 및 평균 +- 표준편차인 2.2, 4.2를 이용하여 상호작용을 시각화해보면 다음과 같다.

A<-c(2.2,3.2,4.2)        
coef<-fit2$coef
labels=as.character(A)
intercept=coef[1]+coef[2]*A
slope=coef[3]+coef[4]*A
intercept
[1] 31.73185 23.51523 15.29860
slope
[1] -0.058836165 -0.030988016 -0.003139868
colour=rainbow(length(A))
df1=data.frame(A,intercept,slope,colour)   
df1
    A intercept        slope    colour
1 2.2  31.73185 -0.058836165 #FF0000FF
2 3.2  23.51523 -0.030988016 #00FF00FF
3 4.2  15.29860 -0.003139868 #0000FFFF
p<-ggplot(data=mtcars,aes(x=hp,y=mpg))+geom_point()
p<-p+geom_abline(data=df1,aes(slope=slope,intercept=intercept,colour=colour))
p<-p+scale_colour_discrete(labels=labels)+labs(colour="wt")
p

위의 상호작용을 시각화해주는 plot을 함수화하면 다음과 같이 할 수 있다.

ggEffect=function(fit,x=1,probs=c(0.25,0.5,0.75),point=TRUE,xvalue=NULL){
        df=fit$model
        coef=fit$coef
        name=colnames(df)
    
        p<-ggplot(data=df,aes_string(x=name[1+x],y=name[1]))
        if(is.null(xvalue)) A=quantile(df[[4-x]],probs,na.rm=TRUE)
        else A=xvalue
        labels=as.character(A)
        intercept=coef[1]+coef[4-x]*A
        slope=coef[1+x]+coef[4]*A
        colour=rainbow(length(A))
        df1=data.frame(A,intercept,slope,colour)
        if(point) p<-p+geom_point()
        p<-p+geom_abline(data=df1,aes(slope=slope,intercept=intercept,colour=colour))
        p<-p+scale_colour_discrete(labels=labels)+labs(colour=name[4-x])
        p
 
}        

위 함수에서 x=1이면 첫번째 설명변수를 x축에 할당하고 x=2이면 두번째 설명변수를 x축에 할당한다. 아무 값도 주지 않으면 c(0.25,0.5,0.75)에 해당하는 quantile값을 구해 회귀직선을 그려준다. xvalue에 값을 주면 그 값으로 회귀직선을 그려준다. 사용예는 다음과 같다.

ggEffect(fit2)