R: fonction glm avec famille = «binôme» et «poids»

14

Je suis très confus avec la façon dont le poids fonctionne en glm avec family = "binomial". À ma connaissance, la probabilité du glm avec family = "binomial" est spécifiée comme suit: où y est la "proportion de succès observée" et n est le nombre connu d'essais.yn

F(y)=(nny)pny(1-p)n(1-y)=exp(n[yJournalp1-p-(-Journal(1-p))]+Journal(nny))
yn

À ma connaissance, la probabilité de succès p est paramétrée avec certains coefficients linéaires β comme p=p(β) et la fonction glm avec family = "binomial" recherche:

argmaxβjeJournalF(yje).
Ensuite, ce problème d'optimisation peut être simplifié comme suit:

argmaxβjeJournalF(yje)=argmaxβjenje[yjeJournalp(β)1-p(β)-(-Journal(1-p(β)))]+Journal(njenjeyje)=argmaxβjenje[yjeJournalp(β)1-p(β)-(-Journal(1-p(β)))]

Par conséquent, si nous laissons nje=njec pour tout je=1,...,N pour une constante c , alors il doit également être vrai que:
argmaxβjeJournalF(yje)=argmaxβjenje[yjeJournalp(β)1-p(β)-(-Journal(1-p(β)))]
De là, je pensais que la mise à l'échelle du nombre d'essais njeavec une constante n'affecte PAS les estimations de vraisemblance maximale de étant donné la proportion de réussiteβyje .

Le fichier d'aide de glm indique:

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

Par conséquent, je m'attendais à ce que la mise à l'échelle du poids n'affecte pas la estimée étant donné la proportion de succès en tant que réponse. Cependant, les deux codes suivants renvoient des valeurs de coefficient différentes:β

 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)

Cela donne:

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

 Coefficients:
 (Intercept)  
      -2.197     

alors que si je multiplie tous les poids par 1000, les coefficients estimés sont différents:

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

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

 Coefficients:
 (Intercept)  
    -3.153e+15  

J'ai vu de nombreux autres exemples comme celui-ci, même avec une mise à l'échelle modérée des poids. Qu'est-ce qui se passe ici?

FairyOnIce
la source
3
Pour ce qu'il vaut, l' weightsargument se retrouve à deux endroits à l'intérieur de la glm.fitfonction (dans glm.R ), ce qui fait le travail dans R: 1) dans les résidus de déviance, par le biais de la fonction C binomial_dev_resids(dans family.c ) et 2) dans l'étape IWLS au moyen de Cdqrls(en lm.c ). Je ne connais pas assez de C pour être plus utile dans le traçage de la logique
shadowtalker
3
Vérifiez les réponses ici .
Stat
@ssdecontrol Je lis par glm.fit dans le lien que vous m'avez donné mais je ne trouve pas où la fonction C "binomial_dev_resids" est appelée dans glm.fit. Cela vous dérangerait-il si vous le montriez?
FairyOnIce
@ssdecontrol Oh, désolé, je pense que je comprends. Chaque "famille" est une liste et l'un des éléments est "dev.resids". Lorsque je tape binomial dans la console R, je vois la définition de l'objet binomial et il a une ligne: dev.resids <- fonction (y, mu, wt) .Call (C_binomial_dev_resids, y, mu, wt)
FairyOnIce

Réponses:

4

Votre exemple provoque simplement une erreur d'arrondi dans R. Les grands poids ne fonctionnent pas bien dans glm. Il est vrai que la mise wà l' échelle de pratiquement n'importe quel nombre plus petit, comme 100, conduit aux mêmes estimations que celles non mises à l'échelle w.

Si vous voulez un comportement plus fiable avec les arguments de poids, essayez d'utiliser la svyglmfonction du surveypackage.

Vois ici:

    > 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
la source
1

glm.fitfamily$initializeglm.fitWXXW

Le $intializecode pertinent est:

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

Voici une version simplifiée glm.fitqui montre mon point

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

Nous pouvons répéter la dernière partie deux fois de plus pour voir que la méthode 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

Cela ne se produit pas si vous commencez weights <- 1:nrow(y)ou dites weights <- 1:nrow(y) * 100.

Notez que vous pouvez éviter la divergence en définissant l' mustartargument. Par exemple faire

> 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
la source
Je pense que les poids affectent plus que les arguments à initialiser. Avec la régression logistique, Newton Raphson estime la probabilité maximale qui existe et qui est unique lorsque les données ne sont pas séparées. La fourniture de différentes valeurs de départ à l'optimiseur n'arrivera pas à des valeurs différentes, mais il faudra peut-être plus de temps pour y arriver.
AdamO
"La fourniture de différentes valeurs de départ à l'optimiseur n'arrivera pas à des valeurs différentes ..." . Eh bien, la méthode Newton ne diverge pas et trouve le maximum unique dans le dernier exemple où j'ai défini les valeurs initiales (voir l'exemple où je fournis l' mustart argument). Cela semble être une question liée à une mauvaise estimation initiale .
Benjamin Christoffersen