glm poisson ou binomiale, pb de contradiction Chi 2/donnees

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

glm poisson ou binomiale, pb de contradiction Chi 2/donnees

Messagepar myriam desanlis » 24 Aoû 2011, 08:05

Bonjour,

je suis de retour avec une nouvelle question. J'ai un modèle qui représente un taux de contamination selon des facteurs de conduite de culture (uniquement des variables qualitatives). Etant donné les valeurs de taux de contamination (varient entre 0 et 1 avec un maximum de valeurs entre 0.5 et 0.6), je m'étais dirigée sur une famille binomiale plutot que poisson.

Code : Tout sélectionner

> mod1=glm(Incidence~Az+Var+Con, family="binomial",weights=Nfeuille)
> summary(mod1)

Call:
glm(formula = Incidence ~ Az + Var + Con, family = "binomial",
    weights = Nfeuille)

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-4.4854  -0.7473  -0.0263   0.6148   5.2458 

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.057587   0.039500   1.458  0.14487   
Az150        0.081697   0.030156   2.709  0.00675 **
VarKe        0.436208   0.048131   9.063  < 2e-16 ***
VarLG        0.008759   0.048322   0.181  0.85617   
VarNK       -0.124632   0.047574  -2.620  0.00880 **
VarSI       -0.133770   0.045557  -2.936  0.00332 **
ConS         0.088388   0.030164   2.930  0.00339 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1193.96  on 599  degrees of freedom
Residual deviance:  991.16  on 593  degrees of freedom
AIC: 3269.2

Number of Fisher Scoring iterations: 4


J'ai fait le test du chi 2 de pearson (j'ai suivi les indications d'un lien de Renaud Lancelot sur la sélection de modèles avec l'AIC) et là je trouve que la famille binomiale n'est pas adaptée !

Code : Tout sélectionner

> x1=sum(residuals(mod1,type="pearson")^2)
> dd1=df.residual(mod1)
> 1-pchisq(x1,dd1)
[1] 0


donc nouvel essai en mettant une loi de poisson

Code : Tout sélectionner

> mod2=glm(Incidence~Az+Var+Con+offset(log(Nfeuille)), family="poisson")
There were 50 or more warnings (use warnings() to see the first 50)
> summary(mod2)

Call:
glm(formula = Incidence ~ Az + Var + Con + offset(log(Nfeuille)),
    family = "poisson")

Deviance Residuals:
     Min        1Q    Median        3Q       Max 
-0.61149  -0.10678  -0.00221   0.08749   0.69975 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -4.05264    0.14817 -27.352   <2e-16 ***
Az150        0.01767    0.11068   0.160    0.873   
VarKe        0.22195    0.16902   1.313    0.189   
VarLG        0.14065    0.17615   0.798    0.425   
VarNK        0.02093    0.17912   0.117    0.907   
VarSI       -0.15508    0.17921  -0.865    0.387   
ConS        -0.01659    0.11068  -0.150    0.881   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 21.564  on 599  degrees of freedom
Residual deviance: 15.919  on 593  degrees of freedom
AIC: Inf

Number of Fisher Scoring iterations: 4


J'ai des warnings qui apparaissent (apparemment ca correspond aux taux de contaminations supérieurs à 0.5) mais le test de pearson a l'air ok

Code : Tout sélectionner

> x2=sum(residuals(mod2,type="pearson")^2)
> dd2=df.residual(mod2)
> 1-pchisq(x2,dd2)
[1] 1


d'où ma question : est-ce que j'ai raté qqch ?

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 24 Aoû 2011, 08:52

Avec votre modèle binomial, vous avez un paramètre de surdispersion de

Code : Tout sélectionner

> 991.16 / 593
[1] 1.671433


ce qui n'est pas monstrueux. J'essaierais d'abord de prendre en compte cette surdispersion avec un modèle béta-binomial, par exemple. Voir la fonction betabin du package aod. Votre modèle devrait s'écrire

Code : Tout sélectionner

library(aod)
mod1 <- betabin(cbind(Incidence * Nfeuille, Nfeuille - Incidence * Nfeuille) ~ Az + Var + Con, random = ~ 1)


