如何解释逻辑回归的系数
How to interpret coefficients of logistic regression
我试图弄清楚具有多项式项的逻辑回归系数与预测的关系。具体来说，我对 x 轴上预测最高的位置感兴趣。下面的示例：
set.seed(42)
# Setup some dummy data
x <- 1:200
y <- rep(0, length(x))
y[51:150] <- rbinom(100, 1, 0.5)
# Fit a model
family <- binomial()
model <- glm(y ~ poly(x, 2), family = family)
# Illustrate model
plot(x, y)
lines(x, family$linkinv(predict(model)), col = 2)
上面的模型给出了这些系数：
coef(model)
#> (Intercept) poly(x, 2)1 poly(x, 2)2
#> -1.990317 -3.867855 -33.299893
由 reprex package (v1.0.0)
于 2021-08-03 创建
poly() 的手册页说明如下：
The orthogonal polynomial is summarized by the coefficients, which can be used to evaluate it via the three-term recursion given in Kennedy & Gentle (1980, pp. 343–4), and used in the predict part of the code.
但是，我无法访问这本书，也无法从 predict.glm S3 方法中辨别这些系数是如何处理的。有没有一种方法可以仅根据系数（即不使用 predict() 来找到最大值）来重建山顶的位置（示例中大约为 100）？
假设您想通过分析找到预测的最大值，对于正交多项式为 1 阶和 2 阶的特定情况，我建议采用以下方法：
摘要
1) 推断多项式系数
这可以通过将线性模型拟合到模型矩阵中包含的相应多项式值来轻松完成。
2) 推导预测表达式w.r.t。 x 并将导数设置为 0
求解由(1)中的多项式拟合推导出的预测表达式中的x，得到预测最大值出现的x的值。
详情
1) 多项式系数
根据您拟合 GLM 模型的行，我们估计 1 阶多项式的系数 p1(x) = a0 + a1*x 和 2 阶多项式的系数 p2(x) = b0 + b1*x + b2*x^2：
X = model.matrix(model)
p1 = X[, "poly(x, 2)1"]
p2 = X[, "poly(x, 2)2"]
p1.lm = lm(p1 ~ x)
a0 = p1.lm$coeff["(Intercept)"]
a1 = p1.lm$coeff["x"]
p2.lm = lm(p2 ~ x + I(x^2))
b0 = p2.lm$coeff["(Intercept)"]
b1 = p2.lm$coeff["x"]
b2 = p2.lm$coeff["I(x^2)"]
这给出：
> c(a0, a1, b0, b1, b2)
(Intercept) x (Intercept) x I(x^2)
-1.2308840e-01 1.2247602e-03 1.6050353e-01 -4.7674315e-03 2.3718565e-05
2) 求最大值的预测导数
预测表达式 z（在反 link 函数之前）是：
z = Intercept + coef1 * p1(x) + coef2 * p2(x)
我们推导这个表达式，将其置0得到：
coef1 * a1 + coef2 * (b1 + 2 * b2 * xmax) = 0
求解 xmax 我们得到：
xmax = - (coef1 * a1 + coef2 * b1) / (2 * coef2 * b2)
在 R 代码中，计算如下：
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
(xmax = - ( coef1 * a1 + coef2 * b1 ) / (2 * coef2 * b2))
给出：
x
97.501114
检查
我们可以通过将最大值作为绿色十字添加到预测曲线来验证最大值：
# Prediction curve computed analytically
Intercept = model$coeff["(Intercept)"]
pred.analytical = family$linkinv( Intercept + coef1 * p1 + coef2 * p2 )
# Find the prediction's maximum analytically
pred.max = family$linkinv( Intercept + coef1 * (a0 + a1 * xmax) +
coef2 * (b0 + b1 * xmax + b2 * xmax^2) )
# Plot
plot(x, y)
# The following two lines should coincide!
lines(x, pred.analytical, col = 3)
lines(x, family$linkinv(predict(model)), col = 2)
# Location of the maximum!
points(xmax, pred.max, pch="x", col="green")
给出：
从正交多项式的理论表达式推导预测最大值的位置
我得到了 Kennedy 和 Gentle（1982 年） 的《统计计算》一书的副本，在 poly 的文档中引用了它，现在分享我的发现正交多项式的计算，以及我们如何使用它们找到最大预测值的位置 x。
书中介绍的正交多项式 (pp. 343-4) 是一元的（即最高阶系数始终为 1）并通过以下递归过程获得：
其中 q 是考虑的正交多项式的数量。
请注意以下上述术语与poly文档的关系：
- 出现在你的问题摘录中的“三项递归”是第三个表达式的RHS，它恰好有三个项。
- 第三个表达式中的rho(j+1)系数称为“中心常数”.
- 第三个表达式中的 gamma(j) 系数在文档中没有名称，但它们与 “归一化常数”直接相关。 =126=]，如下所示。
作为参考，我在这里粘贴了 poly 文档“值”部分的相关摘录：
A matrix with rows corresponding to points in x and columns corresponding to the degree, with attributes "degree" specifying the degrees of the columns and (unless raw = TRUE) "coefs" which contains the centering and normalization constants used in constructing the orthogonal polynomials
回到递归，我们可以推导出参数rho(j+1)和的值gamma(j) 通过在 p(j+1) w.r.t 上施加正交条件来从第三个表达式中得出。 p(j) 和 p(j-1).
（重要的是要注意，正交性条件 不是 积分，而是 n 观察到的 x 点，因此多项式系数取决于数据！——例如 Tchebyshev 正交多项式并非如此）。
参数表达式变为：
对于回归中使用的 1 阶和 2 阶多项式，我们得到以下表达式，已用 R 代码编写：
# First we define the number of observations in the data
n = length(x)
# For p1(x):
# p1(x) = (x - rho1) p0(x) (since p_{-1}(x) = 0)
rho1 = mean(x)
# For p2(x)
# p2(x) = (x - rho2) p1(x) - gamma1
gamma1 = var(x) * (n-1)/n
rho2 = sum( x * (x - mean(x))^2 ) / (n*gamma1)
为此我们得到：
> c(rho1, rho2, gamma1)
[1] 100.50 100.50 3333.25
注意 poly(x,2) 的 coefs 属性是：
> attr(poly(x,2), "coefs")
$alpha
[1] 100.5 100.5
$norm2
[1] 1 200 666650 1777555560
其中 $alpha 包含 居中常数 ，即 rho 值（与我们的值一致——顺便说一句，所有居中当 x 的 分布对于任何 q 是对称的 时，常数等于 x 的平均值！（观察和证明））和 $norm2 包含 归一化常数 （在本例中为 p(-1,x), p(0,x) 、p(1,x)、p(2,x))，即常数c(j) 将递归公式中的多项式归一化——通过将它们除以 sqrt(c(j))——，使生成的多项式 r(j,x) 满足 sum_over_i{ r(j,x_i)^2 } = 1; 注意 r(j,x) 是存储在 poly().
返回的对象中的多项式
从上面已经给出的表达式，我们观察到gamma(j)恰好是两个连续归一化常数的比值，即：gamma(j) = c(j) / c(j-1).
我们可以通过计算来检查我们的 gamma1 值是否与这个比率一致：
gamma1 == attr(poly(x,2), "coefs")$norm2[3] / attr(poly(x,2), "coefs")$norm2[2]
其中 returns TRUE.
回到寻找模型预测值的最大值 的问题，我们可以：
将预测值表示为 r(1,x) 和 r(2,x) 以及逻辑回归系数的函数，即：
pred(x) = beta0 + beta1 * r(1,x) + beta2 * r(2,x)
导出表达式w.r.t。 x，将其设置为0并求解x。
在 R 代码中：
# Get the normalization constants alpha(j) to obtain r(j,x) from p(j,x) as
# r(j,x) = p(j,x) / sqrt( norm(j) ) = p(j,x) / alpha(j)
alpha1 = sqrt( attr(poly(x,2), "coefs")$norm2[3] )
alpha2 = sqrt( attr(poly(x,2), "coefs")$norm2[4] )
# Get the logistic regression coefficients (beta1 and beta2)
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
# Compute the x at which the maximum occurs from the expression that is obtained
# by deriving the predicted expression pred(x) = beta0 + beta1*r(1,x) + beta2*r(2,x)
# w.r.t. x and setting the derivative to 0.
xmax = ( alpha2^-1 * coef2 * (rho1 + rho2) - alpha1^-1 * coef1 ) / (2 * alpha2^-1 * coef2)
给出：
> xmax
[1] 97.501114
即与我的 .
中描述的其他“经验”方法 获得的相同值
获取最大预测值位置 x 的完整代码，从您提供的代码开始，是：
# First we define the number of observations in the data
n = length(x)
# Parameters for p1(x):
# p1(x) = (x - rho1) p0(x) (since p_{-1}(x) = 0)
rho1 = mean(x)
# Parameters for p2(x)
# p2(x) = (x - rho2) p1(x) - gamma1
gamma1 = var(x) * (n-1)/n
rho2 = mean( x * (x - mean(x))^2 ) / gamma1
# Get the normalization constants alpha(j) to obtain r(j,x) from p(j,x) as
# r(j,x) = p(j,x) / sqrt( norm(j) ) = p(j,x) / alpha(j)
alpha1 = sqrt( attr(poly(x,2), "coefs")$norm2[3] )
alpha2 = sqrt( attr(poly(x,2), "coefs")$norm2[4] )
# Get the logistic regression coefficients (beta1 and beta2)
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
# Compute the x at which the maximum occurs from the expression that is obtained
# by deriving the predicted expression pred(x) = beta0 + beta1*r(1,x) + beta2*r(2,x)
# w.r.t. x and setting the derivative to 0.
( xmax = ( alpha2^-1 * coef2 * (rho1 + rho2) - alpha1^-1 * coef1 ) / (2 * alpha2^-1 * coef2) )
我试图弄清楚具有多项式项的逻辑回归系数与预测的关系。具体来说，我对 x 轴上预测最高的位置感兴趣。下面的示例：
set.seed(42)
# Setup some dummy data
x <- 1:200
y <- rep(0, length(x))
y[51:150] <- rbinom(100, 1, 0.5)
# Fit a model
family <- binomial()
model <- glm(y ~ poly(x, 2), family = family)
# Illustrate model
plot(x, y)
lines(x, family$linkinv(predict(model)), col = 2)
上面的模型给出了这些系数：
coef(model)
#> (Intercept) poly(x, 2)1 poly(x, 2)2
#> -1.990317 -3.867855 -33.299893
由 reprex package (v1.0.0)
于 2021-08-03 创建poly() 的手册页说明如下：
The orthogonal polynomial is summarized by the coefficients, which can be used to evaluate it via the three-term recursion given in Kennedy & Gentle (1980, pp. 343–4), and used in the predict part of the code.
但是，我无法访问这本书，也无法从 predict.glm S3 方法中辨别这些系数是如何处理的。有没有一种方法可以仅根据系数（即不使用 predict() 来找到最大值）来重建山顶的位置（示例中大约为 100）？
假设您想通过分析找到预测的最大值，对于正交多项式为 1 阶和 2 阶的特定情况，我建议采用以下方法：
摘要
1) 推断多项式系数
这可以通过将线性模型拟合到模型矩阵中包含的相应多项式值来轻松完成。
2) 推导预测表达式w.r.t。 x 并将导数设置为 0
求解由(1)中的多项式拟合推导出的预测表达式中的x，得到预测最大值出现的x的值。
详情
1) 多项式系数
根据您拟合 GLM 模型的行，我们估计 1 阶多项式的系数 p1(x) = a0 + a1*x 和 2 阶多项式的系数 p2(x) = b0 + b1*x + b2*x^2：
X = model.matrix(model)
p1 = X[, "poly(x, 2)1"]
p2 = X[, "poly(x, 2)2"]
p1.lm = lm(p1 ~ x)
a0 = p1.lm$coeff["(Intercept)"]
a1 = p1.lm$coeff["x"]
p2.lm = lm(p2 ~ x + I(x^2))
b0 = p2.lm$coeff["(Intercept)"]
b1 = p2.lm$coeff["x"]
b2 = p2.lm$coeff["I(x^2)"]
这给出：
> c(a0, a1, b0, b1, b2)
(Intercept) x (Intercept) x I(x^2)
-1.2308840e-01 1.2247602e-03 1.6050353e-01 -4.7674315e-03 2.3718565e-05
2) 求最大值的预测导数
预测表达式 z（在反 link 函数之前）是：
z = Intercept + coef1 * p1(x) + coef2 * p2(x)
我们推导这个表达式，将其置0得到：
coef1 * a1 + coef2 * (b1 + 2 * b2 * xmax) = 0
求解 xmax 我们得到：
xmax = - (coef1 * a1 + coef2 * b1) / (2 * coef2 * b2)
在 R 代码中，计算如下：
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
(xmax = - ( coef1 * a1 + coef2 * b1 ) / (2 * coef2 * b2))
给出：
x
97.501114
检查
我们可以通过将最大值作为绿色十字添加到预测曲线来验证最大值：
# Prediction curve computed analytically
Intercept = model$coeff["(Intercept)"]
pred.analytical = family$linkinv( Intercept + coef1 * p1 + coef2 * p2 )
# Find the prediction's maximum analytically
pred.max = family$linkinv( Intercept + coef1 * (a0 + a1 * xmax) +
coef2 * (b0 + b1 * xmax + b2 * xmax^2) )
# Plot
plot(x, y)
# The following two lines should coincide!
lines(x, pred.analytical, col = 3)
lines(x, family$linkinv(predict(model)), col = 2)
# Location of the maximum!
points(xmax, pred.max, pch="x", col="green")
给出：
从正交多项式的理论表达式推导预测最大值的位置
我得到了 Kennedy 和 Gentle（1982 年） 的《统计计算》一书的副本，在 poly 的文档中引用了它，现在分享我的发现正交多项式的计算，以及我们如何使用它们找到最大预测值的位置 x。
书中介绍的正交多项式 (pp. 343-4) 是一元的（即最高阶系数始终为 1）并通过以下递归过程获得：
其中 q 是考虑的正交多项式的数量。
请注意以下上述术语与poly文档的关系：
- 出现在你的问题摘录中的“三项递归”是第三个表达式的RHS，它恰好有三个项。
- 第三个表达式中的rho(j+1)系数称为“中心常数”.
- 第三个表达式中的 gamma(j) 系数在文档中没有名称，但它们与 “归一化常数”直接相关。 =126=]，如下所示。
作为参考，我在这里粘贴了 poly 文档“值”部分的相关摘录：
A matrix with rows corresponding to points in x and columns corresponding to the degree, with attributes "degree" specifying the degrees of the columns and (unless raw = TRUE) "coefs" which contains the centering and normalization constants used in constructing the orthogonal polynomials
回到递归，我们可以推导出参数rho(j+1)和的值gamma(j) 通过在 p(j+1) w.r.t 上施加正交条件来从第三个表达式中得出。 p(j) 和 p(j-1).
（重要的是要注意，正交性条件 不是 积分，而是 n 观察到的 x 点，因此多项式系数取决于数据！——例如 Tchebyshev 正交多项式并非如此）。
参数表达式变为：
对于回归中使用的 1 阶和 2 阶多项式，我们得到以下表达式，已用 R 代码编写：
# First we define the number of observations in the data
n = length(x)
# For p1(x):
# p1(x) = (x - rho1) p0(x) (since p_{-1}(x) = 0)
rho1 = mean(x)
# For p2(x)
# p2(x) = (x - rho2) p1(x) - gamma1
gamma1 = var(x) * (n-1)/n
rho2 = sum( x * (x - mean(x))^2 ) / (n*gamma1)
为此我们得到：
> c(rho1, rho2, gamma1)
[1] 100.50 100.50 3333.25
注意 poly(x,2) 的 coefs 属性是：
> attr(poly(x,2), "coefs")
$alpha
[1] 100.5 100.5
$norm2
[1] 1 200 666650 1777555560
其中 $alpha 包含 居中常数 ，即 rho 值（与我们的值一致——顺便说一句，所有居中当 x 的 分布对于任何 q 是对称的 时，常数等于 x 的平均值！（观察和证明））和 $norm2 包含 归一化常数 （在本例中为 p(-1,x), p(0,x) 、p(1,x)、p(2,x))，即常数c(j) 将递归公式中的多项式归一化——通过将它们除以 sqrt(c(j))——，使生成的多项式 r(j,x) 满足 sum_over_i{ r(j,x_i)^2 } = 1; 注意 r(j,x) 是存储在 poly().
从上面已经给出的表达式，我们观察到gamma(j)恰好是两个连续归一化常数的比值，即：gamma(j) = c(j) / c(j-1).
我们可以通过计算来检查我们的 gamma1 值是否与这个比率一致：
gamma1 == attr(poly(x,2), "coefs")$norm2[3] / attr(poly(x,2), "coefs")$norm2[2]
其中 returns TRUE.
回到寻找模型预测值的最大值 的问题，我们可以：
将预测值表示为 r(1,x) 和 r(2,x) 以及逻辑回归系数的函数，即：
pred(x) = beta0 + beta1 * r(1,x) + beta2 * r(2,x)
导出表达式w.r.t。 x，将其设置为0并求解x。
在 R 代码中：
# Get the normalization constants alpha(j) to obtain r(j,x) from p(j,x) as
# r(j,x) = p(j,x) / sqrt( norm(j) ) = p(j,x) / alpha(j)
alpha1 = sqrt( attr(poly(x,2), "coefs")$norm2[3] )
alpha2 = sqrt( attr(poly(x,2), "coefs")$norm2[4] )
# Get the logistic regression coefficients (beta1 and beta2)
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
# Compute the x at which the maximum occurs from the expression that is obtained
# by deriving the predicted expression pred(x) = beta0 + beta1*r(1,x) + beta2*r(2,x)
# w.r.t. x and setting the derivative to 0.
xmax = ( alpha2^-1 * coef2 * (rho1 + rho2) - alpha1^-1 * coef1 ) / (2 * alpha2^-1 * coef2)
给出：
> xmax
[1] 97.501114
即与我的
获取最大预测值位置 x 的完整代码，从您提供的代码开始，是：
# First we define the number of observations in the data
n = length(x)
# Parameters for p1(x):
# p1(x) = (x - rho1) p0(x) (since p_{-1}(x) = 0)
rho1 = mean(x)
# Parameters for p2(x)
# p2(x) = (x - rho2) p1(x) - gamma1
gamma1 = var(x) * (n-1)/n
rho2 = mean( x * (x - mean(x))^2 ) / gamma1
# Get the normalization constants alpha(j) to obtain r(j,x) from p(j,x) as
# r(j,x) = p(j,x) / sqrt( norm(j) ) = p(j,x) / alpha(j)
alpha1 = sqrt( attr(poly(x,2), "coefs")$norm2[3] )
alpha2 = sqrt( attr(poly(x,2), "coefs")$norm2[4] )
# Get the logistic regression coefficients (beta1 and beta2)
coef1 = as.numeric( model$coeff["poly(x, 2)1"] )
coef2 = as.numeric( model$coeff["poly(x, 2)2"] )
# Compute the x at which the maximum occurs from the expression that is obtained
# by deriving the predicted expression pred(x) = beta0 + beta1*r(1,x) + beta2*r(2,x)
# w.r.t. x and setting the derivative to 0.
( xmax = ( alpha2^-1 * coef2 * (rho1 + rho2) - alpha1^-1 * coef1 ) / (2 * alpha2^-1 * coef2) )