R help - Mar 2009 - Estimating LC50 from a Weibull distribution

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.

Hi Greg,

you can use the extension package 'drc' from CRAN:


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)

## Loading 'drc'
library(drc)

## Fitting model assuming 100% mortality for high concentrations
## and 0% for concentration 0
d.m1<-drm(ndead/orign~conc, weight=orign, data=d, fct=W1.2(),
type="binomial")
plot(d.m1)

## Fitting model where mortality at conc=0 is estimated
## (ad hoc adjustments should be avoided)
d.m2<-update(d.m1, fct=W1.3u())
plot(d.m2, add=TRUE, type="none", lty=2)

## Calculating LC50
ED(d.m2, 50)


Note that you don't really have time-to-event data as you only have
observations
reflecting the effect of exposure at one particular time point (96h).
Time-to-event data
occur if exposure is monitored repeatedly over time, in several intervals or at
several
time points, e.g. at 24h, 48h, and 96h, possibly subject to censoring. Therefore
the above
analysis is based on binomial distributions that are often appropriate for
quantal data.

Note also, that in ecotoxicology "Weibull model" may refer to one
several related models
depending on which guidelines, papers, or software you refer to.


Christian