¿Puedo interpretar que la inclusión de un término cuadrático en la regresión logística indica un punto de inflexión?

Question

En una Regresión Logística con términos lineales y cuadráticos solamente, si tengo un coeficiente lineal $\beta_1$ y el coeficiente cuadrático $\beta_2$ ¿Puedo decir que hay un punto de inflexión de la probabilidad en $-\beta_1 / (2\beta_2)$ ?

Answer 1

Sí, puedes hacerlo.

---

El modelo es

$$\Pr(Y=1) = \frac{1}{1+ \exp(-[\beta_0 + \beta_1 x + \beta_2 x^2])}.$$

Cuando $\beta_2$ es distinto de cero, tiene un extremo global en $x = -\beta_1 / (2\beta_2)$ .

Regresión logística estimaciones estos coeficientes como $(b_0, b_1, b_2)$ . Como se trata de una estimación de máxima verosimilitud (y las estimaciones ML de las funciones de los parámetros son las mismas funciones de las estimaciones) podemos estimar que la ubicación del extremo está en $-b_1/(2b_2)$ .

Sería interesante disponer de un intervalo de confianza para esa estimación. Para conjuntos de datos lo suficientemente grandes como para que se aplique la teoría asintótica de máxima probabilidad, podemos encontrar los puntos finales de este intervalo reexpresando $\beta_0 + \beta_1 x + \beta_2 x^2$ en la forma

$$\beta_0 + \beta_1 x + \beta_2 x^2 = \beta_0 - \beta_1^2/(4\beta_2) + \beta_2(x + \beta_1/(2\beta_2))^2 = \beta + \beta_2(x + \gamma)^2$$

y encontrar cuánto $\gamma$ puede variar antes de que la probabilidad logarítmica disminuya demasiado. "Demasiado" es, asintóticamente, la mitad del $1-\alpha/2$ cuantil de una distribución chi-cuadrado con un grado de libertad.

Este enfoque funcionará bien siempre que los rangos de $x$ cubren ambos lados del pico y hay un número suficiente de $0$ y $1$ respuestas entre los $y$ para delimitar ese pico. De lo contrario, la ubicación del pico será muy incierta y las estimaciones asintóticas pueden ser poco fiables.

---

R El código para llevarlo a cabo está abajo. Se puede utilizar en una simulación para comprobar que la cobertura de los intervalos de confianza se acerca a la cobertura prevista. Obsérvese cómo el verdadero pico es $-1/2$ y -al observar la fila inferior de histogramas- cómo la mayoría de los límites de confianza inferiores son menores que el valor real y la mayoría de los límites de confianza superiores son mayores que el valor real, como cabría esperar. En este ejemplo, la cobertura prevista era $1 - 2(0.05) = 0.90$ y la cobertura real (descontando cuatro de $500$ casos en los que la regresión logística no convergió) fue $0.91$ , lo que indica que el método funciona bien (para los tipos de datos simulados aquí).

n <- 50            # Number of observations in each trial
beta <- c(-1,2,2)  # Coefficients
x <- seq(from=-3, to=3, length.out=n)
y0 <- cbind(rep(1,length(x)), x, x^2) %*% beta

# Conduct a simulation.
set.seed(17)
sim <- replicate(500, peak(x, rbinom(length(x), 1, logistic(y0)), alpha=0.05))

# Post-process the results to check the actual coverage.
tp <- -beta[2] / (2 * beta[3])
covers <- sim["lcl",] <= tp & tp <= sim["ucl",]
mean(covers, na.rm=TRUE) # Should be close to 1 - 2*alpha

# Plot the distributions of the results.
par(mfrow=c(2,2))
plot(x, logistic(y0), type="l", lwd=2, col="#4040d0", main="Simulated Data",ylim=c(0,1))
points(x, rbinom(length(x), 1, logistic(y0)), pch=19)
hist(sim["peak.x",], main="Estimates"); abline(v=tp, col="Red")
hist(sim["lcl",], main="Lower Confidence Limits"); abline(v=tp, col="Red")
hist(sim["ucl",], main="Upper Confidence Limits"); abline(v=tp, col="Red")

logistic <- function(x) 1 / (1 + exp(-x))

peak <- function(x, y, alpha=0.05) {
  #
  # Estimate the peak of a quadratic logistic fit of y to x
  # and a 1-alpha confidence interval for that peak.
  #
  logL <- function(b) {
    # Log likelihood.
    p <- sapply(cbind(rep(1, length(x)), x, x*x) %*% b, logistic)
    sum(log(p[y==1])) + sum(log(1-p[y==0]))
  }
  f <- function(gamma) {
    # Deviance as a function of offset from the peak.
    b0 <- c(b[1] - b[2]^2/(4*b[3]) + b[3]*gamma^2, -2*b[3]*gamma, b[3])
    -2.0 * logL(b0)
  }
  # Estimation.
  fit <- glm(y ~ x + I(x*x), family=binomial(link = "logit"))
  if (!fit$converged) return(rep(NA,3))

  b <- coef(fit)
  tp <- -b[2] / (2 * b[3])

  # Two-sided confidence interval:
  # Search for where the deviance is at a threshold determined by alpha.
  delta <- qchisq(1-alpha, df=1)
  u <- sd(x)
  while(fit$deviance - f(tp+u) + delta > 0) u <- 2*u # Find an upper bound
  l <- sd(x)
  while(fit$deviance - f(tp-l) + delta > 0) l <- 2*l # Find a lower bound
  upper <- uniroot(function(gamma) fit$deviance - f(gamma) + delta, 
                   interval=c(tp, tp+u))
  lower <- uniroot(function(gamma) fit$deviance - f(gamma) + delta, 
                   interval=c(tp-l, tp))

  # Return a vector of the estimate, lower limit, and upper limit.
  c(peak=tp, lcl=lower$root, ucl=upper$root)
}