문제

디음 plot의 회귀선에 평행한 라벨을 붙이는 방법을 알아보겠습니다.

require(ggplot2)
ggplot(data=iris,aes(x=Sepal.Length,y=Sepal.Width,color=Species))+stat_smooth(method="lm",se=FALSE)

회귀모형에서 회귀식 추출

먼저 iris데이터에서 Sepal.Width를 반응변수로, Sepal.Length와 Species를 설명변수로 하는 회귀모형을 만듭니다. 이때 설명변수 간의 상호작용이 있어야 하므로 Sepal.Length*Species로 사용합니다

fit=lm(Sepal.Width~Sepal.Length*Species,data=iris)
summary(fit)

Call:
lm(formula = Sepal.Width ~ Sepal.Length * Species, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.72394 -0.16327 -0.00289  0.16457  0.60954 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     -0.5694     0.5539  -1.028 0.305622    
Sepal.Length                     0.7985     0.1104   7.235 2.55e-11 ***
Speciesversicolor                1.4416     0.7130   2.022 0.045056 *  
Speciesvirginica                 2.0157     0.6861   2.938 0.003848 ** 
Sepal.Length:Speciesversicolor  -0.4788     0.1337  -3.582 0.000465 ***
Sepal.Length:Speciesvirginica   -0.5666     0.1262  -4.490 1.45e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2723 on 144 degrees of freedom
Multiple R-squared:  0.6227,    Adjusted R-squared:  0.6096 
F-statistic: 47.53 on 5 and 144 DF,  p-value: < 2.2e-16

이 회귀모형에서 Species에 따른 회귀식의 slope와 intercept를 추출합니다. 이를 이용하여 label을 만듭니다.

Species=levels(iris$Species)
intercept=c(fit$coef[1],fit$coef[1]+fit$coef[3],fit$coef[1]+fit$coef[4])
slope=c(fit$coef[2],fit$coef[2]+fit$coef[5],fit$coef[2]+fit$coef[6])
df=data.frame(Species,intercept,slope)
df$label=paste0("italic(y)==",round(df$intercept,2),"+",round(df$slope,2),"*italic(x)")
df
     Species  intercept     slope                          label
1     setosa -0.5694327 0.7985283 italic(y)==-0.57+0.8*italic(x)
2 versicolor  0.8721460 0.3197193 italic(y)==0.87+0.32*italic(x)
3  virginica  1.4463054 0.2318905 italic(y)==1.45+0.23*italic(x)

직선의 slope는 직선의 가지는 각도의 arctan 값입니다. R에서는 atan()합수로 계산됩니다. 단 이렇게 계산한 각도는 radian이라 ggplot에서 쓰려면 degree로 바꾸어주어야 합니다.

df$radian=atan(df$slope)
df$angle=df$radian*180/pi
df
     Species  intercept     slope                          label    radian
1     setosa -0.5694327 0.7985283 italic(y)==-0.57+0.8*italic(x) 0.6738429
2 versicolor  0.8721460 0.3197193 italic(y)==0.87+0.32*italic(x) 0.3094483
3  virginica  1.4463054 0.2318905 italic(y)==1.45+0.23*italic(x) 0.2278632
     angle
1 38.60836
2 17.73008
3 13.05560

또한 라벨을 붙일 x좌표와 y좌표를 계산합니다. 편의상 x축의 5,6,7.2 위치에 라벨을 붙이기 위해 다음과 같은 함수를 만들고 좌표를 계산합니다.

fun=list()
fun[[1]]=function(x){fit$coef[1]+fit$coef[2]*x}
fun[[2]]=function(x){fit$coef[1]+fit$coef[3]+(fit$coef[2]+fit$coef[5])*x}
fun[[3]]=function(x){fit$coef[1]+fit$coef[4]+(fit$coef[2]+fit$coef[6])*x}

df$x=c(5,6,7.2)
df$y=unlist(lapply(1:3,function(i){fun[[i]](df$x[i])}))
df
     Species  intercept     slope                          label    radian
1     setosa -0.5694327 0.7985283 italic(y)==-0.57+0.8*italic(x) 0.6738429
2 versicolor  0.8721460 0.3197193 italic(y)==0.87+0.32*italic(x) 0.3094483
3  virginica  1.4463054 0.2318905 italic(y)==1.45+0.23*italic(x) 0.2278632
     angle   x        y
1 38.60836 5.0 3.423209
2 17.73008 6.0 2.790462
3 13.05560 7.2 3.115917

회귀선 확인

먼저 회귀선이 이 그림과 일치하는지 plot에 직선을 추가해 선이 일치하는지 알아봅니다.

p <-ggplot(data=iris,aes(x=Sepal.Length,y=Sepal.Width,color=Species))+
    stat_smooth(method="lm",se=FALSE)
p+  stat_function(fun=fun[[1]],lty=2)+
    stat_function(fun=fun[[2]],lty=2)+
    stat_function(fun=fun[[3]],lty=2)

당연히 선이 일치합니다. 이제 우리가 계산한 각도로 라벨을 붙여봅니다.

p+ geom_text(data=df,aes(x=x,y=y,label=label,angle=angle,color=Species),parse=TRUE)

회귀선과 라벨의 각도가 맞지 않는 것을 확인할 수 있습니다. 맞지 않는 이유는 이 plot의 x축과 y축의 비율 즉 aspect ratio가 1이 아니기 때문입니다. 만일 다음과 같이 coord_fixed()를 추가하면 정확히 일치합니다.

p+ geom_text(data=df,aes(x=x,y=y,label=label,angle=angle,color=Species),parse=TRUE)+
    coord_fixed()

하지만 coord_fixed()를 사용하니 plot이 너무 납작해집니다. 이를 해결하기 위해 어떻게 해야 할까요? 저는 먼저 ggplot을 그린후 그 ggplot의 aspect ratio를 계산하여 그 ratio로 각도를 변형시켰습니다. 먼저 ggplot의 aspect ratio를 구하는 함수를 하나 만들었습니다.

getAspectRatio=function(p){
    xmin=layer_scales(p)$x$range$range[1]
    xmax=layer_scales(p)$x$range$range[2]
    ymin=layer_scales(p)$y$range$range[1]
    ymax=layer_scales(p)$y$range$range[2]
    cat("xmin=",xmin,"\n")
    cat("xmax=",xmax,"\n")
    cat("ymin=",ymin,"\n")
    cat("ymax=",ymax,"\n")

    (xmax-xmin)/(ymax-ymin)
}

이 plot의 aspect ratio를 구하면 다음과 같습니다.

ratio=getAspectRatio(p)
xmin= 4.3 
xmax= 7.9 
ymin= 2.438771 
ymax= 4.062031 
ratio
[1] 2.217758

이 ratio를 이용해 기존에 계산한 slope를 변형합니다.

df$slope2=df$slope*ratio
df$radian2=atan(df$slope2)
df$angle2=df$radian2*180/pi
df
     Species  intercept     slope                          label    radian
1     setosa -0.5694327 0.7985283 italic(y)==-0.57+0.8*italic(x) 0.6738429
2 versicolor  0.8721460 0.3197193 italic(y)==0.87+0.32*italic(x) 0.3094483
3  virginica  1.4463054 0.2318905 italic(y)==1.45+0.23*italic(x) 0.2278632
     angle   x        y    slope2   radian2   angle2
1 38.60836 5.0 3.423209 1.7709428 1.0567592 60.54784
2 17.73008 6.0 2.790462 0.7090602 0.6167808 35.33894
3 13.05560 7.2 3.115917 0.5142771 0.4750039 27.21572

이제 이렇게 계산한 각도로 라벨을 붙입니다.

p+ geom_text(data=df,aes(x=x,y=y,label=label,angle=angle2,color=Species),parse=TRUE)+
    coord_fixed(ratio=ratio)

마지막으로 직선과 회귀식이 겹치므로 vjust를 이용하여 조절합니다.

df$vjust=c(-0.5,1.5,-0.5)
p+ geom_text(data=df,aes(x=x,y=y,label=label,angle=angle2,color=Species,vjust=vjust),parse=TRUE)+
    coord_fixed(ratio=ratio)

주어진 ggplot의 aspect ratio를 알아내 간단한 삼각함수를 통해 문제를 해결하였습니다.