Który algorytm optymalizacji jest używany w funkcji glm w R?

17

Można wykonać regresję logit w R przy użyciu takiego kodu:

> 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

Wygląda na to, że algorytm optymalizacji jest zbiegnięty - jest informacja o liczbie kroków algorytmu oceniania Fishera:

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

Jestem ciekawy, jaki to optymalny algorytm? Czy to algorytm Newtona-Raphsona (spadek gradientu drugiego rzędu)? Czy mogę ustawić niektóre parametry, aby korzystać z algorytmu Cauchy'ego (spadek gradientu pierwszego rzędu)?

Marcin Kosiński
źródło
5
Czy masz coś przeciwko cytowaniu, gdzie algorytm Newtona-Raphsona jest określany jako „poziom spadku gradientu II”?
Cliff AB
4
Samo ocenianie FIshera jest powiązane, ale różni się od Newtona-Raphsona, w efekcie zastępując Hesjan jego oczekiwaną wartością w modelu.
Glen_b -Reinstate Monica
@CliffAB sory. Mam na myśli, że Newton's methodjest to metoda spadku gradientu drugiego rzędu.
Marcin Kosiński,
1
@ hxd1011, nie należy edytować, aby zmienić czyjeś pytanie. Jedną rzeczą jest edytowanie, gdy myślisz, że wiesz, co one oznaczają, ale ich pytanie jest niejasne (być może angielski nie jest ich językiem ojczystym, np.), Ale nie powinieneś różnicować ich pytania (np. Bardziej ogólnie) niż mieli chciał. Zamiast tego zadaj nowe pytanie z tym, czego chcesz. Cofam edycję.
gung - Przywróć Monikę
1
@ MarcinKosiński Metodą Newtona jest Newton-Raphson, Raphson zbudował jedynie pomysły Newtona na bardziej ogólny przypadek.
AdamO,

Odpowiedzi:

19

Będziesz zainteresowany, aby dowiedzieć się, że dokumentacja, do której glmmożna uzyskać dostęp za pośrednictwem, ?glmzapewnia wiele użytecznych informacji: poniżej methodznajdujemy, że iteracyjnie przeważone najmniejsze kwadraty to metoda domyślna, dla glm.fitktórej jest funkcja konia roboczego glm. Dodatkowo w dokumentacji wspomniano, że tutaj mogą być dostarczone funkcje zdefiniowane przez użytkownika zamiast domyślnych.

Sycorax mówi Przywróć Monikę
źródło
3
Możesz także wpisać nazwę funkcji glmlub fit.glmw Rwierszu polecenia, aby przestudiować kod źródłowy.
Matthew Drury
@MatthewDrury Chyba masz na myśli siłą napędową glm.fit, która nie będzie całkowicie powtarzalne, ponieważ opiera się na kodzie C C_Cdqrls.
AdamO,
Tak, masz rację, często mieszam kolejność w R. Co masz na myśli mówiąc, że nie jest odtwarzalna?
Matthew Drury
16

Zastosowana metoda jest wymieniona w samym danych wyjściowych: jest to punktacja Fishera. W większości przypadków jest to odpowiednik Newtona-Raphsona. Wyjątkiem są sytuacje, w których stosuje się parametry nienaturalne. Względny regres ryzyka jest przykładem takiego scenariusza. Tam oczekiwane i obserwowane informacje są różne. Ogólnie rzecz biorąc, Newton Raphson i Fisher Scoring dają prawie identyczne wyniki.

pp(1-p)Punktacja Fishera. Dodatkowo daje pewną intuicję algorytmowi EM, który jest bardziej ogólną strukturą do szacowania skomplikowanych prawdopodobieństw.

Domyślny ogólny optymalizator w R używa metod numerycznych do oszacowania drugiego momentu, zasadniczo w oparciu o linearyzację (uważaj na przekleństwo wymiarowości). Więc jeśli chciałbyś porównać wydajność i stronniczość, możesz wdrożyć naiwną logistyczną procedurę maksymalnego prawdopodobieństwa z czymś takim jak

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)

daje mi

> ans$estimate
[1] -2.2261225  0.1651472
> coef(glm(y~x, family=binomial))
(Intercept)           x 
 -2.2261215   0.1651474 
AdamO
źródło
Jaka jest różnica w porównaniu z gradientem przyzwoitym / metodą Newtona / BFGS? Myślę, że wyjaśniłeś, ale nie do końca mnie śledzę.
Haitao Du
@ hxd1011 patrz „Metody numeryczne dla nieograniczonej optymalizacji i równań nieliniowych” Dennis, JE i Schnabel, RB
AdamO
1
@ hxd1011, o ile mi wiadomo, Newton Raphson nie wymaga ani nie szacuje Hesji w krokach. Metoda typu Newtona nlmszacuje gradient numerycznie, a następnie stosuje Newtona Raphsona. W BFGS uważam, że gradient jest wymagany, jak w przypadku Newtona Raphsona, ale kolejne kroki są oceniane przy użyciu aproksymacji drugiego rzędu, która wymaga oszacowania Hesji. BFGS jest dobry do wysoce nieliniowych optymalizacji. Ale w przypadku GLM są zazwyczaj bardzo dobrze wychowani.
AdamO,
2
W istniejącej odpowiedzi wspomniano „iteracyjnie przeważone najmniejsze kwadraty”, jak to wchodzi w obraz w porównaniu z algorytmami, o których tu wspomniałeś?
ameba mówi Przywróć Monikę
@AdamO czy mógłbyś odnieść się do komentarzy Amoeby? Dzięki
Haitao Du
1

Dla przypomnienia, prosta implementacja czystego R algorytmu glm R, oparta na punktacji Fishera (iteracyjnie ponownie ważona najmniejszych kwadratów), jak wyjaśniono w drugiej odpowiedzi, daje:

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

Przykład:

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

Dobra dyskusja na temat algorytmów dopasowania GLM, w tym porównanie z Newtonem-Raphsonem (który wykorzystuje obserwowany Hesjan w przeciwieństwie do oczekiwanego Hesji w algorytmie IRLS) i algorytmami hybrydowymi (które zaczynają się od IRLS, ponieważ są one łatwiejsze do zainicjowania, ale potem zakończone dalszą optymalizacją przy użyciu Newtona-Raphsona) można znaleźć w książce „Uogólnione modele liniowe i rozszerzenia” Jamesa W. Hardina i Josepha M. Hilbe .

Tom Wenseleers
źródło