我如何访问 glm() 中的色散参数估计，为什么它似乎没有使用迭代重新加权最小二乘法？
How do I access the dispersion parameter estimate in glm(), and why doesn't it seem to be using Iteratively Reweighted Least Squares?
在 5 年前关于 logLik.lm() 和 glm() 的 / 中，有人指出 R stats 模块中的代码注释表明 lm() 和 glm() 都在内部计算某种尺度或分散参数——大概是描述回归预测的观测值的估计分散的参数。
这自然会引发另一个问题：如果它真的是某个地方的拟合算法正在估计的真实参数（或者即使它只是某种隐式/有效参数），我如何从生成的拟合对象访问这个参数？
我在下面制作了一个 MWE（加上支持的设置/情节代码）：
第 1 部分 构建了一些模拟输入数据，我们将其拟合成一条直线（暗示需要两个拟合参数）。鉴于问题是关于一个隐藏的、内部建模的分散参数，我想确保拟合算法被迫做一些有趣的事情，因此 10% 的点被故意建模为异常值。如果您了解下图中显示的内容，那么您可以跳过阅读这部分代码。
第 2 部分 是 MWE 的主体，说明了我的问题要问的要点：它在输入上运行 glm()数据并检查了一些结果，证明 logLik() 声称三个参数估计值，显然不同意 glm() 似乎给出两个。
第 3 部分 只是根据输入数据和结果生成了一个补充图。包含它只是为了完整性和可重复性；你也可以跳过阅读它。
library(tidyverse)
# ---- Part 1: Simulate sample data, with 10% outliers ----
# Set x values to stepwise 1D grid
n <- 100
xlo <- 0
xhi <- 10
step <- (xhi - xlo) / n
x <- seq(xlo+step/2, xhi-step/2, step)
# Set slope and intercept of a simple line
intercept <- 2
slope <- 0.2
# By default, initially, set 100% of y observations to have standard
# deviation of 0.2, but then override 10% of those at random indices
# so that they become outliers with standard deviation of 5
frac_outlier <- 0.1
sd <- rep(0.2, n)
is_outlier <- rep(FALSE, n)
set.seed(100)
# Generates 10 non-repeating random integer indices between 1 and n
ranidx <- order(runif(n))[1:round(frac_outlier*n)]
sd[ranidx] <- 5
is_outlier[ranidx] <- TRUE
set.seed(200)
# Generate simulated y data, with 90% of points have sd = 0.2, and 10% of
# points having sd = 5 (i.e., they are outliers)
y <- intercept + slope * x + rnorm(n, sd=sd)
# Package simulated data into a tibble
tbl <- tibble(x, y, is_outlier)
# ---- Part 2: Perform regression and demonstrate unexpected results ----
# Estimate fit parameters using glm
fitobj <- glm(formula = y ~ x, data = tbl)
# Verify that the example fit produces two fit coefficients
message("\nEstimated fit parameters are:")
print(fitobj$coefficients)
# Verify that the log-likelihood function says there are effectively three
# fit parameters (i.e., the 3rd estimated parameter presumably being the
# effective "scale" or "dispersion" parameter of the y observation error
# distribution, as described in another Whosebug question)
ll <- logLik(fitobj)
message(paste0("\nlogLik number of estimated fit parameters: ", attr(ll, "df")))
# Verify that the Iteratively Reweighted Least Squares algorithm promised
# by the glm documentation does not actually seem to perform any re-weighting
message(paste0("\nNumber of weights equal to 1: ", sum(fitobj$weights == 1)))
message(paste0("Number of weights not equal to 1: ", sum(fitobj$weights != 1)))
# Print fit object attributes. Crux of this post: where in this list would
# we find the estimated scale parameter (i.e., parameter #3) which
# characterizes the glm algorithm's internal estimate of the y observation
# error?
message("\nFit object attributes:")
print(attributes(fitobj))
# ---- Part 3 (supplementary): Plot results to use as extra visual aid ----
plt <- ggplot(tbl, aes(x=x, y=y)) +
geom_point(aes(color=is_outlier)) +
scale_color_manual(values=c("black", "red")) +
geom_smooth(method = glm, formula = y ~ x)
print(plt)
这是上面代码片段的第 2 部分 生成的控制台输出：
Estimated fit parameters are:
(Intercept) x
2.1012286 0.1869693
logLik number of estimated fit parameters: 3
Number of weights equal to 1: 100
Number of weights not equal to 1: 0
Fit object attributes:
$names
[1] "coefficients" "residuals" "fitted.values" "effects"
[5] "R" "rank" "qr" "family"
[9] "linear.predictors" "deviance" "aic" "null.deviance"
[13] "iter" "weights" "prior.weights" "df.residual"
[17] "df.null" "y" "converged" "boundary"
[21] "model" "call" "formula" "terms"
[25] "data" "offset" "control" "method"
[29] "contrasts" "xlevels"
$class
[1] "glm" "lm"
如果它有帮助，这里有一些模拟数据和拟合结果的可视化：
我的主要问题： 在包含拟合对象的众多 attributes() 和 names() 中，我如何访问被计为的色散参数logLik() 函数的第三个参数？
此外，作为一个附加的附带问题，official R documentation for glm states that "The default method "glm.fit" uses iteratively reweighted least squares (IWLS)." The wikipedia entry for IWLS 指出它是一种稳健的回归技术，旨在最大限度地减少异常值的影响，它通过重新计算与观察相关的后验权重来实现这一点。但是，正如您在控制台输出中看到的那样，所有后验权重都设置为 1；即，glm() 中的 IWLS 功能似乎没有被触发，尽管模拟输入数据包含明显的异常值。
额外的奖励问题：这里发生了什么，为什么 glm() 的行为似乎不像文档建议的那样？
编辑：我问这个问题的目的是寻求进一步澄清 R stats 模块如何处理计算模型中参数的数量，特别是探索它如何估计拟合模型的色散参数。
我最终在评论部分与回答我问题的贡献者进行了有趣的后续交流。根据他的建议，我将那个简短的线程复制粘贴到这里，因为我同意他链接的附加信息显着增加了我对统计包开发人员处理这个问题的更广泛背景的理解。
问：你的表达式sqrt(deviance(obj)/(nobs(obj)-length(coef(obj))))不是一个独立可调参数；即，不是您可以转向以进一步增加 argmax(LL) 的“旋钮”。我断言：1.) attr(logLik(fitobj), "df") 当它说有 3 个参数时是混合概念，因为 2 个是独立可调的，而第 3 个（分散）不是； 2.) 如果模型选择包使用“df”属性来计算 AIC/BIC，如 5 年前链接的原始问题，那么它是错误的，因为 AIC/BIC 应该只计算独立可调的参数。你同意还是不同意？
答：我不同意。 （要进一步讨论，您可能应该在 CrossValidated 上提问（我在那里找不到）。
在高斯 glm() 拟合的情况下，summary() 报告的色散参数是均方误差。如果您使用 lm() 拟合模型，则其在报告摘要中的等价物将是残差标准误差，即它的平方根。
您可以使用
从您的 glm() 对象计算报告的分散度 parameter/MSE
sum((fitobj$residuals^2))/fitobj$df.r
#> [1] 2.336175
sigma(fitobj) 是 generic 从 R 中的拟合模型检索色散参数的方法（适用于 lm 和 glm 对象，以及 lme4::lmer() 等其他人）。这里“分散参数”被定义为残差标准差，因此sigma(fitobj)^2匹配摘要中打印的值（也匹配@Kieran's answer that computes the mean residual sum of直接平方）。- 一般来说，对于 GLM（和 LM），色散参数没有明确估计（这部分是为什么它没有包含在系数向量中）；
stats:::sigma.default() 使用（大约）
sqrt(deviance(object)/(nobs(object)-length(coef(object)))
对于线性模型，偏差是残差平方和。 （也可以使用 Pearson 残差平方和除以 (n-p) 和其他各种方法来估计色散参数……所有这些估计都渐近地收敛到相同的答案，但它们的 属性 有限样本量。）
- 我认为您误解了 IWLS 在 GLM 中的作用。您引用的维基百科文章说
IRLS is used to find the maximum likelihood estimates of a generalized linear model, and in robust regression to find an M-estimator, as a way of mitigating the influence of outliers in an otherwise normally-distributed data set.
最后一个条款（“作为减轻异常值影响的一种方式......”）仅适用于 second 条款，关于稳健回归——它是 不是 IWLS 在 GLM 上下文中应用时的作用。在 GLM 中，IWLS 根据方差-均值关系假设模型下的预期方差 对观察的影响重新加权 （例如，在拟合泊松模型时，我们假设方差等于到预测的平均值）；在这种情况下，不会自适应地处理异常值。在最后一步计算的“权重”可通过 weights(fitobj, "working") 获得（与 先前的权重 相反，后者与例如二项式变量的试验次数相关） .然而，在这种情况下，它们都正好是 1，因为您指定（默认情况下）family="gaussian"，在这种情况下，假定所有观察值都具有 相同 独立于均值的方差 → 所有权重都相同。
- 要得到你想要的，试试这个（使用稳健估计，也使用 IWLS（它称为“IRWLS”），但以不同的方式）：
library(robustbase)
fitobj_rob <- lmrob(y~x, data=tbl)
w <- weights(fitobj_rob, type="robustness")
hist(w, main="")
这表明，如您所料，有几个点的权重非常低（因为您将它们设置为异常值），而许多点的权重相对较高。
在 5 年前关于 logLik.lm() 和 glm() 的 lm() 和 glm() 都在内部计算某种尺度或分散参数——大概是描述回归预测的观测值的估计分散的参数。
这自然会引发另一个问题：如果它真的是某个地方的拟合算法正在估计的真实参数（或者即使它只是某种隐式/有效参数），我如何从生成的拟合对象访问这个参数？
我在下面制作了一个 MWE（加上支持的设置/情节代码）：
第 1 部分 构建了一些模拟输入数据，我们将其拟合成一条直线（暗示需要两个拟合参数）。鉴于问题是关于一个隐藏的、内部建模的分散参数，我想确保拟合算法被迫做一些有趣的事情，因此 10% 的点被故意建模为异常值。如果您了解下图中显示的内容，那么您可以跳过阅读这部分代码。
第 2 部分 是 MWE 的主体，说明了我的问题要问的要点：它在输入上运行
glm()数据并检查了一些结果，证明logLik()声称三个参数估计值，显然不同意glm()似乎给出两个。第 3 部分 只是根据输入数据和结果生成了一个补充图。包含它只是为了完整性和可重复性；你也可以跳过阅读它。
library(tidyverse)
# ---- Part 1: Simulate sample data, with 10% outliers ----
# Set x values to stepwise 1D grid
n <- 100
xlo <- 0
xhi <- 10
step <- (xhi - xlo) / n
x <- seq(xlo+step/2, xhi-step/2, step)
# Set slope and intercept of a simple line
intercept <- 2
slope <- 0.2
# By default, initially, set 100% of y observations to have standard
# deviation of 0.2, but then override 10% of those at random indices
# so that they become outliers with standard deviation of 5
frac_outlier <- 0.1
sd <- rep(0.2, n)
is_outlier <- rep(FALSE, n)
set.seed(100)
# Generates 10 non-repeating random integer indices between 1 and n
ranidx <- order(runif(n))[1:round(frac_outlier*n)]
sd[ranidx] <- 5
is_outlier[ranidx] <- TRUE
set.seed(200)
# Generate simulated y data, with 90% of points have sd = 0.2, and 10% of
# points having sd = 5 (i.e., they are outliers)
y <- intercept + slope * x + rnorm(n, sd=sd)
# Package simulated data into a tibble
tbl <- tibble(x, y, is_outlier)
# ---- Part 2: Perform regression and demonstrate unexpected results ----
# Estimate fit parameters using glm
fitobj <- glm(formula = y ~ x, data = tbl)
# Verify that the example fit produces two fit coefficients
message("\nEstimated fit parameters are:")
print(fitobj$coefficients)
# Verify that the log-likelihood function says there are effectively three
# fit parameters (i.e., the 3rd estimated parameter presumably being the
# effective "scale" or "dispersion" parameter of the y observation error
# distribution, as described in another Whosebug question)
ll <- logLik(fitobj)
message(paste0("\nlogLik number of estimated fit parameters: ", attr(ll, "df")))
# Verify that the Iteratively Reweighted Least Squares algorithm promised
# by the glm documentation does not actually seem to perform any re-weighting
message(paste0("\nNumber of weights equal to 1: ", sum(fitobj$weights == 1)))
message(paste0("Number of weights not equal to 1: ", sum(fitobj$weights != 1)))
# Print fit object attributes. Crux of this post: where in this list would
# we find the estimated scale parameter (i.e., parameter #3) which
# characterizes the glm algorithm's internal estimate of the y observation
# error?
message("\nFit object attributes:")
print(attributes(fitobj))
# ---- Part 3 (supplementary): Plot results to use as extra visual aid ----
plt <- ggplot(tbl, aes(x=x, y=y)) +
geom_point(aes(color=is_outlier)) +
scale_color_manual(values=c("black", "red")) +
geom_smooth(method = glm, formula = y ~ x)
print(plt)
这是上面代码片段的第 2 部分 生成的控制台输出：
Estimated fit parameters are:
(Intercept) x
2.1012286 0.1869693
logLik number of estimated fit parameters: 3
Number of weights equal to 1: 100
Number of weights not equal to 1: 0
Fit object attributes:
$names
[1] "coefficients" "residuals" "fitted.values" "effects"
[5] "R" "rank" "qr" "family"
[9] "linear.predictors" "deviance" "aic" "null.deviance"
[13] "iter" "weights" "prior.weights" "df.residual"
[17] "df.null" "y" "converged" "boundary"
[21] "model" "call" "formula" "terms"
[25] "data" "offset" "control" "method"
[29] "contrasts" "xlevels"
$class
[1] "glm" "lm"
如果它有帮助，这里有一些模拟数据和拟合结果的可视化：
我的主要问题： 在包含拟合对象的众多 attributes() 和 names() 中，我如何访问被计为的色散参数logLik() 函数的第三个参数？
此外，作为一个附加的附带问题，official R documentation for glm states that "The default method "glm.fit" uses iteratively reweighted least squares (IWLS)." The wikipedia entry for IWLS 指出它是一种稳健的回归技术，旨在最大限度地减少异常值的影响，它通过重新计算与观察相关的后验权重来实现这一点。但是，正如您在控制台输出中看到的那样，所有后验权重都设置为 1；即，glm() 中的 IWLS 功能似乎没有被触发，尽管模拟输入数据包含明显的异常值。
额外的奖励问题：这里发生了什么，为什么 glm() 的行为似乎不像文档建议的那样？
编辑：我问这个问题的目的是寻求进一步澄清 R stats 模块如何处理计算模型中参数的数量，特别是探索它如何估计拟合模型的色散参数。
我最终在评论部分与回答我问题的贡献者进行了有趣的后续交流。根据他的建议，我将那个简短的线程复制粘贴到这里，因为我同意他链接的附加信息显着增加了我对统计包开发人员处理这个问题的更广泛背景的理解。
问：你的表达式sqrt(deviance(obj)/(nobs(obj)-length(coef(obj))))不是一个独立可调参数；即，不是您可以转向以进一步增加 argmax(LL) 的“旋钮”。我断言：1.) attr(logLik(fitobj), "df") 当它说有 3 个参数时是混合概念，因为 2 个是独立可调的，而第 3 个（分散）不是； 2.) 如果模型选择包使用“df”属性来计算 AIC/BIC，如 5 年前链接的原始问题，那么它是错误的，因为 AIC/BIC 应该只计算独立可调的参数。你同意还是不同意？
答：我不同意。
在高斯 glm() 拟合的情况下，summary() 报告的色散参数是均方误差。如果您使用 lm() 拟合模型，则其在报告摘要中的等价物将是残差标准误差，即它的平方根。
您可以使用
从您的glm() 对象计算报告的分散度 parameter/MSE
sum((fitobj$residuals^2))/fitobj$df.r
#> [1] 2.336175
sigma(fitobj)是 generic 从 R 中的拟合模型检索色散参数的方法（适用于lm和glm对象，以及lme4::lmer()等其他人）。这里“分散参数”被定义为残差标准差，因此sigma(fitobj)^2匹配摘要中打印的值（也匹配@Kieran's answer that computes the mean residual sum of直接平方）。- 一般来说，对于 GLM（和 LM），色散参数没有明确估计（这部分是为什么它没有包含在系数向量中）；
stats:::sigma.default()使用（大约）
sqrt(deviance(object)/(nobs(object)-length(coef(object)))
对于线性模型，偏差是残差平方和。 （也可以使用 Pearson 残差平方和除以 (n-p) 和其他各种方法来估计色散参数……所有这些估计都渐近地收敛到相同的答案，但它们的 属性 有限样本量。）
- 我认为您误解了 IWLS 在 GLM 中的作用。您引用的维基百科文章说
IRLS is used to find the maximum likelihood estimates of a generalized linear model, and in robust regression to find an M-estimator, as a way of mitigating the influence of outliers in an otherwise normally-distributed data set.
最后一个条款（“作为减轻异常值影响的一种方式......”）仅适用于 second 条款，关于稳健回归——它是 不是 IWLS 在 GLM 上下文中应用时的作用。在 GLM 中，IWLS 根据方差-均值关系假设模型下的预期方差 对观察的影响重新加权 （例如，在拟合泊松模型时，我们假设方差等于到预测的平均值）；在这种情况下，不会自适应地处理异常值。在最后一步计算的“权重”可通过 weights(fitobj, "working") 获得（与 先前的权重 相反，后者与例如二项式变量的试验次数相关） .然而，在这种情况下，它们都正好是 1，因为您指定（默认情况下）family="gaussian"，在这种情况下，假定所有观察值都具有 相同 独立于均值的方差 → 所有权重都相同。
- 要得到你想要的，试试这个（使用稳健估计，也使用 IWLS（它称为“IRWLS”），但以不同的方式）：
library(robustbase)
fitobj_rob <- lmrob(y~x, data=tbl)
w <- weights(fitobj_rob, type="robustness")
hist(w, main="")
这表明，如您所料，有几个点的权重非常低（因为您将它们设置为异常值），而许多点的权重相对较高。