R: función glm con familia = "binomial" y especificación "peso".

Question

Estoy muy confundido con cómo funciona el peso en glm con family="binomial" . A mi entender, la probabilidad de que el glm con family = "binomial" se especifica del siguiente modo: $$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)$$ donde $y$ es la "proporción de éxito observada" y $n$ es el número conocido de ensayos.

A mi entender, la probabilidad de éxito $p$ se parametriza con unos coeficientes lineales $\beta$ como $p=p(\beta)$ y glm con family = "binomial" buscar: $$\textrm{arg}\max_{\beta} \sum_i \log f(y_i).$$ Entonces este problema de optimización se puede simplificar como:

$$\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] \\ $$
Por lo tanto, si dejamos que $n_i^*=n_ic$ para todos $i=1,...,N$ para alguna constante $c$ entonces también debe ser cierto que: $$\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] \\ $$ A partir de esto, pensé que Escalado del número de ensayos $n_i$ con una constante NO afecta a las estimaciones de máxima verosimilitud de $\beta$ dada la proporción de éxito $y_i$ .

El archivo de ayuda de glm dice:

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

Por lo tanto, esperaba que el escalado del peso no afectara a la estimación del $\beta$ dada la proporción de éxito como respuesta. Sin embargo, los dos códigos siguientes devuelven valores de coeficiente diferentes:

 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)

Esto produce:

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

 Coefficients:
 (Intercept)  
      -2.197     

mientras que si multiplico todas las ponderaciones por 1000, los coeficientes estimados son diferentes:

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

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

 Coefficients:
 (Intercept)  
    -3.153e+15  

He visto muchos otros ejemplos como este, incluso con un escalado moderado en los pesos. ¿Qué ocurre aquí?

Answer 1

Su ejemplo no hace más que provocar un error de redondeo en R. Los pesos grandes no funcionan bien en glm . Es cierto que escalar w por prácticamente cualquier número más pequeño, como 100, conduce a las mismas estimaciones que la no escalada w .

Si desea un comportamiento más fiable con los argumentos de los pesos, pruebe a utilizar la función svyglm de la función survey paquete.

Consulte aquí:

    > 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

Answer 2

Creo que se reduce a los valores iniciales que se utiliza en glm.fit del family$initialize que hace que el método divergere. Hasta donde yo sé, glm.fit resolver el problema formando una descomposición QR de $\sqrt{W}X$ donde $X$ es la matriz de diseño y $\sqrt{W}$ es una diagonal con raíces cuadradas de las entradas como se describe aquí . Es decir, utiliza un método Newton-Raphson.

Los $intialize código es:

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

He aquí una versión simplificada de glm.fit lo que demuestra mi punto

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

Podemos repetir la última parte dos veces más para ver que el método Newton-Raphson diverge:

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

Esto no ocurre si se empieza con weights <- 1:nrow(y) o decir weights <- 1:nrow(y) * 100 .

Tenga en cuenta que puede evitar la divergencia estableciendo el parámetro mustart argumento. Por ejemplo

> 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