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)
--
This e-mail and its attachment(s) (if any) may contain confidential and/or
proprietary information and is intended for its addressee(s) only. Any
unauthorized use of the information contained herein (including, but not limited
to, alteration, reproduction, communication, distribution or any other form of
dissemination) is strictly prohibited. If you are not the intended addressee,
please notify the orginator promptly and delete this e-mail and its
attachment(s) (if any) subsequently.
Galapagos nor any of its affiliates shall be liable for direct, special,
indirect or consequential damages arising from alteration of the contents of
this message (by a third party) or as a result of a virus being passed on.
On May 20, 2009, at 11:20 AM, Hans Vermeiren wrote:> 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)Hans, You have three important factors changing here. The version of R, the version of drc and the versions of any relevant drc dependencies (alr3, lattice, magic, MASS, nlme, plotrix). I would first try to install the newer version of drc on the older R system (all else staying the same) and see what you get. Don't run update.packages() here, lest you change other things. Just install the newer version of drc. If you get the same results as the older version, then it might lead you to something in R or one of the package dependencies changing. If you get a different result, then it would lead to something in drc changing. You can also install the old version of drc on your more recent R system to see what you get, which might help to confirm behavior. The old source version of drc would be available from: http://cran.r-project.org/src/contrib/Archive/drc/ I also found a Windows binary of the old package here: http://cran.r-project.org/bin/windows/contrib/2.6/drc_1.4-2.zip I have also copied Christian Ritz, the drc package author/maintainer, so that he may be able to assist you further with the problem. HTH, Marc Schwartz
Hi Hans, I hope I can resolve your problems below (Marc, thank you very much for cc'ing me on your initial response!). Have a look at the following R lines: ## Fitting the model using drm() (from the latest version) m1<- drm(response ~ dose, data = d, fct = LL.4()) summary(m1) plot(m1) ## Checking the fit by using nls() ## (we have very good guesses for the parameter estimates) m2 <- nls(response ~ c + (d - c)/(1 + (dose/e)^b), data=d, start=list(b=-0.95, c=10, d=106, e=1.2745e-06)) summary(m2) The standard errors agree quite well. The minor discrepancies between to two fits are attributable to different numerical approximations of the variance-covariance matrix being used in drm() and nls(). So I would use the latest version of 'drc', especially for datasets with really small doses. One recent change to drm() was to incorporate several layers of scaling prior to estimation (as well as subsequent back scaling after estimation): 1) scaling of parameters with the same scale as the x axis 2) scaling of parameters with the same scale as the y axis 3) scaling of parameters in optim() The effect of scaling is to temporarily "convert" the dataset (and the model) to scales that are more convenient for the estimation procedure. Any feedback on this would be much appreciated. Therefore it should also not be necessary to manually do any scaling prior to using drm() (like what you did). Compare, for instance, your specification of drm() to mine above. Is this explanation useful?! Christian