Puis-je interpréter l'inclusion d'un terme quadratique dans la régression logistique comme indiquant un tournant?

Dans une régression logistique avec des termes linéaires et quadratiques uniquement, si j'ai un coefficient linéaire et un coefficient quadratique , puis-je dire qu'il y a un point tournant de la probabilité à ?β 2 - β 1 / ( 2 β 2$\beta_1$$\beta_2$$-\beta_1 / (2\beta_2)$

Oui, vous pouvez.

Le modèle est

Lorsque $\beta_2$ est différent de zéro, cela a un extremum global à $x = -\beta_1 / (2\beta_2)$ .

La régression logistique estime ces coefficients comme . Comme il s'agit d'une estimation du maximum de vraisemblance (et les estimations ML des fonctions des paramètres sont les mêmes fonctions des estimations), nous pouvons estimer que l'emplacement de l'extremum est à .- b 1 / ( 2 b 2$(b_0, b_1, b_2)$$-b_1/(2b_2)$

Un intervalle de confiance pour cette estimation serait intéressant. Pour les ensembles de données suffisamment volumineux pour que la théorie du maximum de vraisemblance asymptotique s'applique, nous pouvons trouver les points finaux de cet intervalle en ré-exprimant sous la forme$\beta_0 + \beta_1 x + \beta_2 x^2$

et trouver combien peut varier avant que la probabilité logarithmique ne diminue trop. «Trop» est, asymptotiquement, la moitié du quantile d'une distribution chi carré avec un degré de liberté.1 - α /$\gamma$$1-\alpha/2$

Cette approche fonctionnera bien à condition que les plages de couvrent les deux côtés du pic et qu'il y ait un nombre suffisant de réponses et parmi les valeurs pour délimiter ce pic. Sinon, l'emplacement du pic sera très incertain et les estimations asymptotiques pourraient ne pas être fiables.0 1 $x$$0$$1$$y$

Rle code pour effectuer ceci est ci-dessous. Il peut être utilisé dans une simulation pour vérifier que la couverture des intervalles de confiance est proche de la couverture prévue. Remarquez comment le vrai pic est et - en regardant la ligne inférieure des histogrammes - comment la plupart des limites de confiance inférieures sont inférieures à la valeur réelle et la plupart des limites de confiance supérieures sont supérieures à la valeur réelle, comme nous l'espérons. Dans cet exemple, la couverture prévue était de et la couverture réelle (en excluant quatre des cas où la régression logistique n'a pas convergé) était de , ce qui indique que la méthode fonctionne bien (pour les types de données simulées ici).1 - 2 ( 0,05 ) = 0,90 500 $-1/2$$1 - 2(0.05) = 0.90$$500$$0.91$

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