我如何访问 glm() 中的离散参数估计，为什么它似乎没有使用迭代重新加权最小二乘法？

如何解决我如何访问 glm() 中的离散参数估计，为什么它似乎没有使用迭代重新加权最小二乘法？

在这个 5 年前关于 logLik.lm() 和 glm() 的 question / answer 中，有人指出 R 统计模块中的代码注释表明 {{1 }} 和 lm() 都在内部计算某种尺度或离散参数——大概是描述回归预测的观测值的估计离散的参数。

这自然会引出另一个问题：如果它确实是某个地方的拟合算法估计的真实参数（或者即使它只是某种隐式/有效参数），我如何从生成的拟合对象中访问该参数？

我在下面制作了一个 MWE（加上支持的设置/绘图代码）：

• 第 1 部分构建了一些模拟输入数据，我们将这些数据拟合成一条直线（意味着需要两个拟合参数）。鉴于问题是关于一个隐藏的、内部建模的分散参数，我想确保拟合算法被迫做一些有趣的事情，因此 10% 的点被故意建模为异常值。如果您理解下图中显示的内容，那么您可能可以跳过阅读这部分代码。

• 第 2 部分是 MWE 的主体，说明了我的问题要问的点：它对输入数据运行 glm() 并检查一些结果，证明 glm() 声明了三个参数估计，显然与 logLik() 不一致，后者似乎给出了两个。

• 第 3 部分 仅根据输入数据和结果生成一个小补充图。它只是为了完整性和可重复性而包含在内；你也可以跳过阅读。

glm()

这是上面代码片段的第 2 部分生成的控制台输出：

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 stackoverflow 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)

如果有帮助，下面是模拟数据和拟合结果的一些可视化：

我的主要问题：在包含拟合对象的众多 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() 函数？

此外，作为一个额外的附带问题，official R documentation for glm 声明“默认方法“glm.fit”使用迭代重新加权最小二乘法 (IWLS)。” wikipedia entry for IWLS 指出它是一种稳健的回归技术，旨在最小化异常值的影响，它通过重新计算与观察相关的后验权重来实现这一点。但是，正如您在控制台输出中看到的，所有后验权重都设置为等于 1；即，logLik() 中的 IWLS 功能似乎没有被触发，尽管模拟输入数据包含明显的异常值。

额外的问题：这里发生了什么，为什么 glm() 似乎不像文档建议的那样行事？

编辑：我提出这个问题的目的是寻求进一步澄清 R stats 模块如何处理计算模型中参数的数量，特别是探索它如何估计拟合模型的离散参数。

我最终在评论部分与回答我问题的贡献者进行了有趣的后续交流。在他的建议下，我将这个简短的主题复制粘贴到这里，因为我同意他链接的其他信息显着增加了我对 stats 包开发人员处理此问题的更广泛背景的理解。

问：你的表达式 glm() 不是一个独立可调参数；即，不是“旋钮”，您可以转向进一步增加 argmax(LL)。我断言：1.) sqrt(deviance(obj)/(nobs(obj)-length(coef(obj)))) 当它说有 3 个参数时是混合概念，因为 2 个是独立可调的，而第 3 个（分散）不是； 2.) 如果模型选择包使用“df”属性来计算 AIC/BIC，如 5 年前链接的原始问题，那么它是错误的，因为 AIC/BIC 应该只计算独立可调参数。你同意还是不同意？

答：我不同意。 https://stackoverflow.com/a/67134012/3083138（为了进一步讨论，您可能应该在 CrossValidated 上问一个问题（我已经在那里找不到了）。

解决方法

• sigma(fitobj) 是从 R 中的拟合模型中检索离散参数的通用方法（适用于 lm 和 glm 对象，以及其他如 lme4::lmer() 等）。这里“离散参数”被定义为残差标准偏差，所以 sigma(fitobj)^2 匹配摘要中打印的值（也匹配 @Kieran 直接计算平均残差平方和的答案） .
• 一般来说，对于 GLM（和 LM），没有明确估计色散参数（这也是它不包含在系数向量中的部分原因）； stats:::sigma.default() 使用（大约）

sqrt(deviance(object)/(nobs(object)-length(coef(object)))

对于线性模型，偏差是残差平方和。 （也可以使用 Pearson 残差平方和除以 (n-p) 以及各种其他方法来估计色散参数……所有这些估计都渐近收敛到相同的答案，但对于有限样本大小的属性不同。）

• 我认为您误解了 IWLS 在 GLM 中的作用。您引用的维基百科文章说

IRLS 用于寻找广义线性模型的最大似然估计，并在稳健回归中寻找 M 估计量，作为减轻其他正态分布数据集中异常值影响的一种方式。

最后一个子句（“作为减轻异常值影响的一种方式......”）仅适用于关于稳健回归的第二个子句——它不是什么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="")

这表明，如您所料，有些点的权重非常低（因为您将它们设置为异常值），而许多点的权重相对较高。

,

在高斯 glm() 拟合的情况下，summary() 报告的色散参数是均方误差。如果您使用 lm() 拟合模型，则报告摘要中的等效项将是残差标准误差，即其平方根。

您可以使用

从您的 glm() 对象计算报告的分散参数/MSE

sum((fitobj$residuals^2))/fitobj$df.r

#> [1] 2.336175

This Cross-Validated question and answer 可能有用。