¿Qué algoritmo de optimización se usa en la función glm en R?

Se puede realizar una regresión logit en R usando dicho código:

> 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

Parece que el algoritmo de optimización ha convergido: hay información sobre el número de pasos del algoritmo de puntuación de pescador:

Call:
glm(formula = cbind(Menarche, Total - Menarche) ~ Age, family = binomial(logit), 
    data = menarche)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0363  -0.9953  -0.4900   0.7780   1.3675  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -21.22639    0.77068  -27.54   <2e-16 ***
Age           1.63197    0.05895   27.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3693.884  on 24  degrees of freedom
Residual deviance:   26.703  on 23  degrees of freedom
AIC: 114.76

Number of Fisher Scoring iterations: 4

Tengo curiosidad acerca de qué algoritmo optim es? ¿Es el algoritmo Newton-Raphson (descenso de gradiente de segundo orden)? ¿Puedo establecer algunos parámetros para usar el algoritmo Cauchy (descenso de gradiente de primer orden)?

Le interesará saber que la documentación a la que se glmaccede a través de ?glmproporciona muchas ideas útiles: a methodcontinuación, descubrimos que los mínimos cuadrados repesados ​​de forma iterativa son el método predeterminado para glm.fit, que es la función de caballo de batalla glm. Además, la documentación menciona que aquí se pueden proporcionar funciones definidas por el usuario, en lugar de las predeterminadas.

El método utilizado se menciona en el resultado en sí: es Fisher Scoring. Esto es equivalente a Newton-Raphson en la mayoría de los casos. La excepción son las situaciones en las que está utilizando parametrizaciones no naturales. La regresión del riesgo relativo es un ejemplo de tal escenario. Allí, la información esperada y observada es diferente. En general, Newton Raphson y Fisher Scoring dan resultados casi idénticos.

$p$$p(1-p)$Puntuación de pescadores. Además, le da una buena intuición al algoritmo EM, que es un marco más general para estimar probabilidades complicadas.

El optimizador general predeterminado en R usa métodos numéricos para estimar un segundo momento, básicamente basado en una linealización (tenga cuidado con la maldición de la dimensionalidad). Entonces, si estaba interesado en comparar la eficiencia y el sesgo, podría implementar una rutina ingenua de máxima probabilidad logística con algo como

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)

me da

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

Para el registro, una implementación R pura y simple del algoritmo glm de R, basada en la puntuación de Fisher (mínimos cuadrados repesados ​​de forma iterativa), como se explica en la otra respuesta, viene dada por:

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

Ejemplo:

## 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

Una buena discusión sobre los algoritmos de ajuste GLM, incluida una comparación con Newton-Raphson (que usa el Hessian observado en lugar del Hessian esperado en el algoritmo IRLS) y algoritmos híbridos (que comienzan con IRLS, ya que estos son más fáciles de inicializar, pero luego terminar con una mayor optimización con Newton-Raphson) se puede encontrar en el libro "Modelos lineales generalizados y extensiones" de James W. Hardin y Joseph M. Hilbe .