я сказать, что в случае логистической регрессии только с линейными и квадратичными терминами, если у меня есть линейный коэффициент и квадратичный коэффициент , существует переломный момент вероятности в ?β 2 - β 1 / ( 2 β 2 )
я сказать, что в случае логистической регрессии только с линейными и квадратичными терминами, если у меня есть линейный коэффициент и квадратичный коэффициент , существует переломный момент вероятности в ?β 2 - β 1 / ( 2 β 2 )
Да, ты можешь.
Модель
Когда от нуля, это имеет глобальный экстремум в .
Логистическая регрессия оценивает эти коэффициенты как . Поскольку это оценка максимального правдоподобия (и ML-оценки функций параметров являются одинаковыми функциями оценок), мы можем оценить местоположение экстремума в .- b 1 / ( 2 b 2 )
Доверительный интервал для этой оценки будет представлять интерес. Для наборов данных, достаточно больших для применения асимптотической теории максимального правдоподобия, мы можем найти конечные точки этого интервала, повторно выразив в виде
и найти сколько 1 - α / 2 может быть изменено до того, как вероятность регистрации слишком сильно уменьшится. «Слишком много» - это асимптотически половина 2/2 распределения хи-квадрат с одной степенью свободы.
Этот подход будет работать хорошо при условии, что диапазоны охватывают обе стороны пика и имеется достаточное количество ответов и среди значений чтобы очертить этот пик. В противном случае местоположение пика будет очень неопределенным, и асимптотические оценки могут быть ненадежными.
R
код для выполнения этого ниже. Он может использоваться в моделировании для проверки того, что покрытие доверительных интервалов близко к предполагаемому покрытию. Обратите внимание, что истинный пик равен и - глядя на нижний ряд гистограмм - как большинство нижних доверительных границ меньше истинного значения, а большинство верхних доверительных границ больше истинного значения, как мы надеемся. В этом примере предполагаемое покрытие составляло а фактическое покрытие (без учета четырех из случаев, когда логистическая регрессия не сходилась) составляло , что указывает на то, что метод работает хорошо (для типов данных, моделируемых Вот).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)
}