dose response curve with ggplot

library(drc)
library(ggplot2)
library(drcData)
library(doBy)

ggplot package를 활용하면 기본 제공되는 plot() 함수보다 더 다양한 그래프를 그릴 수 있습니다.

ggplot을 사용하기 앞서 준비해야 될 점들이 있습니다.

우선 dose level을 모두 log 로 변환하여야 합니다. 이는 expand.grid()를 사용하기 위함입니다. 만약 dose level에 “0”이 포함되어 있다면 이는 log 변환이 되지 않으므로 아주 작은 수를 더해 주어야 합니다.

그리고 predict() 함수를 통해 얻을 신뢰구간을 넣을 newdata <- expand.grid()를 만들어 주어야 합니다.

ryegrass에 대해 4 parameter drc를 그립니다. 이후 expand.grid() 함수를 이용하여 dataframe을 만들어 줍니다.

expand.grid() 함수 안에는 concentration의 범위를 seq()로 넣어줍니다. 해당 concentration 범위에 대해 predict() 함수를 이용해 신뢰구간을 산출할 것입니다.

predict() 함수를 이용해 신뢰구간을 산출한 뒤 expand.grid() 함수로 산출한 dataframe에 각각의 열을 합칩니다.

그래프 선 위에 실제 관측값을 점 그래프로 표현할 것입니다. 각각의 점들의 평균과 표준편차를 구하기 위해 doBy 패키지를 사용합니다.

data(ryegrass)
ryegrass
##        rootl  conc
## 1  7.5800000  0.00
## 2  8.0000000  0.00
## 3  8.3285714  0.00
## 4  7.2500000  0.00
## 5  7.3750000  0.00
## 6  7.9625000  0.00
## 7  8.3555556  0.94
## 8  6.9142857  0.94
## 9  7.7500000  0.94
## 10 6.8714286  1.88
## 11 6.4500000  1.88
## 12 5.9222222  1.88
## 13 1.9250000  3.75
## 14 2.8857143  3.75
## 15 4.2333333  3.75
## 16 1.1875000  7.50
## 17 0.8571429  7.50
## 18 1.0571429  7.50
## 19 0.6875000 15.00
## 20 0.5250000 15.00
## 21 0.8250000 15.00
## 22 0.2500000 30.00
## 23 0.2200000 30.00
## 24 0.4400000 30.00
summary(ryegrass)
##      rootl             conc       
##  Min.   :0.2200   Min.   : 0.000  
##  1st Qu.:0.8491   1st Qu.: 0.705  
##  Median :5.0778   Median : 2.815  
##  Mean   :4.3272   Mean   : 7.384  
##  3rd Qu.:7.4262   3rd Qu.: 9.375  
##  Max.   :8.3556   Max.   :30.000
ryegrass$conc[ryegrass$conc == 0] <- 0.5
summary(ryegrass)
##      rootl             conc       
##  Min.   :0.2200   Min.   : 0.500  
##  1st Qu.:0.8491   1st Qu.: 0.830  
##  Median :5.0778   Median : 2.815  
##  Mean   :4.3272   Mean   : 7.509  
##  3rd Qu.:7.4262   3rd Qu.: 9.375  
##  Max.   :8.3556   Max.   :30.000
ryegrass.LL.4 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())

newdata <- expand.grid(conc = exp(seq(log(0.5), log(100), length = 100))) # 꼭 log()로 구간을 잡아야 합니다ㄴ
head(newdata)
##        conc
## 1 0.5000000
## 2 0.5274882
## 3 0.5564876
## 4 0.5870812
## 5 0.6193568
## 6 0.6534068
pm <- predict(ryegrass.LL.4, newdata = newdata, interval = "confidence")
head(pm)
##      Prediction    Lower    Upper
## [1,]   7.787993 7.399890 8.176097
## [2,]   7.781914 7.398324 8.165505
## [3,]   7.774812 7.396129 8.153494
## [4,]   7.766515 7.393116 8.139913
## [5,]   7.756827 7.389047 8.124607
## [6,]   7.745519 7.383616 8.107422
newdata$p <- pm[ , 1]
newdata$pmin <- pm[ , 2]
newdata$pmax <- pm[ , 3]
head(newdata)
##        conc        p     pmin     pmax
## 1 0.5000000 7.787993 7.399890 8.176097
## 2 0.5274882 7.781914 7.398324 8.165505
## 3 0.5564876 7.774812 7.396129 8.153494
## 4 0.5870812 7.766515 7.393116 8.139913
## 5 0.6193568 7.756827 7.389047 8.124607
## 6 0.6534068 7.745519 7.383616 8.107422
ryegrass.mean <- summaryBy(rootl ~ conc, data = ryegrass, FUN = mean)
str(ryegrass.mean)
## 'data.frame':    7 obs. of  2 variables:
##  $ conc      : num  0.5 0.94 1.88 3.75 7.5 15 30
##  $ rootl.mean: num  7.75 7.67 6.41 3.01 1.03 ...
ryegrass.sd <- summaryBy(rootl ~ conc, data = ryegrass, FUN = sd)
str(ryegrass.sd)
## 'data.frame':    7 obs. of  2 variables:
##  $ conc    : num  0.5 0.94 1.88 3.75 7.5 15 30
##  $ rootl.sd: num  0.415 0.724 0.476 1.16 0.166 ...
ryegrass.summary <- merge(ryegrass.mean, ryegrass.sd)
ryegrass.summary
##    conc rootl.mean  rootl.sd
## 1  0.50  7.7493452 0.4151924
## 2  0.94  7.6732804 0.7236913
## 3  1.88  6.4145503 0.4755951
## 4  3.75  3.0146825 1.1595582
## 5  7.50  1.0339286 0.1663975
## 6 15.00  0.6791667 0.1501735
## 7 30.00  0.3033333 0.1193035
ggplot(ryegrass.summary, aes(x = conc, y = rootl.mean)) +
  geom_point() +
  geom_errorbar(aes(ymin = rootl.mean - rootl.sd, ymax = rootl.mean + rootl.sd), width = 0.1) +
  geom_line(data = newdata, aes(x = conc, y = p)) +
  geom_ribbon(data = newdata, aes(x = conc, y = p, ymin = pmin, ymax = pmax), alpha = 0.3) + # alpha는 음영의 세기
  coord_trans(x = "log") + #semi-log 그래프를 그립니다. 
  xlab("conc (mM)") +
  ylab("length (cm)")