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?
Example: First I generate a data set.
mydata <- 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 <- function(x, a = 1, K = 1, r = 1) { # Logistic growth function. a = position of turning point
Fx <- (K * exp(r * (x - a))) / (1 + exp(r * (x - a))) # K = carrying capacity
return(Fx) # r = growth rate (larger r -> narrower curve)
}
y <- rep(NA, length = length(mydata$LogPesticide))
y[mydata$LogFood == log(1)] <- growth(x = mydata$LogPesticide[mydata$LogFood == log(1)], a = log(0.1), K = 1, r = 6)
y[mydata$LogFood == log(2)] <- growth(x = mydata$LogPesticide[mydata$LogFood == log(2)], a = log(0.2), K = 1, r = 4)
y[mydata$LogFood == log(4)] <- growth(x = mydata$LogPesticide[mydata$LogFood == log(4)], a = log(0.4), K = 1, r = 3)
y[mydata$LogFood == log(8)] <- growth(x = mydata$LogPesticide[mydata$LogFood == log(8)], a = log(0.8), K = 1, r = 1)
mydata$Dead <- rbinom(n = length(y), size = 20, prob = y)
mydata$Alive <- 20 - mydata$Dead
mydata$Mortality <- cbind(mydata$Dead, mydata$Alive)
Then I fit the full glm. Model diagnostics are ok and all interaction terms are significant.
model <- glm(Mortality ~ LogPesticide * LogFood, family = quasibinomial, data = mydata)
plot(model)
Anova(model)
summary(model)
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 this post. But with two continuous explanatory variables, you get only 1 intercept, 2 slopes and a difference of slopes (for the interaction).
I can estimate the LC50s using effect(), but don’t know how to get the associated CI for LogPesticide.
# Create a set of fitted values.
library(effects)
food.min <- round(min(model$model$LogFood), 3)
food.max <- round(max(model$model$LogFood), 3)
eff <- 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 <- as.data.frame(eff)
# Find fitted values closest to 0.5 when LogFood is minimal and maximal.
fit.min <- which.min(abs(eff2$fit[eff2$LogFood == food.min] - 0.5))
fit.min <- eff2$fit[eff2$LogFood == food.min][fit.min]
fit.max <- which.min(abs(eff2$fit[eff2$LogFood == food.max] - 0.5))
fit.max <- eff2$fit[eff2$LogFood == food.max][fit.max]
# Use those fitted values to predict the LC50s.
lc50.min <- eff2$LogPesticide[eff2$fit == fit.min & eff2$LogFood == food.min]
lc50.max <- eff2$LogPesticide[eff2$fit == fit.max & 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")
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!
