Les poids et le décalage peuvent-ils conduire à des résultats similaires dans la régression du poisson?

8

Dans le "Guide du praticien des modèles linéaires généralisés" au paragraphe 1.83, il est indiqué que:

"Dans le cas particulier d'un GLM multiplicatif de Poisson, il peut être démontré que la modélisation des comptes de sinistres avec un terme de décalage égal au log de l'exposition a produit des résultats identiques à la modélisation des fréquences de revendications avec des poids antérieurs définis pour être égaux à l'exposition de chaque observation. "

Je ne suis pas en mesure de trouver d'autres références de ces résultats, j'ai donc entrepris des tests empiriques dans lesquels je n'ai pas pu trouver la preuve que l'énoncé est correct. Quelqu'un peut-il nous expliquer pourquoi ces résultats peuvent être bons / mauvais?

Pour info, j'ai utilisé le code R suivant pour tester l'hypothèse, dans laquelle je n'ai pas pu obtenir des résultats similaires pour les deux cas mentionnés:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y'=y, 'X'=X, 'offset' = offset)
formula = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula, family = "poisson", df, weights = offset) 
#Third model using poisson model on y/offset as reference
dfNew = df
dfNew$y = dfNew$y/offset
fit3 = glm(formula, family = "poisson", dfNew)

#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, fit3$coefficients, c(intercept,coefs))

Les estimations de coefficient résultant de l'exécution de ce code sont données ci-dessous:

 >  
    (Intercept)       X.1       X.2       X.3        X.4       X.5       X.6
[1,]    1.998277 0.2923091 0.4586666 0.1802960 0.11688860 0.7997154 0.4786655
[2,]    1.588620 0.2708272 0.4540180 0.1901753 0.07284985 0.7928951 0.5100480
[3,]    1.983903 0.2942196 0.4593369 0.1782187 0.11846876 0.8018315 0.4807802
[4,]    2.000000 0.2909240 0.4576965 0.1807591 0.11658183 0.8005451 0.4780123
              X.7       X.8       X.9      X.10
[1,]  0.005772078 0.9154808 0.9078758 0.3512824
[2,] -0.003705015 0.9117014 0.9063845 0.4155601
[3,]  0.007595660 0.9181014 0.9076908 0.3505173
[4,]  0.005881960 0.9150350 0.9084375 0.3511749
> 

et nous pouvons observer que les coefficients ne sont pas identiques.

BDP1
la source
1
Vous ne devriez pas vraiment inclure rm(list=ls() )dans le code R que vous postez ici! Cela pourrait surprendre quelqu'un qui l'exécute, les mettant en colère contre vous. Je l'ai retiré. J'ai également modifié pour inclure les résultats de l'exécution du code. Si vous l'aviez fait à l'origine, vous avez peut-être obtenu une réponse plus rapide, car peu de lecteurs exécuteront le code eux-mêmes.
kjetil b halvorsen
1
@kjetilbhalvorsen Merci pour la rétroaction, fera à l'avenir.
BDP1
@ BDP1 En l'état, votre code R ne teste pas la revendication à laquelle vous faites référence. Je vous suggère d'ajouter le terme de poids pour le fit3 et d'ajouter un addendum directement à la question.
aivanov

Réponses:

7

(avec votre code R, vous pouvez remplacer "poisson" par "quasipoisson" pour éviter tous les avertissements qui sont générés. Rien d'autre d'import ne changera. Voir (*) ci-dessous). Votre référence utilise le terme "glm multiplicatif" qui, je pense, signifie simplement un glm avec un lien de log, car un lien de log peut être considéré comme un modèle multiplicatif. Votre propre exemple montre que l'affirmation est fausse, du moins telle que nous l'avons interprétée (puisque les paramètres estimés ne sont pas égaux). Vous pouvez écrire aux auteurs et leur demander ce qu'ils signifient. Ci-dessous, je vais expliquer pourquoi cette affirmation est fausse.

Laisser λje être le paramètre poisson et ωjeles poids. Laisserηje être le prédicteur linéaire sans décalage, puis ηje+Journal(ωje)être le prédicteur linéaire avec le décalage. La fonction de probabilité du poisson est

F(yje)=e-λjeλjeyje/yje!
Ensuite, la fonction log de vraisemblance pour le modèle avec décalage devient
=-jeωjeeηje+jeyjeηje+jeyjeJournalωje-jeJournalyje!
tandis que la fonction log de vraisemblance pour le modèle avec des poids devient
w=-jeωjeeηje+jeyjeωjeηje-jeωjeJournalyje!
et ce n'est clairement pas la même chose. Donc ce que ces auteurs voulaient dire n'est pas clair pour moi.

