[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