Ich habe ein Quasibinomial-Glm mit zwei kontinuierlichen erklärenden Variablen (sagen wir "LogPesticide" und "LogFood") und einer Interaktion. Ich möchte die LC50 des Pestizids mit Konfidenzintervallen für unterschiedliche Mengen an Lebensmitteln berechnen (z. B. den minimalen und maximalen Lebensmittelwert). Wie kann das erreicht werden?

Beispiel: Zuerst erstelle ich einen Datensatz.

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)

Dann passe ich die volle glm. Die Modelldiagnose ist in Ordnung und alle Interaktionsterme sind von Bedeutung.

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

Ich habe versucht, LC50s mit dose.p () aus dem MASS-Paket zu schätzen. Wenn LogFood ein Faktor wäre, würde dies funktionieren, wenn ich das Modell wie in @ beschrieben neu anpassdieser Beitra. Mit zwei kontinuierlichen erklärenden Variablen erhalten Sie jedoch nur einen Achsenabschnitt, zwei Steigungen und eine Steigungsdifferenz (für die Interaktion).

Ich kann die LC50 mit effect () abschätzen, weiß aber nicht, wie ich das zugehörige CI für LogPesticide abrufen kann.

# 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")

Aus dem Code von dose.p () sehe ich, dass man die VCOV-Matrix verwenden muss. effect () bietet auch eine VCOV-Matrix, aber ich konnte dose.p () nicht ändern, um mit diesen Informationen richtig zu arbeiten. Für Anregungen wäre ich dankbar!