R: função glm com família = especificação "binomial" e "peso"

14

Estou muito confuso com a forma como o peso funciona em glm com family = "binomial". No meu entendimento, a probabilidade do glm com family = "binomial" é especificada da seguinte forma: onde é a "proporção de sucesso observado" e é o número conhecido de tentativas.

f(y)=(nny)pny(1-p)n(1-y)=exp(n[yregistrop1-p-(-registro(1-p))]+registro(nny))
yn

No meu entendimento, a probabilidade de sucesso p é parametrizada com alguns coeficientes lineares β como p=p(β) e função glm com family = "binomial" procure por:

argmaxβEuregistrof(yEu).
Em seguida, esse problema de otimização pode ser simplificado como:

argmaxβEuregistrof(yEu)=argmaxβEunEu[yEuregistrop(β)1-p(β)-(-registro(1-p(β)))]+registro(nEunEuyEu)=argmaxβEunEu[yEuregistrop(β)1-p(β)-(-registro(1-p(β)))]

Portanto, se deixarmos nEu=nEuc para todos Eu=1,...,N para alguma constante c , também deve ser verdade que:
argmaxβEuregistrof(yEu)=argmaxβEunEu[yEuregistrop(β)1-p(β)-(-registro(1-p(β)))]
A partir disso, pensei que o escalonamento do número de tentativas nEucom uma constante NÃO afeta as estimativas de probabilidade máxima de dada a proporção de sucessoβyEu .

O arquivo de ajuda da glm diz:

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

Portanto, eu esperava que a escala de peso não afetasse o valor estimado dada a proporção de sucesso como resposta. No entanto, os dois códigos a seguir retornam valores diferentes de coeficiente:β

 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)

Isso produz:

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

 Coefficients:
 (Intercept)  
      -2.197     

enquanto que se eu multiplicar todos os pesos por 1000, os coeficientes estimados são diferentes:

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

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

 Coefficients:
 (Intercept)  
    -3.153e+15  

Eu vi muitos outros exemplos como esse, mesmo com uma escala moderada de pesos. O que está acontecendo aqui?

FairyOnIce
fonte
3
Pelo que vale, o weightsargumento termina em dois lugares dentro da glm.fitfunção (em glm.R ), que é o que funciona em R: 1) nos resíduos de desvio, por meio da função C binomial_dev_resids(em family.c ) e 2) no IWLS passo a passo de Cdqrls(em lm.c ). Eu não sei o suficiente C a ser de mais ajuda no rastreamento da lógica
shadowtalker
3
Verifique as respostas aqui .
Stat
@ssdecontrol Estou lendo o glm.fit no link que você me forneceu, mas não consigo encontrar onde a função C "binomial_dev_resids" é chamada no glm.fit. Você se importaria se apontasse?
FairyOnIce
@ssdecontrol Oh, desculpe, acho que entendi. Cada "família" é uma lista e um dos elementos é "dev.resids". Quando digito binomial no console do R, vejo a definição do objeto binomial e ele tem uma linha: dev.resids <- function (y, mu, wt) .Call (C_binomial_dev_resids, y, mu, wt)
FairyOnIce

Respostas:

4

Seu exemplo está apenas causando erro de arredondamento em R. Pesos grandes não têm bom desempenho glm. É verdade que o dimensionamento wde praticamente qualquer número menor, como 100, leva às mesmas estimativas que as não dimensionadas w.

Se você deseja um comportamento mais confiável com os argumentos de ponderação, tente usar a svyglmfunção do surveypacote.

Veja aqui:

    > 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
fonte
1

glm.fitfamily$initializeglm.fitWXXW

O $intializecódigo relevante é:

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

Aqui está uma versão simplificada do glm.fitque mostra o meu ponto

> #####
> # 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 a última parte mais duas vezes para ver que o método de 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

Isso não acontece se você começar weights <- 1:nrow(y)ou dizer weights <- 1:nrow(y) * 100.

Observe que você pode evitar divergências configurando o mustartargumento. Eg fazer

> 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
fonte
Eu acho que pesos afeta mais do que argumentos para inicializar. Com a regressão logística, Newton Raphson estima a probabilidade máxima que existe e é única quando os dados não são separados. O fornecimento de valores iniciais diferentes ao otimizador não chegará a valores diferentes, mas talvez demore mais tempo para chegar lá.
precisa saber é o seguinte
"O fornecimento de diferentes valores iniciais ao otimizador não chegará a valores diferentes ..." . Bem, o método Newton não diverge e encontra o máximo único no último exemplo em que eu defino os valores iniciais (veja o exemplo em que forneço o mustart argumento). Parece um assunto relacionado à estimativa inicial ruim .
Benjamin Christoffersen