Quel algorithme d'optimisation est utilisé dans la fonction glm dans R?

17

On peut effectuer une régression logit dans R en utilisant un tel code:

> 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

Il semble que l'algorithme d'optimisation ait convergé - il existe des informations sur le nombre d'étapes de l'algorithme de notation du pêcheur:

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

Je suis curieux de savoir de quel algorithme optim il s'agit? S'agit-il d'un algorithme de Newton-Raphson (descente de gradient de second ordre)? Puis-je définir certains paramètres pour utiliser l'algorithme de Cauchy (descente de gradient de premier ordre)?

Marcin Kosiński
la source
5
Cela vous dérange de citer où un algorithme de Newton-Raphson est appelé "niveau de descente de gradient II"?
Cliff AB du
4
Le score de FIsher lui-même est lié à Newton-Raphson, mais différent de celui-ci, remplaçant en fait le Hessian par sa valeur attendue dans le modèle.
Glen_b -Reinstate Monica
@CliffAB sory. Je mentionne qu'il Newton's methods'agit d'une méthode de descente de gradient de second ordre.
Marcin Kosiński
1
@ hxd1011, vous ne devez pas modifier pour modifier la question de quelqu'un d'autre. C'est une chose à modifier lorsque vous pensez savoir ce qu'ils signifient, mais leur question n'est pas claire (peut-être que l'anglais n'est pas leur langue maternelle, par exemple), mais vous ne devriez pas rendre leur question différente (par exemple, plus générale) que ce qu'ils avaient voulait. Au lieu de cela, posez une nouvelle question avec ce que vous voulez. J'annule l'édition.
gung - Rétablir Monica
1
@ La méthode de MarcinKosiński Newton est Newton-Raphson, Raphson s'est simplement basé sur les idées de Newton pour un cas plus général.
AdamO

Réponses:

19

Vous serez intéressé de savoir que la documentation de glm, accessible via ?glmfournit de nombreuses informations utiles: methodnous trouvons ci- dessous que les moindres carrés itérativement repondérés sont la méthode par défaut pour glm.fit, qui est la fonction de cheval de bataille pour glm. En outre, la documentation mentionne que des fonctions définies par l'utilisateur peuvent être fournies ici, au lieu de la valeur par défaut.

Sycorax dit de réintégrer Monica
la source
3
Vous pouvez également simplement taper le nom de la fonction glmou fit.glmà l' Rinvite pour étudier le code source.
Matthew Drury
@MatthewDrury Je pense que vous voulez dire que le cheval de bataille glm.fitqui ne sera pas entièrement reproductible puisqu'il repose sur le code C C_Cdqrls.
AdamO
Yah, tu as raison, je mélange beaucoup l'ordre dans R. Que voulez-vous dire par non reproductible?
Matthew Drury
16

La méthode utilisée est mentionnée dans la sortie elle-même: c'est le Fisher Scoring. Cela équivaut à Newton-Raphson dans la plupart des cas. L'exception étant les situations où vous utilisez des paramétrisations non naturelles. La régression du risque relatif est un exemple d'un tel scénario. Là, les informations attendues et observées sont différentes. En général, Newton Raphson et Fisher Scoring donnent des résultats presque identiques.

Une autre formulation de la notation de Fisher est celle des moindres carrés itérativement repondérés. Pour donner une certaine intuition, les modèles d'erreur non uniformes ont le modèle des moindres carrés pondérés par la variance inverse comme modèle "optimal" selon le théorème de Gauss Markov. Avec les GLM, il existe une relation moyenne-variance connue. Un exemple est la régression logistique où la moyenne est et la variance estpp(1p)Fisher Scoring. De plus, il donne une belle intuition à l'algorithme EM qui est un cadre plus général pour estimer les probabilités compliquées.

L'optimiseur général par défaut dans R utilise des méthodes numériques pour estimer un second moment, essentiellement basé sur une linéarisation (attention à la malédiction de la dimensionnalité). Donc, si vous vouliez comparer l'efficacité et le biais, vous pourriez mettre en œuvre une routine de naïveté logistique naïve avec quelque chose comme

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)

Donne moi

> ans$estimate
[1] -2.2261225  0.1651472
> coef(glm(y~x, family=binomial))
(Intercept)           x 
 -2.2261215   0.1651474 
AdamO
la source
quelle est la différence par rapport à un gradient décent / méthode de newton / BFGS? Je pense que vous avez expliqué, mais je ne suis pas tout à fait suivre.
Haitao Du
@ hxd1011 voir "Méthodes numériques pour l'optimisation sans contrainte et les équations non linéaires" Dennis, JE et Schnabel, RB
AdamO
1
@ hxd1011 pour autant que je sache, Newton Raphson n'a pas besoin ou estime un Hessian dans les étapes. La méthode de Newton-Type permet d' nlmestimer numériquement le gradient puis applique Newton Raphson. Dans BFGS, je pense que le gradient est requis comme avec Newton Raphson, mais les étapes successives sont évaluées en utilisant une approximation de second ordre qui nécessite une estimation de la Hesse. BFGS est bon pour les optimisations hautement non linéaires. Mais pour les GLM, ils se comportent généralement très bien.
AdamO
2
La réponse existante mentionnait "les moindres carrés repondérés de manière itérative", comment cela entre-t-il dans l'image, par rapport aux algorithmes que vous avez mentionnés ici?
Amoeba dit Reinstate Monica
@AdamO pourriez-vous répondre aux commentaires d'Amoeba? Merci
Haitao Du
1

Pour mémoire, une implémentation pure R simple de l'algorithme glm de R, basée sur le score de Fisher (moindres carrés itérativement repondérés), comme expliqué dans l'autre réponse, est donnée par:

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

Exemple:

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

Une bonne discussion des algorithmes d'ajustement GLM, y compris une comparaison avec Newton-Raphson (qui utilise la Hesse observée par opposition à la Hesse attendue dans l'algorithme IRLS) et les algorithmes hybrides (qui commencent par IRLS, car ils sont plus faciles à initialiser, mais ensuite terminer avec une optimisation supplémentaire en utilisant Newton-Raphson) peut être trouvé dans le livre "Generalized Linear Models and Extensions" par James W. Hardin & Joseph M. Hilbe .

Tom Wenseleers
la source