J'essaie d'adapter les modèles linéaires généralisés à certains ensembles de données de comptage qui pourraient ou non être sur-dispersés. Les deux distributions canoniques qui s'appliquent ici sont le binôme de Poisson et négatif (Negbin), avec EV et variance
qui peuvent être montés dans R en utilisant glm(..,family=poisson)
et glm.nb(...)
, respectivement. Il y a aussi la quasipoisson
famille, qui selon moi est un Poisson ajusté avec le même EV et la même variance
,
c'est-à-dire se situant quelque part entre Poisson et Negbin. Le principal problème avec la famille des quasipoisson est qu'il n'y a pas de probabilité correspondante pour elle, et donc de nombreux tests statistiques et mesures d'ajustement extrêmement utiles (AIC, LR, etc.) ne sont pas disponibles.
Si vous comparez les variances QP et Negbin, vous remarquerez peut-être que vous pouvez les égaliser en mettant . En poursuivant dans cette logique, vous pourriez essayer d'exprimer la distribution quasipoisson comme un cas particulier du Negbin:
,
c'est-à-dire un Negbin avec dépendant linéairement de . J'ai essayé de vérifier cette idée en générant une séquence aléatoire de nombres selon la formule ci-dessus et en l'adaptant avec :μglm
#fix parameters
phi = 3
a = 1/50
b = 3
x = 1:100
#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison
mu = exp(a*x+b)
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator
#fit a generalized linear model y = f(x)
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial
> glmQP
Call: glm(formula = y ~ x, family = quasipoisson)
Coefficients:
(Intercept) x
3.11257 0.01854
(Dispersion parameter for quasipoisson family taken to be 3.613573)
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 2097
Residual Deviance: 356.8 AIC: NA
> glmNB
Call: glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)
Coefficients:
(Intercept) x
3.10182 0.01873
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 578.1
Residual Deviance: 107.8 AIC: 824.7
Les deux ajustements reproduisent les paramètres et le quasipoisson donne une estimation «raisonnable» pour . On peut maintenant aussi définir une valeur AIC pour le quasipoisson:
df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values
#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329
(J'ai dû copier manuellement la valeur ajustée depuis , car je ne pouvais pas la trouver dans l' objet)summary(glmQP)
glmQP
Puisque , cela indiquerait que le quasipoisson est, sans surprise, le meilleur ajustement; donc au moins fait ce qu'il devrait faire, et donc ce pourrait être une définition raisonnable pour l'AIC (et par extension, la vraisemblance) d'un quasipoisson. Les grandes questions qui me restent sont alors A I C Q P
- Cette idée a-t-elle un sens? Ma vérification est-elle basée sur un raisonnement circulaire?
- La principale question pour quiconque «invente» quelque chose qui semble manquer dans un sujet bien établi: si cette idée a du sens, pourquoi n'est-elle pas déjà mise en œuvre
glm
?
Modifier: figure ajoutée
la source
Réponses:
Le quasi-Poisson n'est pas un modèle à maximum de vraisemblance maximale (ML) mais un modèle quasi-ML. Vous utilisez simplement la fonction d'estimation (ou la fonction de score) du modèle de Poisson pour estimer les coefficients, puis utilisez une certaine fonction de variance pour obtenir des erreurs standard appropriées (ou plutôt une matrice de covariance complète) pour effectuer l'inférence. Par conséquent,
glm()
ne fournit pas etlogLik()
ouAIC()
ici etc.size
S'il n'y a pas régresseurs (juste une interception) la paramétrisation NB1 et la paramétrisation NB2 employé par
MASS
de »glm.nb()
la coïncidence. Avec les régresseurs, ils diffèrent. Dans la littérature statistique, le paramétrage NB2 est plus fréquemment utilisé mais certains progiciels proposent également la version NB1. Par exemple, dans R, vous pouvez utiliser legamlss
package pour le fairegamlss(y ~ x, family = NBII)
. Notez que quelque peu confusegamlss
utiliseNBI
pour le paramétrage NB2 etNBII
pour NB1. (Mais le jargon et la terminologie ne sont pas unifiés dans toutes les communautés.)Vous pourriez alors vous demander, bien sûr, pourquoi utiliser du quasi-Poisson si NB1 est disponible? Il y a encore une différence subtile: le premier utilise du quasi-ML et obtient l'estimation à partir de la dispersion à partir des résidus de déviance au carré (ou Pearson). Ce dernier utilise le ML complet. Dans la pratique, la différence n'est souvent pas grande, mais les motivations à utiliser l'un ou l'autre modèle sont légèrement différentes.
la source
gamlss
maintenant et il semble que c'est exactement ce dont j'avais besoin. Pourriez-vous nous expliquer les motivations de l'utilisation de la quasi-vraisemblance par rapport au ML complet?vignette("countreg", package = "pscl")
.