Konfidenzintervalle für die letale Dosis (LD) für die logistische Regression in R
Ich will Lethal Dose finden LD50
) mit seinem Konfidenzintervall inR
. Andere Softwarelinien von Minitab, SPSS und SAS bieten drei verschiedene Versionen solcher Konfidenzintervalle. Ich konnte solche Intervalle in keinem Paket in @ findR
(Ich habe auch @ verwendefindFn
Funktion vonsos
package).
Wie finde ich solche Intervalle? Ich habe für einen Intervalltyp basierend auf der Delta-Methode codiert (da ich nicht sicher bin, ob er korrekt ist), möchte aber eine beliebige etablierte Funktion von @ verwendeR
package. Vielen Dan
MWE:
dose <- c(10.2, 7.7, 5.1, 3.8, 2.6, 0)
total <- c(50, 49, 46, 48, 50, 49)
affected <- c(44, 42, 24, 16, 6, 0)
finney71 <- data.frame(dose, total, affected)
fm1 <- glm(cbind(affected, total-affected) ~ log(dose),
family=binomial(link = logit), data=finney71[finney71$dose != 0, ])
summary(fm1)$coef
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.886912 0.6429272 -7.601035 2.937717e-14
log(dose) 3.103545 0.3877178 8.004650 1.198070e-15
library(MASS)
xp <- dose.p(fm1, p=c(0.50, 0.90, 0.95)) # from MASS
xp.ci <- xp + attr(xp, "SE") %*% matrix(qnorm(1 - 0.05/2)*c(-1,1), nrow=1)
zp.est <- exp(cbind(xp, attr(xp, "SE"), xp.ci[,1], xp.ci[,2]))
dimnames(zp.est)[[2]] <- c("LD", "SE", "LCL","UCL")
zp.est
LD SE LCL UCL
p = 0.50: 4.828918 1.053044 4.363708 5.343724
p = 0.90: 9.802082 1.104050 8.073495 11.900771
p = 0.95: 12.470382 1.133880 9.748334 15.952512