# using drc package, calculating LD50 of 2-mercaptobenzimidazole in rats # written by Minato Nakazawa on 23 June 2014. # data source: http://dra4.nihs.go.jp/mhlw_data/home/paper/paper583-39-1a.html # The report used probit analysis or Behrens-Karber method, and the resulted # LD50s were 208 for Females, 218 (95%CI: 166-287) for Males # The log-logit method [fct=LL.2()] has been applied here. if (require(drc)==FALSE) { install.packages("drc"); library(drc) } rats <- data.frame( sex = factor(c(rep(1, 7), rep(2, 7)), labels=c("M", "F")), dose = rep(c(0, 79, 119, 178, 267, 400, 600), 2), ndead = c(0, 0, 0, 1, 4, 5, 5, 0, 0, 0, 1, 5, 5, 5), ntotal = rep(5, 14)) mx <- drm(ndead/ntotal ~ dose, curveid=sex, weights=ntotal, data=rats, fct=LL.2(), type="binomial") summary(mx) ED(mx, 50, interval="delta") # Estimate Std. Error Lower Upper # F:50 186.636 60.373 68.306 304.97 # M:50 218.013 22.809 173.308 262.72 plot(mx) # additional information on 20170623 # using glm() and car::deltaMethod(), it's also possible # mlogit <- glm(ndead/ntotal ~ sex + dose, weights=ntotal, data=rats, family=binomial) summary(mlogit) mprobit <- glm(ndead/ntotal ~ sex + dose, weights=ntotal, data=rats, family=binomial(probit)) summary(mprobit) library(car) deltaMethod(mlogit, "-(Intercept)/dose") # male's LD50 using logit deltaMethod(mlogit, "-(Intercept+sexF)/dose") # female's LD50 using logit deltaMethod(mprobit, "-(Intercept)/dose") # male's LD50 using probit deltaMethod(mprobit, "-(Intercept+sexF)/dose") # male's LD50 using probit # Estimate SE 2.5 % 97.5 % # -(Intercept)/dose 223.5115 21.15981 182.039 264.9839 # Estimate SE 2.5 % 97.5 % # -(Intercept + sexF)/dose 204.3379 19.99189 165.1546 243.5213 # Estimate SE 2.5 % 97.5 % # -(Intercept)/dose 222.9661 19.19045 185.3535 260.5787 # Estimate SE 2.5 % 97.5 % # -(Intercept + sexF)/dose 204.479 18.70611 167.8157 241.1423 # # If you want to apply log-logit model (same as LL.2()), ... # rats$logdose <- ifelse(rats$dose<=0, NA, log(rats$dose)) mloglogit <- glm(ndead/ntotal ~ sex + logdose, weight=ntotal, data=rats, family=binomial) summary(mloglogit) exp(deltaMethod(mloglogit, "-(Intercept)/logdose")) # males exp(deltaMethod(mloglogit, "-(Intercept+sexF)/logdose")) # females # Estimate SE 2.5 % 97.5 % # -(Intercept)/logdose 218.0088 1.103206 179.8326 264.2895 # Estimate SE 2.5 % 97.5 % # -(Intercept + sexF)/logdose 199.41 1.09814 165.9812 239.5713