drc results differ for different versions

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

drc results differ for different versions

Hans Vermeiren
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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results differ for different versions

Marc Schwartz-3
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results differ for different versions

Hans Vermeiren
Hello
Thanks a lot Marc, for the suggestion to explore the issue a bit more
systematically
So I did and the conclusion is that with the old drc 1.4-2, I get a
SE=0.003, with the new drc 1.5-2, I get a SE=0.4
irrespective of the R version or the version of the packages drc depends
on

I hope somebody can help me further, I still have the feeling that a SE
of 0.003 on an IC50 of 1.2, for a reasonable good fit is more realistic
than a SE of 0.4 on an IC50 of 1.2

Regards,
Hans

Here are the results of the following test script:
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:
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.5-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  tools_2.7.0

RESULT: ec50:(Intercept)=1.27447   SE=0.43541
CONCLUSION: R 2.7.0 with recent drc 1.5-2 (older dependencies) gives
SE=0.43541
========================================================================
=======
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.4-2       plotrix_2.5-5   nlme_3.1-90     MASS_7.2-46    
[5] lattice_0.17-22 alr3_1.1.7    

loaded via a namespace (and not attached):
[1] grid_2.9.0  tools_2.9.0
RESULT: ec50:(Intercept)    SE=1.2039e+00
CONCLUSION: R 2.9.0 with old drc 1.4-2 (newer dependencies) gives
SE=0.003


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

--
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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results differ for different versions

Christian Ritz-3
In reply to this post by Hans Vermeiren
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results differ for different versions

Hans Vermeiren
Yes, thanks that's very useful
Apart from checking the fit with nls() as you suggested, I've also used Prism, which gave the following results

Equation 1
Best-fit values
     BOTTOM 10.96
     TOP 106.4
     LOGEC50 -5.897
     HILLSLOPE 0.9501
     EC50 1.2670e-006
Std. Error
     BOTTOM 2.196
     TOP 9.337
     LOGEC50 0.1439
     HILLSLOPE 0.2270
95% Confidence Intervals
     BOTTOM 6.301 to 15.61
     TOP 86.62 to 126.2
     LOGEC50 -6.202 to -5.592
     HILLSLOPE 0.4689 to 1.431
     EC50 6.2750e-007 to 2.5560e-006
Goodness of Fit
     Degrees of Freedom 16
     R² 0.9622
     Absolute Sum of Squares 787.5
     Sy.x 7.015
Data
     Number of X values 20
     Number of Y replicates 1
     Total number of values 20
     Number of missing values 0

In other words: also in line with the drc 1.6-3 and nls() results
As for the scaling: yes this is useful because I can't predict whether concentrations are in molar, micromolar,..., right now I indeed scaled dose-values "manually", it's better/ more robust when the drm-function takes care of that
I suppose this also means I don't have to do the log transformation anymore?
Thanks (both of you) for your swift feedback

Hans

-----Original Message-----
From: Christian Ritz [mailto:[hidden email]]
Sent: vrijdag 22 mei 2009 11:30
To: Hans Vermeiren
Cc: [hidden email]; [hidden email]
Subject: Re: [R] drc results differ for different versions

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

--
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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results differ for different versions

Christian Ritz-3
Yes, you're right: taking logarithms is no longer needed!

Christian

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: drc results for binomial data

PatGauthier
I've been having trouble making sense of the drc results for my binomial response toxicity data. Firstly, the standard errors are far too large for how well the data fit the log-logistic model, particularly compared to other methods of LC50 estimation (e.g. probit and trimmed spearman-karber). Secondly, the LC50 estimates produced do not even fit with the DRC curve (see example below). I've tried all sort of log-logistic functions with varying fixed parameters (i.e. upper and lower limits for mortality data) and still fail to produce reliable estimates. Furthermore, my data response curves exhibit asymmetry so I feel the need to use the 5 parm function. There seems to be some inconsistency here which is causing me a serious headache.


#---------------data 1------------------#
n<-c(20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20)
x<-c(256,256,256,847,847,847,2208,2208,2208,3565,3565,3565,5756,5756,5756,7313,7313,7313,15000,15000,15000)
r<-c(0,0,0,0,2,0,1,3,4,7,6,9,17,13,18,18,17,20,20,20,20)
p<-r/n
d<-data.frame(x=x,n=n,r=r,p=p)

fct <- LL2.5(fixed = c(NA,0,1,NA,NA),names = c("Hillslope", "Lower Limit", "Upper Limit", "LC50", "Asymmetry"))
drc <- drm(p~log(x), weight=n, data=d, fct=fct, type="binomial", logDose=exp(1))
summary(drc)
plot(drc)
est1<-8.67
se1<-0.14
points(est1,0.5)
arrows(est1,0.5,est1-se1,angle=90,code=2,length=0.025)
arrows(est1,0.5,est1+se1,angle=90,code=2,length=0.025)

#---------------data 2------------------#
n<-c(20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20)
x<-c(798,798,798,6770,6770,6770,12130,12130,12130,21736,21736,21736,38948,38948,38948,69791,69791,69791)
r<-c(0,0,0,4,3,3,5,8,7,12,12,7,17,15,17,20,20,19)
p<-r/n
d2<-data.frame(x=x,n=n,r=r)
d2

fct <- LL2.5(fixed = c(NA,0,1,NA,NA),names = c("Hillslope", "Lower Limit", "Upper Limit", "LC50", "Asymmetry"))
drc2 <- drm(p~log(x), weight=n, data=d2, fct=fct, type="binomial", logDose=exp(1))
summary(drc2)
plot(drc2)
est2<-10.56
se2<-0.271
points(est2,0.5)
arrows(est2,0.5,est2-se2,angle=90,code=2,length=0.025)
arrows(est2,0.5,est2+se2,angle=90,code=2,length=0.025)

The idea is to use the drc package for these types of data to produce isoboles. Has anyone else encountered this issue before and can help me? OR I noticed Dr. Ritz has responded to this thread, do you have an suggestions? As far as I know, I have the latest version of every package for R.

Sincerely,

Patrick


Reply | Threaded
Open this post in threaded view
|

Re: drc results for binomial data

PatGauthier
P.S.
Here are the results I get for;

#---------------data 1------------------#

Model fitted: Generalised log-logistic (log(ED50) as parameter) (3 parms)

Parameter estimates:

                               Estimate   Std. Error    t-value      p-value
Hillslope:(Intercept)   -6.05453      2.18455    -2.77152    0.0056
LC50:(Intercept)         8.67707      0.14016   61.90663   0.0000
Asymmetry:(Intercept)  0.32382   0.16608   1.94978    0.0512


#---------------data 2------------------#
Model fitted: Generalised log-logistic (log(ED50) as parameter) (3 parms)

Parameter estimates:

                                 Estimate      Std. Error    t-value     p-value
Hillslope:(Intercept)      -3.98218     1.91035     -2.08453   0.0371
LC50:(Intercept)           10.55736    0.27120     38.92901  0.0000
Asymmetry:(Intercept)  0.26840     0.18244      1.47120   0.1412

Pat