Cómo estimar los intervalos de confianza para la CL50

Question

Esta es mi primera pregunta, así que espero que la pregunta esté bien hecha (mis disculpas si no lo está...)

Estoy utilizando un modelo GLM binomial (logit) para algunos datos de toxicología que investigan los efectos de un pesticida en algunos organismos. Se expuso la misma cantidad de organismos a diferentes concentraciones del pesticida y se evaluó la supervivencia al final del experimento. Los datos y el modelo son los siguientes:

## data:
Concentration <- as.numeric(c(0, 100, 200, 300, 400, 500, 600, 700, 800))
Survival <- as.integer(c(20, 20, 20, 18, 13, 8, 3, 0, 0))
Total <- as.integer(c(20, 20, 20, 20, 20, 20, 20, 20, 20))
Data <- data.frame(Concentration, Survival, Total)
Data$Dead <- Data$Total-Data$Survival
Data$Proport_Surv <- (Data$Total - Data$Dead) / Data$Total

mod <- glm(cbind(Survival, Dead) ~ Concentration, Data, family=binomial(link = "logit"))

Esta es la salida resumida:

    Call:
    glm(formula = cbind(Survival, Dead) ~ Concentration, family = binomial(link = "logit"), data = Data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0156  -0.4775   0.1917   0.4356   0.8721  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    6.992127   1.121066   6.237 4.46e-10 ***
Concentration -0.015196   0.002366  -6.423 1.34e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 163.5934  on 8  degrees of freedom
Residual deviance:   3.2464  on 7  degrees of freedom
AIC: 19.401

Number of Fisher Scoring iterations: 5

Un resultado importante en toxicología para inferir sobre la toxicidad de cualquier sustancia química, es estimar el valor LC50, que es la concentración estimada por el modelo que provoca el 50% de mortalidad. He utilizado la función dose.p() del MASS. Para este conjunto de datos:

library(MASS)
dose.p(mod, p=0.5)
             Dose       SE
p = 0.5: 460.1353 18.16643

La CL50 para este pesticida sería de 460,1353 pero este valor debería ir acompañado de intervalos de confianza. Mi duda es ¿Cómo estimar los intervalos de confianza para un valor X? Trazar el modelo ayuda a entender cuál debería ser mi resultado:

Cálculo de los intervalos de confianza del modelo:

NewData1 <- expand.grid(Conc  = seq(0, 800, length = 9))

P_logit <- predict(mod, newdata = NewData1, se = TRUE, type = "link")
NewData1$P_logit <- exp(P_logit$fit) / (1 + exp(P_logit$fit))
NewData1$SeUp <- exp(P_logit$fit + 1.96*P_logit$se.fit) / (1 + exp(P_logit$fit + 1.96*P_logit$se.fit))
NewData1$SeLo <- exp(P_logit$fit - 1.96*P_logit$se.fit) / (1 + exp(P_logit$fit - 1.96*P_logit$se.fit))

Trama:

library(ggplot2)

plot <- ggplot()
plot <- plot + geom_point(data = Data, aes(y = Proport_Surv, x = Concentration), shape = 1, size = 2.5)
plot <- plot + xlab("Concentration") + ylab("Proportion of surviving organisms")
plot <- plot + theme(text = element_text(size=15)) + theme_bw()
plot <- plot + geom_line(data = NewData1, aes(x = Concentration, y = P_logit), colour = "black")
plot <- plot + geom_ribbon(data = NewData1, aes(x = Concentration, ymax = SeUp, ymin = SeLo), alpha = 0.2)
plot <- plot + geom_hline(yintercept = 0.5, linetype = "dashed")
plot <- plot + annotate("point", x = 460.1353, y = 0.5, size = 3.25, colour = "blue")
plot

La anotación del punto azul es la CL50 estimada por el modelo. Cómo puedo estimar los valores donde la línea discontinua intercepta los intervalos de confianza inferior y superior del modelo?

Answer 1

La forma más funcional es utilizar el ajuste de la función dose.p que se puede encontrar en Venables, W. N. y Ripley, B. D. (2002) Modern Applied Statistics with S. Springer.

El código siguiente utiliza el estadístico de Wald para el IC del 95%

ec = 0.5
library(VGAM)
eta <- logit(ec) 
beta <- coef(mod)[1:2] 
ecx <- (eta - beta[1])/beta[2] 
pd <- -cbind(1, ecx)/beta[2] 
ff = as.matrix(vcov(mod)[1:2,1:2])
se <- sqrt(((pd %*% ff )* pd) %*% c(1, 1))
upper = (ecx+se*1.96)
lower = (ecx-se*1.96)
df1 = data.frame(ecx, lower, upper)

plot <- plot + geom_vline(xintercept = df1$ecx, linetype = "dashed")
plot <- plot + geom_vline(xintercept = df1$lower, linetype = "dashed")
plot <- plot + geom_vline(xintercept = df1$upper, linetype = "dashed")
plot

Answer 2

Volvamos a los primeros principios para ver lo que está intentando estimar. El MLG con la familia binomial y la función de enlace logit es simplemente el modelo de regresión logística por lo que podemos utilizar la ecuación de regresión para ese modelo. Supongamos que $p \equiv \mathbb{P}(Y=1)$ sea la probabilidad de un resultado de respuesta positivo (aquí la muerte del organismo), su modelo se basa en la ecuación del modelo:

$$\log \Big( \frac{p}{1-p} \Big) = \beta_0 + \beta_1 x.$$

Configuración $p = \tfrac{1}{2}$ da el valor explicativo correspondiente (en este caso, la concentración LC50):

$$x_* \equiv \frac{\log(\tfrac{1}{2}/\tfrac{1}{2}) - \beta_0}{\beta_1} = \frac{\log(1) - \beta_0}{\beta_1} = - \frac{\beta_0}{\beta_1}.$$

Por lo tanto, en realidad está intentando encontrar un intervalo de confianza para la relación de los coeficientes verdaderos en la regresión logística. Normalmente aproximamos la distribución conjunta de los coeficientes estimados mediante una distribución normal (apelando al teorema del límite central para grandes $n$ ) y, por lo tanto, estamos intentando encontrar un intervalo de confianza para la proporción verdadera basado en estimadores observados que se supone que se distribuyen normalmente de forma conjunta alrededor de los valores verdaderos (con correlación distinta de cero). (La matriz de varianza-covarianza estimada para los estimadores se puede encontrar utilizando la función vcov en su modelo).

La otra respuesta aquí ya muestra cómo se puede utilizar bootstrapping para obtener un intervalo de confianza en este caso, pero hay una serie de aproximaciones analíticas que puede utilizar si lo prefiere; para algunos ejemplos de formas analíticas para este tipo de intervalo de confianza, véase, por ejemplo, Malley (1982) , Shanmugalingam (1982) y Dunlap y Silver (1986) .

Answer 3

Si tiene suficiente poder de cálculo, un intervalo de confianza bootstrap podría ser una opción.