Hello,
We use drc to fit dose-response curves, recently we discovered that
there are quite different standard error values returned for the same
dataset depending on the drc-version / R-version that was used (not
clear which factor is important)
On R 2.9.0 using drc_1.6-3 we get an IC50 of 1.27447 and a standard
error on the IC50 of 0.43540
Whereas on R 2.7.0 using drc_1.4-2 the IC50 is 1.2039e+00 and the
standard error is 3.7752e-03
Normally I would use the most recent version (both R and drc library)
but it seems to me that a standard error of 0.4 on a mean of 1.2 is
too
big, so I trust the values we get with the older versions more
Has anyone suggestions on
- how to solve these discrepancies, if possible
- how to calculate which one of the 2 solutions is the correct one?
Thanks a lot,
Hans Vermeiren
Demo (on a windows machine, while the issue was actually discovered on
our ubuntu linux server):
1)
sessionInfo()
R version 2.7.0 (2008-04-22)
i386-pc-mingw32
locale:
LC_COLLATE=Dutch_Belgium.1252;LC_CTYPE=Dutch_Belgium.
1252;LC_MONETARY=Du
tch_Belgium.1252;LC_NUMERIC=C;LC_TIME=Dutch_Belgium.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] drc_1.4-2 plotrix_2.4-2 nlme_3.1-89 MASS_7.2-41
lattice_0.17-6
[6] alr3_1.1.7
loaded via a namespace (and not attached):
[1] grid_2.7.0
d<-data.frame(dose=c(2.00e-05, 4.00e-06, 8.00e-07, 1.60e-07, 3.20e-08,
6.40e-09, 1.28e-09, 2.56e-10, 5.10e-11, 1.00e-11, 2.00e-05, 4.00e-06,
8.00e-07, 1.60e-07, 3.20e-08, 6.40e-09, 1.28e-09, 2.56e-10, 5.10e-11,
1.00e-11),
response=c(97.202,81.670,47.292,16.924, 16.832, 6.832, 11.118,
1.319, 5.495, -3.352, 102.464, 83.114, 50.631, 22.792, 18.348,
19.066, 27.794, 14.682, 11.992, 12.868))
m<- drm(response ~ (log10(dose*1e6)), data = d, fct = l4(fixed =
c(NA,NA,NA,NA), names = c("hs", "bottom", "top", "ec50")), logDose =
10,
control = drmc(useD = T))
summary(m)
results in:
Model fitted: Log-logistic (ED50 as parameter) (4 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
hs:(Intercept) -9.8065e-01 2.5821e-03 -3.7979e+02 2.248e-33
bottom:(Intercept) 1.0955e+01 2.2546e-02 4.8591e+02 4.364e-35
top:(Intercept) 1.0502e+02 9.0935e-02 1.1549e+03 4.210e-41
ec50:(Intercept) 1.2039e+00 3.7752e-03 3.1890e+02 3.681e-32
Residual standard error: 7.026655 (16 degrees of freedom)
=
=
======================================================================
===========================================
2)
sessionInfo()
R version 2.9.0 (2009-04-17)
i386-pc-mingw32
locale:
LC_COLLATE=Dutch_Belgium.1252;LC_CTYPE=Dutch_Belgium.
1252;LC_MONETARY=Du
tch_Belgium.1252;LC_NUMERIC=C;LC_TIME=Dutch_Belgium.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] drc_1.6-3 plotrix_2.5-5 nlme_3.1-90 MASS_7.2-46
magic_1.4-4 abind_1.1-0 lattice_0.17-22 alr3_1.1.7
loaded via a namespace (and not attached):
[1] grid_2.9.0 tools_2.9.0
d<-data.frame(dose=c(2.00e-05, 4.00e-06, 8.00e-07, 1.60e-07, 3.20e-08,
6.40e-09, 1.28e-09, 2.56e-10, 5.10e-11, 1.00e-11, 2.00e-05, 4.00e-06,
8.00e-07, 1.60e-07, 3.20e-08, 6.40e-09, 1.28e-09, 2.56e-10, 5.10e-11,
1.00e-11),
response=c(97.202,81.670,47.292,16.924, 16.832, 6.832, 11.118,
1.319, 5.495, -3.352, 102.464, 83.114, 50.631, 22.792, 18.348,
19.066, 27.794, 14.682, 11.992, 12.868))
m<- drm(response ~ (log10(dose*1e6)), data = d, fct = l4(fixed =
c(NA,NA,NA,NA), names = c("hs", "bottom", "top", "ec50")), logDose =
10,
control = drmc(useD = T))
summary(m)
gives:
Model fitted: Log-logistic (ED50 as parameter) (4 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
hs:(Intercept) -0.95266 0.25778 -3.69564 0.0020
bottom:(Intercept) 10.97437 2.24421 4.89009 0.0002
top:(Intercept) 106.38373 9.98378 10.65565 1.127e-08
ec50:(Intercept) 1.27447 0.43540 2.92712 0.0099
Residual standard error:
7.020175 (16 degrees of freedom)