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
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