Multiple Dose Response Curve Analysis
Cato. LEE
2020 6 26
library(drc)
library(devtools)
library(drcData) # devtools library 설치 후 install_github("DoseResponse/drcData") 로 설치
library(boot)
library(lmtest)
library(metafor)
library(sandwich)1.2 Analysis of multiple dose-response curves
두 개의 약물을 처리하여 반응 시킨 뒤 효능을 비교하기 위해서는 두 개의 그래프를 비교하는 과정이 필요합니다.
이 과정에서 두 그래프의 EC50를 비교하는 것도 좋은 예시가 될 수 있습니다.
우리가 진행할 첫 번째 예제에서는 두 종류의 약물에 대한 반복 실험을 1회 진행한 실험에 대한 그래프를 살펴볼 겁니다. 두 그래프에서 random variation을 분리할 수는 없습니다. 한 번만 실험을 했기 때문입니다. 이 상황에서는 두 약물 처리결과에 대해 차이점을 강하게 주장하기는 어렵습니다.
두 번째 예제에서는 두 종류의 약물에 대한 반복 실험을 2회 진행한 실험에 대해 살펴봅니다. 이 경우는 첫 번째 실험보다 차이점을 강하게 주장할 수 있습니다.
1.2.1 Effect of an herbicide mixture on Galium aparine
head(G.aparine)## dose drymatter treatment
## 1 0 1146 1
## 2 0 1005 1
## 3 0 756 1
## 4 0 1108 1
## 5 0 956 1
## 6 0 989 1summary(G.aparine)## dose drymatter treatment
## Min. : 0.00 Min. : 123.0 1:140
## 1st Qu.: 5.53 1st Qu.: 495.0 2:100
## Median : 38.01 Median : 773.0
## Mean : 270.00 Mean : 715.8
## 3rd Qu.: 253.38 3rd Qu.: 937.2
## Max. :1621.60 Max. :1377.0str(G.aparine)## 'data.frame': 240 obs. of 3 variables:
## $ dose : num 0 0 0 0 0 0 0 0 0 0 ...
## $ drymatter: int 1146 1005 756 1108 956 989 1109 867 864 997 ...
## $ treatment: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...G.aparine은 두 종류의 제초제를 농도별로 처리하였을 때 식물의 중량이 어떻게 변하는지 측정한 자료입니다.
일단, 일반적인 4 parametric concentration curve를 그려봅시다.
G.aparine.LL.4 <- drm(drymatter ~ dose, curveid = treatment, data = G.aparine, fct = LL.4(),
pmodels = list( b = ~ treatment, c = ~ treatment, d = ~ 1, e = ~ treatment))
G.aparine.LL.4.d <- drm(drymatter ~ dose, curveid = treatment, data = G.aparine, fct = LL.4(),
pmodels = data.frame(b = treatment, c = treatment, d = 1, e = treatment))
# curveid = curveid는 numeric vector나 factor를 통해 그룹을 나누는 기준이 됩니다.
# pmodels = b : hillslope / c : lower limit / d : upper limit / e : EC50 에 대해 curveid 에서 나눈 그룹 별 수치를 보여주는 여부를 정합니다. ~ 1을 작성할 경우 해당 parameter는 모든 커브 간에 서로 공유됩니다.
# pmoels을 list로 받을 경우 function 형태로 받습니다. dataframe 형태로 받을 때는 vector 형태로 받습니다.
# dataframe을 사용하면 각 그래프의 parameter 값을 알 수 있고, list로 받으면 두 그래프의 차이를 알 수 있습니다.
# 어떤 경우든 1을 넣으면 global share로 인정되고 curveid로 나눈 column 명을 입력하면 각각의 결과를 산출합니다.
plot(G.aparine.LL.4, broken = T, xlab = "Dose", ylab = "Mass")
summary() 함수를 이용해서 용량반응곡선의 parameter 값들을 요약할 수 있습니다.
summary(G.aparine.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 1.61291 0.33330 4.8392 2.373e-06 ***
## b:treatment2 0.13809 0.37301 0.3702 0.7115667
## c:(Intercept) 509.50242 23.25903 21.9056 < 2.2e-16 ***
## c:treatment2 -357.58391 34.71802 -10.2997 < 2.2e-16 ***
## d:(Intercept) 984.88794 12.63336 77.9593 < 2.2e-16 ***
## e:(Intercept) 50.80033 7.87857 6.4479 6.468e-10 ***
## e:treatment2 42.64599 10.94359 3.8969 0.0001274 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 111.403 (233 degrees of freedom)결과를 보면 (parameter)는 첫 번째 약물에 대한 값이고, (parameter)2는 첫 번째 약물 대비 두 번째 약물에 대한 결과입니다. 순서는 알파벳 순서나 숫자의 증가에 따라 매겨집니다.
parameter b를 살펴보면 1번 곡선의 hill slope는 1.61 입니다. 1번 곡선 대비 2번 곡선 기울기의 차이는 0.14로 t-test 결과 유의미한 차이가 없습니다. 따라서 두 그래프의 hill slope는 같다고 볼 수 있습니다.
반면, parameter c, parameter e는 두 곡선 간에 유의한 차이를 보입니다. parameter c lower limit는 1번곡선 대비 2번 곡선이 358만큼 낮고 이는 유의한 차이입니다. parameter e EC50는 1번곡선 대비 2번 곡선이 42만큼 높고 이 또한 유의한 차이입니다.
보통 약물의 효능 차이를 볼 때 EC50를 많이 봅니다, 이 그래프에서는 1번이 2번보다 EC50 값이 낮으므로 더 효능이 우수한 제초제라고 볼 수도 있겠지만 식물을 더 많이 죽이는 제초제는 2번입니다. 따라서 생물학 데이터를 볼 때는 EC50 외에도 약물이 가지는 최대 효능도 같이 비교해야 합니다.
위에서 언급하였지만, 두 곡선의 parameter 간의 차이가 아닌 각각의 값을 보고 싶다면 pmodels = data.frame() 으로 받으면 됩니다.
summary(G.aparine.LL.4.d)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:1 1.61291 0.33330 4.8392 2.372e-06 ***
## b:2 1.75100 0.20392 8.5869 1.311e-15 ***
## c:1 509.50367 23.25885 21.9058 < 2.2e-16 ***
## c:2 151.91840 26.00899 5.8410 1.734e-08 ***
## d:(Intercept) 984.88779 12.63335 77.9594 < 2.2e-16 ***
## e:1 50.80009 7.87851 6.4479 6.467e-10 ***
## e:2 93.44626 8.11091 11.5211 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 111.403 (233 degrees of freedom)1.2.2 Glyphostate and bentazone treatment of Sinapis alba
head(S.alba.comp)## exp herbicide dose drymatter Tf area Fo Fm
## 1 ben1 bentazone 0 4.1 200 31200 278 1662
## 2 ben1 bentazone 0 3.4 230 30600 278 1670
## 3 ben1 bentazone 0 2.6 210 27400 299 1646
## 4 ben1 bentazone 0 3.5 260 34600 288 1715
## 5 ben1 bentazone 0 4.3 200 31000 272 1651
## 6 ben1 bentazone 0 4.2 240 31400 286 1681str(S.alba.comp)## 'data.frame': 141 obs. of 8 variables:
## $ exp : Factor w/ 4 levels "ben1","ben2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ herbicide: Factor w/ 2 levels "bentazone","glyphosate": 1 1 1 1 1 1 1 1 1 1 ...
## $ dose : int 0 0 0 0 0 0 0 0 10 10 ...
## $ drymatter: num 4.1 3.4 2.6 3.5 4.3 4.2 4.1 4.5 3.5 3.8 ...
## $ Tf : int 200 230 210 260 200 240 250 210 1000 160 ...
## $ area : int 31200 30600 27400 34600 31000 31400 27000 29600 42800 26600 ...
## $ Fo : int 278 278 299 288 272 286 294 280 301 289 ...
## $ Fm : int 1662 1670 1646 1715 1651 1681 1620 1659 1510 1564 ...summary(S.alba.comp)## exp herbicide dose drymatter Tf
## ben1:36 bentazone :72 Min. : 0.0 Min. :0.500 Min. : 17.0
## ben2:36 glyphosate:69 1st Qu.: 10.0 1st Qu.:0.900 1st Qu.: 160.0
## gly1:36 Median : 40.0 Median :1.600 Median : 200.0
## gly2:33 Mean : 174.3 Mean :2.171 Mean : 256.4
## 3rd Qu.: 160.0 3rd Qu.:3.700 3rd Qu.: 250.0
## Max. :1280.0 Max. :4.800 Max. :1200.0
## area Fo Fm
## Min. : 0 Min. :259.0 Min. : 378
## 1st Qu.: 1400 1st Qu.:288.0 1st Qu.:1217
## Median :23400 Median :303.0 Median :1511
## Mean :18765 Mean :370.4 Mean :1371
## 3rd Qu.:28200 3rd Qu.:418.0 3rd Qu.:1653
## Max. :52800 Max. :857.0 Max. :1858exp는 각각의 실험군을 의미합니다. ben1, ben2, gly1, gly2가 있는데 개별적인 실험 batch라고 생각을 하면 됩니다.
하나의 제초제 당 2개의 실험이 존재하는데 우리는 제초제에 따라 실험을 하나로 모아 살펴볼 필요가 있습니다.
그렇게 하면 random variation의 크기를 줄일 수 있고 통계적 이해에도 도움이 됩니다.
1.2.2.1 A joint dose-response model
일단 두 종류의 제초제에 따라서 농도의존곡선을 도출합니다.
S.alba.LL.4 <- drm(drymatter ~ dose, curveid = herbicide, data = S.alba.comp, fct = LL.4())
plot(S.alba.LL.4, xlab = "dose", ylab = "drymatter", broken = T)
summary(S.alba.LL.4)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:bentazone 4.421705 1.376978 3.2112 0.001658 **
## b:glyphosate 1.606429 0.322669 4.9786 1.949e-06 ***
## c:bentazone 0.679894 0.092323 7.3643 1.673e-11 ***
## c:glyphosate 0.894469 0.138322 6.4666 1.756e-09 ***
## d:bentazone 4.018895 0.107033 37.5481 < 2.2e-16 ***
## d:glyphosate 3.874354 0.122103 31.7302 < 2.2e-16 ***
## e:bentazone 21.293475 1.348860 15.7863 < 2.2e-16 ***
## e:glyphosate 46.245932 6.379062 7.2496 3.074e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.4952302 (133 degrees of freedom)약물마다의 parameter 확인이 가능합니다. 앞에서 언급한 pmodels = list()를 넣으면 각 parameter마다 유의적인 차이가 존재하는지 확인할 수 있습니다.
S.alba.LL.4.p <- drm(drymatter ~ dose, curveid = herbicide, data = S.alba.comp, fct = LL.4(),
pmodels = list(b = ~ herbicide, c = ~ herbicide, d = ~ herbicide, e = ~ herbicide))
summary(S.alba.LL.4.p)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 4.419951 1.374743 3.2151 0.0016369 **
## b:herbicideglyphosate -2.812699 1.412168 -1.9918 0.0484457 *
## c:(Intercept) 0.679807 0.092317 7.3638 1.677e-11 ***
## c:herbicideglyphosate 0.214970 0.166245 1.2931 0.1982206
## d:(Intercept) 4.018973 0.107036 37.5480 < 2.2e-16 ***
## d:herbicideglyphosate -0.144456 0.162387 -0.8896 0.3752984
## e:(Intercept) 21.294644 1.348997 15.7855 < 2.2e-16 ***
## e:herbicideglyphosate 24.931083 6.514514 3.8270 0.0001988 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.4952302 (133 degrees of freedom)coef(summary(S.alba.LL.4.p))## Estimate Std. Error t-value p-value
## b:(Intercept) 4.4199508 1.3747428 3.2151112 1.636888e-03
## b:herbicideglyphosate -2.8126987 1.4121683 -1.9917589 4.844570e-02
## c:(Intercept) 0.6798066 0.0923174 7.3637974 1.677151e-11
## c:herbicideglyphosate 0.2149701 0.1662450 1.2930927 1.982206e-01
## d:(Intercept) 4.0189726 0.1070357 37.5479689 5.889097e-73
## d:herbicideglyphosate -0.1444557 0.1623866 -0.8895786 3.752984e-01
## e:(Intercept) 21.2946442 1.3489974 15.7855339 1.394524e-32
## e:herbicideglyphosate 24.9310827 6.5145140 3.8270058 1.987790e-04나눠진 집단 간 특정 parameter를 비교하고 싶을 땐 comParm()을 활용합니다.
compParm(S.alba.LL.4, strVal = "e", operator = "-")##
## Comparison of parameter 'e'
##
## Estimate Std. Error t-value p-value
## bentazone-glyphosate -24.9525 6.5201 -3.827 0.0001988 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1compParm(S.alba.LL.4, strVal = "e", operator = "/")##
## Comparison of parameter 'e'
##
## Estimate Std. Error t-value p-value
## bentazone/glyphosate 0.460440 0.069889 -7.7202 2.473e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1strVal은 어떤 parameter를 보고 싶은지 지정해주고 operator는 “-”, “/”를 통해 차이 또는 비율을 산출합니다. 그리고 둘 사이에 유의한 차이가 있을 경우 p value에서 표기해 줍니다.
둘의 차이를 구한 경우와 비를 구한 경우의 p value 값이 항상 같지는 않습니다. 비율에 대한 표준오차를 구하는데 있어 delta method를 사용하기 때문입니다. 이는 비선형적인 변화에 대한 추정치를 제공하는 방법입니다.
1.2.2.2 Fitting separate dose-response models
two-step 접근법은 농도-반응 분석에서 nonlinear mixed effect model에 대한 대안이 될 수 있습니다.
첫 번째 스텝으로 각각의 sub-experiment에 대한 dose-response curve를 그립니다. 해당 curve 모델들 사이에서 모든 parameter가 동일할 필요는 없습니다. 추정된 parameter는 표준오차와 함께 새로운 dataset에 처리 조건과 함께 결합됩니다.
두 번째 스텝에서 추정치들은 고정 효과와 랜덤 효과를 설명하기 위해 설계 된 실험 디자인이 바탕이 된 메타분석의 반응변수로 사용됩니다.
일반적인 linear mixed model과 다르게, residual standard error는 메타분석 랜덤 효과 모델에서 추정되지 않습니다. 대신 이에 해당하는 standard error로 추정됩니다.
S.alba.comp data는 bentazone에 대한 실험이 ben1, ben2로 나뉘어 있고 hlyphostate에 대한 실험 역시 gly1, gly2로 나뉘어 있습니다. 각각의 sub-experiments에 대해 4 parametric model을 만듭니다.
S.alba.LL.4.b1 <- drm(drymatter ~ dose, data = subset(S.alba.comp, exp == "ben1"), fct = LL.4())
S.alba.LL.4.b2 <- drm(drymatter ~ dose, data = subset(S.alba.comp, exp == "ben2"), fct = LL.4())
S.alba.LL.4.g1 <- drm(drymatter ~ dose, data = subset(S.alba.comp, exp == "gly1"), fct = LL.4())
S.alba.LL.4.g2 <- drm(drymatter ~ dose, data = subset(S.alba.comp, exp == "gly2"), fct = LL.4())다음으로 각각의 모델에서 ED50와 ED50에 해당하는 표준오차를 뽑아냅니다.
ed50.estimates <- matrix(NA, 4, 2) #빈 행렬을 만들어 줍니다
ed50.estimates## [,1] [,2]
## [1,] NA NA
## [2,] NA NA
## [3,] NA NA
## [4,] NA NAsummary(S.alba.LL.4.b1)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 5.125377 1.042018 4.9187 2.517e-05 ***
## c:(Intercept) 0.681752 0.087984 7.7486 7.778e-09 ***
## d:(Intercept) 3.806283 0.102051 37.2977 < 2.2e-16 ***
## e:(Intercept) 29.264785 2.065501 14.1684 2.444e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.3450769 (32 degrees of freedom)# ED50는 parameter e이므로 4행의 1, 2열이 측정값과 표준오차입니다.
# coef() 함수로 summary()함수를 묶어주면 summary 결과의 계수들을 dataframe 형식으로 대체해서 뽑아낼 수 있습니다.
ed50.estimates[1, ] <- coef(summary(S.alba.LL.4.b1))[4, 1:2]
ed50.estimates[2, ] <- coef(summary(S.alba.LL.4.b2))[4, 1:2]
ed50.estimates[3, ] <- coef(summary(S.alba.LL.4.g1))[4, 1:2]
ed50.estimates[4, ] <- coef(summary(S.alba.LL.4.g2))[4, 1:2]
ed50.estimates## [,1] [,2]
## [1,] 29.26479 2.065501
## [2,] 18.85085 3.967966
## [3,] 25.56620 5.793800
## [4,] 62.73827 6.390966최종적으로 matrix를 dataframe으로 변경하여 column 이름을 부여하고 treatment를 같이 기록하게끔 컬럼을 추가합니다.
ed50.estimates <- as.data.frame(ed50.estimates)
names(ed50.estimates) <- c("est", "se")
ed50.estimates <- cbind(ed50.estimates, treatment = c("bentazone", "bentazone", "glyphosate", "glyphosate"))
ed50.estimates## est se treatment
## 1 29.26479 2.065501 bentazone
## 2 18.85085 3.967966 bentazone
## 3 25.56620 5.793800 glyphosate
## 4 62.73827 6.390966 glyphosatemetafor package의 rma() 함수를 활용하여 메타분석을 할 수 있습니다. 위에서 만든 ed50.estimates dataframe을 사용하겠습니다.
rma()의 첫 번째 요소에 ED50 값에 해당하는 “est”를 넣고, 두 번째 요소에는 “est”에 해당하는 분산인 “se”^2을 넣습니다.
표준오차가 작을수록 해당 est에 가중치를 부여합니다.
S.alba.ED50 <- rma(est, se^2, mods = ~ treatment, data = ed50.estimates)
summary(S.alba.ED50)##
## Mixed-Effects Model (k = 4; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -8.7257 17.4514 23.4514 19.5308 47.4514
##
## tau^2 (estimated amount of residual heterogeneity): 325.2213 (SE = 348.0325)
## tau (square root of estimated tau^2 value): 18.0339
## I^2 (residual heterogeneity / unaccounted variability): 95.38%
## H^2 (unaccounted variability / sampling variability): 21.62
## R^2 (amount of heterogeneity accounted for): 3.50%
##
## Test for Residual Heterogeneity:
## QE(df = 2) = 23.9885, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.1262, p-val = 0.2886
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 24.1470 12.9447 1.8654 0.0621 -1.2241 49.5180 .
## treatmentglyphosate 19.8187 18.6751 1.0612 0.2886 -16.7838 56.4213
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1모델 결과를 살펴보면 glyphostate의 ed50 추정 치는 bentazone 보다 19.8만큼 더 큽니다. 하지만 이러한 차이는 유의적이지는 않습니다.
joint model과 seperate model은 때때로 다른 결과 값을 도출합니다. sub-experiments 값을 더 많이 포함하는 것은 적절한 model의 방법이긴 하지만 데이터의 변동이 더 많이 포함되므로 유의한 차이를 선언하기 어려워지곤 합니다.
적절한 model에 대해 수학적 타당성을 입증하거나 해당 실험에 대한 도메인 지식을 바탕으로 모델을 선정하는 것이 필요합니다.