Cela suffit souvent à faire revenir le test de la qualité globale de l'ajustement dans de meilleures dispositions.

NB:

1) mettez vos variables dans un data.frame et utilisez l'argument data, que ce soit pour glm ou betabin. Cela vous évitera des ennuis.

2) il vaut mieux repartir des comptages initiaux pour éviter les problèmes d'arrondis: R va se plaindre si le numérateur et le dénominateur de la proportion ne sont pas entiers.

Si vous voulez utiliser un modèle de Poisson, il faut que votre réponse soit un comptage et non une proportion. Ce serait donc plutôt Nfeuille * Incidence, en prenant en compte la remarque 2 ci-dessus. Les résultats devraient être proches de la régression binomiale. La surdispersion se juge de la même manière que pour un modèle binomial, et une manière simple de la prendre en compte le cas échéant est l'utilisation d'une distribution binomiale négative: fonction glm.nb de MASS (plutôt que de negbin dans aod).
Renaud

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

Messagepar myriam desanlis » 24 Aoû 2011, 09:21

Bonjour,

un grand merci pour la réponse!!! Une petite question concernant le paramètre de surdispersion : à partir de quelle valeur on considère qu'il y a surdispersion?>1?

Je viens de tester la fonction betabin en utilisant mon nombre de contamination (Nconta) à la place de l'incidence. Tout à l'air ok. J'ai juste une question sur l'interprétation de mes données (je n'ai pas l'habitude d'utiliser glm) : aucun facteur n'a l'air d'être significatif même si par exemple pour le facteur Azote (Az), j'ai plus de contamination avec la modalité 150, est-ce bien cela?

Code : Tout sélectionner

> mod1 <- betabin(cbind(Nconta, Nfeuille - Nconta) ~ Az + Var + Con, random = ~ 1, data=donnee)
> summary(mod1)
Beta-binomial model
-------------------
betabin(formula = cbind(Nconta, Nfeuille - Nconta) ~ Az + Var +
    Con, random = ~1, data = donnee)

Convergence was obtained after 1859 iterations.

Fixed-effect coefficients:
              Estimate Std. Error    z value Pr(> |z|)
(Intercept)  5.477e-02  4.964e-02  1.103e+00 2.699e-01
Az150        8.934e-02  3.797e-02  2.353e+00 1.863e-02
VarKe        4.435e-01  6.064e-02  7.315e+00 2.580e-13
VarLG        1.376e-02  6.016e-02  2.286e-01 8.192e-01
VarNK       -1.270e-01  5.964e-02 -2.129e+00 3.327e-02
VarSI       -1.320e-01  5.811e-02 -2.272e+00 2.310e-02
ConS         8.652e-02  3.797e-02  2.279e+00 2.269e-02

Overdispersion coefficients:
                 Estimate Std. Error   z value   Pr(> z)
phi.(Intercept) 2.009e-02  3.073e-03 6.537e+00 3.132e-11

Log-likelihood statistics
   Log-lik      nbpar    df res.   Deviance        AIC       AICc
-1.590e+03          8        592  9.163e+02  3.196e+03  3.197e+03
> x1=sum(residuals(mod1,type="pearson")^2)
> dd1=df.residual(mod1)
> 1-pchisq(x1,dd1)
[1] 0.4935208

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 24 Aoû 2011, 09:49

myriam desanlis a écrit :Bonjour,

un grand merci pour la réponse!!! Une petite question concernant le paramètre de surdispersion : à partir de quelle valeur on considère qu'il y a surdispersion?>1?


Oui... Le test sur phi dans la sortie de betabin vous donne une réponse approximative à la question.

Je viens de tester la fonction betabin en utilisant mon nombre de contamination (Nconta) à la place de l'incidence. Tout à l'air ok. J'ai juste une question sur l'interprétation de mes données (je n'ai pas l'habitude d'utiliser glm) : aucun facteur n'a l'air d'être significatif même si par exemple pour le facteur Azote (Az), j'ai plus de contamination avec la modalité 150, est-ce bien cela?


Les valeurs de p indiquées dans les sorties sont des tests de Wald testant H0 = le coefficient a une valeur de 0. La plupart des coefs sont donc signficativement différents de 0 au seuil alpha = 0.05 (tous sauf varLG).

Si vous voulez tester l'effet global d'une variable, ajuster le modèle avec et sans cette variable, et faites un test du rapport de vraisemblance entre les deux modèles:

Code : Tout sélectionner

mod1 <- betabin(cbind(Nconta, Nfeuille - Nconta) ~ Az + Var + Con, random = ~ 1, data=donnee)
mod2 <- betabin(cbind(Nconta, Nfeuille - Nconta) ~ Var + Con, random = ~ 1, data=donnee)
anova(mod1, mod2)


Comme votre variable Az n'a que deux modalités, ce test devrait donner le même résultat que le test de Wald: p = 1.863e-02 (< alpha = 0.05).
Renaud

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

Messagepar myriam desanlis » 24 Aoû 2011, 09:58

Un grand merci pour les précisions

Dernière question pour vérifier que j'ai tout compris: si je teste les interactions et que j'obtiens:

Code : Tout sélectionner

Az150:VarNK:ConS  3.879e-01  2.295e-01  1.690e+00 9.102e-02


ca veut bien dire que j'ai une interaction significative pour la modalité Az150, la conduite S et la variété NK ?

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 24 Aoû 2011, 10:35

Non: 9.102e-2 = 0.09102 (> alpha = 0.05).

Je n'aime pas ces tests partiels (le danger des comparaisons multiples rôde !). Faire d'abord le test global pour l'interaction azote * conduite * variété.
Renaud

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

Messagepar myriam desanlis » 24 Aoû 2011, 11:01

oups j'avais mal décalé la virgule dans ma tête...dsl

en fait cette ligne vient du test du modèle :

Code : Tout sélectionner

mod2 <- betabin(cbind(Nconta, Nfeuille - Nconta) ~ Az * Var * Con, random = ~ 1, data=donnee)


donc a priori pas de danger de comparaisons multiples

J'ai aussi testé si j'avais des différences entre le modèle de base sans interaction et celui avec interaction :

Code : Tout sélectionner

> anova(mod1,mod2,test="Chisq")
Analysis of Deviance Table (beta-binomial models)

mod1: fixed = cbind(Nconta, Nfeuille - Nconta) ~ Az + Var + Con; random = ~1
mod2: fixed = cbind(Nconta, Nfeuille - Nconta) ~ Az * Var * Con; random = ~1

      logL  k  AIC AICc  BIC Resid. dev. Resid. Df      Test Deviance Df P(> Chi2)
mod1 -1590  8 3196 3197 3232       916.3       592                               
mod2 -1562 21 3167 3169 3259       860.9       579 mod1-mod2    55.37 13 3.476e-07


pour moi, ils sont différents et d'après les AIC celui avec interaction est meilleur que celui sans interaction (même si ca n'a pas l'air d'être énorme comme différence)

Code : Tout sélectionner

> res=AIC(mod1,mod2)
> res
     df      AIC     AICc
mod1  8 3196.350 3196.593
mod2 21 3166.977 3168.576

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 24 Aoû 2011, 15:08

Si, la différence est énorme même avec AICc: Delta(AICc) de 28 ==> le modèle avec interaction est de très loin meilleur, ce qui est cohérent avec le test du rapport des vraisemblances.
Renaud

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

Messagepar myriam desanlis » 25 Aoû 2011, 12:04

Bonjour

tout d'abord encore merci pour toutes vos indications!!
Pour l'AIC, à partir de quand considère-t-on (approximativement) que la différence est grande?
J'ai un autre jeu de données que je viens de traiter et j'aurais besoin de savoir si j'interprète bien les résultats.
J'ai des données de comptage de symptômes que je ne peux pas cette fois ci rapporter à ma population de plantes puisque les attaques peuvent avoir lieu partout. J'ai donc choisi une loi de Poisson. Le paramètre de surdispersion est également > 1 donc j'ai utilisé glm.nb comme vous me l'aviez suggérer :

Code : Tout sélectionner

> mod1 <- glm.nb(Nconta ~ Dens * Var * Con, data=donnee)   
> summary(mod1)

Call:
glm.nb(formula = Nconta ~ Dens * Var * Con, data = donnee, init.theta = 6.283334928,
    link = log)

Deviance Residuals:
    Min       1Q   Median       3Q      Max 
-2.2939  -0.8253  -0.2002   0.5271   4.1566 

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)   
(Intercept)        1.780586   0.104513  17.037  < 2e-16 ***
DensD2            -0.969656   0.202738  -4.783 1.73e-06 ***
VarKe              0.049325   0.146887   0.336  0.73702   
VarLG              0.208704   0.145355   1.436  0.15105   
VarNK              1.279685   0.133381   9.594  < 2e-16 ***
VarSI              1.319506   0.133152   9.910  < 2e-16 ***
ConS              -0.927096   0.190686  -4.862 1.16e-06 ***
DensD2:VarKe       0.199136   0.266852   0.746  0.45552   
DensD2:VarLG       0.714967   0.249990   2.860  0.00424 **
DensD2:VarNK       0.121316   0.238590   0.508  0.61112   
DensD2:VarSI      -0.006181   0.239162  -0.026  0.97938   
DensD2:ConS        0.809313   0.324522   2.494  0.01264 * 
VarKe:ConS         0.500009   0.245999   2.033  0.04210 * 
VarLG:ConS         0.478251   0.242928   1.969  0.04899 * 
VarNK:ConS        -0.231067   0.231320  -0.999  0.31784   
VarSI:ConS        -0.158093   0.230037  -0.687  0.49193   
DensD2:VarKe:ConS -0.150633   0.410679  -0.367  0.71378   
DensD2:VarLG:ConS -0.822103   0.393496  -2.089  0.03669 * 
DensD2:VarNK:ConS -0.039773   0.379143  -0.105  0.91645   
DensD2:VarSI:ConS -0.271465   0.381825  -0.711  0.47710   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(6.2833) family taken to be 1)

    Null deviance: 1292.2  on 545  degrees of freedom
