Eu tenho um glm quasibinomial com duas variáveis explicativas contínuas (digamos "LogPesticide" e "LogFood") e uma interação. Eu gostaria de calcular a CL50 do pesticida com intervalos de confiança em diferentes quantidades de alimento (por exemplo, o valor mínimo e máximo de alimento). Como isso pode ser alcançado?

Exemplo: Primeiro eu giro um conjunto de dados.

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)

Então eu me encaixo no glm completo. O diagnóstico do modelo está correto e todos os termos de interação são significativos.

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

Tentei estimar LC50s com dose.p () do pacote MASS. Se o LogFood fosse um fator, isso funcionaria quando eu reajustasse o modelo, conforme discutido emesta postagem. Mas com duas variáveis explicativas contínuas, você obtém apenas 1 interceptação, 2 inclinações e uma diferença de inclinações (para a interação).

Posso estimar os LC50s usando effect (), mas não sei como obter o IC associado ao 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")

A partir do código de dose.p (), vejo que é preciso usar a matriz vcov. effect () também fornece uma matriz vcov, mas não pude modificar dose.p () para trabalhar com essas informações corretamente. Ficaria muito grato por todas as idéias!