R: funkcja glm z rodziną = specyfikacja „dwumianowa” i „waga”

Jestem bardzo zdezorientowany, jak waga działa w glm z rodziną = „dwumianowy”. W moim rozumieniu prawdopodobieństwo glm z rodziną = „dwumianowy” jest określone w następujący sposób:

f (y) = (\binom{n}{n y}) p^{n y} (1 - p)^{n (1 - y)} = \exp (n [y \log \frac{p}{1 - p} - (- \log (1 - p))] + \log (\binom{n}{n y}))

$f(y) = {n\choose{ny}} p^{ny} (1-p)^{n(1-y)} = \exp \left(n \left[ y \log \frac{p}{1-p} - \left(-\log (1-p)\right) \right] + \log {n \choose ny}\right)$ gdzie

y

$y$ to „odsetek zaobserwowanego sukcesu”, a

n

$n$ to znana liczba prób.

W moim rozumieniu prawdopodobieństwo sukcesu $p$ jest sparametryzowane za pomocą niektórych współczynników liniowych $\beta$ jako $p=p(\beta)$ i funkcji glm z wyszukiwaniem rodziny = „dwumianowe” dla:

arg max_{β} \sum_{i} \log f (y_{i}) .

$\textrm{arg}\max_{\beta} \sum_i \log f(y_i).$ Następnie ten problem optymalizacji można uprościć:

arg max_{β} \sum_{i} \log f (y_{i}) = arg max_{β} \sum_{i} n_{i} [y_{i} \log \frac{p (β)}{1 - p (β)} - (- \log (1 - p (β)))] + \log (\binom{n_{i}}{n_{i} y_{i}}) = arg max_{β} \sum_{i} n_{i} [y_{i} \log \frac{p (β)}{1 - p (β)} - (- \log (1 - p (β)))]

$\textrm{arg}\max_{\beta} \sum_i \log f(y_i)= \textrm{arg}\max_{\beta} \sum_i n_i \left[ y_i \log \frac{p(\beta)}{1-p(\beta)} - \left(-\log (1-p(\beta))\right) \right] + \log {n_i \choose n_iy_i}\\ = \textrm{arg}\max_{\beta} \sum_i n_i \left[ y_i \log \frac{p(\beta)}{1-p(\beta)} - \left(-\log (1-p(\beta))\right) \right] \\$
Dlatego jeśli pozwolimy

n_{i}^{*} = n_{i} c

$n_i^*=n_ic$ dla wszystkich

i = 1, . . ., N

$i=1,...,N$ dla jakiejś stałej

c

$c$ , to musi być również prawdą, że:

arg max_{β} \sum_{i} \log f (y_{i}) = arg max_{β} \sum_{i} n_{i}^{*} [y_{i} \log \frac{p (β)}{1 - p (β)} - (- \log (1 - p (β)))]

$\textrm{arg}\max_{\beta} \sum_i \log f(y_i) = \textrm{arg}\max_{\beta} \sum_i n^*_i \left[ y_i \log \frac{p(\beta)}{1-p(\beta)} - \left(-\log (1-p(\beta))\right) \right] \\$ Z tego pomyślałem, że Skalowanie liczby prób $n_i$ ze stałą NIE wpływa na oszacowania maksymalnego prawdopodobieństwa $\beta$ biorąc pod uwagę odsetek powodzenia $y_i$ .

Plik pomocy glm mówi:

 "For a binomial GLM prior weights are used to give the number of trials 
  when the response is the proportion of successes"

Dlatego spodziewałem się, że skalowanie wagi nie wpłynie na szacowany $\beta$ biorąc pod uwagę odsetek sukcesu jako odpowiedzi. Jednak następujące dwa kody zwracają różne wartości współczynników:

 Y <- c(1,0,0,0) ## proportion of observed success
 w <- 1:length(Y) ## weight= the number of trials
 glm(Y~1,weights=w,family=binomial)

