Tengo una película cuasibinomial con dos variables explicativas continuas (digamos "LogPesticide" y "LogFood") y una interacción. Me gustaría calcular la LC50 del pesticida con intervalos de confianza en diferentes cantidades de alimentos (por ejemplo, el valor mínimo y máximo de los alimentos). ¿Cómo se puede lograr esto?

Ejemplo: Primero genero un conjunto de datos.

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)

Luego me ajusto a la película completa. El diagnóstico del modelo está bien y todos los términos de interacción son significativos.

model <- glm(Mortality ~ LogPesticide * LogFood, family = quasibinomial, data = mydata)
plot(model)
Anova(model)
summary(model)

Traté de estimar LC50s con dose.p () del paquete MASS. Si LogFood fuera un factor, esto funcionaría cuando vuelva a ajustar el modelo como se describe enesta publicación. Pero con dos variables explicativas continuas, obtienes solo 1 intersección, 2 pendientes y una diferencia de pendientes (para la interacción).

Puedo estimar las LC50 utilizando el efecto (), pero no sé cómo obtener el IC asociado para 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")

Del código de dose.p () veo que uno debe usar la matriz vcov. effect () también proporciona una matriz vcov pero no pude modificar dose.p () para trabajar con esa información correctamente. Le agradecería cualquier idea!