Single Dose Response Curve Analysis
Cato.L
2020. 06. 24
Dose response analysis using R 교재의 내용을 정리합니다. 정확한 지식은 해당 책을 구매하여 읽어보시길 권장합니다.
사용하는 library list
library(drc)
library(devtools)
library(drcData) # devtools library 설치 후 install_github("DoseResponse/drcData") 로 설치치
library(boot)
library(lmtest)
library(metafor)
library(sandwich)1. 연속형 데이터
농도, 용량과 같은 연속형 데이터에서 반응 값인 \(y\)는 농도 \(x\)와 오차 \(\epsilon\)으로 나눌 수 있습니다.
\[ yi = f(xi, \beta) + \epsilon i \]
약물 dose (\(xi\)) 의 증가에 따라 결과 값인 (\(yi\)) 의 값이 일정하게 변하나 random error인 \(\epsilon i\)가 항상 영향을 주는 함수식을 만들어볼 수 있습니다.
같은 약물을 같은 농도로 반복하여 처리할 시, 약물에 의한 변화 외에도 랜덤오차가 발생함으로써 항상 같은 결과값을 얻을 수는 없습니다. 일반적으로 랜덤오차는 \(\epsilon\) ~ \(N(0, \sigma^2)\) 을 따르며 모든 농도에서 등분산이라고 가정합니다.
1.1 반복이 없는 DRC
1.1.1 Inhibitory effect of secalonic acid
secalonic## dose rootl
## 1 0.000 5.5
## 2 0.010 5.7
## 3 0.019 5.4
## 4 0.038 4.6
## 5 0.075 3.3
## 6 0.150 0.7
## 7 0.300 0.4secalonic acid (mM) 가 plant root (cm) 의 길이에 주는 영향을 나타낸 데이터입니다.
1.1.1.1 Fitting the model
secalonic.LL.4 <- drm(rootl ~ dose, data = secalonic, fct = LL.4(), type = "continuous")
# LL.4() = 4 parametric logistic model
# type = continuous가 기본값
names(secalonic.LL.4)## [1] "varParm" "fit" "curve" "summary" "start"
## [6] "parNames" "predres" "call" "data" "parmMat"
## [11] "fct" "robust" "estMethod" "df.residual" "sumList"
## [16] "scaleFct" "pmFct" "pfFct" "type" "indexMat"
## [21] "logDose" "cm" "deriv1" "curveVarNam" "origData"
## [26] "weights" "dataList" "coefficients" "boxcox" "indexMat2"# drm 내 결과 값들의 종류를 확인할 수 있음plot(secalonic.LL.4, bp = 1e-3, broken = T, ylim = c(0, 7), xlab = "Dose (mM)", ylab = "Root length (cm)")
broken = 농도가 0인 경우 log 값 변환이 불가하므로 그래프에 broken을 넣어줍니다.
bp = 그래프에서 break point를 잡아줍니다. 위에서는 0.001을 break point로 설정하였습니다.
summary(secalonic.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 2.6542086 0.6962333 3.8122 0.0317398 *
## c:(Intercept) 0.0917852 0.3747246 0.2449 0.8223012
## d:(Intercept) 5.5297495 0.2010300 27.5071 0.0001055 ***
## e:(Intercept) 0.0803547 0.0078829 10.1935 0.0020121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.2957497 (3 degrees of freedom)# b : hill slope
# c : Emin
# d : Enax
# e : EC50EC50가 0.08 mM이 산출되었고 secalonic data를 보면 거의 비슷한 것을 직관적으로 알 수 있습니다.
t-value는 Estimate / Std.Error의 값입니다. t-분포를 기반으로 하고 이에 대해 p-value가 주어져 있습니다.
해당 검정의 H0은 “parameter = 0” 입니다. 반드시 p-value가 유의해야만 해당 factor를 채택하는 것은 아닙니다.
Residual standard error는 자주 사용하는 값은 아니며 해당 실험의 분산에 대한 참고치로 생각하면 됩니다.
1.1.1.2 Estimation of arbitrary ED values
ED(secalonic.LL.4, c(10, 25, 50, 75, 90))##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:10 0.0351149 0.0078689
## e:1:25 0.0531192 0.0071527
## e:1:50 0.0803547 0.0078829
## e:1:75 0.1215547 0.0191012
## e:1:90 0.1838789 0.0462775ED() 함수는 EC50를 구하는데 사용합니다. 위의 예와 같이 여러 x 값에 해당하는 y 값을 산출할 수 있습니다.
ED 값은 model parmter 중 holl slope (b) 값을 따르므로 simulation 하는 것이 어렵지 않습니다.
predict(secalonic.LL.4, data.frame(dose = ED(secalonic.LL.4, c(50), display = F)))## Prediction
## 2.810767predict() 함수를 이용하면 역으로 EC50 값에서의 y 값을 산출할 수 있습니다. display = F 옵션을 넣지 않으면 ED50 값이 중복으로 표기되니 F 지정을 추천합니다.
1.1.2 Data from a fish test in ecotoxicology
head(O.mykiss) ; tail(O.mykiss)## conc weight
## 1 0 2.8
## 2 0 3.0
## 3 0 2.7
## 4 0 3.9
## 5 0 3.1
## 6 0 1.8## conc weight
## 65 46 NA
## 66 46 NA
## 67 46 NA
## 68 46 NA
## 69 46 NA
## 70 46 NAsummary(O.mykiss)## conc weight
## Min. : 0.00 Min. :0.900
## 1st Qu.: 1.00 1st Qu.:2.200
## Median : 4.60 Median :2.700
## Mean :12.26 Mean :2.641
## 3rd Qu.:22.00 3rd Qu.:3.000
## Max. :46.00 Max. :4.000
## NA's :9독성시험 데이터이며 몇 개의 missing value를 포함하고 있습니다. 한 농도당 10반복의 실험을 수행했습니다.
missing value는 원하지 않는 bias 결과를 유도할 수 있으므로 제거합니다.
해당 분석에서는 2-parmeter exponential decay model 을 사용하며 이는 OECD의 권고사항입니다.
O.mykiss.EXD.2 <- drm(weight ~ conc, data = O.mykiss, fct = EXD.2(), na.action = na.omit)
summary(O.mykiss.EXD.2)##
## Model fitted: Exponential decay with lower limit at 0 (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## d:(Intercept) 2.846794 0.092526 30.7674 < 2.2e-16 ***
## e:(Intercept) 111.738614 33.196876 3.3659 0.001347 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.5598508 (59 degrees of freedom)plot(O.mykiss.EXD.2, broken = T, type = "all", xlab = "Conc. (mg/l)", ylab = "weight (g)",
xlim = c(0, 500), ylim = c(0, 4))
만약 평균에 해당하는 점만 찍고 싶다면 type = “average”로 설정하면 됩니다.
해당 data가 위에서 언급한 잔차들의 랜덤성, 정규성을 따르는지 확인해 봅시다.
plot(x = fitted(O.mykiss.EXD.2), y = residuals(O.mykiss.EXD.2))
abline(h = 0, lty = 2)
잔차는 0을 기준으로 경향성 없이 고르게 퍼져 있습니다.
qqnorm(residuals(O.mykiss.EXD.2))
abline(a = 0, b = sd(residuals(O.mykiss.EXD.2)))
qq-plot을 그리고 intercept = 0, slope는 잔차들의 SD를 반영하여 그립니다.
실제 값과 기대값이 매칭되어 linear한 그래프가 형성되므로 정규성을 따른다고 볼 수 있습니다.
1.1.3 Ferulic acid as an herbicide
head(ryegrass) ; tail(ryegrass)## rootl conc
## 1 7.580000 0
## 2 8.000000 0
## 3 8.328571 0
## 4 7.250000 0
## 5 7.375000 0
## 6 7.962500 0## rootl conc
## 19 0.6875 15
## 20 0.5250 15
## 21 0.8250 15
## 22 0.2500 30
## 23 0.2200 30
## 24 0.4400 30summary(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천연성분 제초제인 효능을 평가하기 위한 실험입니다.conc는 ferulic acid 의 농도(mM), rootl은 뿌리길이(cm) 값입니다.
7개의 농도구간이며 0 농도에서는 6반복, 다른 농도에서는 3반복 실험을 진행했습니다.
ryegrass.LL.4 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
summary(ryegrass.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 2.98222 0.46506 6.4125 2.960e-06 ***
## c:(Intercept) 0.48141 0.21219 2.2688 0.03451 *
## d:(Intercept) 7.79296 0.18857 41.3272 < 2.2e-16 ***
## e:(Intercept) 3.05795 0.18573 16.4644 4.268e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.5196256 (20 degrees of freedom)par(mfrow = c(1, 3))
plot(ryegrass.LL.4, broken = T, type = "all", xlab = "mM", ylab = "cm")
plot(x = fitted(ryegrass.LL.4), y = residuals(ryegrass.LL.4))
abline(h = 0, lty = 2)
qqnorm(residuals(ryegrass.LL.4))
abline(a = 0, b = sd(residuals(ryegrass.LL.4)))
par(mfrow = c(1, 1))농도의존 곡선을 확인해보면 중간 농도에서 값의 분산이 가장 커지고 농도가 높아질수록 분산이 작아집니다. 따라서 등분산이 성립하는지 의심을 해보아야 합니다.
residual plot에서도 중간 농도에서 잔차가 커지는 것을 볼 수 있습니다. 하지만 위 아래의 퍼짐정도가 비슷하므로 해당 모델은 적절하게 average를 설명하고 있다고 볼 수 있습니다.
따라서 해당 모델은 등분산 가정을 만족하지 않는다고 보기 어렵습니다. 이와 같은 패턴은 생물학 실험에서 흔히 관찰되는 경우이므로 적절한 배경지식을 바탕으로 한 해석이 중요합니다.
qq-plot에서도 정규성을 벗어난 경향이 없습니다.
1.1.4 Glyphosate in barley
head(barley) ; tail(barley)## Dose weight
## 1 0.00000 57.2
## 2 0.00000 49.8
## 3 21.09375 62.2
## 4 21.09375 30.6
## 5 42.18750 40.9
## 6 42.18750 70.9## Dose weight
## 13 675 16.0
## 14 675 10.4
## 15 1350 12.8
## 16 1350 9.8
## 17 2700 10.4
## 18 2700 7.4summary(barley)## Dose weight
## Min. : 0.00 Min. : 7.40
## 1st Qu.: 42.19 1st Qu.:13.60
## Median : 168.75 Median :35.75
## Mean : 597.66 Mean :34.60
## 3rd Qu.: 675.00 3rd Qu.:52.70
## Max. :2700.00 Max. :70.90약물 처리에 따라 작물이 자라는 것을 확인한 데이터입니다.
barley.LL.4 <- drm(weight ~ Dose, data = barley, fct = LL.4(), na.action = na.omit)
summary(barley.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 9.7084 42.9166 0.2262 0.82430
## c:(Intercept) 11.1275 3.7803 2.9435 0.01068 *
## d:(Intercept) 52.0478 3.2487 16.0212 2.123e-10 ***
## e:(Intercept) 286.2600 209.6374 1.3655 0.19364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 9.241941 (14 degrees of freedom)plot(barley.LL.4, broken = T, xlab = "Dose", ylab = "Weight", type = "all")
plot(x = fitted(barley.LL.4), y = residuals(barley.LL.4))
abline(h = 0, lty = 2)
qqnorm(residuals(barley.LL.4))
abline(a = 0, b = sd(residuals(barley.LL.4)))
정규성에는 이상이 없지만 잔차의 크기가 굉장히 큰 값들이 존재합니다. 몇몇의 경우에는 data의 형태를 바꾸는 것이 등분산을 가정하는데 도움이 됩니다.
dose-response model의 가정들을 유지하기 위해서는 response와 predicted value를 model에 근거하여 변환시켜주어야 합니다.
여기서는 Box-Cox transformation을 사용해보도록 하겠습니다.
단, 변수의 변환은 dose-response curve에 부적절한 결과를 가져올 수 있으니 일반적 방법으로 고려해서는 안됩니다.
barley.LL.4.bc <- boxcox(barley.LL.4, method = "anova", plotit = F)
summary(barley.LL.4.bc)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 4.6567 4.3207 1.0778 0.299362
## c:(Intercept) 10.5303 1.1542 9.1231 2.875e-07 ***
## d:(Intercept) 51.4583 5.4983 9.3590 2.109e-07 ***
## e:(Intercept) 250.8884 52.3177 4.7955 0.000285 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.1121892 (14 degrees of freedom)
##
## Non-normality/heterogeneity adjustment through optimal Box-Cox transformation
##
## Estimated lambda: -0.25
## Confidence interval for lambda: [-0.576, 0.370]plot(x = fitted(barley.LL.4.bc), y = residuals(barley.LL.4.bc))
abline(h = 0, lty = 2)
\(\lambda\)는 Box-Cox 변환에서 optimal power exponent를 나타냅니다. summary() 결과를 보면 -0.25, 신뢰구간 -0.576 ~ 0.370을 출력했습니다.
신뢰구간이 1을 포함하지 않았고 따라서 \(\lambda\) 값은 1과는 유의적으로 다르나 0과는 다르지 않습니다.
변환 후 결과를 보면 변환 전과 b (hill slope) 값을 제외하고 큰 변화가 없습니다.
다음과 같이 drm 함수 안에 bcVal 항목을 추가함으로써 transformation을 할 수도 있습니다.
barley.LL.4.log <- drm(weight ~ Dose, data = barley, fct = LL.4(), na.action = na.omit, bcVal = 0)
# bcVal = lambda in Box-Cox transformation / lambda = 0 인 경우 Box-Cox transformation은 log transfomration 이다.
# 만약 lambda = 2라면 Box-Cox transformation은 root transformation이다.
summary(barley.LL.4.log)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 6.6523 8.8818 0.7490 0.466271
## c:(Intercept) 10.7738 1.1385 9.4628 1.843e-07 ***
## d:(Intercept) 51.1127 4.4105 11.5890 1.461e-08 ***
## e:(Intercept) 269.2678 76.8634 3.5032 0.003513 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.2457367 (14 degrees of freedom)
##
## Non-normality/heterogeneity adjustment through optimal Box-Cox transformation
##
## Specified lambda: 0plot(x = fitted(barley.LL.4.log), y = residuals(barley.LL.4.log))
abline(h = 0, lty = 2)
transformation이 되기 전과 후의 plot을 겹쳐서 그려봅시다.
plot(barley.LL.4, type = "all")
plot(barley.LL.4.bc, add = T, lty = 2, type = "none")
plot(barley.LL.4.log, add = T, lty = 2, col = "red", type = "none")
legend(x = "topright", legend = c("not transfomration", "Box-cox transformation", "log transformation"),
col = c("black", "black", "red"), lty = c(1, 2, 2))
transformation 방법에 따라 그래프의 fitting이 조금씩 다르나 큰 차이는 없어보입니다.
따라서, transformation을 하지 않은 원 그래프를 채택해도 문제는 없습니다. 다만 standard error는 sandwich method를 채택하는 것이 더 좋아보입니다.
해당 data는 그래프가 꺽어는 부분에서의 point가 너무 적어 많은 결과를 추론에 의존할 수 밖에 없습니다. 이런 경우, 적합한 model을 선정하기 어렵습니다.
1.1.5 Lower limits for dose-response data
DRC에서 basal 값이 0이 아닌 경우는 자주 등장합니다. 시험계의 특징일수도 있고 실험오차에 의해 0이 아닌 값이 나타날 수도 있습니다.
head(ryegrass2) ; tail(ryegrass2)## dose biomass day
## 1 0 77.5 0
## 2 0 78.0 0
## 3 0 75.0 0
## 4 0 214.4 15
## 5 0 215.4 15
## 6 0 227.5 15## dose biomass day
## 22 50 80.0 15
## 23 50 93.3 15
## 24 50 84.0 15
## 25 100 78.8 15
## 26 100 84.6 15
## 27 100 71.5 15summary(ryegrass2)## dose biomass day
## Min. : 0.000 Min. : 71.5 0 : 3
## 1st Qu.: 1.562 1st Qu.: 84.3 15:24
## Median : 6.250 Median :139.8
## Mean : 22.049 Mean :143.5
## 3rd Qu.: 25.000 3rd Qu.:207.2
## Max. :100.000 Max. :234.3str(ryegrass2)## 'data.frame': 27 obs. of 3 variables:
## $ dose : num 0 0 0 0 0 ...
## $ biomass: num 77.5 78 75 214.4 215.4 ...
## $ day : Factor w/ 2 levels "0","15": 1 1 1 2 2 2 2 2 2 2 ...ryegrass2 data는 제초제를 뿌린 날과 뿌리고 15일 지난 날의 biomass 차이를 기록한 것입니다.
ryegrass2.LL.4 <- drm(biomass ~ dose, data = subset(ryegrass2, day == "15"), fct = LL.4())
summary(ryegrass2.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 1.64515 0.21750 7.5639 2.741e-07 ***
## c:(Intercept) 76.89206 5.65271 13.6027 1.440e-11 ***
## d:(Intercept) 223.24892 4.44878 50.1820 < 2.2e-16 ***
## e:(Intercept) 9.34112 0.91776 10.1782 2.349e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 9.391495 (20 degrees of freedom)plot(ryegrass2.LL.4, type = "all", xlab = "dose", ylab = "biomass", ylim = c(0, 250), broken = T)
abline(h = 50, lty = 2)
baseline이 0보다 위에 있습니다.
head(red.fescue) ; tail(red.fescue)## dose day biomass
## 1 0 0 45.0
## 2 0 0 69.0
## 3 0 16 137.0
## 4 0 16 102.0
## 5 0 16 101.4
## 6 87 16 139.7## dose day biomass
## 21 2805 16 19.5
## 22 2805 16 33.8
## 23 2805 16 33.2
## 24 5610 16 13.9
## 25 5610 16 18.9
## 26 5610 16 25.9summary(red.fescue)## dose day biomass
## Min. : 0 0 : 2 Min. : 13.90
## 1st Qu.: 87 16:24 1st Qu.: 36.25
## Median : 350 Median : 60.45
## Mean :1284 Mean : 69.76
## 3rd Qu.:1402 3rd Qu.: 98.78
## Max. :5610 Max. :145.20str(red.fescue)## 'data.frame': 26 obs. of 3 variables:
## $ dose : int 0 0 0 0 0 87 87 87 175 175 ...
## $ day : Factor w/ 2 levels "0","16": 1 1 2 2 2 2 2 2 2 2 ...
## $ biomass: num 45 69 137 102 101 ...red.fescue data는 제초제의 효능을 0일과 16일에 비교한 data 입니다.
red.fescue.LL.4 <- drm(biomass ~ dose, data = subset(red.fescue, day == "16"), fct = LL.4())
summary(red.fescue.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 1.93197 0.62841 3.0744 0.0059823 **
## c:(Intercept) 26.54978 8.63284 3.0754 0.0059679 **
## d:(Intercept) 127.37245 8.08135 15.7613 9.626e-13 ***
## e:(Intercept) 347.66456 88.34741 3.9352 0.0008186 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 17.97319 (20 degrees of freedom)plot(red.fescue.LL.4, broken = T, xlab = "dose", ylab = "biomass", ylim = c(0, 150))
abline(h = coef(red.fescue.LL.4)["c:(Intercept)"], col = "blue", lty = 2)
abline(h = mean(subset(red.fescue, day == "0")[["biomass"]]), lty = 2, col = "red")
# data.frame 내 변수를 [x] 로 추출 시 data.frame 형태로 뽑히고, [[x]]로 추출 시 vector 형태로 뽑힙니다. 만약 추출 변수로 평균 등의 연산을 할 경우 [[x]] 형태로 뽑아야 합니다.baseline이 day 0에서의 mean 보다 낮은 경우를 볼 수 있습니다. 그래프의 baseline은 해당 data의 산출과정에 대한 이해를 바탕으로 분석해야 합니다.
1.1.6 A hormesis effect on lettuce growth
head(lettuce) ; tail(lettuce)## conc weight
## 1 0.00 1.126
## 2 0.00 0.833
## 3 0.32 1.096
## 4 0.32 1.106
## 5 1.00 1.163
## 6 1.00 1.336## conc weight
## 9 10 0.716
## 10 10 0.683
## 11 32 0.560
## 12 32 0.488
## 13 100 0.375
## 14 100 0.344str(lettuce)## 'data.frame': 14 obs. of 2 variables:
## $ conc : num 0 0 0.32 0.32 1 1 3.2 3.2 10 10 ...
## $ weight: num 1.126 0.833 1.096 1.106 1.163 ...summary(lettuce)## conc weight
## Min. : 0.00 Min. :0.3440
## 1st Qu.: 0.49 1st Qu.:0.5907
## Median : 3.20 Median :0.7935
## Mean : 20.93 Mean :0.8261
## 3rd Qu.: 26.50 3rd Qu.:1.1035
## Max. :100.00 Max. :1.3360hormesis effect는 원래 독성 물질이지만 낮은 농도에서는 긍정적 효과가 있는 것을 말합니다. 암세포에게 항암제를 처리하면 아무 효과가 없다가 (NOEL, No Oveserved Effect Level), 농도가 약간 증가함에 따라 암세포의 생장이 오히려 증가하고 (Effect), 농도가 계속 증가하면 다시 아무 효과가 없다가 (NOAEL, NOAdverseEL), 암세포가 농도의존적으로 사멸하게 됩니다.
주로 농약, 항생제 등의 효능을 측정할 때 고려하는 factor입니다.
lettuce data는 배양액에 isobutul alcohol이 포함된 정도에 따라 생장을 비교한 data 입니다.
우선 scatterplot을 그려봅시다.
plot(y = lettuce$weight, x = lettuce$conc, log = "x", xlab = "conc", ylab = "weight")## Warning in xy.coords(x, y, xlabel, ylabel, log): 2 x values <= 0 omitted from
## logarithmic plot
x가 0인 경우 log 값을 취할 수 없기 때문에 제거되었습니다.
그래프의 형태는 J를 90도 대칭한 모습이고, 농도가 증가할수록 무게가 감소되는 경향을 보입니다.
Brain-Cousens hormesis model을 적용해 봅시다. 최저값은 0으로 고정해봅니다.
lettuce.BC.4 <- drm(weight ~ conc, data = lettuce, fct = BC.4())
summary(lettuce.BC.4)##
## Model fitted: Brain-Cousens (hormesis) with lower limit fixed at 0 (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 1.282812 0.049346 25.9964 1.632e-10 ***
## d:(Intercept) 0.967302 0.077123 12.5423 1.926e-07 ***
## e:(Intercept) 0.847633 0.436093 1.9437 0.08059 .
## f:(Intercept) 1.620703 0.979711 1.6543 0.12908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.1117922 (10 degrees of freedom)# b : e에서의 hill slope
# d : Emax
# e : EC50
# f : hormesis effect 매개변수
plot(lettuce.BC.4, broken = T)
f는 0 이상의 값을 가지며 0일 경우 hormesis 효과는 없습니다. 해당 결과에서 f는 1.62므로 hormesis effect가 존재합니다.
hormesis effect를 무시할 경우 EC00의 값은 실제보다 낮게 측정됩니다. hormesis effect에 관심이 없더라도, 그래프의 개형을 보고 hormesis effect가 존재하는지 고려하여야 합니다.
1.1.7 Nonlinear calibration
DRC를 위해서는 실험 결과를 non-linear regression model에 적합시킵니다. 이후 역회귀를 통해 n%의 효능을 가지기 위해서 필요한 약물의 농도를 계산하기도 합니다.
cell assay에서 하나의 target에 대한 약물을 처리할 시 DRC는 s자 모양의 커브를 가집니다. 동물이나 조직에서는 linear한 경우도 있습니다만 흔하지 않고 사실 동물실험부터는 DRC를 정확히 구하기 쉽지 않습니다.
head(nasturtium) ; tail(nasturtium)## conc rep wt
## 1 0.000 1 920
## 2 0.025 1 919
## 3 0.075 1 870
## 4 0.250 1 880
## 5 0.750 1 693
## 6 2.000 1 429## conc rep wt
## 37 0.025 6 850
## 38 0.075 6 875
## 39 0.250 6 810
## 40 0.750 6 591
## 41 2.000 6 257
## 42 4.000 6 221summary(nasturtium)## conc rep wt
## Min. :0.000 Min. :1.0 Min. : 128.0
## 1st Qu.:0.025 1st Qu.:2.0 1st Qu.: 430.5
## Median :0.250 Median :3.5 Median : 817.5
## Mean :1.014 Mean :3.5 Mean : 676.9
## 3rd Qu.:2.000 3rd Qu.:5.0 3rd Qu.: 873.8
## Max. :4.000 Max. :6.0 Max. :1017.0str(nasturtium)## 'data.frame': 42 obs. of 3 variables:
## $ conc: num 0 0.025 0.075 0.25 0.75 2 4 0 0.025 0.075 ...
## $ rep : int 1 1 1 1 1 1 1 2 2 2 ...
## $ wt : int 920 919 870 880 693 429 200 889 878 825 ...plot(nasturtium)
nasturtium.LL.3 <- drm(wt ~ conc, data = nasturtium, fct = LL.3())
summary(nasturtium.LL.3)##
## Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 1.350256 0.113557 11.891 1.518e-14 ***
## d:(Intercept) 897.862776 13.844884 64.852 < 2.2e-16 ***
## e:(Intercept) 1.576200 0.095694 16.471 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 55.55935 (39 degrees of freedom)plot(nasturtium.LL.3, xlim = c(0, 100))
# 3 parameter는 Emin을 0이라고 놓습니다.
# normalization을 거치지 않았다면 3 paramter를 쓸 이유는 없습니다.
# 여기서는 책에서 시키는대로 합니다.적합된 model에서 실제 x값과 모델에서 예측하는 x 값을 비교해보려면 backfit() 함수를 쓰면 됩니다.
backfit(nasturtium.LL.3)## Warning in log(exp(-tempVal/parmVec[5]) - 1): NaN이 생성되었습니다## dose Estimate
## [1,] 0.000 NaN
## [2,] 0.025 0.1152911
## [3,] 0.075 0.1522799
## [4,] 0.250 0.2418827
## [5,] 0.750 0.7209820
## [6,] 2.000 2.0729493
## [7,] 4.000 3.8933864# dose 0은 log0 값이 없으므로 여기서는 무시해버립니다.결과에서 보면 dose와 estimated dose가 비슷한 것을 알 수 있습니다.
이번에는 반응변수가 690, 693, 722인 경우를 생각해 봅니다.
각각의 변수에서 예측을 하면 ED()를 이용합니다.
ED(nasturtium.LL.3, 690, type = "absolute", interval = "delta")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:690 0.648194 0.070708 0.505174 0.791214# type = "absolute" : y값에 해당하는 x값을 추정합니다.
# type = "relative" : x값에 해당하는 y값을 추정합니다.아니면 3 포인트의 평균을 내서 진행해도 됩니다.
ED(nasturtium.LL.3, (690 + 693 + 722) / 3, type = "absolute", interval = "delta")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:1:701.666666666667 0.613385 0.069337 0.473137 0.753633샘플에 대해서 부트스트랩을 진행해 볼 수 있습니다.
nasturtium.boot.res1 <- boot((690 + 693 + 722) / 3,
statistic = function(simyVal) {
ED(nasturtium.LL.3, simyVal, type = "absolute", display = F)[1]
},
ran.gen = function(yVal, mle) {
rnorm(1, yVal, 55.55935 / sqrt(3))
},
sim = "parametric",
R = 1000)
# boot()
# 1번째 항목은 data
# statistic은 data 값을 측정된 농도 값으로 변환. 해당 값은 변화되는 값임
# ran.gen은 정규분포 가정 하에 추출된 샘플을 반환하는 함수
# sim은 ran.gen 사용시 필요함
# R은 반복 횟수, 일반적으로 1000이 쓰임
summary(nasturtium.boot.res1)## R original bootBias bootSE bootMed
## 1 1000 0.61338 -0.0016472 0.096441 0.6098895% CI를 구할 수 있습니다.
boot.ci(nasturtium.boot.res1, type = "basic")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = nasturtium.boot.res1, type = "basic")
##
## Intervals :
## Level Basic
## 95% ( 0.4157, 0.7982 )
## Calculations and Intervals on Original Scale해당 결과 도출된 95% CI 값은 직접 data points를 입력한 95% CI 보다 넓은 추정치를 보여줍니다. 부트스트랩 결과를 개선하기 위해서는 특정 point에서의 mean 값으로 진행하기보다 전체 데이터를 넣는 것이 당연히 좋습니다. nasturtium은 총 42개의 데이터로 이루어져 있는데 기존 42개에 우리가 원하는 3개의 포인트를 넣어 부트스트랩 해봅시다.
nasturtium.boot.res2 <- boot(c(nasturtium[["wt"]], 690, 693, 722),
statistic = function(yValues) {
ED(drm(head(yValues, -3) ~ conc, data = nasturtium, fct = LL.3()),
mean(tail(yValues, 3)),
type = "absolute", display = F)[1]
},
ran.gen = function(yVal, mle) {
rnorm(42 + 3, c(fitted(nasturtium.LL.3), 690, 693, 722), 55.55935)
},
sim = "parametric", R = 1000)
boot.ci(nasturtium.boot.res2, type = "basic")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = nasturtium.boot.res2, type = "basic")
##
## Intervals :
## Level Basic
## 95% ( 0.3718, 0.8289 )
## Calculations and Intervals on Original Scale위에서 진행한 결과 대비 95% 신뢰구간이 아주 살짝 넓어졌습니다. 이것은 nonlinear calibration 상에서 추가적인 분산의 발생이 거의 없었음을 의미합니다.