Welcher Optimierungsalgorithmus wird in glm-Funktion in R verwendet?

Mit folgendem Code kann eine logit-Regression in R durchgeführt werden:

> 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

Es sieht so aus, als ob der Optimierungsalgorithmus konvergiert hat - es gibt Informationen über die Schrittanzahl des Fisher-Scoring-Algorithmus:

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

Ich bin gespannt, um welchen Optimierungsalgorithmus es sich handelt. Ist es ein Newton-Raphson-Algorithmus (Gradientenabstieg zweiter Ordnung)? Kann ich einige Parameter für die Verwendung des Cauchy-Algorithmus einstellen (Gradientenabstieg erster Ordnung)?

Es wird Sie interessieren, ob die Dokumentation für glm, auf die über zugegriffen ?glmwird, viele nützliche Einblicke bietet: Unter methodfinden wir, dass iterativ gewichtete kleinste Quadrate die Standardmethode für ist glm.fit, für die die Arbeitspferdefunktion ist glm. Darüber hinaus wird in der Dokumentation erwähnt, dass benutzerdefinierte Funktionen anstelle der Standardfunktionen bereitgestellt werden können.

Die verwendete Methode wird in der Ausgabe selbst erwähnt: Es handelt sich um Fisher Scoring. Dies entspricht in den meisten Fällen Newton-Raphson. Die Ausnahme bilden Situationen, in denen Sie nicht-natürliche Parametrisierungen verwenden. Die relative Risikorückführung ist ein Beispiel für ein solches Szenario. Dort sind die erwarteten und beobachteten Informationen unterschiedlich. Im Allgemeinen liefern Newton Raphson und Fisher Scoring nahezu identische Ergebnisse.

$p$$p(1-p)$Fisher Scoring. Darüber hinaus wird der EM-Algorithmus, der einen allgemeineren Rahmen für die Schätzung komplizierter Wahrscheinlichkeiten darstellt, ein wenig näher erläutert.

Der standardmäßige allgemeine Optimierer in R verwendet numerische Methoden zum Schätzen eines zweiten Moments, die im Wesentlichen auf einer Linearisierung basieren (Vorsicht vor dem Fluch der Dimensionalität). Wenn Sie also die Effizienz und Verzerrung vergleichen möchten, können Sie eine naive logistische Maximum-Likelihood-Routine mit so etwas wie implementieren

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)

gibt mir

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

Für die Aufzeichnung ist eine einfache reine R-Implementierung des glm-Algorithmus von R basierend auf der Fisher-Wertung (iterativ gewichtete kleinste Quadrate) wie in der anderen Antwort erklärt gegeben durch:

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

Beispiel:

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

Eine gute Diskussion der GLM-Anpassungsalgorithmen, einschließlich eines Vergleichs mit Newton-Raphson (der das beobachtete Hessische im Gegensatz zum erwarteten Hessischen im IRLS-Algorithmus verwendet) und hybriden Algorithmen (die mit IRLS beginnen, da diese dann leichter zu initialisieren sind) Finishing mit weiterer Optimierung unter Verwendung von Newton-Raphson) finden Sie im Buch "Generalized Linear Models and Extensions" von James W. Hardin & Joseph M. Hilbe .