Daje to:

 Call:  glm(formula = Y ~ 1, family = "binomial", weights = w)

 Coefficients:
 (Intercept)  
      -2.197

podczas gdy jeśli pomnożę wszystkie wagi przez 1000, szacowane współczynniki są różne:

 glm(Y~1,weights=w*1000,family=binomial)

 Call:  glm(formula = Y ~ 1, family = binomial, weights = w * 1000)

 Coefficients:
 (Intercept)  
    -3.153e+15

Widziałem wiele innych takich przykładów, nawet z umiarkowanym skalowaniem wag. Co tu się dzieje?

— Wróżka
źródło

Dla tego, co jest warte, weightsargument kończy się w dwóch miejscach wewnątrz glm.fitfunkcji (w glm.R ), czyli co działa w R: 1) binomial_dev_residsw resztkach odchylenia, za pomocą funkcji C (w rodzinie. C ) i 2) w kroku IWLS za pomocą Cdqrls(w lm.c ). Nie znam wystarczającej ilości C, aby pomóc w śledzeniu logiki

— shadowtalker

Sprawdź odpowiedzi tutaj .

— Stat.

@ssdecontrol Czytam przez glm.fit w linku, który mi dałeś, ale nie mogę znaleźć, gdzie w glm.fit jest wywoływana funkcja C „binomial_dev_resids”. Czy miałbyś coś przeciwko, gdybyś to wskazał?

— FairyOnIce

@ssdecontrol Och, przepraszam, myślę, że rozumiem. Każda „rodzina” jest listą, a jednym z elementów jest „dev.resids”. Kiedy wpisuję dwumianowy w konsoli R, widzę definicję obiektu dwumianowego i ma on linię: dev.resids <- funkcja (y, mu, wt). Call (C_binomial_dev_resids, y, mu, wt)

— FairyOnIce

Odpowiedzi:

Twój przykład powoduje jedynie błąd zaokrąglania w R. Duże ciężary nie działają dobrze w glm. Prawdą jest, że skalowanie wo praktycznie każdą mniejszą liczbę, na przykład 100, prowadzi do takich samych oszacowań, co nieskalowane w.

Jeśli chcesz bardziej niezawodnego zachowania z argumentami wag, spróbuj użyć svyglmfunkcji z surveypakietu.

Spójrz tutaj:

    > svyglm(Y~1, design=svydesign(ids=~1, weights=~w, data=data.frame(w=w*1000, Y=Y)), family=binomial)
Independent Sampling design (with replacement)
svydesign(ids = ~1, weights = ~w, data = data.frame(w = w * 1000, 
    Y = Y))

Call:  svyglm(formula = Y ~ 1, design = svydesign(ids = ~1, weights = ~w2, 
    data = data.frame(w2 = w * 1000, Y = Y)), family = binomial)

Coefficients:
(Intercept)  
     -2.197  

Degrees of Freedom: 3 Total (i.e. Null);  3 Residual
Null Deviance:      2.601 
Residual Deviance: 2.601    AIC: 2.843

— AdamO
źródło

Myślę, że sprowadza się to do początkowych wartości użytych w glm.fitstosunku do tego, family$initializeco sprawia, że metoda się różni. O ile wiem, glm.fitrozwiąż problem, tworząc rozkład QR z gdzie jest macierzą projektową, a to przekątna z pierwiastkami kwadratowymi wpisów, jak opisano tutaj . Oznacza to, że wykorzystuje metodę Newtona-Raphsona. $\sqrt{W}X$ $X$ $\sqrt{W}$

Odpowiedni $intializekod to:

if (NCOL(y) == 1) {
    if (is.factor(y)) 
        y <- y != levels(y)[1L]
    n <- rep.int(1, nobs)
    y[weights == 0] <- 0
    if (any(y < 0 | y > 1)) 
        stop("y values must be 0 <= y <= 1")
    mustart <- (weights * y + 0.5)/(weights + 1)
    m <- weights * y
    if (any(abs(m - round(m)) > 0.001)) 
        warning("non-integer #successes in a binomial glm!")
}

