J'essaye de mettre en application la modélisation de poisson et j'ai besoin de votre expertise.
J'ai un dataframe contenant le nombre de pathologies. Ce data frame a été obtenu par un dénombrement sur un fichier d'individu.
voici la méthode de dénombrement:
Code : Tout sélectionner
nbT <- ddply(dtPrt, .variables = .(annee,sexe,classeage,statut,origine,activite), nrow)
Ici, dtPrt correspond au jeu de données initial au format csv que j'ai préalablement ouvert par un read.csv et nrow va correspondre à la variable de dénombrement (V1) dans mon dataframe nbT.
1/ Identification des écarts par rapport à une distribution poissonnienne
L'approche graphique montre bien la forme d'un L. Cependant je n'ai pas de zéro dans cette variable V1:
Code : Tout sélectionner
summary(nbT$V1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 1.000 2.000 4.228 5.000 39.000
moy <- paste("Moyenne = ",mean(V1))
var <- paste("Variance = ",var(V1))
disp <- paste("Dispersion= ",var(V1)/mean(V1))
list=c(moy,var,disp)
print(cbind(list))
[1,] "Moyenne = 4.22842377260982"
[2,] "Variance = 28.4214253839227"
[3,] "Dispersion= 6.72151773623691"
A priori on a une surdispersion
2/ Définition du modèle en fonction du rapport moyenne/variance
On choisit ici deux modèles : quasi poisson et négatif binomial
Code : Tout sélectionner
library(MASS)
mod1 <- glm(V1~classeage+sexe+statut+annee+origine+activite,family=quasipoisson(link="log"), data=nbT)
mod2 <- glm.nb(mod1)
Code : Tout sélectionner
mod1$deviance/mod1$df.residual # modèle 1 (quasi poisson)
[1] 1.652043
mod2$deviance/mod2$df.residual # modèle 2 (binomial negative)
[1] 0.8528109
3/ Qualité des modèles
Test des paramètres du modèle basé sur le khi-deux :
Comparaison de la déviance résiduelle à une distribution du khi-deux
Code : Tout sélectionner
> pchisq(mod1$deviance,mod1$df.residual, lower=F)
[1] 1.710149e-64
Calcul de la qualité de l'ajustement des modèles avec les résidus de pearson (Ref : Lancelot, Lesnoff, Sélection de modèle avec AIC)
Code : Tout sélectionner
> x2<-sum(residuals(mod1,type="pearson")^2)
> dd2<-df.residual(mod1)
> 1-pchisq(x2,dd2)
[1] 0
> x2<-sum(residuals(mod2,type="pearson")^2)
> dd2<-df.residual(mod2)
> 1-pchisq(x2,dd2)
[1] 0.9937055
Peut-on conclure que le modèle 1 ajuste mal les données et que par contre le modèle 2 les ajuste bien?
La valeur du 1-pchisq(x2,dd2) correspond bien à un p-value ?
4/ comparaison des deux modèles :
Peut-on utiliser l’analyse de la déviance par :
Code : Tout sélectionner
> 1-pchisq(deviance(mod1)-deviance(mod2),df.residual(mod1)-df.residual(mod2))
[1] 0
Il me semble également qu’on peut comparer avec une anova(mod1,mod2, test= ‘Chisq’)
Cependant, il faut que les deux modèles ajustent bien les données.
Je vous remercie par avance.