Residual deviance:  533.4  on 526  degrees of freedom
AIC: 2840.2

Number of Fisher Scoring iterations: 1


              Theta:  6.283
          Std. Err.:  0.750

 2 x log-likelihood:  -2798.244
> x1=sum(residuals(mod1,type="pearson")^2)
> dd1=df.residual(mod1)
> 1-pchisq(x1,dd1)
[1] 0.107817


Est-ce le même système d'interprétation qu'avec betabin, à savoir que si la valeur de p est <0.05 le facteur est significatif et que la première colonne permet de classer les modalités les unes par rapport aux autres?

myriam desanlis
Messages : 27
Enregistré le : 10 Avr 2009, 08:19

Messagepar myriam desanlis » 25 Aoû 2011, 14:55

Nouveau problème à l'horizon:

avec un de mes jeux de données de comptage, j'ai beaucoup de zéros (environs 2000 contre 2074 données au total), le test du Chi 2 de Pearson donne 2.213337e-07 (modèle avec glm.nb).
J'ai fait une recherche sur les modèles de poisson avec excès de zéro et j'ai testé zeroinfl (pck pscl) qui donne un Chi 2 de 0. Avec ce modèle, plus aucun facteur n'est significatif contrairement au modèle précédent.

Est-ce qu'il existe d'autre type de modèle plus adapté à ce cas particulier

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 26 Aoû 2011, 10:25

Le fait d'avoir bcp de zéros n'est pas forcément une indication qu'il faut utiliser un modèle avec excès de zéros: si la moyenne lambda d'une loi de Poisson est faible, normal qu'l y ait bcp de zéros:

Code : Tout sélectionner

> x <- rpois(lambda = .01, n = 2000)
> table(x)
x
   0    1
1979   21


Il faut juger de l'excès de zéro par rapport à une loi de Poisson conditionnellement aux effets fixes: pour chaque combinaison, calculer lambda et comparer prévisions et observations.

Pour les modèles avec excès de zéros, j'aime bien la fonction hurdle du package pscl car les deux processus (génération des zéros et comptages censurés) et leurs coefs s'interprètent facilement. Voir l'article JSS accompagnant le package: http://www.jstatsoft.org/v27/i08
Renaud


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité

cron