Oto uproszczona wersja, glm.fitktóra pokazuje mój punkt widzenia

> #####
> # setup
> y <- matrix(c(1,0,0,0), ncol = 1)
> weights <- 1:nrow(y) * 1000
> nobs <- length(y)
> family <- binomial()
> X <- matrix(rep(1, nobs), ncol = 1) # design matrix used later
> 
> # set mu start as with family$initialize
> if (NCOL(y) == 1) {
+   n <- rep.int(1, nobs)
+   y[weights == 0] <- 0
+   mustart <- (weights * y + 0.5)/(weights + 1)
+   m <- weights * y
+   if (any(abs(m - round(m)) > 0.001)) 
+     warning("non-integer #successes in a binomial glm!")
+ }
> 
> mustart # starting value
             [,1]
[1,] 0.9995004995
[2,] 0.0002498751
[3,] 0.0001666111
[4,] 0.0001249688
> (eta <- family$linkfun(mustart))
          [,1]
[1,]  7.601402
[2,] -8.294300
[3,] -8.699681
[4,] -8.987322
> 
> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] -5.098297
> (eta <- .coef * X)
          [,1]
[1,] -5.098297
[2,] -5.098297
[3,] -5.098297
[4,] -5.098297
> 
> # repeat a few times from "start loop to fit"

Możemy powtórzyć ostatnią część jeszcze dwa razy, aby zobaczyć, że metoda Newtona-Raphsona różni się:

> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] 10.47049
> (eta <- .coef * X)
         [,1]
[1,] 10.47049
[2,] 10.47049
[3,] 10.47049
[4,] 10.47049
> 
> 
> #####
> # Start loop to fit
> mu <- family$linkinv(eta)
> mu_eta <- family$mu.eta(eta)
> z <- drop(eta + (y - mu) / mu_eta)
> w <- drop(sqrt(weights * mu_eta^2 / family$variance(mu = mu)))
> 
> # code is simpler here as (X^T W X) is a scalar
> X_w <- X * w
> (.coef <- drop(crossprod(X_w)^-1 * ((w * z) %*% X_w)))
[1] -31723.76
> (eta <- .coef * X)
          [,1]
[1,] -31723.76
[2,] -31723.76
[3,] -31723.76
[4,] -31723.76

Nie dzieje się tak, jeśli zaczniesz od weights <- 1:nrow(y)lub powiesz weights <- 1:nrow(y) * 100.

Zauważ, że można uniknąć rozbieżności, ustawiając mustartargument. Np

> glm(Y ~ 1,weights = w * 1000, family = binomial, mustart = rep(0.5, 4))

Call:  glm(formula = Y ~ 1, family = binomial, weights = w * 1000, mustart = rep(0.5, 
    4))

Coefficients:
(Intercept)  
     -2.197  

Degrees of Freedom: 3 Total (i.e. Null);  3 Residual
Null Deviance:      6502 
Residual Deviance: 6502     AIC: 6504

— Benjamin Christoffersen
źródło

Myślę, że wagi mają wpływ na więcej niż argumenty do inicjalizacji. Dzięki regresji logistycznej Newton Raphson ocenia maksymalne prawdopodobieństwo, które istnieje i jest unikalne, gdy dane nie są rozdzielane. Podanie optymalizatora różnych wartości początkowych nie doprowadzi do uzyskania różnych wartości, ale osiągnięcie go może potrwać dłużej.

— AdamO,

„Podanie różnych wartości początkowych do optymalizatora nie doprowadzi do różnych wartości ...” . No cóż, metoda Newtona nie różni się i znajduje unikalne maksimum w ostatnim przykładzie, w którym ustawiam wartości początkowe (patrz przykład, w którym podam mustart argument). Wydaje się, że jest to kwestia związana ze złym początkowym oszacowaniem .

— Benjamin Christoffersen