Vraissemblance glm(quasi poisson)

Bonjour,

Comme le titre l'indique, je voudrai savoir quelles méthodes sont applicables pour estimer la pertinence d'un GLM où la fonction de lien est de la famille quasi-poisson. Ma variables explicative p est un ratio de comptage et est "injectée" dans le modèle par log(p/1-p) avec en offset le log du dénominateur du ratio p.

D'avance merci!
Cdlt.

Alex

Salut,

Je ne suis pas sûr de bien comprendre ta démarche. Une famille quasi-poisson n'est pas une fonction de lien mais une distribution d'erreur. Le lien est en général la fonction log pour les données de comptage.
Je ne vois pas trop non plus le pourquoi du comment tu intègres une variable transformée en logit. Pour l'offset je comprends à peu près l'idée bien que je ne sois pas bien sûr du log mais bon j'ai pas trop l'habitude des offsets.

Je suppose que tu prends une quasi à cause de problème de surdispersion ? Un modèle binomial négatif fera tout aussi bien l'affaire avec une base statistique plus robuste et bien plus d'outils pour évaluer ton modèle. Le must pour moi étant un modèle beta permettant de modélisant indépendamment la moyenne et la variance.

Pour enfin répondre à ta question, McCullagh & Nelder (1989) préconisent de tester un modèle quasi via un test du rapport des vraisemblances qui suit alors une loi de Fisher.

dans R:

Code:

anova (mod0,mod1, test="F")

Nik

quasi poisson est souvent utilisé pour prendre en compte de la surdispersion dans les données tout comme quasibinomial.
Pour le reste je partage l'avis de Nik, tout ça n'est pas clair !

Merci d'avoir pris le temps de me répondre.

Clairement, ma question trahi mon inexpérience dans le domaine et ma maîtrise imparfaite du vocabulaire. Je reprends donc pour clarifier en essayant de rester succinct. J'ai des données de comptages correspondant au dénombrement de parasites tués par 3 produits testés, appelons ça x. Un traitement de contrôle est réalisé après coup pour avoir une efficacité relative et pour lequel un 2e comptage est fait, disons y. Le ratio sur lequel je travail est : p= x/(x+y). J'ai pu voir que la méthode préconisée pour faire un glm (le but étant purement descriptif) sur de telles données était bien de faire une transformation log de la façon suivante :
log(p/(1-p))
Dans le cadre de ce glm, il est également préconisé de prendre une loi de poisson pour la distribution de l'erreur (et non en fonction de lien) et log comme fonction de lien, tout en plaçant log(y) en offset. J'ai une surdispersion qui m'impose de passer par du quasipoisson (le facteur de dispersion est supérieur à 5 dans le modèle suivant). La formule du modèle est donc :

glm(log(p/(1-p))~site+type+S0, offset = log(y), family = "quasipoisson") # Vous l'avez compris, je travaille sous R.

"S0" est une covariable dont je veux tester l'effet, "type" est le type de traitement testé et "site" renseigne sur la provenance des données. J'ai également d'autres covariables dont j'aimerai tester l'effet dont des données thermiques ainsi que des variables qui peuvent éventuellement être impacté par les traitements testés.

Maintenant que je pense avoir clarifié la situation j'en profite pour vous poser les questions sur lesquelles je n'arrive pas à trancher malgré de nombreuses recherches :
1) Les différents log présents dans le modèle me donnent des valeurs négatives, ce qui pose problème pour le recours à une quasipoisson. J'ai donc ajouté une constante à la sortie de ces log (+2, n'apparait pas dans le code ci-dessus). Conscient que c'est du "bricolage", je pense que cela ne change pas la donne au niveau de la variabilité des données et de la régression, toutefois je m'interroge sur le bien-fondé de cette pratique.

