LC50 / LD50 confidence intervals from multiple regression glm with interaction
<p>I have a quasibinomial glm with two continuous explanatory variables (let's say "LogPesticide" and "LogFood") and an interaction. I would like to calculate the LC50 of the pesticide with confidence intervals at different amounts of food (e. g. the minimum and maximum food value). How can this be achieved?</p>

<p>Example: First I generate a data set.</p>

<pre><code>mydata &lt;- data.frame(
            LogPesticide = rep(log(c(0, 0.1, 0.2, 0.4, 0.8, 1.6) + 0.05), 4),
            LogFood = rep(log(c(1, 2, 4, 8)), each = 6)
          )

set.seed(seed=16) 

growth &lt;- function(x, a = 1, K = 1, r = 1) {            # Logistic growth function. a = position of turning point
  Fx &lt;- (K * exp(r * (x - a))) / (1 + exp(r * (x - a))) # K = carrying capacity
  return(Fx)                                            # r = growth rate (larger r -&gt; narrower curve)
}

y &lt;- rep(NA, length = length(mydata$LogPesticide))
y[mydata$LogFood == log(1)] &lt;- growth(x = mydata$LogPesticide[mydata$LogFood == log(1)], a = log(0.1), K = 1, r = 6)
y[mydata$LogFood == log(2)] &lt;- growth(x = mydata$LogPesticide[mydata$LogFood == log(2)], a = log(0.2), K = 1, r = 4)
y[mydata$LogFood == log(4)] &lt;- growth(x = mydata$LogPesticide[mydata$LogFood == log(4)], a = log(0.4), K = 1, r = 3)
y[mydata$LogFood == log(8)] &lt;- growth(x = mydata$LogPesticide[mydata$LogFood == log(8)], a = log(0.8), K = 1, r = 1)
mydata$Dead &lt;- rbinom(n = length(y), size = 20, prob = y)
mydata$Alive &lt;- 20 - mydata$Dead
mydata$Mortality &lt;- cbind(mydata$Dead, mydata$Alive)
</code></pre>

<p>Then I fit the full glm. Model diagnostics are ok and all interaction terms are significant.</p>

<pre><code>model &lt;- glm(Mortality ~ LogPesticide * LogFood, family = quasibinomial, data = mydata)
plot(model)
Anova(model)
summary(model)
</code></pre>

<p>I tried to estimate LC50s with dose.p() from the MASS package. If LogFood was a factor, this would work when I re-fit the model as discussed in <a href="https://stat.ethz.ch/pipermail/r-help/2008-February/154177.html" rel="nofollow">this post</a>. But with two continuous explanatory variables, you get only 1 intercept, 2 slopes and a difference of slopes (for the interaction).</p>

<p>I can estimate the LC50s using effect(), but don’t know how to get the associated CI for LogPesticide. </p>

<pre><code># Create a set of fitted values.

library(effects)

food.min &lt;- round(min(model$model$LogFood), 3)
food.max &lt;- round(max(model$model$LogFood), 3)

eff &lt;- effect("LogPesticide:LogFood", model, 
              xlevels = list(LogPesticide = seq(min(model$model$LogPesticide), max(model$model$LogPesticide), length = 100),
                             LogFood = c(food.min, food.max)))

eff2 &lt;- as.data.frame(eff)


# Find fitted values closest to 0.5 when LogFood is minimal and maximal.

fit.min &lt;- which.min(abs(eff2$fit[eff2$LogFood == food.min] - 0.5))
fit.min &lt;- eff2$fit[eff2$LogFood == food.min][fit.min]

fit.max &lt;- which.min(abs(eff2$fit[eff2$LogFood == food.max] - 0.5))
fit.max &lt;- eff2$fit[eff2$LogFood == food.max][fit.max]


# Use those fitted values to predict the LC50s.

lc50.min &lt;- eff2$LogPesticide[eff2$fit == fit.min &amp; eff2$LogFood == food.min]
lc50.max &lt;- eff2$LogPesticide[eff2$fit == fit.max &amp; eff2$LogFood == food.max]


# Plot the results.

plot(fit ~ LogPesticide, data = eff2[eff2$LogFood == round(min(model$model$LogFood), 3),], type = "l")
lines(fit ~ LogPesticide, data = eff2[eff2$LogFood == round(max(model$model$LogFood), 3),], col = "red")
points(y = 0.5, x = lc50.min, pch = 19)
points(y = 0.5, x = lc50.max, pch = 19, col = "red")
</code></pre>

<p>From the code of dose.p() I see one must use the vcov matrix. effect() also provides a vcov matrix but I could not modify dose.p() to work with that information correctly. I would be grateful for any ideas!</p>