有没有办法在glm函数中获取优化算法每一步的系数？

Question

当在R中执行logit回归时，可以使用

coefficients()

函数获得优化算法收敛（或不收敛）后的系数：

library(MASS)
data(menarche)
glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
               family=binomial(logit), data=menarche)
coefficients(glm.out)
## (Intercept)         Age 
## -21.226395    1.631968

有没有办法获取优化算法每一步的系数来追踪其步骤？

Answer 1

glm.fit 的内部结构已更改（请参阅@John 的评论），因此请使用它。它不依赖于内部的行位置，而是拦截 glm.fit 中 cat 的每个实例，并向迭代消息添加一条消息，因此尽管它仍然依赖于内部，但它应该不那么脆弱。这在 R 4.1 和 4.2 中对我有用。

library(MASS)
data(menarche)

trace(glm.fit, quote(cat <- function(...) {
  base::cat(...)
  if (...length() >= 3 && identical(..3, " Iterations - ")) print(coefold)
}))
glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
                     family=binomial(logit), data=menarche,
                     control = glm.control(trace = TRUE))
untrace(glm.fit)

---

之前的解决方案

带有所示值的

control=

参数会导致打印偏差，而

trace

语句将导致打印系数值：

trace(glm.fit, quote(print(coefold)), at = list(c(22, 4, 8, 4, 19, 3)))
glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
                     family=binomial(logit), data=menarche,
                     control = glm.control(trace = TRUE))

输出将如下所示：

Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
NULL
Deviance = 27.23412 Iterations - 1
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] -20.673652   1.589536
Deviance = 26.7041 Iterations - 2
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] -21.206854   1.630468
Deviance = 26.70345 Iterations - 3
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  .... step 22,4,8,4,19,3 
[1] -21.226370   1.631966
Deviance = 26.70345 Iterations - 4

要删除痕迹，请使用：

untrace(glm.fit)

请注意，在

trace

调用中，

coefold

是

glm.fit

源代码内部使用的变量名称，使用的数字指的是源代码中的语句编号，因此如果

glm.fit 则可能需要更改

来源变更。我正在使用“R 版本 3.2.2 已修补 (2015-10-19 r69550)”。

Answer 2

我建议采用手动方法。将参数修改为

glm.control

并在

while

循环中提供拟合值。如下：

library(MASS)
data(menarche)
converged <- F
coeftrace <- matrix(0, 25, 2)
i <- 1
mu <- NULL

while(!converged & i <= 25) {
  glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
                mustart = mu,
                family=binomial(logit), data=menarche, 
                control=glm.control(maxit=1))
  mu <- glm.out$fitted.values
  coeftrace[i,] <- coef(glm.out)
  i <- i+1
  converged<- glm.out$converged
}
coefficients(glm.out)

给予：

> coeftrace
           [,1]     [,2]
 [1,] -20.67365 1.589536
 [2,] -21.20685 1.630468
 [3,] -21.22637 1.631966
 [4,] -21.22639 1.631968
 [5,]   0.00000 0.000000
 [6,]   0.00000 0.000000
 [7,]   0.00000 0.000000
 [8,]   0.00000 0.000000

好处是这些值可以用于分析，例如绘图。请注意，

glm

提供的初始条件并未在此答案或 @GGrothendieck 的答案中捕获。也就是说，

glm.fit

行为是在

binomial()$initialize

环境中使用

glm.fit

来设置初始条件，并对应于系数。这种方法比提供狂野的起始参数更有效。如果我更改上面的

mu

，我可能会添加 4 个额外的迭代。

library(MASS)
data(menarche)
converged <- F
coeftrace <- matrix(0, 25, 2)
i <- 1
mu <- rep(0.1, 25)

while(!converged & i <=25) {
  glm.out = glm(cbind(Menarche, Total-Menarche) ~ Age,
                mustart = mu,
                family=binomial(logit), data=menarche, 
                control=glm.control(maxit=1))
  mu <- glm.out$fitted.values
  coeftrace[i,] <- coef(glm.out)
  i <- i+1
  converged<- glm.out$converged
}

           [,1]      [,2]
 [1,] -18.18140 1.5441246
 [2,] -10.04138 0.6936207
 [3,] -11.95066 0.9018941
 [4,] -16.27171 1.2450011
 [5,] -19.71743 1.5147288
 [6,] -21.07972 1.6206137
 [7,] -21.22499 1.6318598
 [8,] -21.22639 1.6319683
 [9,] -21.22639 1.6319683
[10,]   0.00000 0.0000000
[11,]   0.00000 0.0000000