Groupe des utilisateurs du logiciel R Index du Forum Groupe des utilisateurs du logiciel R
Un forum francophone d'échange autour du logiciel de calcul statistique R dans le domaine de la recherche agronomique tropicale
 
 FAQFAQ   RechercherRechercher   Liste des MembresListe des Membres   Groupes d'utilisateursGroupes d'utilisateurs   S'enregistrerS'enregistrer 
 ProfilProfil   Se connecter pour vérifier ses messages privésSe connecter pour vérifier ses messages privés   ConnexionConnexion 

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

 
Poster un nouveau sujet   Répondre au sujet    Groupe des utilisateurs du logiciel R Index du Forum -> Questions en cours
Voir le sujet précédent :: Voir le sujet suivant  
Auteur Message
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Mer Aoû 24, 2011 8:05 am    Sujet du message: glm poisson ou binomiale, pb de contradiction Chi 2/donnees Répondre en citant

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:
> 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:
> 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:
> 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:
> 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 ?
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
Renaud Lancelot



Inscrit le: 16 Déc 2004
Messages: 2459
Localisation: CIRAD, Montpellier

MessagePosté le: Mer Aoû 24, 2011 8:52 am    Sujet du message: Répondre en citant

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

Code:
> 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:
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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé Envoyer un e-mail Visiter le site web de l'utilisateur
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Mer Aoû 24, 2011 9:21 am    Sujet du message: Répondre en citant

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:
> 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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
Renaud Lancelot



Inscrit le: 16 Déc 2004
Messages: 2459
Localisation: CIRAD, Montpellier

MessagePosté le: Mer Aoû 24, 2011 9:49 am    Sujet du message: Répondre en citant

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.

Citation:
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:
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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé Envoyer un e-mail Visiter le site web de l'utilisateur
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Mer Aoû 24, 2011 9:58 am    Sujet du message: Répondre en citant

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:
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 ?
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
Renaud Lancelot



Inscrit le: 16 Déc 2004
Messages: 2459
Localisation: CIRAD, Montpellier

MessagePosté le: Mer Aoû 24, 2011 10:35 am    Sujet du message: Répondre en citant

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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé Envoyer un e-mail Visiter le site web de l'utilisateur
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Mer Aoû 24, 2011 11:01 am    Sujet du message: Répondre en citant

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

en fait cette ligne vient du test du modèle :
Code:
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:
> 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:
> res=AIC(mod1,mod2)
> res
     df      AIC     AICc
mod1  8 3196.350 3196.593
mod2 21 3166.977 3168.576
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
Renaud Lancelot



Inscrit le: 16 Déc 2004
Messages: 2459
Localisation: CIRAD, Montpellier

MessagePosté le: Mer Aoû 24, 2011 3:08 pm    Sujet du message: Répondre en citant

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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé Envoyer un e-mail Visiter le site web de l'utilisateur
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Jeu Aoû 25, 2011 12:04 pm    Sujet du message: Répondre en citant

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:
> 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?
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
myriam desanlis



Inscrit le: 10 Avr 2009
Messages: 27
Localisation: Toulouse

MessagePosté le: Jeu Aoû 25, 2011 2:55 pm    Sujet du message: Répondre en citant

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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé
Renaud Lancelot



Inscrit le: 16 Déc 2004
Messages: 2459
Localisation: CIRAD, Montpellier

MessagePosté le: Ven Aoû 26, 2011 10:25 am    Sujet du message: Répondre en citant

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:
> 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
Revenir en haut de page
Voir le profil de l'utilisateur Envoyer un message privé Envoyer un e-mail Visiter le site web de l'utilisateur
Montrer les messages depuis:   
Poster un nouveau sujet   Répondre au sujet    Groupe des utilisateurs du logiciel R Index du Forum -> Questions en cours Toutes les heures sont au format GMT
Page 1 sur 1

 
Sauter vers:  
Vous ne pouvez pas poster de nouveaux sujets dans ce forum
Vous ne pouvez pas répondre aux sujets dans ce forum
Vous pouvez éditer vos messages dans ce forum
Vous ne pouvez pas supprimer vos messages dans ce forum
Vous ne pouvez pas voter dans les sondages de ce forum


Powered by phpBB © 2001, 2005 phpBB Group
Traduction par : phpBB-fr.com

Anti Bot Question MOD - phpBB MOD against Spam Bots
Inscriptions bloqués / messages: 80144 / 702