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

12

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 )β1β2β1/(2β2)

FZo
la source

Réponses:

17

Oui, vous pouvez.


Le modèle est

Pr(Y=1)=11+exp([β0+β1x+β2x2]).

Lorsque β2 est différent de zéro, cela a un extremum global à x=β1/(2β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 )(b0,b1,b2)b1/(2b2)

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β0+β1x+β2x2

β0+β1x+β2x2=β0β12/(4β2)+β2(x+β1/(2β2))2=β+β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 - α / 2γ1α/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 yx01y


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 0,911/212(0.05)=0.905000.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")

Résultats

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)
}
whuber
la source
+1, excellente réponse. Vous mentionnez quelques mises en garde concernant l'approche asymptotique; que pensez-vous de l'amorçage du CI dans des cas comme celui-ci? J'ai fait cela une fois pour montrer que le pic d'une courbe quadratique pour un groupe était supérieur à celui de l'autre groupe.
gung - Rétablir Monica
Cela pourrait fonctionner, @gung, mais la théorie du bootstrap s'applique également aux grands échantillons. Dans votre application, un test de permutation pourrait peut-être être justifié.
whuber
Cool. Mais le tournant ne serait-il pas en dehors de la plage de données? Et puis il serait dangereux d'extrapoler.
Peter Flom - Réintègre Monica
@Peter C'est vrai, c'est pourquoi j'ai commenté que "Cette approche fonctionnera bien à condition que les plages de x couvrent les deux côtés du pic."
whuber
@whuber Oups, j'ai raté ça. Désolé
Peter Flom - Réintègre Monica