您将有兴趣知道，glm通过访问的文档?glm提供了许多有用的见解：method我们发现迭代重新加权最小二乘法是 的默认方法glm.fit，它是 的主力函数glm。此外，文档提到这里可能会提供用户定义的函数，而不是默认的。

使用的方法在输出本身中提到：它是 Fisher Scoring。在大多数情况下，这等效于 Newton-Raphson。使用非自然参数化的情况除外。相对风险回归就是这种情况的一个例子。在那里，预期和观察到的信息是不同的。一般来说，Newton Raphson 和 Fisher Scoring 给出几乎相同的结果。

Fisher 评分的另一个公式是迭代重加权最小二乘法。直观地说，根据高斯马尔可夫定理，非均匀误差模型将逆方差加权最小二乘模型作为“最优”模型。对于 GLM，存在已知的均值-方差关系。一个例子是逻辑回归，其中平均值为$p$方差是$p(1-p)$. 因此，一种算法是通过在朴素模型中估计均值来构建的，从预测的均值创建权重，然后使用更精细的精度重新估计均值，直到收敛。事实证明，这就是费舍尔得分。此外，它为 EM 算法提供了一些很好的直觉，这是一种用于估计复杂可能性的更通用的框架。

R 中默认的通用优化器使用数值方法来估计二阶矩，基本上基于线性化（小心维度灾难）。因此，如果您对比较效率和偏差感兴趣，您可以实现一个简单的逻辑最大似然例程，例如

set.seed(1234)
x <- rnorm(1000)
y <- rbinom(1000, 1, exp(-2.3 + 0.1*x)/(1+exp(-2.3 + 0.1*x)))
f <- function(b) {
  p <- exp(b[1] + b[2]*x)/(1+exp(b[1] + b[2]*x))
  -sum(dbinom(y, 1, p, log=TRUE))
}
ans <- nlm(f, p=0:1, hessian=TRUE)

给我

> ans$estimate
[1] -2.2261225  0.1651472
> coef(glm(y~x, family=binomial))
(Intercept)           x 
 -2.2261215   0.1651474 

作为记录，R的glm算法的简单纯R实现，基于Fisher评分（迭代重新加权最小二乘法），如另一个答案中所解释的：

glm_irls = function(X, y, weights=rep(1,nrow(X)), family=poisson(log), maxit=25, tol=1e-16) {
    if (!is(family, "family")) family = family()
    variance = family$variance
    linkinv = family$linkinv
    mu.eta = family$mu.eta
    etastart = NULL

    nobs = nrow(X)    # needed by the initialize expression below
    nvars = ncol(X)   # needed by the initialize expression below
    eval(family$initialize) # initializes n and fitted values mustart
    eta = family$linkfun(mustart) # we then initialize eta with this
    dev.resids = family$dev.resids
    dev = sum(dev.resids(y, linkinv(eta), weights))
    devold = 0
    beta_old = rep(1, nvars)

    for(j in 1:maxit)
    {
      mu = linkinv(eta) 
      varg = variance(mu)
      gprime = mu.eta(eta)
      z = eta + (y - mu) / gprime # potentially -offset if you would have an offset argument as well
      W = weights * as.vector(gprime^2 / varg)
      beta = solve(crossprod(X,W*X), crossprod(X,W*z), tol=2*.Machine$double.eps)
      eta = X %*% beta # potentially +offset if you would have an offset argument as well
      dev = sum(dev.resids(y, mu, weights))
      if (abs(dev - devold) / (0.1 + abs(dev)) < tol) break
      devold = dev
      beta_old = beta
    }
    list(coefficients=t(beta), iterations=j)
}

例子：

## Dobson (1990) Page 93: Randomized Controlled Trial :
y <- counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
X <- model.matrix(counts ~ outcome + treatment)

coef(glm.fit(x=X, y=y, family = poisson(log))) 
  (Intercept)      outcome2      outcome3    treatment2    treatment3 
 3.044522e+00 -4.542553e-01 -2.929871e-01 -7.635479e-16 -9.532452e-16

coef(glm_irls(X=X, y=y, family=poisson(log)))
     (Intercept)   outcome2   outcome3    treatment2   treatment3
[1,]    3.044522 -0.4542553 -0.2929871 -3.151689e-16 -8.24099e-16

对 GLM 拟合算法的一个很好的讨论，包括与 Newton-Raphson（使用观察到的 Hessian，而不是 IRLS 算法中的预期 Hessian）和混合算法（从 IRLS 开始，因为这些更容易初始化，但随后使用 Newton-Raphson 进行进一步优化）可以在 James W. Hardin 和 Joseph M. Hilbe 的“广义线性模型和扩展”一书中找到。