2) Ma première question portait sur la vraisemblance du modèle. Comme je suis dans une logique descriptive (et non prédictive), le but était de faire un tri sur les covariables (donc de ne retenir que la combinaison donnant le modèle le mieux "fitté"). Il me semble que c'est la démarche indiquée et c'est pq j'étais un peu désemparé de voir les AIC et autre likelihood estimates non disponible (et j'ai conscience que c'est pour une bonne raison). Je me suis donc simplement basé sur la variance résiduelle pour le moment, les écarts étant souvent important d'un modèle à l'autre, ainsi que sur l'absence de significativité des covariables pour faire mon tri. Est-ce suffisant? J'ai vu lors de mes recherches que certains utilisent les résultats obtenus par un glm "poisson" pour faire un indice approximatif (mais considéré comme acceptable) pour ce type de sélection (voir: https://stat.ethz.ch/pipermail/r-help/2008-October/178537.html). Dans mon cas, j'ai un message d'erreur lorsque j'utilise la loi de poisson : "In dpois(y, mu, log = TRUE) : non-integer x = 6.31..." J'avoue que je ne le comprends pas... Pouvez vous m'éclairer sur ce point? Et donc pour résumer : Ai-je vraiment besoin d'un indice de vraisemblance dans ma démarche? Si oui, quel recours sont disponibles? Je suis en train de me renseigner sur la fonction dredge() du pkg MuMIn, peut être est-ce une solution : http://cran.r-project.org/web/packages/bbmle/vignettes/quasi.pdf

3) La variance résiduelle étant, dans le meilleur des cas, d'environ 50 % de la variance nulle (Null deviance), puis-je tout de même faire une interprétation des p.valeurs significatives à l'issu du modèle? Les résultats sont cohérents par rapport à ce que j'ai pu voir lors des analyses préliminaires...

4) Conscient des problèmes de sur-ajustement, j'ai tout de même intégré les covariables potentiellement co-évolutives par curiosité. Elles se sont révélées non-explicatives et donc la variance résiduelle n'a pas fortement diminué. Puis-je me servir de ces résultats pour rejeter leur intérêt dans les expérimentations futures?

Voilà les questions de bases que je me pose, vous l'aurez compris je débute dans ce type d'analyse mais après moult recherches j'ai le sentiment d'être dans un cas de figure assez particulier (bien que répandu) et les docs dispo sur le net portent généralement plus sur la théorie que la pratique et sur la réalisation que sur l’interprétation des sorties de modèle, d'où ma présence sur ce forum. Si vous avez de la doc à me conseiller je suis donc preneur!

Je vous remercie encore de l'attention et du temps que vous voudrez bien me consacrer.
Cordialement

