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?

r logistic generalized-linear-model binomial weighted-regression 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