Article intéressant. Voir en particulier la figure 4: un article publié d'abord soumis à ArXiv est beaucoup, *beaucoup* plus cité qu'un article directement publié. Les physiciens se font plus chier: la publication permet de "référencer" une étude bien proprement, mais en pratique, toute la discussion scientifique se passe autour des articles sur ArXiv...

Belle vidéo

Une discussion intéressante sur la modélisation des données de comptage.

Un nouveau type de données apparaît: le "medium" data: 215 Go de données!!! Putain, on en est pas là encore en écologie.
Chose intéressante, ils ont mis le code source sur github. Ça peut être intéressant de voir comment ils gèrent cette quantité de données. J'ai regardé un peu, c'est surtout du PostgreSQL, avec des requêtes depuis R. Je n'ai pas l'impression qu'il y ait des outils propres à cette taille de donnée (à part postgresql, qui certes est rarement utilisé pour les "small" data, ou les nano data, je sais pas comment on doit appeler nos jeux de données de 50 lignes, mais qui est quand même plus fréquemment utilisé pour les BD plus grandes comme les données réseau, etc.).

Intéressant. Je résume. Un écologue reviewant un papier devra:

1. Be Honest about What You Know and What You Don’t
2. Old Doesn’t Mean Bad (putain celui là on l'aura entendu! "cette méthode n'est plus utilisée, aujourd'hui ya mieux!")
3. New Doesn’t Mean Good (Daniel Chessel avait souligné dans son HDR les points 2 et 3 comme réaction fréquente de la part des biologistes face aux nouvelles méthodes).
4. Understand Assumptions (Suggesting unnecessarily complicated methodology makes things harder for the authors, the editors, for you as the reviewer (...), and finally, (...) the readers.)
5. Match the Method to the Assumptions
6. There is Rarely a Single Correct Method
7. p-values are not the Work of the Devil
8. Beware of Uninformed Pronunciamentos* (le comportement qui consiste à dire qu'une méthode ou une pratique -- e.g. calculer des P-values, s'intéresser à l'autocorrélation spatiale, etc. -- est toujours mauvaise, quel que soit le contexte).
9. Recommend a Statistician Look at the Paper (i.e. indiquer son incompétence en stats, et demander un autre reviewer sur ce thème, à relier au point 1).
10. Read up on the Statistical Methods

J'avais déjà vu passer ça...

Ben putain, c'est bon à savoir:
some-command <(another-command)
va exécuter "another-command" -- il faut que cette autre commande écrive le résultat dans stdout -- et va passer le résultat à some-command (avec some-command une commande qui prend un nom de fichier comme argument). Donc pas besoin de stocker dans les fichiers, on passe directement le résultat du l'une à l'autre. Comme un pipe, sauf que (i) on peut substituer autant de commande qu'on veut, par exemple:
diff <(curl http://somesite/file1) <(curl http://somesite/file2)
va substituer les deux curl et les passer à diff; et (ii) les commandes substituées sont lancées en parallèle et pas à la suite.

J'adore ce genre de conneries...

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.

Ah ben tiens, je connaissais pas: pour se créer des adresses bidons, jetables, et personnalisées.
Peut être pratique, je me garde ça sous le coude.