(*) Remarque à l'aide de la glmfonction de R :

Les «poids» non «NULS» peuvent être utilisés pour indiquer que différentes observations ont des dispersions différentes (les valeurs en «poids» étant inversement proportionnelles aux dispersions); ou de manière équivalente, lorsque les éléments de "poids" sont des entiers positifs w_i, que chaque réponse y_i est la moyenne des observations de poids unitaire w_i. Pour un GLM binomial, des poids antérieurs sont utilisés pour donner le nombre d'essais lorsque la réponse est la proportion de succès: ils seraient rarement utilisés pour un GLM de Poisson.

L'examen de la signification des arguments de poids explique cela, il donne peu de sens à la fonction de la famille poisson, qui suppose un paramètre d'échelle constant ϕ=1 tandis que les arguments de poids modifient ϕ. Cela donne plus de sens à la fonction de la famille des quasipossons. Voir la réponse à l' entrée "poids" dans les fonctions glm et lm dans R La réponse donnée ici permet également de comprendre pourquoi la probabilité dans le cas pondéré prend la forme donnée ci-dessus.

La réponse donnée ici pourrait être pertinente: comment une régression du taux de Poisson est-elle égale à une régression de Poisson avec le terme de décalage correspondant? et est très intéressant.

kjetil b halvorsen
la source
Merci d'avoir répondu. Aller à travers les probabilités de contributions clarifie beaucoup. En recherchant un peu plus en réponse à votre réponse, j'ai trouvé la question suivante: stats.stackexchange.com/questions/66791/… Dans laquelle il est démontré qu'en divisant la réponse d'origine par le décalage et en définissant la variable de décalage comme poids, les mêmes résultats peuvent être obtenus que le modèle de base, où le log (offset) entre dans le prédicteur linéaire avec un coefficient constant de 1.
BDP1
J'ai également tenté de déduire que les contributions probables de ce nouveau modèle sont égales aux contributions du modèle nouvellement mentionné, mais je n'ai pas pu le faire, étant donné la (yje/wje)!terme qui apparaît pour ce dernier.
BDP1
Au moins dans l'expérience du script R fourni, cette affirmation semble être vraie. cela peut être facilement vu en ajoutant un argument weights = offset dans la ligne fit3 = glm (...)
BDP1
Je ne comprends pas ce que vous dites ici, si vous pensez que votre expérience corrobore, pourquoi ne pas modifier ce changement dans l'expérience du PO et expliquer pourquoi vous pensez que cela corrobore la revendication. J'ai essayé d'expliquer pourquoi ce n'est pas vrai, quel est le problème avec mon argument?
kjetil b halvorsen
1
La citation du guide du praticien est en fait correcte, elle n'a tout simplement pas été correctement mise en œuvre par OP. Cela a été souligné par Alan Chalk dans une autre réponse ici, par moi stats.stackexchange.com/questions/430001 et aussi par Cokes stats.stackexchange.com/questions/264071 (bien que la conclusion correcte soit enterrée dans la dernière ligne de cette réponse).
Gordon Smyth
4

Désolé de ne pas simplement ajouter un commentaire ci-dessus, mais je n'ai pas assez de représentants.

L'allégation originale - mais un peu modifiée - est en fait correcte.

Les deux modèles suivants donnent exactement la même réponse dans R en utilisant un poisson glm avec log-link:

  • modèle y, utiliser un décalage étant log (offset)
  • modélisez y / offset, utilisez des poids égaux à offset

ajuster légèrement votre code d'origine montre des réponses identiques:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y' = y,
                'y_over_offset' = y/offset,
                'X' = X,
                'offset' = offset)

formula_offset = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))
formula_weights = paste("y_over_offset ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula_offset, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula_weights, family = "poisson", df, weights = offset) 


#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, c(intercept,coefs))

J'espère que cela devrait donner des réponses identiques.

Il est possible de montrer que les deux modèles sont statistiquement équivalents (il y a un article CAS quelque part qui le montre - je posterai un lien si j'ai le temps).

Soit dit en passant, si vous effectuez une régression pénalisée, la façon dont différents packages tels que glmnet et H2o mesurent la déviance pour les deux façons différentes de définir un modèle peut conduire à des résultats différents.

Alan Chalk
la source
1
Question rapide, vous pouvez voir qu'avec l'ajustement2 vous n'avez pas d'AIC et les tracés entre les deux ajustements sont légèrement différents (ajustement1); plot (fit2) - savez-vous aussi pourquoi cela se produit?]
Charl Francois Marais