Bon, rapidement. Pour commencer il semble nécessaire que tu lises des choses sur les glms. Je te conseille vivmement, extending the minear model de Faraway dispo ici : Faraway. Quand tu utilises un glm, la relation linéaire entre les variables explicatives et la variable a expliquer ne se fait pas directement, pour cela tu passes par une fonction de lien. Autrement dit, au lieu d'avoir Y = aX+b, tu as f(Y) = aX+b. Par exemple si le lien est un log alors tu as : f=log et log(Y) = aX+b.
Quand dans ton modèle tu mets log(p/(1-p) ~ X, avec comme argument family = quasipoisson, ça veut dire que le logiciel va faire : log(log(p/(1-p))) ~ aX + b. Autrement dit ça n'a pas de sens.
Pour modéliser une proportion tu as deux possibilités, tu utlises la distribution binomial avec par exemple un lien logit, ce qui veut dire que tu vas estimer les coefficients de sorte que : log(p/(1-p)) = aX + b. Dans ce cas p est nomalemet la probabilité de "succès", ici ce sera ta proportion. Pour faire ça dans R :
glm(cbind(n,N-n)~X, family=binomial), avec n qui est le nombre de succès et N le nombre d'observations. Si tu as une proportion de 0.6 pour 100 individus ça donne : glm(cbind(60,40)~X, etc.).
Sinon tu peux aussi utiliser une famille de poisson en intégrant pour offset log(N). Un offset est une variable dont on fixe le coefficient à 1. Dans ces cas la voila ce qui se passe : si le lien est un log, alors tu auras : log(n) ~ aX+b + log(N), soit log(n)-log(N) = aX+b et donc log(n/N) = aX+b.
Sous R ça donne : glm(n ~ X+log(N), family=poisson).
Le choix de la distribution de poisson ou de quasi poisson se fait sur une éventuelle sur-dispersion des données.

Si le but est de tester les effets des différentes variables alors tu n'as pas de sélection de modèle à faire. Tu fais le modèle complet et tu testes les différents effets.

La distribution de poisson est faite pour des données de comptage, pas des données continues. Donc si tu lui mets autre chose que des entiers la fonction foire.

Dans le cas de glm (or gaussien), on ne résonne pas sur les variances résiduelles, mais sur les déviances.
HTH

Merci pour la référence, celle-ci m'avait échappé. Merci également pour toutes ces précisions. Le modèle que j'ai présenté ici est issu d'une discussion téléphonique avec un statisticien "puriste" (malheureusement je n'ai eu qu'un contact avec lui et il est indisponible pdt un long moment), et il y a eu confusion (de ma part) au sujet de la transformation logarithmique. J'avais bien compris ce qu'était la fonction de lien et comment elle fonctionnait (du moins, après coup) et cette question log(log(x)) m'a alors interpellé, j'aurais du être plus vigilant et m'y intéresser d'avantage. Lorsqu'il m'a conseillé de faire du log, je n'avais donc pas compris que c'était implicitement fait par le glm. Mea Culpa! Finalement je m'aperçois qu'il faisait soit référence au lien logistique utilisé lors du recours à la famille binomiale, soit au fait d'intégrer le offset dans le poisson. Merci également pour la traduction du offset dans la formule, maintenant je comprends de fait à quoi cela correspond (au delà de son utilité, vu par ex page 11 de cette doc http://www.casact.org/pubs/forum/07wforum/07w263.pdf).

Je crois avoir bien compris les modèles que tu proposes. Je vais essayer de rester sur du poisson (puisque mon approche était dans cet esprit) pour le moment, je verrai le cas binomial plus tard (la notion de "succès" me dérange un peu dans mon cas)... Je reviens ici en cas de besoin, sans vouloir abuser  

Merci vraiment beaucoup pour ton temps et pour la précision de tes réponses. Il est parfois difficile d'avoir du concret par soit même.
Cdlt.

Bonjour,
J'ai donc opté pour le modèle: glm(n~X, offset=N, family="poisson")
Et R me retourne le message suivant : "Error: no valid set of coefficients has been found: please supply starting values"

Il s'agit bien de valeur servant de point de départ pour finalement obtenir les paramètres qui "fit" le mieux le modèle après x itération? J'ai essayé de lui fournir un peu le tout venant (coef issus d'un quasipoisson, plus souple etc...) mais en vain. En creusant un peu, notamment là http://stackoverflow.com/questions/8212063/r-glm-starting-values-not-accepted-log-link et ailleurs, la présence de 0 dans l'offset serai une source éventuelle du problème, mais dans mon cas je n'ai pas de 0 nul part (min(N)=57).
Cela dit, au regard de ce que tu m'as expliqué dans les précédentes réponses, j'ai essayé de formuler de cette manière le modèle :
glm(n ~ X+log(N), family="quasipoisson") # On est bien d'accord c'est log(N) que j'introduis?
Et là tout va bien. Je compte donc travailler sur ce modèle mais je me demandais si l'erreur mentionnée plus haut viens de la gestion de l'offset par R ou bien d'autres subtilités théoriques sur l'utilisation d'un offset en glm qui impliquerai que ce n'est pas strictement la même chose de mettre N en offset ou sous la forme +log(N) dans la formule...

J'espère que la lecture du livre de Julian Faraway m'éclairera de manière générale mais je ne peux malheureusement pas mettre ce travail en stand-by pour le moment. Et donc merci encore de ton aide.
Bien cordialement.

On est bien d'accord c'est log(N) que j'introduis?

oui sinon ça n'a pas de sens. Regarde le calcul indiqué par Droopy sur la manière dont l'offset dans la relation linéaire avec la fonction de lien.

Je crois pas ailleurs qu'une confirmation t'a été donnée par Renaud Lancelot sur un autre forum bien connu de certains ici

Damned! Je suis repéré!  
Plus sérieusement, j'avais oublié le log qd j'ai fait l'offset, c'est pour ça que ça ne marchais pas. Cela dit, ça ne reviens pas exactement au même sous R puisque le coef est alors recalculé et on a un effet de log(N) dans la sortie du modèle. Dans l'analyse résiduelle il emporte donc clairement l'explication de la variance résiduelle, ce qui n'est pas étonnant puisque "N" inclus "n". Bref, j'ai bien compris l'idée sur le papier grâce aux indications de Droopy, mais je doute pouvoir m'appuyer sur un tel modèle pour R, l'argument offset existe pour cette raison, j'imagine. Je ne peux retourner sur le forum GuR pour le moment (je suis en congé aujourd'hui et j'ai oublié mon mdp...) mais si j'ai bonne mémoire Renaud passe par predict() et je n'ai pas eu le temps d'étudier sa réponse. Au final, maintenant la modélisation fonctionne très bien avec glm(n~X, offset=log(N), family="quasipoisson") et je vais rester là dessus...
Par rapport à la discussion des résultats (voir mon autre post mais je groupe ici), finalement, un modèle nous dit si il existe un lien significatif de nature logarithmique (dans mon cas) entre la variable réponse (dont la distribution suit une loi de Poisson, comme le sera celle des erreurs du modèle) et les variables explicatives, et le "sens" (>0 ou <0 en qq sortes) de ce lien. J'ai tout bon ou je suis à coté de la plaque?
Merci encore pour votre patience
Cdlt
Alex

un modèle nous dit si il existe un lien significatif de nature logarithmique (dans mon cas) entre la variable réponse (dont la distribution suit une loi de Poisson, comme le sera celle des erreurs du modèle) et les variables explicatives

Si on veut être complètement rigoureux sur la démarche, le modèle fait l'estimation des paramètres selon un critère choisi (maximum de vraisemblance) puis quantifie la variabilité de la réponse associée aux diverses variables explicatives. L'aspect test et donc significativité est extérieur au modèle et repose sur une stratégie d'analyse définie par l'utilisateur.
Donc attention, t'es pas à côté de la plaque mais plutôt "border line" . Un modèle porte bien son nom : il modélise un processus par rapport aux contraintes que tu lui impose par tes choix (type de modèle, transfo des données etc...etc...). Il n'est pas en mesure de quantifier l'importance des effets (valeurs des paramètres).

Nik

PS : oui le plus simple est de faire via l'argument offset.

D'accord... On a tout de même in fine une quantification de l'information portée par les variables explicatives utilisée pour appréhender l'effet de ces variables sur la réponse... J'avais cru comprendre que c'était le postulat ou la démarche lorsqu'on utilise les modèles dans une logique descriptive... On est donc forcément border line dans cette optique?

Si on veut être complètement rigoureux sur la démarche

J'y tiens absolument, c'est ma plus grande qualité et mon plus grand problème!

J'ai fini par comprendre ce qu'était un offset, avec le recul mes questions n'avaient pas beaucoup de sens, mais elles m'ont permis d'apprendre

oui mon message n'était pas clair. Le modèle quantifie bien mais il n'est pas en mesure de juger si c'est important ou pas. Les quelques tests stats utilisés n'apportent d'ailleurs pas bcp d'info à ce sujet en général.