Comment obtenir la distribution de Poisson "lambda" à partir des coefficients R glm()

Question

Mon R-script produit des coeffs glm() ci-dessous. Qu'est-ce que le lambda de Poisson, alors? Ca devrait être ~ 3.0 puisque c'est ce que j'ai utilisé pour créer la distribution.

Call: 
glm(formula = h_counts ~ ., family = poisson(link = log), data = pois_ideal_data) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-22.726 -12.726 -8.624 6.405 18.515 

Coefficients: 
      Estimate Std. Error z value Pr(>|z|)  
(Intercept) 8.222532 0.015100 544.53 <2e-16 *** 
h_mids  -0.363560 0.004393 -82.75 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1) 

Null deviance: 11451.0 on 10 degrees of freedom 
Residual deviance: 1975.5 on 9 degrees of freedom 
AIC: 2059 

Number of Fisher Scoring iterations: 5 


random_pois = rpois(10000,3) 
h=hist(random_pois, breaks = 10) 
mean(random_pois) #verifying that the mean is close to 3. 
h_mids = h$mids 
h_counts = h$counts 
pois_ideal_data <- data.frame(h_mids, h_counts) 
pois_ideal_model <- glm(h_counts ~ ., pois_ideal_data, family=poisson(link=log)) 
summary_ideal=summary(pois_ideal_model) 
summary_ideal 

Answer 1

Que faites-vous ici ??? Vous avez utilisé un glm pour adapter une distribution ???

Eh bien, il est impossible de le faire, mais il se fait par ceci:

set.seed(0) 
x <- rpois(10000,3) 
fit <- glm(x ~ 1, family = poisson()) 

à savoir, nous ajuster les données avec un modèle de régression ordonnée à l'origine uniquement.

fit$fitted[1] 
# 3.005 

Ceci est la même chose que:

mean(x) 
# 3.005 

Answer 2

On dirait que vous essayez de faire un Poisson ajusté aux données agrégées ou mis en cellule; ce n'est pas ce que fait glm. J'ai jeté un rapide coup d'œil aux méthodes d'ajustement des distributions aux données en conserve, mais je n'ai pas pu en trouver un; il semble que les versions antérieures du package bda l'aient proposé, mais pas maintenant.

À la racine, ce que vous devez faire est de mettre en place une fonction de log-vraisemblance négative qui calcule (# counts)*prob(count|lambda) et le minimise en utilisant optim(); la solution donnée ci-dessous à l'aide du paquet bbmle est un peu plus d'avance complexe, mais vous donne des avantages supplémentaires comme calculer facilement les intervalles de confiance etc ..

Définissez des données:

set.seed(101) 
random_pois <- rpois(10000,3) 
tt <- table(random_pois) 
dd <- data.frame(counts=unname(c(tt)), 
       val=as.numeric(names(tt))) 

J'utilise ici table plutôt que hist parce histogrammes sur des données discrètes sont difficiles (ayant des points de coupe entiers rend souvent les choses confuses parce que vous devez faire attention à droite par rapport à la fermeture de gauche)

Régler la fonction de densité pour les données regroupées par casiers Poisson (pour travailler avec bbmle ' s forme interface ula, le premier argument doit être appelé x, et il doit avoir un argument log).

dpoisbin <- function(x,val,lambda,log=FALSE) { 
    probs <- dpois(val,lambda,log=TRUE) 
    r <- sum(x*probs) 
    if (log) r else exp(r) 
} 

Fit lambda (log lien permet d'éviter les problèmes numériques/avertissements de valeurs lambda négatives):

library(bbmle) 
m1 <- mle2(counts~dpoisbin(val,exp(loglambda)), 
     data=dd, 
     start=list(loglambda=0)) 
all.equal(unname(exp(coef(m1))),mean(random_pois),tol=1e-6) ## TRUE 
exp(confint(m1)) 
## 2.5 % 97.5 % 
## 2.972047 3.040009