## Fitting a four-parameter Weibull (type 1) model
terbuthylazin.m1 <- drm(rgr~dose, data = terbuthylazin, fct = W1.4())
summary(terbuthylazin.m1)
## Fitting a first-order multistage model
## to data from BMDS by EPA
bmds.ex1 <- data.frame(ad.dose=c(0,50,100), dose=c(0, 2.83, 5.67), num=c(6,10,19), total=c(50,49,50))
bmds.ex1.m1<-drm(num/total~dose, weights=total, data=bmds.ex1, fct=W2.4(fixed=c(1,NA,1,NA)), type="binomial")
plot(bmds.ex1.m1, ylim = c(0.05, 0.4), log = "")
modelFit(bmds.ex1.m1) # same as in BMDS
summary(bmds.ex1.m1) # same background estimate as in BMDS
logLik(bmds.ex1.m1)
## BMD estimate identical to BMDS result
## BMDL estimate differs from BMDS result (different method)
ED(bmds.ex1.m1, 10, ci="delta")
## Better fit
bmds.ex1.m2<-drm(num/total~dose, weights=total, data=bmds.ex1, fct=W1.4(fixed=c(-1,NA,1,NA)), type="binomial")
modelFit(bmds.ex1.m2)
summary(bmds.ex1.m2)
plot(bmds.ex1.m2, ylim = c(0.05, 0.4), log = "", add = TRUE, lty = 2)
ED(bmds.ex1.m2, 50, ci = "delta")
Run the code above in your browser using DataLab