[R] Estimating LC50 from a Weibull distribution

Greg chaoborid at gmail.com
Sun Mar 22 07:19:51 CET 2009


I am attempting to estimate LC50 (analogous to LD50, but uses exposure
concentration rather than dose) by fitting a Weibull model; but I
can't seem to get it to work.  From what I can gather, I should be
using survreg() from the survival package.  The survreg() function
relies on time-to-event data; my data result from 96 h exposures
(i.e., dead or alive after a fixed period; 96 h).  I've tried the
following (doesn't work):

> conc <- c(10.3, 10.8, 11.6, 13.2, 15.8, 20.1) # Exposure concentrations
> orign <- c(76, 79, 77, 76, 78, 77) # Original number of subjects
> ndead <- c(16, 22, 40, 69, 78, 77) # Number dead after 96 h
> d <- data.frame(conc=conc, orign=orign, ndead=ndead)
> d$prop <- d$ndead/d$orign  # Calculate proportion dead after 96 h
# Adjust for 100% mortalities
> d$prop[d$prop==1.00] <- 1-(1/(2*d$orign[d$prop==1.00]))
> > fit <- survreg(Surv(d$prop) ~ d$conc, dist="weibull")
> summary(fit)

Call:
survreg(formula = Surv(d$prop) ~ d$conc, dist = "weibull")
             Value Std. Error     z        p
(Intercept) -2.254     0.9506 -2.37 0.017737
d$conc       0.135     0.0686  1.97 0.048532
Log(scale)  -1.061     0.3203 -3.31 0.000927

Scale= 0.346

Weibull distribution
Loglik(model)= 0.6   Loglik(intercept only)= -1.6
	Chisq= 4.56 on 1 degrees of freedom, p= 0.033
Number of Newton-Raphson Iterations: 5
n= 6

Estimating the LC50 from these coefficients yields an unreasonable
answer:

LC50 = (0.5 + 2.254)/0.135 = 20.4; i.e., higher than the highest
exposure concentration.

I'm sure I'm doing something silly--I just don't know what.

Essentially, I'm trying to do this, but with a Weibull model:

> library(MASS)
> resp <- cbind(d$ndead, nalive = d$orign - d$ndead)
> mod <- glm(resp ~ d$conc, family = binomial(link = "probit"))
> result <- dose.p(mod, p = 0.5)
> result
             Dose        SE
p = 0.5: 11.49053 0.1069564

Any help would be greatly appreciated.

Sincerely,

Greg.




More information about the R-help mailing list