Bisakah saya mengartikan dimasukkannya istilah kuadrat dalam regresi logistik sebagai indikasi titik balik?

Dalam Regresi Logistik dengan istilah linier dan kuadrat saja, jika saya memiliki koefisien linier dan koefisien kuadratik , dapatkah saya mengatakan bahwa ada titik balik probabilitas di ?β 2 - β 1 / ( $\beta_1$$\beta_2$$-\beta_1 / (2\beta_2)$

Ya kamu bisa.

Modelnya adalah

Ketika $\beta_2$ bukan nol, ini memiliki ekstrem global di $x = -\beta_1 / (2\beta_2)$ .

Regresi logistik memperkirakan koefisien ini sebagai . Karena ini adalah perkiraan kemungkinan maksimum (dan ML perkiraan fungsi parameter adalah fungsi yang sama dengan perkiraan), kami dapat memperkirakan lokasi ekstrem berada di .- b 1 / ( 2 b 2$(b_0, b_1, b_2)$$-b_1/(2b_2)$

Interval kepercayaan untuk perkiraan itu akan menarik. Untuk kumpulan data yang cukup besar untuk diterapkan teori kemungkinan maksimum asimptotik, kita dapat menemukan titik akhir interval ini dengan menyatakan kembali dalam formulir$\beta_0 + \beta_1 x + \beta_2 x^2$

dan menemukan berapa banyak dapat divariasikan sebelum kemungkinan log berkurang terlalu banyak. "Terlalu banyak" adalah, asimtotik, satu-setengah kuantil dari distribusi chi-kuadrat dengan satu derajat kebebasan.1 - α /$\gamma$$1-\alpha/2$

Pendekatan ini akan bekerja dengan baik asalkan rentang mencakup kedua sisi puncak dan ada cukup dan tanggapan di antara nilai-nilai untuk menggambarkan puncak itu. Jika tidak, lokasi puncak akan sangat tidak pasti dan perkiraan asimptotik mungkin tidak dapat diandalkan.0 1 $x$$0$$1$$y$

Rkode untuk melakukan ini di bawah. Ini dapat digunakan dalam simulasi untuk memeriksa bahwa cakupan interval kepercayaan dekat dengan cakupan yang dimaksud. Perhatikan bagaimana puncak sebenarnya adalah dan - dengan melihat baris bawah histogram - bagaimana sebagian besar batas kepercayaan yang lebih rendah kurang dari nilai sebenarnya dan sebagian besar batas kepercayaan atas lebih besar dari nilai sebenarnya, seperti yang kita harapkan. Dalam contoh ini cakupan yang dimaksudkan adalah dan cakupan aktual (diskon empat dari kasus di mana regresi logistik tidak bertemu) adalah , menunjukkan metode ini bekerja dengan baik (untuk jenis data yang disimulasikan sini).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)
}