3
$\begingroup$

This is my first question, so I hope the question is properly done (my apologies if it's not...)

I am using a binomial GLM model (logit) for some toxicology data investigating the effects of a pesticide to some organisms. The same amount of organisms were exposed to different concentrations of the pesticide and survival was assessed at the end of the experiment. The data and model go as follows:

## data:
Concentration <- as.numeric(c(0, 100, 200, 300, 400, 500, 600, 700, 800))
Survival <- as.integer(c(20, 20, 20, 18, 13, 8, 3, 0, 0))
Total <- as.integer(c(20, 20, 20, 20, 20, 20, 20, 20, 20))
Data <- data.frame(Concentration, Survival, Total)
Data$Dead <- Data$Total-Data$Survival
Data$Proport_Surv <- (Data$Total - Data$Dead) / Data$Total

mod <- glm(cbind(Survival, Dead) ~ Concentration, Data, family=binomial(link = "logit"))

This is the summary output:

    Call:
    glm(formula = cbind(Survival, Dead) ~ Concentration, family = binomial(link = "logit"), data = Data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0156  -0.4775   0.1917   0.4356   0.8721  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    6.992127   1.121066   6.237 4.46e-10 ***
Concentration -0.015196   0.002366  -6.423 1.34e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 163.5934  on 8  degrees of freedom
Residual deviance:   3.2464  on 7  degrees of freedom
AIC: 19.401

Number of Fisher Scoring iterations: 5

An important output in toxicology to infer about the toxicity of any chemical, is to estimate the LC50 value, which is the concentration estimated by the model that elicits 50% mortality. I used the dose.p() function from MASS. For this dataset:

library(MASS)
dose.p(mod, p=0.5)
             Dose       SE
p = 0.5: 460.1353 18.16643

The LC50 for this pesticide would be 460.1353 but this value should be accompanied by confidence intervals. My doubt is: How to estimate confidence intervals for an X value? Plotting the model helps to understand what should be my output:

Calculating confidence intervals for the model:

NewData1 <- expand.grid(Conc  = seq(0, 800, length = 9))

P_logit <- predict(mod, newdata = NewData1, se = TRUE, type = "link")
NewData1$P_logit <- exp(P_logit$fit) / (1 + exp(P_logit$fit))
NewData1$SeUp <- exp(P_logit$fit + 1.96*P_logit$se.fit) / (1 + exp(P_logit$fit + 1.96*P_logit$se.fit))
NewData1$SeLo <- exp(P_logit$fit - 1.96*P_logit$se.fit) / (1 + exp(P_logit$fit - 1.96*P_logit$se.fit))

Plot:

library(ggplot2)

plot <- ggplot()
plot <- plot + geom_point(data = Data, aes(y = Proport_Surv, x = Concentration), shape = 1, size = 2.5)
plot <- plot + xlab("Concentration") + ylab("Proportion of surviving organisms")
plot <- plot + theme(text = element_text(size=15)) + theme_bw()
plot <- plot + geom_line(data = NewData1, aes(x = Concentration, y = P_logit), colour = "black")
plot <- plot + geom_ribbon(data = NewData1, aes(x = Concentration, ymax = SeUp, ymin = SeLo), alpha = 0.2)
plot <- plot + geom_hline(yintercept = 0.5, linetype = "dashed")
plot <- plot + annotate("point", x = 460.1353, y = 0.5, size = 3.25, colour = "blue")
plot

The blue point annotation is the LC50 estimated by the model. How can I estimate the values where the dashed line intercepts the lower and upper confidence intervals of the model?

$\endgroup$
0
$\begingroup$

If you have enough computational power, a bootstrap confidence interval might be an option.

$\endgroup$

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

  • $\begingroup$ Does the model give better results in estimating the LC50 with both the "zero concentration" and "zero survival" data points removed? I ask as it seems that those two data points are not informative - the "zero survival" concentrations can be very large, and the "zero concentration" is not making a mortality test at all. I'm thinking these two data points might be distorting the regression analysis. $\endgroup$ – James Phillips Mar 21 at 22:05
  • $\begingroup$ Thanks for the answers. @DavidMoseler, to be honest I am not familiar with bootstrap confidence intervals, so I am still trying to understand how would that allow me to calculate the confidence intervals for the LC50 (concentration that elicits 50% effect, which is in the X axis). Could you please show an example? $\endgroup$ – Mr Krinkle Mar 27 at 11:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.