Les bookmarks de Clem
2353 shaares
7 results
tagged
poisson
Équivalence entre un modèle de survie à risques proportionnels et un modèle de Poisson. Le truc est très rigolo, et va grandement me simplifier la vie.
Au passage, il y a des liens vers ses cours, qui sont super intéressants...
Au passage, il y a des liens vers ses cours, qui sont super intéressants...
Ver Hoef et Boveng, 2007. Très intéressant. Résumé:
Les auteurs comparent l'approche quasi-poisson et binomiale négative pour prendre en compte la surdispersion. Sur un plan théorique et sur un plan pratique. Les deux principales différences entre ces approches sont:
* Dans la relation entre moyenne et variance: pour la quasi-Poisson, on a (Var = theta*mu) et pour la binomiale négative, on a (Var = mu + kappa*mu). Pour savoir lequel des deux est meilleur, il est recommandé d'ajuster les deux modèles, puis de représenter les carrés des résidus (y-mu)^2, qui représentent la variance, en fonction de mu. Comme ces graphes sont en général assez bordéliques, les auteurs recommandent de découper en catégories de mu et de calculer la moyenne des carrés des résidus (donc la variance) dans chaque catégorie. La relation entre les deux est-elle linéaire ou quadratique?
* Dans les poids pris par les observations lors de l'ajustement. En général, on utilise l'IRLS pour ajuster ces modèles. C'est un moindre carré dans lequel on utilise une matrice de poids particulière pour les observations. La seule différence entre quasi-poisson et binomiale négative tient dans ces poids (le reste est identique entre les deux approche). On voit alors que
- Pour la quasi-Poisson, le poids de l'observation i est mu_i/theta (avec theta le coef de surdispersion)
- Pour la binomiale négative, le poids est (mu_i / (1+kappa*mu_i))
Donc, quand mu_i devient grand, le poids de l'observation i devient grand en quasi-poisson, alors qu'il tends vers 1/kappa avec la BN. Il faut alors se poser la question du comportement le plus désirable en fonction de l'objectif. Dans celui présenté par les auteurs, i.e. estimer l'effectif de phoques, le comportement de la BN est problématique: "Our goal is to estimate overall abundance, which is dominated by the larger sites, and we prefer to have adjustments dominated by the effects at those larger sites". En plus, le graphe suggéré au premier point ci-dessus tend à favoriser la quasi-poisson.
J'aime bien la conclusion: "an important way to choose an appropriate model is based on sound scientific reasoning rather than a data-driven method". Toujours bon à rappeler. J'aime bien ces auteurs.
Les auteurs comparent l'approche quasi-poisson et binomiale négative pour prendre en compte la surdispersion. Sur un plan théorique et sur un plan pratique. Les deux principales différences entre ces approches sont:
* Dans la relation entre moyenne et variance: pour la quasi-Poisson, on a (Var = theta*mu) et pour la binomiale négative, on a (Var = mu + kappa*mu). Pour savoir lequel des deux est meilleur, il est recommandé d'ajuster les deux modèles, puis de représenter les carrés des résidus (y-mu)^2, qui représentent la variance, en fonction de mu. Comme ces graphes sont en général assez bordéliques, les auteurs recommandent de découper en catégories de mu et de calculer la moyenne des carrés des résidus (donc la variance) dans chaque catégorie. La relation entre les deux est-elle linéaire ou quadratique?
* Dans les poids pris par les observations lors de l'ajustement. En général, on utilise l'IRLS pour ajuster ces modèles. C'est un moindre carré dans lequel on utilise une matrice de poids particulière pour les observations. La seule différence entre quasi-poisson et binomiale négative tient dans ces poids (le reste est identique entre les deux approche). On voit alors que
- Pour la quasi-Poisson, le poids de l'observation i est mu_i/theta (avec theta le coef de surdispersion)
- Pour la binomiale négative, le poids est (mu_i / (1+kappa*mu_i))
Donc, quand mu_i devient grand, le poids de l'observation i devient grand en quasi-poisson, alors qu'il tends vers 1/kappa avec la BN. Il faut alors se poser la question du comportement le plus désirable en fonction de l'objectif. Dans celui présenté par les auteurs, i.e. estimer l'effectif de phoques, le comportement de la BN est problématique: "Our goal is to estimate overall abundance, which is dominated by the larger sites, and we prefer to have adjustments dominated by the effects at those larger sites". En plus, le graphe suggéré au premier point ci-dessus tend à favoriser la quasi-poisson.
J'aime bien la conclusion: "an important way to choose an appropriate model is based on sound scientific reasoning rather than a data-driven method". Toujours bon à rappeler. J'aime bien ces auteurs.
Une discussion intéressante sur la modélisation des données de comptage.
Bon, j'y vois un peu plus clair sur cette histoire de surdispersion et de modèle log-linéaire sur table de contingence (décrit ici: http://caloine.ouvaton.org/shaarli/?voVBrw). Comme j'avais besoin d'écrire pour me fixer les idées sur un sujet, j'ai mis à plat mes cogitations dans un fichier texte au fur et à mesure de la progression dans ma recherche biblio, et j'ai synthétisé ma review pour y voir plus clair. Je copie donc le résultat de mes cogitations ici, si ça peut servir à d'autres. Pis moi, comme ça, je sais où le trouver. Je décris alors:
* Le contexte: la surdispersion dans le cas d'un modèle log-linéaire d'une variable distribuée selon une loi de Poisson.
* Le fond du problème qui m'intéresse: la surdispersion dans le cas d'un modèle supposé poissonnien ajusté à une table de contingence.
* Une remarque finale sur le diagnostic de la surdispersion dans le contexte "table de contingence"
EDIT: j'ai fait quelques edits pour préciser quelques points, et une nouvelle section sur le test du rapport de vraisemblance.
==========================================
Contexte: Un modèle log-linéaire d'une variable distribuée selon une loi de Poisson.
La notion de surdispersion décrit le cas où la variance de la variable réponse est supérieure à la variance nominale, i.e. celle qui est attendue sous l'hypothèse du modèle supposé pour cette variable (e.g. dans le cas Poissonnien, lorsque cette variance est supérieure à l'espérance de la variable). C'est la définition donnée par McCullagh & Nelder (1989).
Lorsque l'on ajuste un modèle log-linéaire reliant une variable poissonienne Y à un ensemble de variables explicatives X1, X2, ..., Xp, cette surdispersion peut être diagnostiquée en comparant la déviance résiduelle et le nombre de degrés de liberté résiduels du modèle. Asymptotiquement, la déviance résiduelle devrait être égale au nombre de degrés de libertés résiduels. Asymptotiquement. On *diagnostique* une surdispersion lorsque la déviance est bien supérieure au nombre de degrés de libertés. Après, il y a toute une discussion relative à "supérieure de combien?"; Lindsey (1999) indique qu'on doit la prendre en compte lorsque la déviance est plus de deux fois plus grande que le nombre de dll, ce qui est contesté -- à juste titre à mon avis -- par Ripley, cf. e.g. http://caloine.ouvaton.org/shaarli/?dGq1kA et le lien correspondant. En outre, Venables & Ripley (2002) mettent quand même fortement en garde contre une interprétation trop littérale du rapport déviance/ddl ("This can be seriously misleading"), car cette théorie est asymptotique, et ne s'applique que pour des espérances larges dans le cas Poisson. McCullagh et Nelder (1989, p. 36) indiquent d'ailleurs "In general, however, the \chi^2 approximations for the deviance are not very good even as n \rightarrow \infty".
Bon bref, ce n'est pas ce qui m'intéresse le plus ici. En fait, la notion de *diagnostic* est ici importante, car on verra plus loin dans le cas des tableaux de contingence que l'on peut passer à côté d'une surdispersion réelle d'une variable réponse avec ce simple diagnostic. Autrement dit, il peut y avoir une surdispersion de la variable réponse non diagnostiquée par le modèle (cf. deuxième partie de cette revue, j'y reviendrai).
Laissons de côté pour le moment la question du diagnostic de la surdispersion pour nous intéresser plus à celle de la solution à lui apporter. Et pour solutionner le problème de la surdispersion, il faut en connaître la cause. Qu'est ce qui cause la surdispersion dans un modèle supposé poissonnien? Cameron et Trivedi (1998) comme Agresti (2002) indiquent que la source de la surdispersion est l'unobserved heterogeneity, ce qui est presque une tautologie (basiquement, ça revient à dire qu'on a surdispersion parce que l'on a une variabilité plus grande).
Agresti (2002, p. 130) donne un exemple qui permet d'y voir plus clair: "A common cause of overdispersion is subject heterogeneity. For instance, suppose that width, weight, color, and spine condition are the four predictors that affect a female crab's number of satellites. Suppose that Y has a Poisson distribution at each fixed combination of those predictors. Our model uses width alone as predictor. Crabs having a certain width are then a mixture of crabs of various weights, colors, and spine conditions. Thus, the population of crabs having that width is a mixture of several Poisson populations, each having its own mean for the response. This heterogeneity results in an overall response distribution at that width having a greater variation than the Poisson predicts". De façon très intéressante, Agresti s'appuie sur cet exemple pour développer une solution possible au cas du crabe, en modélisant la réponse non pas comme une réponse Poissonnienne, mais comme une réponse binomiale négative.
C'est très intéressant, parce qu'à travers cet exemple, on voit les deux points de vue possibles sur la surdispersion: (i) un prédicteur linéaire mal spécifié (i.e. on est passé à côté d'une ou plusieurs variables explicatives importantes), (ii) un modèle d'erreur mal spécifié (i.e. si je considère l'exemple d'Agresti, conditionnellement à la largeur du crabe, les autres variables explicatives inconnues sont elles-même aléatoires, ce qui ajoute à la variabilité de la réponse: la réponse n'est plus une loi de Poisson, mais une distribution plus variable). Et l'on voit que ces deux causes possibles sont, dans une certaine mesure, équivalentes.
Ainsi, de nombreux auteurs distinguent ces deux causes comme source possible de surdispersion. Quand on diagnostique une surdispersion importante, il y a donc, potentiellement, deux familles de solutions possibles:
1. Modifier le prédicteur linéaire en ajoutant les variables essentielles (ici, voir les commentaires de la cellule d'appui stat de l'univ. de Cornell: http://caloine.ouvaton.org/shaarli/?voVBrw)
2. Modifier le modèle d'erreur: soit en définissant pour la variable réponse une distribution statistique qui va prendre en compte cette surdispersion, comme une distribution modèle binomial négative, soit en abandonnant la procédure d'estimation par le maximum de vraisemblance pour une procédure de type quasi-vraisemblance (i.e. en ne supposant aucune distribution précise pour la variable réponse, mais en supposant qu'il existe une distribution de la famille exponentielle à un seul paramètre pour laquelle la relation entre moyenne et variance est de la forme Var(mu) = c*mu, avec c un paramètre de surdispersion). Voir Ver Hoef et Boveng (2007) pour une comparaison des deux solutions binomialNégative vs. quasi-vraisemblance (je sors un peu du scope de ma review là).
Ces deux solutions sont également notées par de nombreux auteurs, comme Crawley (2007, qui en discute longuement), Cameron et Trivedi (1998), ou Williams (2002, dans un contexte binomial et pas poissonnien, mais c'est la même logique). Ce dernier indique explicitement: "even when all available explanatory variables have been fitted, the residual variation may be greater than can be attributed to the binomial sampling variation assumed by the model. In this event we can either seek additional explanatory variables, or postulate a source of extra-binomial random variation between observations".
Quelle approche préférer? la 1 ou la 2? Anderson et al. (1994) donnent une idée de la philosophie dans un contexte de modèles de capture-recapture: "Under the CJS [Cormack-Jolly-Seber] model theory, c=1; however, with real data we expect c>1, but we do not expect c to exceed 4 (see Eberhardt 1978). Substantially larger values of c (say, 6-10) are usually caused partly by a model structure that is inadequate, that is the fitted model does not actually represent all the explainable variation in the data. Quasi-likelihood methods of variance inflation are appropriate only after the structural adequacy of the model has been achieved". Je n'ai pas bien compris d'où Anderson et al. tirent leur valeur de 4 citée ici (je n'ai pas trouvé ça dans Eberhardt 1978; et Anderson et al. indiquent dans leur discussion que ce seuil de 5 est déduit des expériences dans le domaine de la CMR), mais ce n'est pas important ici. L'important est que si la surdispersion est modérée, on peut la prendre en compte par quasi-vraisemblance. Si elle est plus importante on est à côté de la plaque, et il faut laisser tomber le modèle construit.
Voir aussi la bonne synthèse de Carruthers et al. (2008) trouvée ici: http://www.mun.ca/biology/dschneider/b7932/B7932Final4Mar2008.pdf
En s'appuyant sur Anderson et al. et sur leur expérience, ils fixent entre 5 et 10 le "seuil" de surdispersion à partir duquel on peut raisonnablement commencer à taxer de mauvaise foi le modélisateur bourrin qui chercherait à prendre en compte la surdispersion par quasi-vraisemblance.
==========================================
Le fond du problème qui m'intéresse: le cas d'un modèle supposé poissonnien ajusté à une table de contingence.
Mon problème est le suivant. J'ai une table de contingence croisant trois facteurs A, B, et C. Il peut y avoir une trentaine de niveaux au facteur A (ce sont des espèces), une centaine au facteur B (ce sont des sites), et il y en aura toujours 2 pour le facteur C (ce sont des jeux de données).
Je note N ma variable réponse (effectifs de chaque espèce dans chaque site et chaque jeu de données, dans les cases du tableau). Des considérations théoriques m'amènent à ajuster le modèle suivant:
log(N) ~ Intercept + A + B + C + A:B + A:C + B:C (modèle 1)
Ce modèle ne sort pas du chapeau, il est le résultat d'un développement mathématique, et s'appuie sur des hypothèses que je suis prêt à poser sur mon processus biologique et sur le processus de collecte des données. Et surtout, sous ces hypothèses, ce modèle me permet d'estimer mes quantités d'intérêt (des effectifs) qui, sous les hypothèses biologiques que nous fixons sont les paramètres (Intercept+A+B+C+A:B). Je ne rentre pas dans le détail, mais disons que dans le cas présent, ce modèle a une justification théorique intéressante. Du coup, c'est ce modèle-là -- et aucun autre -- qui m'intéresse.
Lorsque je l'ajuste, j'ai une déviance, et un nombre de ddl résiduels. Du coup, il est tentant de comparer les deux pour diagnostiquer une surdispersion (Bon, dans mon cas précis, je n'ai pas surdispersion, mais dans un cas plus général, e.g. une autre application de la méthode, je pourrais avoir une grande déviance résiduelle associée à un tel modèle). Imaginons que j'ai une déviance importante associée à ce modèle. Que faire?
*À première vue*, il semble que l'on soit dans un cas particulier puisque les effectifs modélisés sont stockés dans un tableau de contingence à trois entrées. On connaît l'ensemble de facteurs qui permettraient de rendre nulle la déviance: les interactions de deuxième ordre A:B:C. Ajuster un modèle (1) avec ces interactions nous conduit à un modèle saturé, avec une déviance nulle. Donc, on sait comment faire pour supprimer cette surdispersion. Mais est-ce une démarche pertinente? En effet, si l'on optait pour cette solution:
1. ça servirait à quoi? Un modèle saturé permettrait en effet de supprimer la surdispersion, mais ne serait pas interprétable, et n'aurait aucun intérêt biologique, puisque ce qui nous intéresse, c'est une certaine combinaison des paramètres d'un modèle très précis.
2. Et surtout: il faut avoir conscience que ces interactions impliquent un grand très grand nombre de variables: Il peut y avoir une trentaine de niveaux au facteur A (ce sont des espèces), une centaine au facteur B (ce sont des sites), et il y en aura toujours 2 pour le facteur C. Donc, les interactions A:B:C pourront impliquer d'ajouter plusieurs centaines, voire des milliers, de variables à un modèle déjà compliqué!!!
Aitkin et al. (1989) discutent très en détail de ce cas particulier de la surdispersion dans le contexte des modèles log-linéaires sur tables de contingence (sa discussion se place dans un contexte binomial, mais est tout-à-fait valide dans un contexte Poissonnien): "When a contingency table is analysed, it may be found that high-order interactions are inexplicably significant, and sometimes even that no model other than the saturated model can provide a satisfactory representation of the data. This phenomenon frequently occurs with very large samples where the contingency table is not classified by some factors relevant to the response: the number of 'successes' in each cell of the table then has a /mixed/ binomial (or binomial mixture) distribution. That is, within a given cell the probability p of a success is not constant, but varies systematically with other factors which have not been recorded or included in the model. Since these factors are nt identified, the success probability p behaves like a random variable with a probability distribution".
Ces auteurs proposent alors les solutions classiques, notamment le principe de l'ajustement par quasi-vraisemblance, mais notent "This procedure is probably satisfactory for small amounts of overdispersion, but it is not a substitute for correct model specification, and it is impossible to test the goodness of fit of the model". Anderson et al. (1994) indiquent dans la discussion "In a treatise on analysis of count data, Cox and Snell (1989) assert that the simple approach, such as we have studied here, should often be adequate, as opposed to the much more arduous task of seeking an explanatory model for the overdispersion. This comment of Cox and Snell is supported by the results of Liang and McCullagh (1993), who found that causal modeling of overdispersion was clearly better than use of a single overdisperson parameter, c, in only one of five case examined".
En fait, comme on le voit ci-dessus, on n'est pas vraiment dans un cas particulier en termes de *solutions* à apporter à la surdispersion (mais voir plus bas concernant le problème du *diagnostic* de la surdispersion, qui lui, est particulier dans le cas des tableaux de contingence). Et c'est ça qui m'a pris du temps à comprendre. C'est pour cette raison que ce cas de figure, pourtant fréquent, n'est pas si souvent évoqué dans les bouquins. Autrement dit: on ne doit pas parler de l'alternative "ajustement du modèle saturé" vs. "modélisation de la surdispersion par quasi-vraisemblance" comme d'un choix cornélien dans ce contexte. On est dans le même cas de figure que les régressions Poisson classiques décrites ci-dessus, et la solution est de la même nature: le modèle (1) est le modèle d'intérêt. Si on a une déviance modérément importante, on peut corriger ça en ajustant le modèle par quasi-vraisemblance. Si on a une surdispersion très importante, c'est signe que le modèle est structurellement mauvais, et que les hypothèses sur lesquelles il repose sont probablement violées. Le modèle doit alors être évité, et une autre solution doit-être recherchée.
======================
Une remarque finale quand même, concernant le *diagnostic* de la surdispersion dans le contexte d'une analyse sur tableau de contingence.
Aitkin et al. (1989) notent que la surdispersion ne sera pas nécessairement *visible* dans un tableau de contingence. Par exemple, si une variable N est influencée essentiellement par un facteur X, mais que l'on somme ces valeurs de N sur tous les niveaux de deux autres facteurs A et B pour former une table de contingence, alors on peut parfaitement ajuster un modèle saturé N~A*B. Et comme le modèle est saturé, la déviance résiduelle sera nulle. On ne diagnostiquera pas la surdispersion occasionnée par le "collapsing" de N sur deux facteurs sans intérêt. Et pourtant, le modèle sera faux. Aitkin et al. donnent un exemple de ça, et notent "inappropriate collapsing of the table over important variable leads to distorted conclusions about the importance of other explanatory variables". La solution à ça consiste alors à inclure dès le départ un maximum de variables dépendantes dans la construction du tableau de contingence (et éviter les collapsings inappropriés). "Paradoxically, the higher the dimension of the table, the more likely we are to obtain a simple final model 'more is less'".
======================
Test du rapport de vraisemblance et Quasi-vraisemblance:
Un dernier point pour la route. Notre éditeur suggérait l'utilisation d'un test du rapport de vraisemblance entre le modèle (1) décrit plus haut et le modèle saturé comme test de qualité d'ajustement (goodness of fit ou GOF). Mais ça me posait un problème. En effet, ce test est valide sous l'hypothèse que le modèle d'erreur est correctement spécifié. Si l'on suppose un modèle de Poisson, alors la déviance va suivre, asymptotiquement un chi-deux (au passage attention: asymptotiquement, et beaucoup d'auteurs mettent en garde contre ça, cf. précédemment. McCullagh et Nelder (1989, p. 36) notent d'ailleurs que la table d'analyse de déviance ne doit être considérée, au mieux, que comme un "screening device", et qu'aucune tentative ne devrait être faite pour assigner des niveaux de significativité précis).
Pourtant, si l'on compare le modèle (1) au modèle saturé avec un test du rapport de vraisemblance, ce test peut apparaître significatif (indiquant donc un mauvais ajustement), et pourtant la solution ne serait pas de corriger ça par l'inclusion des interactions de deuxième ordre A:B:C manquantes. En effet, on peut parfaitement avoir une surdispersion modérée, mais qui apparaîtrait significative. Par exemple, imaginons que le modèle d'intérêt pour nous soit caractérisé par une déviance résiduelle de 100 sur 60 ddl. L'interaction de deuxième ordre A:B:C serait alors très hautement significative. Pourtant, cette surdispersion pourrait parfaitement être prise en compte à l'aide d'un ajustement par quasi-vraisemblance. Le paramètre de surdispersion, de phi = 100/60 = 1.67, ne serait même pas très élevé en comparaison de ce que l'on peut voir dans ce type d'étude.
On en revient toujours au même problème: le test du rapport de vraisemblance suggéré par notre éditeur indique une surdispersion, donc que le modèle de Poisson est mauvais. La solution à ce défaut peut être de modifier le modèle d'erreur en intégrant une surdispersion modérée, ou alors de modifier le prédicteur linéaire en intégrant les interactions. Oui, je sais je radote, je suis en train de redire ici ce que j'ai dit plus haut.
Mais le point intéressant ici, c'est ce qu'on trouve dans la littérature sur les tests de GOF dans le cas de la quasi-vraisemblance, à savoir pas grand chose, et le peu qu'on trouve n'est pas clair.
Par exemple, McCullagh (1983, Eq. 11) indique dans un papier très théorique que la différence de log-quasi-vraisemblance -- appelée déviance [et curieusement, pas quasi-déviance] -- entre deux modèles emboités m1 et m2 caractérisés par la même surdispersion phi, suit asymptotiquement un chi-deux multiplié par phi. Rien n'est dit dans ce papier théorique sur les conséquences pratiques de cette propriété (et sur le fait que de façon générale, phi ne sera jamais connue, donc sera estimée, et donc sur les conséquences de cette estimation sur la disttribution suivie par D/phi, où D est la déviance).
Mais Breslow (1990), dans l'intro de son article, cite la première édition du McCullagh et Nelder (de 1983) qui reporte cette dernière propriété, et indique en s'appuyant dessus que lorsque l'on cherche à comparer deux modèles, le test du (quasi) rapport de vraisemblance est disponible. Dans des simus, cet auteur utilise ce test en calculant pour chaque modèle une scaled deviance (deviance/phi, avec phi estimé classiquement à partir des données). Il en conclut que si la surdispersion est modérée, un tel test reste acceptable.
Pourtant, de façon intéressante, la propriété notée par McCullagh (1983) disparaît du McCullagh & Nelder, édition de 1989. Et cette dernière édition ne cite McCullagh 1983 que noyée au milieu d'autres refs en la référençant comme traitant de l'efficacité, de l'optimalité et de la robustesse (ce qui est vrai). Pas un mot sur le test... Cherchent-ils à éviter de cautionner cette pratique? Considèrent-ils ce résultat comme peu intéressant?
======================
Donc, en résumé, dans le contexte de la modélisation par un modèle log-linéaire d'une table de contingence:
* Le diagnostic de la surdispersion ne peut se faire que si l'on dispose des variables pertinentes comme entrées dans ce tableau (i.e. éviter de fusionner des dimensions pertinentes du tableau avant modélisation)
* La philosophie sous-tendant la correction des problèmes causés par la surdispersion est la même que dans les autres modèles log-linéaires: (i) si la surdispersion n'est pas très importante et que l'on ne tient pas à modifier la structure du modèle (soit parce qu'on n'a pas les variables qui en sont la cause, soit parce qu'on tient à cette structure pour d'autres raisons, notamment théoriques), on modifie le modèle d'erreur pour la réponse (e.g. ajustement par quasi-vraisemblance, définition d'un modèle binomial négatif); (ii) si la surdispersion est trop importante, le modèle est mauvais, il faut l'abandonner.
======================
Références:
Agresti, A. 2002. Categorical Data Analysis. Second Edition. Wiley - Interscience.
Aitkin, M.; Anderson, D.; Francis, B. & Hinde, J. 1989. Statistical modelling in GLIM. Clarendon Press.
Anderson, D.; Burnham, K. & White, G. 1994. AIC model selection in overdispersed capture-recapture data. Ecology 75, 1780-1793.
Breslow, N. 1990. Tests of hypotheses in overdispersed Poisson regression and other quasi-likelihood models. Journal of the American Statistical Association,85, 565-571.
Cameron, A. & Trivedi, P. 1998. Regression analysis of count data. Econometric society Monographs, No 30.
Crawley, M. 2007. The R book. Wiley.
Eberhardt, L. 1978. Appraising variability in population studies. The Journal of Wildlife Management, 207-238
Lindsey, J. 1999. On the use of corrections for overdispersion. Journal of the Royal Statistical Society: Series C (Applied Statistics), 48, 553-561.
McCullagh, P. 1983. Quasi-Likelihood Functions. The Annals of Statistics 11, 59-67.
McCullagh, P. & Nelder, J. 1989. Generalized linear models. Second Edition. Chapman & Hall.
Venables, W. & Ripley, B. 2002. Modern applied statistics with S-Plus. Fourth Edition. Springer.
ver Hoef, J. & Boveng, P. 2007. Quasi-Poisson vs. negative binomial regression: how should we model overdispersed count data? Ecology 88, 2766-2772
Williams, D. A. 1982. Extra-binomial variation in logistic linear models. Applied statistics, 144-148.
* Le contexte: la surdispersion dans le cas d'un modèle log-linéaire d'une variable distribuée selon une loi de Poisson.
* Le fond du problème qui m'intéresse: la surdispersion dans le cas d'un modèle supposé poissonnien ajusté à une table de contingence.
* Une remarque finale sur le diagnostic de la surdispersion dans le contexte "table de contingence"
EDIT: j'ai fait quelques edits pour préciser quelques points, et une nouvelle section sur le test du rapport de vraisemblance.
==========================================
Contexte: Un modèle log-linéaire d'une variable distribuée selon une loi de Poisson.
La notion de surdispersion décrit le cas où la variance de la variable réponse est supérieure à la variance nominale, i.e. celle qui est attendue sous l'hypothèse du modèle supposé pour cette variable (e.g. dans le cas Poissonnien, lorsque cette variance est supérieure à l'espérance de la variable). C'est la définition donnée par McCullagh & Nelder (1989).
Lorsque l'on ajuste un modèle log-linéaire reliant une variable poissonienne Y à un ensemble de variables explicatives X1, X2, ..., Xp, cette surdispersion peut être diagnostiquée en comparant la déviance résiduelle et le nombre de degrés de liberté résiduels du modèle. Asymptotiquement, la déviance résiduelle devrait être égale au nombre de degrés de libertés résiduels. Asymptotiquement. On *diagnostique* une surdispersion lorsque la déviance est bien supérieure au nombre de degrés de libertés. Après, il y a toute une discussion relative à "supérieure de combien?"; Lindsey (1999) indique qu'on doit la prendre en compte lorsque la déviance est plus de deux fois plus grande que le nombre de dll, ce qui est contesté -- à juste titre à mon avis -- par Ripley, cf. e.g. http://caloine.ouvaton.org/shaarli/?dGq1kA et le lien correspondant. En outre, Venables & Ripley (2002) mettent quand même fortement en garde contre une interprétation trop littérale du rapport déviance/ddl ("This can be seriously misleading"), car cette théorie est asymptotique, et ne s'applique que pour des espérances larges dans le cas Poisson. McCullagh et Nelder (1989, p. 36) indiquent d'ailleurs "In general, however, the \chi^2 approximations for the deviance are not very good even as n \rightarrow \infty".
Bon bref, ce n'est pas ce qui m'intéresse le plus ici. En fait, la notion de *diagnostic* est ici importante, car on verra plus loin dans le cas des tableaux de contingence que l'on peut passer à côté d'une surdispersion réelle d'une variable réponse avec ce simple diagnostic. Autrement dit, il peut y avoir une surdispersion de la variable réponse non diagnostiquée par le modèle (cf. deuxième partie de cette revue, j'y reviendrai).
Laissons de côté pour le moment la question du diagnostic de la surdispersion pour nous intéresser plus à celle de la solution à lui apporter. Et pour solutionner le problème de la surdispersion, il faut en connaître la cause. Qu'est ce qui cause la surdispersion dans un modèle supposé poissonnien? Cameron et Trivedi (1998) comme Agresti (2002) indiquent que la source de la surdispersion est l'unobserved heterogeneity, ce qui est presque une tautologie (basiquement, ça revient à dire qu'on a surdispersion parce que l'on a une variabilité plus grande).
Agresti (2002, p. 130) donne un exemple qui permet d'y voir plus clair: "A common cause of overdispersion is subject heterogeneity. For instance, suppose that width, weight, color, and spine condition are the four predictors that affect a female crab's number of satellites. Suppose that Y has a Poisson distribution at each fixed combination of those predictors. Our model uses width alone as predictor. Crabs having a certain width are then a mixture of crabs of various weights, colors, and spine conditions. Thus, the population of crabs having that width is a mixture of several Poisson populations, each having its own mean for the response. This heterogeneity results in an overall response distribution at that width having a greater variation than the Poisson predicts". De façon très intéressante, Agresti s'appuie sur cet exemple pour développer une solution possible au cas du crabe, en modélisant la réponse non pas comme une réponse Poissonnienne, mais comme une réponse binomiale négative.
C'est très intéressant, parce qu'à travers cet exemple, on voit les deux points de vue possibles sur la surdispersion: (i) un prédicteur linéaire mal spécifié (i.e. on est passé à côté d'une ou plusieurs variables explicatives importantes), (ii) un modèle d'erreur mal spécifié (i.e. si je considère l'exemple d'Agresti, conditionnellement à la largeur du crabe, les autres variables explicatives inconnues sont elles-même aléatoires, ce qui ajoute à la variabilité de la réponse: la réponse n'est plus une loi de Poisson, mais une distribution plus variable). Et l'on voit que ces deux causes possibles sont, dans une certaine mesure, équivalentes.
Ainsi, de nombreux auteurs distinguent ces deux causes comme source possible de surdispersion. Quand on diagnostique une surdispersion importante, il y a donc, potentiellement, deux familles de solutions possibles:
1. Modifier le prédicteur linéaire en ajoutant les variables essentielles (ici, voir les commentaires de la cellule d'appui stat de l'univ. de Cornell: http://caloine.ouvaton.org/shaarli/?voVBrw)
2. Modifier le modèle d'erreur: soit en définissant pour la variable réponse une distribution statistique qui va prendre en compte cette surdispersion, comme une distribution modèle binomial négative, soit en abandonnant la procédure d'estimation par le maximum de vraisemblance pour une procédure de type quasi-vraisemblance (i.e. en ne supposant aucune distribution précise pour la variable réponse, mais en supposant qu'il existe une distribution de la famille exponentielle à un seul paramètre pour laquelle la relation entre moyenne et variance est de la forme Var(mu) = c*mu, avec c un paramètre de surdispersion). Voir Ver Hoef et Boveng (2007) pour une comparaison des deux solutions binomialNégative vs. quasi-vraisemblance (je sors un peu du scope de ma review là).
Ces deux solutions sont également notées par de nombreux auteurs, comme Crawley (2007, qui en discute longuement), Cameron et Trivedi (1998), ou Williams (2002, dans un contexte binomial et pas poissonnien, mais c'est la même logique). Ce dernier indique explicitement: "even when all available explanatory variables have been fitted, the residual variation may be greater than can be attributed to the binomial sampling variation assumed by the model. In this event we can either seek additional explanatory variables, or postulate a source of extra-binomial random variation between observations".
Quelle approche préférer? la 1 ou la 2? Anderson et al. (1994) donnent une idée de la philosophie dans un contexte de modèles de capture-recapture: "Under the CJS [Cormack-Jolly-Seber] model theory, c=1; however, with real data we expect c>1, but we do not expect c to exceed 4 (see Eberhardt 1978). Substantially larger values of c (say, 6-10) are usually caused partly by a model structure that is inadequate, that is the fitted model does not actually represent all the explainable variation in the data. Quasi-likelihood methods of variance inflation are appropriate only after the structural adequacy of the model has been achieved". Je n'ai pas bien compris d'où Anderson et al. tirent leur valeur de 4 citée ici (je n'ai pas trouvé ça dans Eberhardt 1978; et Anderson et al. indiquent dans leur discussion que ce seuil de 5 est déduit des expériences dans le domaine de la CMR), mais ce n'est pas important ici. L'important est que si la surdispersion est modérée, on peut la prendre en compte par quasi-vraisemblance. Si elle est plus importante on est à côté de la plaque, et il faut laisser tomber le modèle construit.
Voir aussi la bonne synthèse de Carruthers et al. (2008) trouvée ici: http://www.mun.ca/biology/dschneider/b7932/B7932Final4Mar2008.pdf
En s'appuyant sur Anderson et al. et sur leur expérience, ils fixent entre 5 et 10 le "seuil" de surdispersion à partir duquel on peut raisonnablement commencer à taxer de mauvaise foi le modélisateur bourrin qui chercherait à prendre en compte la surdispersion par quasi-vraisemblance.
==========================================
Le fond du problème qui m'intéresse: le cas d'un modèle supposé poissonnien ajusté à une table de contingence.
Mon problème est le suivant. J'ai une table de contingence croisant trois facteurs A, B, et C. Il peut y avoir une trentaine de niveaux au facteur A (ce sont des espèces), une centaine au facteur B (ce sont des sites), et il y en aura toujours 2 pour le facteur C (ce sont des jeux de données).
Je note N ma variable réponse (effectifs de chaque espèce dans chaque site et chaque jeu de données, dans les cases du tableau). Des considérations théoriques m'amènent à ajuster le modèle suivant:
log(N) ~ Intercept + A + B + C + A:B + A:C + B:C (modèle 1)
Ce modèle ne sort pas du chapeau, il est le résultat d'un développement mathématique, et s'appuie sur des hypothèses que je suis prêt à poser sur mon processus biologique et sur le processus de collecte des données. Et surtout, sous ces hypothèses, ce modèle me permet d'estimer mes quantités d'intérêt (des effectifs) qui, sous les hypothèses biologiques que nous fixons sont les paramètres (Intercept+A+B+C+A:B). Je ne rentre pas dans le détail, mais disons que dans le cas présent, ce modèle a une justification théorique intéressante. Du coup, c'est ce modèle-là -- et aucun autre -- qui m'intéresse.
Lorsque je l'ajuste, j'ai une déviance, et un nombre de ddl résiduels. Du coup, il est tentant de comparer les deux pour diagnostiquer une surdispersion (Bon, dans mon cas précis, je n'ai pas surdispersion, mais dans un cas plus général, e.g. une autre application de la méthode, je pourrais avoir une grande déviance résiduelle associée à un tel modèle). Imaginons que j'ai une déviance importante associée à ce modèle. Que faire?
*À première vue*, il semble que l'on soit dans un cas particulier puisque les effectifs modélisés sont stockés dans un tableau de contingence à trois entrées. On connaît l'ensemble de facteurs qui permettraient de rendre nulle la déviance: les interactions de deuxième ordre A:B:C. Ajuster un modèle (1) avec ces interactions nous conduit à un modèle saturé, avec une déviance nulle. Donc, on sait comment faire pour supprimer cette surdispersion. Mais est-ce une démarche pertinente? En effet, si l'on optait pour cette solution:
1. ça servirait à quoi? Un modèle saturé permettrait en effet de supprimer la surdispersion, mais ne serait pas interprétable, et n'aurait aucun intérêt biologique, puisque ce qui nous intéresse, c'est une certaine combinaison des paramètres d'un modèle très précis.
2. Et surtout: il faut avoir conscience que ces interactions impliquent un grand très grand nombre de variables: Il peut y avoir une trentaine de niveaux au facteur A (ce sont des espèces), une centaine au facteur B (ce sont des sites), et il y en aura toujours 2 pour le facteur C. Donc, les interactions A:B:C pourront impliquer d'ajouter plusieurs centaines, voire des milliers, de variables à un modèle déjà compliqué!!!
Aitkin et al. (1989) discutent très en détail de ce cas particulier de la surdispersion dans le contexte des modèles log-linéaires sur tables de contingence (sa discussion se place dans un contexte binomial, mais est tout-à-fait valide dans un contexte Poissonnien): "When a contingency table is analysed, it may be found that high-order interactions are inexplicably significant, and sometimes even that no model other than the saturated model can provide a satisfactory representation of the data. This phenomenon frequently occurs with very large samples where the contingency table is not classified by some factors relevant to the response: the number of 'successes' in each cell of the table then has a /mixed/ binomial (or binomial mixture) distribution. That is, within a given cell the probability p of a success is not constant, but varies systematically with other factors which have not been recorded or included in the model. Since these factors are nt identified, the success probability p behaves like a random variable with a probability distribution".
Ces auteurs proposent alors les solutions classiques, notamment le principe de l'ajustement par quasi-vraisemblance, mais notent "This procedure is probably satisfactory for small amounts of overdispersion, but it is not a substitute for correct model specification, and it is impossible to test the goodness of fit of the model". Anderson et al. (1994) indiquent dans la discussion "In a treatise on analysis of count data, Cox and Snell (1989) assert that the simple approach, such as we have studied here, should often be adequate, as opposed to the much more arduous task of seeking an explanatory model for the overdispersion. This comment of Cox and Snell is supported by the results of Liang and McCullagh (1993), who found that causal modeling of overdispersion was clearly better than use of a single overdisperson parameter, c, in only one of five case examined".
En fait, comme on le voit ci-dessus, on n'est pas vraiment dans un cas particulier en termes de *solutions* à apporter à la surdispersion (mais voir plus bas concernant le problème du *diagnostic* de la surdispersion, qui lui, est particulier dans le cas des tableaux de contingence). Et c'est ça qui m'a pris du temps à comprendre. C'est pour cette raison que ce cas de figure, pourtant fréquent, n'est pas si souvent évoqué dans les bouquins. Autrement dit: on ne doit pas parler de l'alternative "ajustement du modèle saturé" vs. "modélisation de la surdispersion par quasi-vraisemblance" comme d'un choix cornélien dans ce contexte. On est dans le même cas de figure que les régressions Poisson classiques décrites ci-dessus, et la solution est de la même nature: le modèle (1) est le modèle d'intérêt. Si on a une déviance modérément importante, on peut corriger ça en ajustant le modèle par quasi-vraisemblance. Si on a une surdispersion très importante, c'est signe que le modèle est structurellement mauvais, et que les hypothèses sur lesquelles il repose sont probablement violées. Le modèle doit alors être évité, et une autre solution doit-être recherchée.
======================
Une remarque finale quand même, concernant le *diagnostic* de la surdispersion dans le contexte d'une analyse sur tableau de contingence.
Aitkin et al. (1989) notent que la surdispersion ne sera pas nécessairement *visible* dans un tableau de contingence. Par exemple, si une variable N est influencée essentiellement par un facteur X, mais que l'on somme ces valeurs de N sur tous les niveaux de deux autres facteurs A et B pour former une table de contingence, alors on peut parfaitement ajuster un modèle saturé N~A*B. Et comme le modèle est saturé, la déviance résiduelle sera nulle. On ne diagnostiquera pas la surdispersion occasionnée par le "collapsing" de N sur deux facteurs sans intérêt. Et pourtant, le modèle sera faux. Aitkin et al. donnent un exemple de ça, et notent "inappropriate collapsing of the table over important variable leads to distorted conclusions about the importance of other explanatory variables". La solution à ça consiste alors à inclure dès le départ un maximum de variables dépendantes dans la construction du tableau de contingence (et éviter les collapsings inappropriés). "Paradoxically, the higher the dimension of the table, the more likely we are to obtain a simple final model 'more is less'".
======================
Test du rapport de vraisemblance et Quasi-vraisemblance:
Un dernier point pour la route. Notre éditeur suggérait l'utilisation d'un test du rapport de vraisemblance entre le modèle (1) décrit plus haut et le modèle saturé comme test de qualité d'ajustement (goodness of fit ou GOF). Mais ça me posait un problème. En effet, ce test est valide sous l'hypothèse que le modèle d'erreur est correctement spécifié. Si l'on suppose un modèle de Poisson, alors la déviance va suivre, asymptotiquement un chi-deux (au passage attention: asymptotiquement, et beaucoup d'auteurs mettent en garde contre ça, cf. précédemment. McCullagh et Nelder (1989, p. 36) notent d'ailleurs que la table d'analyse de déviance ne doit être considérée, au mieux, que comme un "screening device", et qu'aucune tentative ne devrait être faite pour assigner des niveaux de significativité précis).
Pourtant, si l'on compare le modèle (1) au modèle saturé avec un test du rapport de vraisemblance, ce test peut apparaître significatif (indiquant donc un mauvais ajustement), et pourtant la solution ne serait pas de corriger ça par l'inclusion des interactions de deuxième ordre A:B:C manquantes. En effet, on peut parfaitement avoir une surdispersion modérée, mais qui apparaîtrait significative. Par exemple, imaginons que le modèle d'intérêt pour nous soit caractérisé par une déviance résiduelle de 100 sur 60 ddl. L'interaction de deuxième ordre A:B:C serait alors très hautement significative. Pourtant, cette surdispersion pourrait parfaitement être prise en compte à l'aide d'un ajustement par quasi-vraisemblance. Le paramètre de surdispersion, de phi = 100/60 = 1.67, ne serait même pas très élevé en comparaison de ce que l'on peut voir dans ce type d'étude.
On en revient toujours au même problème: le test du rapport de vraisemblance suggéré par notre éditeur indique une surdispersion, donc que le modèle de Poisson est mauvais. La solution à ce défaut peut être de modifier le modèle d'erreur en intégrant une surdispersion modérée, ou alors de modifier le prédicteur linéaire en intégrant les interactions. Oui, je sais je radote, je suis en train de redire ici ce que j'ai dit plus haut.
Mais le point intéressant ici, c'est ce qu'on trouve dans la littérature sur les tests de GOF dans le cas de la quasi-vraisemblance, à savoir pas grand chose, et le peu qu'on trouve n'est pas clair.
Par exemple, McCullagh (1983, Eq. 11) indique dans un papier très théorique que la différence de log-quasi-vraisemblance -- appelée déviance [et curieusement, pas quasi-déviance] -- entre deux modèles emboités m1 et m2 caractérisés par la même surdispersion phi, suit asymptotiquement un chi-deux multiplié par phi. Rien n'est dit dans ce papier théorique sur les conséquences pratiques de cette propriété (et sur le fait que de façon générale, phi ne sera jamais connue, donc sera estimée, et donc sur les conséquences de cette estimation sur la disttribution suivie par D/phi, où D est la déviance).
Mais Breslow (1990), dans l'intro de son article, cite la première édition du McCullagh et Nelder (de 1983) qui reporte cette dernière propriété, et indique en s'appuyant dessus que lorsque l'on cherche à comparer deux modèles, le test du (quasi) rapport de vraisemblance est disponible. Dans des simus, cet auteur utilise ce test en calculant pour chaque modèle une scaled deviance (deviance/phi, avec phi estimé classiquement à partir des données). Il en conclut que si la surdispersion est modérée, un tel test reste acceptable.
Pourtant, de façon intéressante, la propriété notée par McCullagh (1983) disparaît du McCullagh & Nelder, édition de 1989. Et cette dernière édition ne cite McCullagh 1983 que noyée au milieu d'autres refs en la référençant comme traitant de l'efficacité, de l'optimalité et de la robustesse (ce qui est vrai). Pas un mot sur le test... Cherchent-ils à éviter de cautionner cette pratique? Considèrent-ils ce résultat comme peu intéressant?
======================
Donc, en résumé, dans le contexte de la modélisation par un modèle log-linéaire d'une table de contingence:
* Le diagnostic de la surdispersion ne peut se faire que si l'on dispose des variables pertinentes comme entrées dans ce tableau (i.e. éviter de fusionner des dimensions pertinentes du tableau avant modélisation)
* La philosophie sous-tendant la correction des problèmes causés par la surdispersion est la même que dans les autres modèles log-linéaires: (i) si la surdispersion n'est pas très importante et que l'on ne tient pas à modifier la structure du modèle (soit parce qu'on n'a pas les variables qui en sont la cause, soit parce qu'on tient à cette structure pour d'autres raisons, notamment théoriques), on modifie le modèle d'erreur pour la réponse (e.g. ajustement par quasi-vraisemblance, définition d'un modèle binomial négatif); (ii) si la surdispersion est trop importante, le modèle est mauvais, il faut l'abandonner.
======================
Références:
Agresti, A. 2002. Categorical Data Analysis. Second Edition. Wiley - Interscience.
Aitkin, M.; Anderson, D.; Francis, B. & Hinde, J. 1989. Statistical modelling in GLIM. Clarendon Press.
Anderson, D.; Burnham, K. & White, G. 1994. AIC model selection in overdispersed capture-recapture data. Ecology 75, 1780-1793.
Breslow, N. 1990. Tests of hypotheses in overdispersed Poisson regression and other quasi-likelihood models. Journal of the American Statistical Association,85, 565-571.
Cameron, A. & Trivedi, P. 1998. Regression analysis of count data. Econometric society Monographs, No 30.
Crawley, M. 2007. The R book. Wiley.
Eberhardt, L. 1978. Appraising variability in population studies. The Journal of Wildlife Management, 207-238
Lindsey, J. 1999. On the use of corrections for overdispersion. Journal of the Royal Statistical Society: Series C (Applied Statistics), 48, 553-561.
McCullagh, P. 1983. Quasi-Likelihood Functions. The Annals of Statistics 11, 59-67.
McCullagh, P. & Nelder, J. 1989. Generalized linear models. Second Edition. Chapman & Hall.
Venables, W. & Ripley, B. 2002. Modern applied statistics with S-Plus. Fourth Edition. Springer.
ver Hoef, J. & Boveng, P. 2007. Quasi-Poisson vs. negative binomial regression: how should we model overdispersed count data? Ecology 88, 2766-2772
Williams, D. A. 1982. Extra-binomial variation in logistic linear models. Applied statistics, 144-148.
J'aime bien la première réponse de Ripley, au sujet des tests sur la surdispersion: "There are, but like formal tests for outliers I would not advise using them, as you may get misleading inferences before they are significant, and they can reject when the inferences from the small model are perfectly adequate".
The moment estimator \phi of over-dispersion (calculé par le rapport Chi-2/ddl) gives you an indication of the likely effects on your inferences: e.g. estimated standard errors are proportional to \sqrt(\phi). having standard errors which need inflating by 40% seems to indicate that the rule you quote is too optimistic (even when its estimator is reliable).
Edit: remarque de Ripley concernant l'article de Lindsey: "And I can add 'it is helpful to know whose authority is being invoked', since (e.g.) some authors are not at all careful." Dommage que la critique ne soit pas plus détaillée.
The moment estimator \phi of over-dispersion (calculé par le rapport Chi-2/ddl) gives you an indication of the likely effects on your inferences: e.g. estimated standard errors are proportional to \sqrt(\phi). having standard errors which need inflating by 40% seems to indicate that the rule you quote is too optimistic (even when its estimator is reliable).
Edit: remarque de Ripley concernant l'article de Lindsey: "And I can add 'it is helpful to know whose authority is being invoked', since (e.g.) some authors are not at all careful." Dommage que la critique ne soit pas plus détaillée.
Une autre démonstration simple et intéressante du fait que quand on a une variable aléatoire X distribuée selon une loi binomiale de paramètres N et p, dont le paramètre N est tiré d'une loi de Poisson de paramètre lambda, la variable aléatoire X est tirée d'une loi de Poisson de paramètre lambda*p. Cette démonstration utilise l'identité
\exp(x) = \sum_{n=0}^\infty (x^n)/(n!)
\exp(x) = \sum_{n=0}^\infty (x^n)/(n!)
Théorème de Rao-Rubin qui permet la caractérisation de la distribution de Poisson. Récupéré.
Une conséquence intéressante: quand X suit une distribution de Poisson de paramètre lambda, et que Y suit une binomiale de paramètre n, p, et que n est une réalisation de X, alors Y est une loi de Poisson de paramètre lambda*p.
Une conséquence intéressante: quand X suit une distribution de Poisson de paramètre lambda, et que Y suit une binomiale de paramètre n, p, et que n est une réalisation de X, alors Y est une loi de Poisson de paramètre lambda*p.