Kann ich die Einbeziehung eines quadratischen Terms in die logistische Regression als Hinweis auf einen Wendepunkt interpretieren?

Kann ich in einer logistischen Regression mit nur linearen und quadratischen Termen, wenn ich einen linearen Koeffizienten $\beta_1$ und einen quadratischen Koeffizienten $\beta_2$ , sagen, dass es einen Wendepunkt der Wahrscheinlichkeit bei $-\beta_1 / (2\beta_2)$ ?

Ja, du kannst.

Das Modell ist

Wenn ungleich Null ist, hat dies ein globales Extremum bei x = - β 1 / ( 2 β 2 ) .$\beta_2$$x = -\beta_1 / (2\beta_2)$

Die logistische Regression schätzt diese Koeffizienten als . Da dies eine Maximum-Likelihood-Schätzung ist (und ML-Schätzungen von Funktionen von Parametern die gleichen Funktionen der Schätzungen sind), können wir die Position des Extremums bei - b 1 / ( 2 b 2 ) schätzen .$(b_0, b_1, b_2)$$-b_1/(2b_2)$

Ein Konfidenzintervall für diese Schätzung wäre von Interesse. Für Datensätze, die groß genug sind, um die asymptotische Maximum-Likelihood-Theorie anzuwenden, können wir die Endpunkte dieses Intervalls durch erneutes Ausdrücken von in der Form$\beta_0 + \beta_1 x + \beta_2 x^2$

und Finden, wie viel & variiert werden kann, bevor die logarithmische Wahrscheinlichkeit zu stark abnimmt. "Zu viel" ist asymptotisch die Hälfte des 1 - α / 2- Quantils einer Chi-Quadrat-Verteilung mit einem Freiheitsgrad.$\gamma$$1-\alpha/2$

Dieser Ansatz funktioniert gut, vorausgesetzt, die Bereiche von decken beide Seiten des Peaks ab und es gibt eine ausreichende Anzahl von 0 und 1 Antworten unter den y- Werten, um diesen Peak abzugrenzen. Andernfalls ist der Ort des Peaks sehr unsicher und die asymptotischen Schätzungen können unzuverlässig sein.$x$$0$$1$$y$

RCode, um dies auszuführen, ist unten. Es kann in einer Simulation verwendet werden, um zu überprüfen, ob die Abdeckung der Konfidenzintervalle nahe an der beabsichtigten Abdeckung liegt. Beachten Sie, wie die wahre Spitze ist und - durch in der unteren Reihe von Histogrammen suchen - wie die meisten der unteren Vertrauensgrenzen weniger sind als der wahre Wert und die meisten der oberen Vertrauensgrenzen sind größer als der wahre Wert, wie wir hoffen würden. In diesem Beispiel betrug die beabsichtigte Abdeckung 1 - 2 ( 0,05 ) = 0,90 und die tatsächliche Abdeckung (abzüglich vier von 500 Fällen, in denen die logistische Regression nicht konvergierte) $-1/2$$1 - 2(0.05) = 0.90$$500$$0.91$Dies zeigt an, dass die Methode gut funktioniert (für die hier simulierten Datentypen).

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)
}