Une introduction à la distribution binomiale négative

L’analyse des phénomènes discrets constitue l’un des piliers méthodologiques les plus cruciaux des sciences contemporaines, qu’il s’agisse de la psychologie quantitative, de l’épidémiologie, de l’économétrie ou des sciences cognitives. Très souvent, les chercheurs sont confrontés à des données de comptage représentant la fréquence d’apparition d’un événement au cours d’une période donnée, ou le nombre d’essais nécessaires avant d’atteindre un objectif prédéterminé. Si la distribution de Poisson a historiquement occupé une place centrale dans la modélisation de tels événements en raison de son élégance et de sa simplicité mathématique, son postulat fondamental d’équidispersion — imposant l’égalité stricte entre la moyenne et la variance — s’avère fréquemment violé dans les données empiriques réelles. Cette violation, caractérisée par une variabilité observée largement supérieure à celle prédite par le modèle théorique, porte le nom de surdispersion.

Face aux limites manifestes du modèle de Poisson, la distribution binomiale négative s’impose comme une alternative statistique majeure et robuste. Issue initialement de la formalisation des temps d’attente dans les processus de Bernoulli répétés, cette distribution a connu une extension théorique remarquable pour devenir le modèle par excellence de l’hétérogénéité non observée à travers les mélanges continus de lois. Elle permet de concilier la nature discrète des variables de comptage avec la flexibilité requise pour modéliser des populations composées d’individus présentant des propensions intrinsèques hétérogènes face à l’événement étudié.

Cet article propose une immersion exhaustive au sein de la distribution binomiale négative, depuis ses fondements axiomatiques ancrés dans la théorie des probabilités jusqu’à ses applications avancées dans la modélisation statistique moderne. À travers une progression méthodique, nous explorerons les déductions analytiques de sa fonction de masse, ses propriétés distributionnelles et ses moments d’ordre supérieur, avant d’aborder les défis computationnels liés à l’estimation des paramètres, aux tests de diagnostic et aux architectures complexes telles que les modèles avec excès de zéros. Ce parcours permettra au lecteur de saisir toute la portée théorique et la pertinence pratique de cet outil d’analyse indispensable.

1. Fondements conceptuels et émergence à partir des épreuves de Bernoulli

1.1 Définition axiomatique d’un processus de Bernoulli

À la racine de la distribution binomiale négative se trouve l’une des structures stochastiques les plus élémentaires et pourtant les plus fécondes des mathématiques : l’épreuve de Bernoulli. Axiomatiquement, une épreuve de Bernoulli est une expérience aléatoire ne comportant que deux issues possibles, mutuellement exclusives et conjointement exhaustives, conventionnellement désignées sous les vocables de « succès » et d’« échec ». L’espace fondamental associé à cette épreuve unitaire se résume ainsi à l’ensemble binaire {0, 1}, où le chiffre 1 symbolise la réalisation de l’événement cible et 0 son absence. Cette dichotomie simplificatrice constitue l’atome probabiliste sur lequel reposent des architectures bien plus élaborées.

Le déploiement de ces épreuves sous forme de processus stochastique séquentiel requiert l’adhésion stricte à deux postulats d’invariance. Le premier est l’invariance temporelle de la probabilité de succès, traditionnellement notée p (avec 0 < p < 1), ce qui implique que le mécanisme sous-jacent générant l’événement demeure rigoureusement stationnaire au fil des itérations. Le second postulat est la condition d’indépendance stochastique mutuelle entre les essais successifs. Selon cette hypothèse, l’issue de l’essai de rang n ne conserve aucune mémoire des réalisations antérieures et n’exerce aucune influence causale sur l’essai de rang n + 1. La probabilité conditionnelle d’obtenir un succès sachant l’historique des tirages passés demeure donc invariablement égale à p.

C’est précisément dans la gestion des conditions d’arrêt que s’opère la césure épistémologique fondamentale entre les grandes lois discrètes. Dans un schéma binomial classique, le chercheur fixe par avance une variable déterministe, à savoir le nombre total d’essais n, et s’intéresse à la variable aléatoire mesurant le nombre de succès cumulés. À l’inverse, dans le cadre de l’échantillonnage inverse ou séquentiel menant à la loi binomiale négative, le statut des grandeurs est inversé : le nombre de succès cibles, noté r, est une constante prédéterminée, tandis que le nombre d’essais nécessaires pour atteindre ce quota devient la variable aléatoire régie par les lois du hasard.

1.2 Mécanisme de passage de la loi géométrique à la loi binomiale négative

Pour appréhender le mécanisme de constitution de la loi binomiale négative, il convient d’examiner son précurseur direct : la loi géométrique. Cette dernière modélise le temps d’attente discret, exprimé en nombre d’échecs survenant avant l’apparition d’un tout premier succès. Si l’on désigne par X cette variable, la probabilité d’observer exactement k échecs avant le premier succès s’obtient simplement par l’application de la règle de multiplication des probabilités indépendantes, donnant (1 − p)^kp. La loi géométrique illustre de manière épurée la notion d’échantillonnage d’arrêt à seuil unitaire.

La distribution binomiale négative émerge alors naturellement comme la généralisation directe de ce processus d’attente lorsqu’on élargit le critère d’arrêt à un nombre entier r de succès, où r ≥ 1. Dès lors, la variable aléatoire représentant le nombre total d’échecs cumulés avant d’obtenir le r-ième succès peut être rigoureusement formalisée comme la somme de r variables aléatoires géométriques indépendantes et identiquement distribuées. Chacune de ces variables géométriques capture le nombre d’échecs intercalés entre deux succès consécutifs. Cette propriété de convolution confère à la distribution une assise théorique solide, démontrant que l’accumulation d’échecs sous condition d’un quota de réussites obéit à une structure additive stable.

Cette approche par échantillonnage inverse modifie profondément les propriétés statistiques par rapport à un échantillonnage à taille fixe. Alors que dans un protocole d’échantillonnage fixe, le support de la variable aléatoire est borné par le nombre maximal d’essais n, le protocole inverse ouvre l’espace des états à un support dénombrable mais potentiellement infini. En effet, il est théoriquement possible, bien que d’une probabilité asymptotiquement nulle, d’observer une suite indéfinie d’échecs sans jamais parvenir à réunir les r succès exigés pour clore l’expérience. Cette ouverture du support engendre des distributions à queue lourde, particulièrement adaptées pour capter des dynamiques d’attente prolongées.

1.3 Pertinence épistémologique dans les sciences comportementales

Dans l’univers des sciences comportementales, de la psychologie clinique et des neurosciences cognitives, la modélisation des processus d’apprentissage et de régulation cognitive se heurte souvent à la nature discrète et asymétrique des comportements humains. De nombreux paradigmes expérimentaux n’imposent pas un nombre arbitraire d’essais, mais évaluent au contraire la persévérance d’un sujet jusqu’à ce qu’il atteigne un critère de maîtrise ou de performance clinique. C’est le cas, par exemple, des protocoles d’apprentissage par essais et erreurs où l’expérimentateur comptabilise les erreurs commises avant d’observer une série ininterrompue de réponses correctes.

Sur le plan psychologique, la binomiale négative capture de manière élégante des construits fondamentaux tels que la tolérance à la frustration, l’inflexibilité cognitive ou la fatigabilité. Un individu présentant un déficit au niveau des fonctions exécutives aura tendance à répéter des choix non optimaux, accumulant ainsi un nombre d’échecs disproportionné avant d’atteindre le niveau de réussite imposé par le protocole expérimental. La distribution permet d’assigner une probabilité formelle à ces profils atypiques sans les exclure en tant que simples valeurs aberrantes, offrant une modélisation mathématique respectueuse de la variabilité clinique.

L’utilisation de distributions continues, telles que la distribution gaussienne ou log-normale, pour analyser ces variables de comptage discrètes engendre des distorsions statistiques majeures. D’une part, les données de comptage sont intrinsèquement non négatives et restreintes aux entiers naturels, rendant les densités continues conceptuellement inadaptées, en particulier au voisinage de zéro. D’autre part, la variance de ces processus comportementaux croît presque systématiquement avec la moyenne, brisant l’hypothèse d’homoscédasticité requise par les modèles linéaires classiques. La loi binomiale négative apporte une solution rigoureuse à ces apories méthodologiques en fournissant un cadre génératif fidèle à la réalité empirique du laboratoire.

2. Définition mathématique et formulations de la loi binomiale négative

2.1 Dérivation analytique de la fonction de masse de probabilité

Pour établir formellement la fonction de masse de probabilité d’une variable aléatoire X suivant une loi binomiale négative, notée traditionnellement NB(r, p), considérons l’événement précis suivant : observer exactement k échecs avant l’obtention du r-ième succès. La réalisation de cet événement impose une contrainte déterministe incontournable sur la chronologie des tirages : le dernier essai — c’est-à-dire l’essai numéro (k + r) — doit impérativement être un succès. Si tel n’était pas le cas, le processus se serait soit arrêté prématurément, soit poursuivi au-delà.

Dès lors, les (r − 1) succès restants et les k échecs peuvent survenir dans n’importe quel ordre le long des (k + r − 1) premiers essais. Le nombre d’ordonnancements distincts permettant de distribuer ces k échecs parmi les (k + r − 1) positions disponibles est régi par le coefficient binomial usuel :

C(k + r − 1, k) = (k + r − 1)! / [k! × (r − 1)!]

En vertu de l’indépendance stochastique des épreuves, chaque séquence élémentaire particulière constituée de k échecs et de r succès se voit attribuer la probabilité multiplicative p^r(1 − p)^k. En combinant le terme combinatoire et la probabilité élémentaire, nous aboutissons à la fonction de masse classique :

P(X = k) = [ (k + r − 1)! / (k! (r − 1)!) ] × p^r(1 − p)^k, pour tout k ∈ {0, 1, 2, …}

Pour démontrer que la somme de ces probabilités sur l’ensemble du support dénombrable est rigoureusement égale à l’unité, on fait appel au développement en série de Taylor de la fonction (1 − u)^−r au voisinage de zéro. En substituant u par (1 − p), on vérifie analytiquement que la sommation infinie reconstitue l’expression p^r × (1 − (1 − p))^−r = p^r × p^−r = 1, confirmant la validité axiomatique de la loi.

2.2 Formulation alternative : nombre total d’essais nécessaires

Dans la littérature statistique et les sciences appliquées, il existe une dualité de paramétrisation qu’il importe de clarifier afin d’éviter toute confusion conceptuelle. Alors que la première formulation s’intéresse exclusivement au nombre d’échecs k, la formulation alternative définit la variable aléatoire d’intérêt, notée Y, comme le nombre total d’essais nécessaires pour parvenir à recueillir les r succès cibles. Il s’ensuit la relation arithmétique fondamentale Y = X + r.

Cette transformation modifie le support de la distribution ainsi que l’expression analytique de sa fonction de masse. Le support de Y ne débute plus à zéro, mais est translaté pour s’étendre sur l’intervalle discret {r, r + 1, r + 2, …}. En procédant au changement de variable adéquat au sein de la fonction de masse originelle, nous obtenons pour tout entier n ≥ r :

P(Y = n) = [ (n − 1)! / ((r − 1)! (n − r)!) ] × p^r(1 − p)^{n − r}

Le choix entre ces deux formulations dépend des finalités de l’investigateur et de la nature de l’appareillage expérimental. Les ingénieurs de fiabilité et les psychométriciens privilégient fréquemment la seconde formulation lorsqu’ils cherchent à quantifier la durée totale d’un test ou le coût énergétique global d’une tâche. En revanche, les économètres et les biostatisticiens recourent presque exclusivement à la première formulation (axée sur k) car celle-ci permet un traitement direct des données de comptage commençant à zéro, rendant possible la comparaison immédiate avec d’autres lois de comptage comme la loi de Poisson.

2.3 Extension aux valeurs non entières du paramètre r

L’une des avancées les plus remarquables de l’analyse stochastique réside dans l’extension de la loi binomiale négative au-delà de la contrainte d’un paramètre r entier. En s’affranchissant du cadre rigide des tirages séquentiels discrets, il devient possible d’admettre que r soit un nombre réel strictement positif quelconque (r ∈ ℜ⁺). Cette généralisation analytique est rendue opérationnelle grâce à l’introduction de la fonction Gamma d’Euler, définie pour tout complexe ou réel de partie réelle strictement positive par l’intégrale :

Γ(z) = ∫₀^∞ t^{z − 1} e^−t dt

Sachant que la fonction Gamma vérifie l’identité fondamentale Γ(z + 1) = zΓ(z), elle permet d’étendre la notion de factorielle aux grandeurs continues selon l’égalité n! = Γ(n + 1). Le coefficient combinatoire s’écrit alors au moyen de rapports de fonctions Gamma :

P(X = k) = [ Γ(k + r) / (k! × Γ(r)) ] × p^r(1 − p)^k

Cette formulation analytique justifie pleinement l’appellation de « binomiale négative ». En effet, le terme combinatoire généralisé peut être réécrit à l’aide du symbole factoriel décroissant étendu aux coefficients négatifs :

(−1)^k × [ (−r)(−r − 1)…(−r − k + 1) / k! ] = C(−r, k)

Cette extension offre une continuité mathématique essentielle. Le paramètre r perd son statut de « quota de succès » pour se muer en un paramètre de forme continu, gouvernant la dispersion de la distribution, conférant à la loi une souplesse d’ajustement empirique remarquable pour des phénomènes où les quotas discrets n’ont pas de sens physique direct.

3. Propriétés statistiques fondamentales et moments

3.1 Espérance mathématique et interprétation centrale

L’étude des moments d’une loi de probabilité permet de caractériser sa géométrie distributionnelle et son comportement asymptotique. Pour la loi binomiale négative sous sa forme centrée sur le nombre d’échecs X ~ NB(r, p), l’espérance mathématique, ou moment ordinaire d’ordre 1, s’obtient en exploitant la linéarité de l’opérateur d’espérance appliquée à la somme des r variables géométriques indépendantes qui la composent. Sachant que pour une épreuve géométrique unitaire, le nombre moyen d’échecs avant un succès est égal à (1 − p) / p, nous obtenons immédiatement par additivité :

E(X) = μ = r(1 − p) / p

Cette relation offre une lecture intuitive lumineuse : l’espérance est le produit direct du nombre de succès visés par le ratio de cotes (odds) d’échec unitaire (1 − p) / p. Si l’événement cible est hautement probable (p proche de 1), le dénominateur compense largement le numérateur, entraînant une espérance tendant vers zéro, ce qui signifie que l’accumulation d’échecs est extrêmement brève. Inversement, si le succès s’avère un événement rare (p tendant vers 0), l’espérance croît de façon hyperbolique vers l’infini, témoignant de la nécessité d’endurer un nombre considérable d’échecs pour accumuler les r réussites.

La sensibilité de l’espérance mathématique aux variations marginales de p peut être mesurée par sa dérivée partielle par rapport à ce paramètre : dμ / dp = −r / p². Cette dérivée, systématiquement négative sur l’intervalle ouvert ]0, 1[, démontre le caractère strictement décroissant de l’espérance vis-à-vis de l’accessibilité du succès. L’accélération de cette décroissance à l’approche de zéro illustre la grande réactivité du modèle face aux fluctuations d’efficacité dans les régimes d’événements rares.

3.2 Variance et phénomène de surdispersion

La variance théorique de la distribution binomiale négative constitue sans doute sa propriété la plus séduisante pour les praticiens des données statistiques. Toujours en exploitant la propriété d’indépendance mutuelle unissant les r réalisations géométriques sous-jacentes, la variance s’établit comme la somme directe des variances géométriques individuelles, soit :

Var(X) = σ² = r(1 − p) / p²

Il est particulièrement instructif de réécrire cette variance en fonction directe de l’espérance mathématique μ précédemment calculée. En procédant à une décomposition arithmétique, nous pouvons exprimer la variance sous la forme fondamentale :

Var(X) = [ r(1 − p) / p ] × (1 / p) = μ × (1 / p) = μ + (1 − p)μ / p = μ + μ² / r

Cette décomposition met en lumière le phénomène central de surdispersion. Étant donné que la probabilité p est strictement inférieure à 1 dans tout cadre stochastique non dégénéré, le facteur d’échelle 1 / p est strictement supérieur à 1, ce qui conduit inévitablement à la relation d’inégalité :

Var(X) > E(X)

Cette supériorité systématique de la variance sur la moyenne contraste radicalement avec la distribution de Poisson, dont l’identité constitutive est l’équidispersion (Var(X) = E(X) = λ). La binomiale négative intègre ainsi de manière endogène un terme quadratique additionnel (μ² / r), désigné sous le nom d’effet de dispersion. Lorsque r prend des valeurs infiniment grandes, ce terme correctif s’annule asymptotiquement, restaurant l’équidispersion poissonnienne. À l’inverse, plus r est faible, plus la variance explose par rapport à la moyenne, conférant au modèle la capacité de capturer des hétérogénéités empiriques considérables.

3.3 Moments d’ordre supérieur : Asymétrie et aplatissement

Pour parfaire la caractérisation géométrique de la densité discrète, il convient de dériver ses moments centrés d’ordre supérieur, spécifiquement les coefficients d’asymétrie (skewness) et d’aplatissement (kurtosis). Le calcul formel du moment centré d’ordre 3 conduit à l’expression de l’asymétrie normalisée de Fisher γ₁ :

γ₁ = E[(X − μ)³] / σ³ = (2 − p) / √[r(1 − p)]

Puisque la probabilité p appartient à l’intervalle ]0, 1[, le numérateur (2 − p) demeure strictement positif en toute circonstance. Il en résulte que le coefficient d’asymétrie de la distribution binomiale négative est invariablement positif (γ₁ > 0). La distribution présente donc une queue étirée de manière caractéristique vers la droite (right-skewed), reflétant une concentration des masses de probabilité sur les petites valeurs d’échecs associée à une décroissance lente vers les valeurs élevées.

L’évaluation du moment centré d’ordre 4 permet de quantifier le coefficient d’aplatissement de Pearson γ₂, ou excès d’aplatissement, défini par rapport à la référence de la distribution normale :

γ₂ = { E[(X − μ)⁴] / σ⁴ } − 3 = [ 6 + (1 − p)² / (r(1 − p)) ] / r = [ 6 / r ] + [ p² / (r(1 − p)) ]

Cet excès d’aplatissement étant strictement supérieur à zéro pour toute combinaison finie des paramètres, la loi binomiale négative est formellement leptokurtique. Ses queues de distribution sont nettement plus épaisses que celles d’une loi normale de même variance, conférant une probabilité non négligeable à des réalisations extrêmes. Néanmoins, lorsque le paramètre de forme r tend vers l’infini tout en maintenant la probabilité p fixe, tant l’asymétrie γ₁ que l’excès de kurtosis γ₂ convergent asymptotiquement vers zéro à un taux d’ordre O(r^−1/2) et O(r⁻¹). En vertu du théorème central limite, la distribution standardisée (X − μ) / σ converge alors en loi vers la distribution normale standard N(0, 1).

4. Comparaison critique avec les lois de probabilité connexes

4.1 Divergences structurelles avec la distribution binomiale classique

Bien que la distribution binomiale négative partage avec la distribution binomiale classique une racine étymologique commune et une filiation directe avec le processus de Bernoulli, ces deux structures probabilistes incarnent des logiques opératoires diamétralement opposées. La première divergence réside dans la nature du cadre d’échantillonnage. La loi binomiale classique découle d’un protocole à horizon fini où le nombre total d’expériences n est fixé de manière exogène, le nombre de succès observés constituant l’unique aléa. À l’opposé, la binomiale négative fige l’objectif de succès r et abandonne l’horizon expérimental à la contingence stochastique.

Cette divergence méthodologique a des répercussions directes sur l’espace d’échantillonnage. Le support mathématique de la distribution binomiale est rigoureusement borné par l’entier déterministe n, son ensemble d’états étant {0, 1, …, n}. À l’inverse, le support de la binomiale négative est infini dénombrable, transcendant toute limite supérieure prédéterminée. Cette propriété rend la binomiale classique intrinsèquement incapable de modéliser des temps d’attente ou des processus d’accumulation non contraints dans le temps.

Enfin, leur comportement vis-à-vis de la dispersion statistique met en évidence une rupture structurelle majeure. La variance d’une loi binomiale classique s’exprime par la relation np(1 − p). Sachant que le facteur multiplicatif (1 − p) est strictement inférieur à 1, il s’ensuit que la variance d’une distribution binomiale est systématiquement inférieure à son espérance mathématique np, un phénomène désigné sous le terme de sous-dispersion. La loi binomiale négative résout le problème opposé : par sa construction même, elle engendre une surdispersion structurelle, se positionnant à l’autre extrémité du spectre de variabilité par rapport à la binomiale usuelle.

4.2 Convergence vers la distribution de Poisson

La distribution de Poisson occupe une position charnière au sein de la théorie des probabilités en tant que loi limite universelle pour le comptage d’événements rares et indépendants. Il existe un pont analytique d’une remarquable élégance reliant la distribution binomiale négative à celle de Poisson. Ce passage à la limite s’opère lorsque le paramètre de succès cible r tend vers l’infini, tandis que la probabilité p tend vers 1, de telle manière que le produit r(1 − p) demeure constant et converge vers une valeur λ strictement positive :

lim_{r → ∞, p → 1} [ r(1 − p) ] = λ

Sous ces conditions de convergence, examinons le comportement de l’espérance et de la variance de la loi binomiale négative. L’espérance μ = r(1 − p) / p converge directement vers λ / 1 = λ. Simultanément, le terme quadratique de surdispersion présent dans la décomposition de la variance s’évanouit :

lim_{r → ∞} [ μ² / r ] = λ² / ∞ = 0

En conséquence, la variance Var(X) converge elle aussi vers λ, restaurant à la limite stricte l’état d’équidispersion caractéristique de la loi de Poisson. Sur le plan analytique, la fonction de masse de la binomiale négative converge ponctuellement vers e^−λλ^k / k! pour chaque entier k. Ce résultat fournit un critère de décision statistique essentiel : lorsque l’estimation empirique de r aboutit à des grandeurs extrêmement élevées, le recours à la binomiale négative cesse d’être justifié, et le principe de parcimonie scientifique enjoint alors d’adopter le modèle de Poisson, moins complexe en termes de paramètres libres.

4.3 Relations avec la loi de Pascal et la loi de Pólya

L’histoire des sciences mathématiques a souvent vu se superposer des terminologies distinctes pour désigner des objets probabilistes similaires ou étroitement intriqués. C’est tout particulièrement le cas de la relation entre la loi binomiale négative et la distribution de Pascal. Historiquement, le vocable de distribution de Pascal a été réservé par plusieurs auteurs, notamment dans la tradition probabiliste francophone, au cas spécifique où le paramètre d’arrêt r est strictement restreint à un entier positif (r ∈ ℕ^*), correspondant à la stricte réalité physique du cumul de succès discrets. L’appellation de binomiale négative est alors réservée au cadre élargi dans lequel r est libéré de cette contrainte arithmétique pour devenir un réel positif continu.

Une généralisation encore plus puissante relie la binomiale négative aux processus d’urnes théorisés par George Pólya. Dans le modèle d’urne de Pólya avec contagion, chaque fois qu’une boule d’une couleur déterminée est tirée au sort, elle est remise dans l’urne accompagnée d’un nombre déterminé de boules de la même couleur. Ce mécanisme crée une corrélation temporelle positive : la survenue d’un événement accroît la probabilité que ce même événement se reproduise ultérieurement, générant un phénomène d’auto-renforcement.

Cette dynamique de contagion capture avec acuité les phénomènes d’apprentissage, d’accoutumance ou d’amplification comportementale en psychologie. Lorsqu’un individu subit un échec, sa probabilité d’échec subséquente peut se trouver temporairement accrue par des mécanismes d’anxiété de performance ou d’impuissance acquise. La loi de Pólya fournit une base théorique démontrant que l’intégration de tels effets de rétroaction positive au sein d’un processus de comptage engendre précisément une distribution marginale équivalente à la binomiale négative. Cela prouve que la surdispersion n’est pas un simple artifice mathématique, mais l’empreinte tangible de dynamiques d’interaction temporelle au sein du système observé.

5. Modélisation de la surdispersion et mélange de lois

5.1 La binomiale négative comme mélange continu Poisson-Gamma

Au-delà de son émergence axiomatique à partir des processus de Bernoulli inverses, la distribution binomiale négative admet une seconde genèse théorique, extraordinairement féconde pour les sciences empiriques : celle d’un mélange continu de distributions. Considérons une population d’individus où chacun émet des événements discrets conformément à un processus de Poisson classique de taux d’occurrence λ. Si la population était strictement homogène, la distribution agrégée des événements suivrait fidèlement une loi de Poisson. Or, dans toute population biologique, clinique ou sociale, le taux de base λ fluctue d’un sujet à l’autre en fonction de caractéristiques inobservées ou non mesurables.

Pour modéliser mathématiquement cette variabilité inter-individuelle, faisons l’hypothèse que le paramètre λ est lui-même une variable aléatoire continue suivant une loi Gamma, caractérisée par un paramètre de forme α et un paramètre d’échelle β. La fonction de masse marginale de la variable de comptage observée Y s’obtient alors en intégrant la distribution conditionnelle de Poisson sur l’ensemble du domaine de support de la distribution a priori de λ :

P(Y = y) = ∫₀^∞ P(Y = y | λ) × f_Gamma(λ; α, β) dλ

En remplaçant les fonctions par leurs densités explicites respectives, cette intégrale prend la forme suivante :

P(Y = y) = ∫₀^∞ [ (λ^y e^−λ / y!) ] × [ (β^α / Γ(α)) λ^{α − 1} e^−βλ ] dλ

En regroupant les puissances de λ et les exponentielles, on extrait les constantes vis-à-vis de la variable d’intégration :

P(Y = y) = [ β^α / (y! Γ(α)) ] ∫₀^∞ λ^{y + α − 1} e^{−(β + 1)λ} dλ

L’intégrale restante se résout analytiquement par identification avec la fonction Gamma standard, valant Γ(y + α) / (β + 1)^{y + α}. En réinjectant ce résultat, nous parvenons à la forme close :

P(Y = y) = [ Γ(y + α) / (y! Γ(α)) ] × [ β / (β + 1) ]^α × [ 1 / (β + 1) ]^y

Cette expression est rigoureusement identique à la fonction de masse de la distribution binomiale négative sous son extension continue, en posant formellement le paramètre r = α et la probabilité p = β / (β + 1). Ce théorème de mélange continu, initialement mis en lumière par Greenwood et Yule en 1920 dans leurs travaux pionniers sur l’accidentologie, éclaire d’une lumière nouvelle la structure de la loi : la binomiale négative est l’empreinte mathématique laissée par l’hétérogénéité individuelle non observée — souvent qualifiée de « fragilité » (frailty) — se superposant à un processus stochastique sous-jacent de Poisson.

5.2 Typologies de la surdispersion dans les données empiriques

L’existence d’une variance supérieure à la moyenne dans un jeu de données empiriques constitue la signature universelle de la surdispersion. Toutefois, sur le plan méthodologique, il est impératif de distinguer deux formes ontologiques de surdispersion : la surdispersion apparente et la surdispersion réelle. La surdispersion apparente n’est pas intrinsèque au processus générateur, mais résulte d’erreurs de spécification commises par le statisticien. Elle survient typiquement lors de l’omission de covariables explicatives de premier ordre, de la mauvaise transformation d’un prédicteur continu ou de l’inadéquation de la fonction de lien reliant le prédicteur linéaire à la réponse.

À l’opposé, la surdispersion réelle ou intrinsèque procède de la nature stochastique profonde du phénomène étudié. Elle se manifeste selon deux vecteurs fondamentaux : la contagion temporelle et la stratification populationnelle. La contagion renvoie aux situations où l’occurrence d’un premier événement modifie organiquement le taux d’occurrence des événements futurs pour ce même individu, invalidant de facto l’hypothèse de stationnarité poissonnienne. La stratification correspond quant à elle à l’agrégation non identifiée de sous-groupes cliniques ou démographiques présentant des niveaux de vulnérabilité intrinsèquement divergents face au risque étudié.

De surcroît, les effets de grappe spatiaux ou organisationnels constituent une source récurrente de surdispersion en milieu empirique. Lorsque des sujets sont regroupés au sein d’unités écologiques partagées — telles que des classes scolaires, des services hospitaliers ou des aires géographiques — ils partagent des déterminants contextuels communs non mesurés. Cette dépendance intra-grappe viole l’indépendance conditionnelle des résidus et induit un étalement de la distribution des comptages totaux que le modèle binomial négatif parvient à absorber efficacement grâce à son paramètre d’hétérogénéité.

5.3 Risques statistiques de l’ignorance de la surdispersion

L’un des écueils les plus pernicieux de l’économétrie et de la psychométrie appliquées consiste à estimer un modèle de régression de Poisson en présence de données présentant une surdispersion non corrigée. Sur le plan théorique, l’impact d’une telle omission sur les estimateurs des coefficients de régression (β) peut sembler bénin à première vue : si la fonction de moyenne conditionnelle est correctement spécifiée, les estimateurs obtenus par le maximum de vraisemblance de Poisson conservent leur convergence et restent asymptotiquement sans biais.

En revanche, les conséquences sur l’inférence statistique sont dramatiques. Le modèle de Poisson impose arbitrairement que la matrice d’information de Fisher soit calculée sous l’hypothèse restrictive de l’équidispersion. Lorsque la variance réelle surpasse la moyenne, l’algorithme de Poisson sous-estime systématiquement et parfois massivement les erreurs-types associées aux coefficients de régression. Dans de nombreux cas pratiques où le facteur de surdispersion atteint des valeurs de 3 à 5, les erreurs-types réelles sont de deux à trois fois plus larges que celles rapportées par le logiciel.

Cette atrophie artificielle des erreurs-types entraîne une inflation incontrôlée du taux d’erreur de première espèce (α). Le statisticien se retrouve à rejeter l’hypothèse nulle d’absence d’effet avec des valeurs de p hautement significatives (p < 0,001) pour des prédicteurs qui, en réalité, n’exercent aucune influence au-delà du hasard. Les intervalles de confiance calculés s’avèrent fallacieusement étroits, menant à des conclusions scientifiques erronées et non reproductibles. Le recours systématique à la distribution binomiale négative restaure la validité de l’inférence en ajustant les erreurs-types à la variabilité empirique réelle de l’échantillon.

6. Paramétrisations alternatives dans l’écosystème statistique

6.1 Paramétrisation en moyenne et dispersion (Modèles NB1 et NB2)

Bien que la formalisation canonique fondée sur les paramètres (r, p) soit idéale pour l’exposition probabiliste des épreuves de Bernoulli, elle s’avère malcommode pour la régression statistique, où l’objectif premier est de relier une espérance conditionnelle μ à une matrice de covariables explicatives. Dans leur traité de référence sur l’analyse des données de comptage, Colin Cameron et Pravin Trivedi ont systématisé la reformulation de la loi binomiale négative sous la forme d’un modèle linéaire généralisé à deux paramètres d’intérêt : la moyenne conditionnelle μ et un paramètre de surdispersion scalaire, souvent noté α ou θ selon les traditions.

Cette approche a donné naissance à deux grandes familles de spécifications fonctionnelles pour la variance conditionnelle. La première est la formulation à variance linéaire, désignée sous l’acronyme NB1, où la structure de variabilité s’exprime sous la forme proportionnelle :

Var(Y | X) = μ + αμ = μ(1 + α)

Dans ce modèle, le coefficient d’amplification de la variance par rapport à la moyenne est une constante indépendante du niveau d’espérance lui-même. Bien que séduisante par sa simplicité, la formulation NB1 est moins souvent privilégiée dans les logiciels contemporains que la seconde variante.

La seconde formulation, universellement reconnue sous le nom de modèle standard NB2, adopte une structure de variance quadratique :

Var(Y | X) = μ + αμ²

Cette spécification correspond rigoureusement à la dérivation issue du mélange Poisson-Gamma présentée précédemment, avec α = 1 / r. Le terme αμ² implique que plus l’espérance conditionnelle d’un sujet est élevée, plus l’impact de la surdispersion sur sa variabilité absolue s’accroît de façon exponentielle. Le modèle NB2 constitue aujourd’hui la spécification canonique par défaut dans la quasi-totalité des progiciels de biostatistique et d’économétrie.

6.2 Variations selon les plateformes logicielles

L’un des défis opérationnels majeurs pour le chercheur réside dans la variabilité déconcertante des notations et des conventions algorithmiques adoptées d’un logiciel statistique à l’autre. Il est impératif de maîtriser ces équivalences conceptuelles sous peine de commettre des contresens interprétatifs majeurs lors du passage d’un environnement de calcul à un autre.

Dans l’environnement de programmation statistique R, la fonction fondamentale de densité du noyau dnbinom(x, size, prob, mu) accepte deux paramétrisations alternatives : soit le couple classique (size = r, prob = p), soit le couple moyenne-dispersion (size = r, mu = μ). Lorsque l’on ajuste un modèle de régression par le package tiers MASS à l’aide de la fonction glm.nb(), R estime un paramètre étiqueté Theta. Ce Theta logiciel correspond rigoureusement à l’inverse du paramètre de dispersion de Cameron et Trivedi, c’est-à-dire θ_R = r = 1 / α_NB2. Une valeur de Theta tendant vers l’infini dans R indique donc une absence totale de surdispersion, signifiant la convergence pure vers la loi de Poisson.

À l’inverse, dans les progiciels tels que Stata (via la commande nbreg) ou SAS (au sein de la procédure PROC GENMOD), la modélisation standard rapporte directement le paramètre α. Dans la sortie de résultats de ces logiciels, le paramètre est souvent désigné sous les étiquettes alpha ou dispersion. La variance y est expressément formulée comme μ + αμ². Enfin, des logiciels spécialisés en modélisation par équations structurelles comme Mplus manipulent des fonctions de lien spécifiques qui paramétrisent l’échelle de variance par rapport à la précision log-inclinée. Le tableau suivant synthétise les équivalences analytiques fondamentales régissant ces environnements :

• R (MASS::glm.nb) : Rapporte Theta, où Var(Y) = μ + μ² / Theta.
• Stata (nbreg) : Rapporte /alpha, où Var(Y) = μ + alpha × μ². Équivalence : alpha = 1 / Theta.
• Python (statsmodels) : Paramètre alpha dans NegativeBinomial(link='log'), équivalent direct de la notation α de Cameron et Trivedi.
• Formulation canonique : r = Theta = 1 / α ; p = r / (μ + r).

6.3 Implications pour l’interprétation des coefficients estimés

Dans le cadre de la régression binomiale négative formulée selon le modèle linéaire généralisé, la fonction de lien adoptée de manière quasi universelle est le logarithme népérien :

ln(μ_i) = β₀ + β₁X_i1 + β₂X_i2 + … + β_pX_ip

Le recours à ce lien logarithmique garantit mathématiquement la non-négativité de l’espérance prédite μ_i = exp(X_iβ), une condition sine qua non pour toute variable de comptage. Cette non-linéarité fondamentale transforme la nature de l’interprétation des coefficients β. Contrairement aux régressions linéaires ordinaires où le coefficient β représente la variation additive absolue de la réponse pour une augmentation d’une unité du prédicteur, le coefficient d’une régression binomiale négative opère sur une échelle multiplicative.

L’exponentiation du coefficient estimé, soit exp(β_j), fournit ce que la littérature biomédicale et épidémiologique nomme le ratio de taux d’incidence (ou Incidence Rate Ratio, IRR). Si exp(β_j) est égal à 1,25, cela signifie qu’une augmentation unitaire du prédicteur X_j, toutes les autres covariables du modèle étant maintenues strictement constantes, est associée à un accroissement multiplicatif de 25 % du taux attendu d’événements. Réciproquement, un IRR de 0,80 traduit une contraction de 20 % du nombre moyen d’occurrences.

Pour des applications cliniques ciblées, la restitution des effets sur une échelle purement multiplicative peut obscurcir la portée concrète des résultats auprès des praticiens de santé. Il est alors recommandé d’évaluer les effets marginaux sur l’échelle des comptages bruts, définis par la dérivée partielle ∂μ / ∂X_j = β_j exp(Xβ). Puisque cette dérivée dépend de l’ensemble du vecteur de régression X, elle est usuellement calculée pour un profil clinique représentatif (effets marginaux à la moyenne) ou moyennée sur l’ensemble des observations de l’échantillon (effets marginaux moyens), offrant ainsi une mesure directement tangible en termes de nombre de consultations ou d’épisodes de rechute évités.

7. Applications directes en psychologie quantitative et sciences cognitives

7.1 Modélisation des comportements addictifs et impulsifs

L’addictologie comportementale et l’étude des troubles du contrôle de l’impulsivité constituent des domaines d’application privilégiés pour la régression binomiale négative. Qu’il s’agisse de mesurer le nombre de consommations compulsives de substances au cours d’un mois, le nombre d’épisodes de jeu pathologique chez des parieurs ou les crises d’hyperphagie boulimique, les distributions observées violent systématiquement les postulats de la normalité et de l’équidispersion de Poisson. Ces données sont caractérisées par une asymétrie droite aiguë : une proportion substantielle d’individus présente une fréquence d’actes modérée ou nulle, tandis qu’une minorité clinique présente une fréquence d’actes vertigineuse.

Cette variabilité extrême reflète une hétérogénéité individuelle profonde quant aux facteurs de vulnérabilité biologique et environnementale, tels que les polymorphismes génétiques des récepteurs dopaminergiques, les antécédents de traumatismes infantiles ou le niveau de détresse psychologique actuelle. Dans cette configuration, le modèle binomial négatif capture cette vulnérabilité non observée par le truchement de sa composante de mélange continu. Tenter d’ajuster une régression de Poisson sur de tels comptages aboutirait à une sous-estimation critique des intervalles de confiance entourant l’efficacité des protocoles thérapeutiques évalués.

De surcroît, le phénomène de dépendance pharmacologique induit intrinsèquement un processus de contagion temporelle. La survenue d’une prise de substance altère les circuits neurobiologiques de la récompense et affaiblit les capacités d’inhibition du cortex préfrontal, augmentant substantiellement la probabilité d’une prise ultérieure dans un intervalle rapproché. Comme nous l’avons explicité à travers le modèle de Pólya, cette violation de l’indépendance inter-essais trouve dans la binomiale négative une formalisation mathématique naturelle respectant la réalité biologique de l’addiction.

7.2 Psychométrie et protocoles d’apprentissage séquentiel

En neuropsychologie expérimentale, de nombreux outils d’évaluation standardisés reposent sur des tâches d’apprentissage séquentiel par critères de réussite. Considérons le célèbre test de classement de cartes du Wisconsin (Wisconsin Card Sorting Test, WCST), largement utilisé pour mesurer la flexibilité cognitive et l’inhibition de la persévération. Dans ce paradigme, les sujets doivent identifier une règle de tri changeante (couleur, forme, nombre) sur la base de rétroactions binaires fournies par le clinicien (succès ou échec). Le test comptabilise le nombre d’erreurs commises par le participant avant d’atteindre un bloc stable de dix réponses consécutives exactes.

Ce protocole expérimental s’aligne rigoureusement sur le schéma axiomatique de la loi binomiale négative : le quota de réussites (r = 10) est déterministe et fixé par l’examinateur, tandis que le nombre d’erreurs accumulées (k) est la variable aléatoire sous observation. L’estimation du paramètre p chez un patient fournit une mesure quantitative pure de son efficacité d’apprentissage, épurée des biais induits par des seuils de terminaison arbitraires. La modélisation binomiale négative permet de comparer rigoureusement des populations cliniques — par exemple des patients atteints de lésions préfrontales versus des témoins sains — en testant formellement l’égalité des paramètres d’échelle et de dispersion.

Par ailleurs, cette approche offre une solution élégante pour modéliser les effets de fatigue cognitive survenant au fil d’une batterie neuropsychologique prolongée. En introduisant le rang de l’épreuve comme covariable au sein d’une régression binomiale négative, l’analyste peut dissocier la vitesse d’apprentissage intrinsèque du taux de détérioration de la performance lié à l’épuisement attentionnel du sujet, apportant ainsi une finesse d’interprétation inaccessible aux méthodes traditionnelles d’analyse de variance.

7.3 Accidents, agressions et événements rares en milieu clinique

La gestion de la sécurité au sein des unités psychiatriques fermées et des institutions médico-sociales requiert une quantification rigoureuse des comportements hautement perturbateurs, tels que les passages à l’acte hétéro-agressifs, les gestes auto-mutilatoires ou les tentatives d’évasion. Ces événements partagent la double propriété d’être rares à l’échelle d’une journée donnée, mais d’une gravité institutionnelle considérable. Les séries temporelles de comptage de ces actes révèlent presque systématiquement des variances qui dépassent largement leur espérance moyenne, invalidant l’usage du modèle de Poisson.

Une pratique encore trop répandue consiste à appliquer une transformation mathématique déformante sur les données brutes — telle que la transformation logarithmique ln(Y + 1) ou la racine carrée — afin de tenter de normaliser artificiellement la distribution pour pouvoir appliquer des modèles de régression linéaire gaussienne. Cette démarche est vivement critiquée par les méthodologistes modernes pour de multiples raisons statistiques fondamentales :

• Elle introduit un biais systématique sur l’estimation de l’espérance d’origine en raison de l’inégalité de Jensen, faussant les prédictions réelles.
• Le choix de la constante arbitraire ajoutée aux zéros (par exemple 1, 0,5 ou 0,1) modifie substantiellement les valeurs des coefficients estimés et les seuils de significativité statistique.
• Elle échoue à stabiliser la variance lorsque la surdispersion est hétérogène à travers les sous-groupes cliniques observés.
• Elle produit des prédictions qui, une fois retransformées sur l’échelle naturelle, peuvent s’avérer incohérentes avec la réalité empirique du comptage.

La régression binomiale négative traite directement les comptages sur leur support originel discret et non négatif. En intégrant formellement la queue lourde de distribution propre à ces événements violents, elle préserve l’intégrité de l’échelle de mesure clinique tout en modélisant avec justesse le risque cumulé au sein des unités de soins.

8. Démonstrations empiriques et calculs pas à pas

8.1 Cas pratique standard : Expérience de choix binaire

Afin de concrétiser le fonctionnement analytique de la distribution, examinons une démonstration numérique pas à pas issue d’un protocole expérimental de conditionnement opérant en neurobiologie. Soit un modèle animal placé dans un dispositif expérimental où il doit explorer deux leviers. L’animal doit obtenir exactement r = 4 succès (distribution de récompenses nutritionnelles) pour valider une phase d’apprentissage. Supposons que l’animal opère ses choix de façon purement aléatoire et non biaisée, fixant la probabilité unitaire de succès à chaque essai à p = 0,5. Nous nous posons la question statistique suivante : quelle est la probabilité exacte que cet animal commette précisément k = 6 erreurs avant de finaliser son quota de 4 succès ?

Pour résoudre ce problème, nous mobilisons la fonction de masse de la loi binomiale négative. La première étape consiste à déterminer le coefficient combinatoire approprié régissant l’ordonnancement des événements. L’expérience s’interrompt obligatoirement au quatrième succès. Par conséquent, les 6 erreurs et les 3 premiers succès se sont déployés de façon libre le long des (6 + 4 − 1) = 9 premiers essais. Le coefficient combinatoire se calcule par l’évaluation de 9 combiné 6 :

C(9, 6) = 9! / [6! × (9 − 6)!] = 9! / [6! × 3!] = (9 × 8 × 7) / (3 × 2 × 1) = 504 / 6 = 84

Il existe donc exactement 84 permutations séquentielles d’essais distinctes conduisant à cette configuration globale. L’étape suivante consiste à calculer la probabilité multiplicative d’une configuration élémentaire individuelle composée de 4 succès et de 6 échecs. Sous l’hypothèse d’indépendance stochastique :

P_{élémentaire} = p⁴ × (1 − p)⁶ = (0,5)⁴ × (0,5)⁶ = (0,5)¹⁰ = 1 / 1024 ≈ 0,0009765625

Enfin, la probabilité totale s’obtient par le produit direct du nombre de configurations possibles par la probabilité d’une séquence élémentaire :

P(X = 6) = 84 × (1 / 1024) = 84 / 1024 = 21 / 256 ≈ 0,08203125

La probabilité théorique que l’animal commette rigoureusement 6 erreurs avant d’obtenir ses 4 récompenses est d’environ 8,20 %. Ce calcul manuel démontre sans ambiguïté la décomposition structurelle de la loi : le terme combinatoire pondère la décroissance exponentielle géométrique pour créer une masse de probabilité maximale autour des valeurs centrales.

8.2 Calcul de probabilités cumulées et seuils de décision

Dans la pratique de l’évaluation diagnostique, l’analyste ne s’intéresse que rarement à une probabilité ponctuelle isolée. Il cherche bien plus fréquemment à déterminer une probabilité cumulée, notée P(X ≤ K), représentant la probabilité que le nombre d’échecs n’excède pas un certain seuil critique K, ce qui permet d’attester d’une performance minimale ou de décider de l’interruption précoce d’une expérimentation clinique.

La fonction de répartition cumulative d’une variable binomiale négative s’exprime par la sommation finie de ses masses de probabilités élémentaires. Toutefois, cette sommation peut être avantageusement reliée à une fonction mathématique spéciale de premier ordre : la fonction bêta incomplète régularisée, conventionnellement notée I_x(a, b). L’identité analytique fondamentale s’établit comme suit :

P(X ≤ k) = I_p(r, k + 1) = 1 − I_{1 − p}(k + 1, r)

Où la fonction bêta incomplète s’exprime par le ratio d’intégrales :

I_p(r, k + 1) = [ 1 / B(r, k + 1) ] ∫₀^p t^{r − 1} (1 − t)^k dt

Cette relation analytique est d’une portée pratique fondamentale pour les calculs computationnels. Elle évite d’exécuter des sommations itératives coûteuses en temps de calcul lorsque k atteint des valeurs très élevées. De surcroît, elle permet de déterminer avec une précision rigoureuse les quantiles de la distribution nécessaires à l’établissement de bornes de contrôle de qualité ou de critères d’arrêt précoce dans les protocoles d’essais cliniques adaptatifs :

Seuil d’échecs (K)	P(X = K) ponctuelle	P(X ≤ K) cumulée	Décision diagnostique (α = 0,05)
K = 0	0,0625	0,0625	Performance non exceptionnelle
K = 1	0,1250	0,1875	Zone de performance normale
K = 2	0,1563	0,3438	Zone de performance normale
K = 4	0,1367	0,6367	Médiane dépassée
K = 10	0,0269	0,9531	Déficit d’inhibition avéré (> 95e percentile)

Le tableau ci-dessus illustre la façon dont l’intégration séquentielle permet de fixer un seuil clinique normatif : observer plus de 10 erreurs avant de parvenir à 4 succès survient avec une probabilité cumulée de 1 − 0,9531 = 0,0469 (soit moins de 5 % du temps sous l’hypothèse de hasard pur), permettant au clinicien de conclure à l’existence d’une persévération déficitaire statistiquement significative.

8.3 Analyse de sensibilité face aux variations des paramètres

La morphologie de la courbe de probabilité de la loi binomiale négative fait preuve d’une grande sensibilité aux modifications structurelles de ses paramètres générateurs r et p. Pour saisir pleinement cette plasticité, il est indispensable de modéliser les conséquences géométriques de la translation de ces paramètres sur l’aplatissement et la dispersion du système.

Considérons tout d’abord l’impact d’une contraction de la probabilité de succès unitaire p lorsque le paramètre de quota r est maintenu constant. Lorsque p décroît, l’espérance μ = r(1 − p) / p glisse inexorablement vers la droite, étirant le support effectif des comptages. Simultanément, la variance subit une amplification d’ordre quadratique vis-à-vis de 1 / p. La distribution, initialement pointue et compacte pour des valeurs de p > 0,7, s’aplatit de manière prononcée, formant un plateau étendu caractérisé par une très longue traîne droite. La probabilité de masse associée au mode s’effondre, signifiant une imprévisibilité croissante du temps d’attente exact.

Examinons désormais l’effet d’une variation conjointe des deux paramètres maintenant l’espérance μ rigoureusement constante. Supposons une espérance fixée à μ = 10. Nous pouvons générer cette même moyenne par une infinité de couples (r, p) : par exemple un modèle hautement dispersé où r = 1 (auquel cas p = 1/11 ≈ 0,0909) ou un modèle modérément dispersé où r = 20 (auquel cas p = 20/30 ≈ 0,6667). Dans le premier cas (r = 1), la distribution est une loi géométrique pure : son mode absolu se situe obligatoirement à zéro (P(X = 0) ≈ 0,091) et la variance s’élève à 10 + 10² / 1 = 110. Dans le second cas (r = 20), la distribution adopte une forme en cloche unimodale bien plus symétrique, dont le mode se décale au voisinage de la moyenne (μ = 10) tandis que la variance est comprimée à 10 + 10² / 20 = 15. Cette comparaison démontre à quel point la flexibilité de la binomiale négative lui permet de représenter des profils phénoménologiques radicalement distincts sous une même espérance globale.

9. Méthodes d’estimation des paramètres

9.1 Estimation par le maximum de vraisemblance (EMV)

L’estimation des paramètres d’un modèle binomial négatif repose principalement sur la méthode du maximum de vraisemblance. Soit un échantillon de n observations indépendantes et identiquement distribuées y = (y₁, y₂, …, y_n) issues d’une loi binomiale négative paramétrée par sa moyenne μ et son paramètre de dispersion α (selon la formulation NB2). La fonction de log-vraisemblance de l’échantillon global s’écrit comme la sommation des log-vraisemblances individuelles :

ln L(μ, α; y) = ∑_i=1ⁿ { ln Γ(y_i + α⁻¹) − ln Γ(α⁻¹) − ln(y_i!) − α⁻¹ ln(1 + αμ_i) + y_i ln(αμ_i) − y_i ln(1 + αμ_i) }

Pour maximiser cette fonction par rapport au vecteur des coefficients de régression β (sachant que μ_i = exp(X_iβ)), on calcule le vecteur des dérivées partielles de premier ordre, désigné sous le nom d’équations de score. Le gradient par rapport à β prend la forme compacte :

∂ ln L / ∂β = ∑_i=1ⁿ [ (y_i − μ_i) / (1 + αμ_i) ] X_i = 0

Cette équation met en évidence le mécanisme de pondération intrinsèque : chaque écart empirique (y_i − μ_i) est inversement pondéré par la variance conditionnelle théorique de l’observation (1 + αμ_i). Cependant, la dérivation par rapport au paramètre de surdispersion α fait intervenir la fonction digamma ψ(z) = d ln Γ(z) / dz, introduisant une non-linéarité sévère qui interdit l’obtention d’une solution analytique exacte sous forme close :

∂ ln L / ∂α = ∑_i=1ⁿ { α⁻² [ ln(1 + αμ_i) − ∑_j=0^y_i−1 (j + α⁻¹)⁻¹ ] + (y_i − μ_i) / [ α(1 + αμ_i) ] } = 0

La résolution conjointe de ce système d’équations non linéaires exige le déploiement d’algorithmes d’optimisation numérique itératifs, tels que la méthode de Newton-Raphson ou l’algorithme de Fisher Scoring. Ces algorithmes alternent séquentiellement l’estimation du vecteur des coefficients β conditionnellement à une valeur fixée de α, puis l’optimisation unidimensionnelle de α conditionnellement au vecteur de coefficients mis à jour, jusqu’à la convergence globale des gradients vers zéro. La matrice hessienne inversée obtenue à l’ultime itération fournit l’estimateur de la matrice de variance-covariance asymptotique des paramètres, permettant de calculer les erreurs-types de Wald.

9.2 Méthode des moments et méthodes robustes

Bien que l’estimateur du maximum de vraisemblance possède des propriétés asymptotiques d’efficacité théorique inégalées, il s’avère extrêmement gourmand en ressources computationnelles et peut échouer à converger dans des échantillons de taille modeste ou hautement stratifiés. La méthode des moments offre alors une alternative computationnelle directe et séduisante. Elle consiste à égaliser les moments théoriques de la distribution avec les moments empiriques calculés directement sur l’échantillon d’observations.

En désignant par m la moyenne arithmétique empirique de l’échantillon et par s² sa variance empirique sans biais, le système d’égalisation s’établit simplement :

m = r(1 − p) / p   et   s² = r(1 − p) / p²

La résolution algébrique de ce système de deux équations à deux inconnues conduit directement aux estimateurs des moments pour p et r :

p̂_moments = m / s²   et   r̂_moments = m² / (s² − m)

Ces estimateurs possèdent l’avantage indiscutable de s’évaluer de façon instantanée sans recourir à la moindre boucle d’optimisation numérique. Néanmoins, ils souffrent de faiblesses statistiques majeures. D’une part, leur efficacité asymptotique est nettement inférieure à celle du maximum de vraisemblance, générant des erreurs-types empiriques sensiblement plus amples. D’autre part, ils sont exposés à une instabilité structurelle : si, par inadvertance d’échantillonnage au sein d’un sous-groupe, la variance empirique s² s’avère accidentellement inférieure ou égale à la moyenne m, le dénominateur de r̂ devient nul ou négatif, projetant l’estimateur en dehors de son domaine de définition mathématique admissible.

Face à ces limitations, les statisticiens ont développé des estimateurs semi-paramétriques fondés sur la quasi-vraisemblance. Dans cette approche popularisée par McCullagh et Nelder, le chercheur s’abstient de spécifier la densité exacte de la distribution de probabilité pour ne stipuler que la relation liant la moyenne à la variance : Var(Y) = φV(μ). Les paramètres de régression sont alors estimés par des équations d’estimation généralisées (GEE), tandis que les erreurs-types sont corrigées a posteriori par un estimateur sandwich « robuste » de Huber-White, conférant au modèle une résistance remarquable face à d’éventuelles violations de la distribution Gamma de mélange.

9.3 Approche bayésienne de la distribution binomiale négative

L’inférence bayésienne fournit un paradigme méthodologique particulièrement puissant pour appréhender la distribution binomiale négative, en particulier dans les contextes où l’expérimentateur dispose d’informations préliminaires issues d’études antérieures ou lorsque la taille de l’échantillon est trop réduite pour assurer la validité des théorèmes asymptotiques sous-tendant le maximum de vraisemblance. Dans le cadre bayésien, les paramètres r et p (ou μ et α) cessent d’être considérés comme des constantes inconnues figées, pour devenir eux-mêmes des variables aléatoires dotées de lois de probabilité a priori.

Un défi computationnel classique réside dans le fait qu’il n’existe pas de loi a priori conjuguée canonique bivariée simple pour les deux paramètres d’une binomiale négative sous sa forme générale. Pour contourner cette complexité, les modélisateurs ont recours à des structures hiérarchiques exploitant la représentation en mélange Poisson-Gamma :

• Niveau de comptage : Y_i ~ Poisson(λ_i)
• Niveau individuel : λ_i ~ Gamma(α, β) ou λ_i = exp(X_iβ + ε_i) avec exp(ε_i) ~ Gamma(α⁻¹, α⁻¹)
• Hyper-a priori : β ~ Normale(0, σ²_prior) et α ~ Demi-Cauchy(0, γ) ou Gamma(1, 1)

L’estimation conjointe des densités a posteriori ne pouvant être résolue par simple dérivation analytique, elle est exécutée au moyen d’algorithmes d’échantillonnage par chaîne de Markov Monte Carlo (MCMC), tels que l’échantillonneur de Gibbs ou, plus récemment, l’échantillonneur hamiltonien (HMC) implémenté dans le moteur probabiliste Stan. Cette démarche offre un avantage décisif en recherche neuropsychologique : elle fournit des distributions a posteriori complètes pour chaque paramètre, permettant de dériver des intervalles de crédibilité d’une rigueur absolue sans faire reposer l’inférence clinique sur des hypothèses d’approximation gaussienne asymptotique des résidus.

10. Diagnostic statistique et adéquation du modèle

10.1 Tests formels de surdispersion

Avant d’entériner la supériorité opérationnelle d’un modèle binomial négatif sur un modèle de Poisson plus économe, il est rigoureusement nécessaire de soumettre les données à des tests formels d’adéquation diagnostique visant à tester l’hypothèse nulle d’équidispersion H₀ : α = 0 contre l’hypothèse alternative unilatérale de surdispersion H₁ : α > 0.

L’une des approches les plus élégantes et simples à mettre en œuvre sans avoir à estimer le modèle binomial négatif complet est le test du score développé par Cameron et Trivedi. Ce test évalue la pertinence de la structure de variance en effectuant une régression auxiliaire par moindres carrés ordinaires sans constante. L’analyste ajuste préalablement un modèle de régression de Poisson standard, en extrait les prédictions μ̂_i et construit la variable dépendante standardisée :

z_i = [ (y_i − μ̂_i)² − y_i ] / [ √(2) × μ̂_i ]

On régresse ensuite la variable auxiliaire z_i sur la fonction de variance sélectionnée (par exemple μ̂_i pour le modèle NB2). La statistique de test t associée au coefficient de pente suit asymptotiquement une loi normale centrée réduite sous l’hypothèse nulle. Une valeur de t significativement supérieure à 1,96 (au seuil unilatéral α = 0,025) commande le rejet catégorique du modèle de Poisson en faveur du modèle surdispersé.

Une démarche alternative repose sur le test du rapport de vraisemblance (Likelihood Ratio Test, LRT), dont la statistique est définie par LR = 2(ln L_NB − ln L_Poisson). Cependant, une subtilité méthodologique cruciale doit impérativement être respectée ici. Étant donné que la valeur testée sous l’hypothèse nulle (α = 0) se situe exactement à la frontière physique du domaine admissible des paramètres (α ≥ 0), la théorie asymptotique classique de Wilks ne s’applique plus directement. Comme l’ont démontré Chernoff et Self-Gelfand, la statistique du LRT ne suit pas une loi du Chi-deux à un degré de liberté χ²(1), mais une distribution composite de mélange 0,5χ²(0) + 0,5χ²(1). En pratique, cela impose de diviser par deux la valeur de p classique issue de la distribution χ²(1) standard pour obtenir le seuil d’inférence nominal exact.

10.2 Critères d’information et comparaison de modèles concurrents

L’arbitrage entre différentes spécifications distributionnelles s’articule également autour des critères d’information pénalisée, qui formalisent le compromis fondamental entre la qualité d’ajustement empirique et la parcimonie paramétrique du modèle. Les deux métriques d’information les plus universellement rapportées sont le critère d’information d’Akaike (AIC) et le critère d’information bayésien de Schwarz (BIC), régis par les expressions formelles suivantes :

AIC = −2 ln L + 2k   et   BIC = −2 ln L + k ln(n)

Où ln L représente la log-vraisemblance maximisée, k le nombre total de paramètres libres estimés et n la taille de l’échantillon d’observations. Dans la mesure où le modèle binomial négatif introduit exactement un paramètre scalaire supplémentaire par rapport au modèle de Poisson (le paramètre de dispersion), le modèle binomial négatif sera jugé significativement supérieur dès lors que le gain en log-vraisemblance compense largement la pénalité de deux unités imposée par l’AIC (ΔAIC > 2, une baisse supérieure à 10 unités apportant un soutien empirique indiscutable).

Au-delà des critères scalaires globaux, le diagnostic moderne de la qualité d’ajustement fait appel à des inspections visuelles avancées, au premier rang desquelles figurent les rootograms suspendus (hanging rootograms) théorisés par John Tukey et formalisés pour les modèles de comptage par Kleiber et Zeileis. Contrairement aux histogrammes résiduels classiques qui masquent les déviations locales dans les queues de distribution, le rootogram suspendu trace la racine carrée des fréquences observées directement suspendues à partir des valeurs théoriques prédites par la courbe du modèle. Si le modèle est adéquat, le bas de chaque barre verticale fluctue aléatoirement le long de l’axe horizontal zéro. Des barres s’enfonçant sous l’axe traduisent une sous-estimation systématique des fréquences par le modèle, tandis que des barres flottant au-dessus identifient des zones de sur-estimation, révélant immédiatement les défauts d’ajustement structurels.

10.3 Vérification des hypothèses distributionnelles et robustesse

La validation d’un modèle binomial négatif ne saurait se limiter au constat statistique de la supériorité de son log-vraisemblance globale sur celle de Poisson. Il est impératif d’évaluer la robustesse interne du modèle face à la violation potentielle de ses hypothèses distributionnelles et d’examiner l’influence démesurée que certaines observations atypiques peuvent exercer sur l’estimation du paramètre de surdispersion.

L’analyse des résidus requiert une attention méthodologique particulière. En raison de l’asymétrie marquée des distributions de comptage, les résidus bruts (y_i − μ̂_i) sont foncièrement hétéroscédastiques. Les statisticiens ont donc recours aux résidus de déviance standardisés ou aux résidus de Pearson, définis pour chaque point d’observation par :

r_{Pearson, i} = (y_i − μ̂_i) / √[ μ̂_i(1 + α̂μ̂_i) ]

Ces résidus doivent idéalement présenter une espérance nulle et une variance approximativement unitaire. La présence d’observations individuelles associées à des résidus de Pearson excédant des valeurs absolues de 3 ou 4 doit immédiatement alerter le chercheur sur l’existence potentielle de points aberrants ou de sous-populations cliniques non modélisées.

De surcroît, le paramètre de dispersion α (ou son inverse r) est réputé pour sa sensibilité extrême aux valeurs de levier élevées situées dans la queue supérieure de distribution. L’introduction d’un unique comptage artificiellement gonflé au sein d’un échantillon moyen peut diviser par deux l’estimation du paramètre r, induisant un gonflement généralisé des erreurs-types sur l’ensemble des prédicteurs du modèle. Il est dès lors vivement préconisé de réaliser des analyses de sensibilité systématiques, consistant à réestimer les paramètres par des procédures de rééchantillonnage de type bootstrap ou par validation croisée à k blocs, garantissant que l’inférence clinique ne repose pas sur les artefacts idiosyncrasiques d’une fraction marginale de données.

11. Implémentation pratique et programmation statistique

11.1 Application dans le langage R

L’écosystème open-source R offre une collection d’outils statistiques de pointe pour la manipulation probabiliste et l’ajustement de régressions binomiales négatives. Les fonctions probabilistes fondamentales intégrées au package de base stats permettent de manipuler directement la densité (dnbinom), la distribution cumulée (pnbinom), les fonctions quantiles (qnbinom) et la génération de nombres pseudo-aléatoires (rnbinom). Ces commandes s’exécutent avec souplesse en spécifiant soit le vecteur de probabilités prob, soit l’espérance mu conjointe au paramètre de forme size.

Pour la modélisation multivariée, l’outil canonique incontournable est fourni par la fonction glm.nb() issue de la bibliothèque MASS développée par Venables et Ripley. Contrairement à la fonction standard glm() qui exige la fixation manuelle et arbitraire d’un paramètre d’échelle, glm.nb() intègre un processus d’optimisation conjointe par maximum de vraisemblance alterné, convergeant vers l’estimation simultanée du vecteur des coefficients de régression et du paramètre de forme Theta. L’accès aux ratios de taux d’incidence (IRR) et à leurs intervalles de confiance s’obtient directement par exponentiation des coefficients et de la matrice de covariance :

IRR = exp(coef(modele))   et   IC_95% = exp(confint(modele))

Pour le diagnostic graphique et la vérification des hypothèses, l’extension countreg — hébergée sur le référentiel de développement R-Forge et issue des travaux de Zeileis, Kleiber et Jackman — permet de produire des rootograms suspendus d’une précision graphique absolue grâce à l’instruction générique rootogram(modele, style = "hanging"), permettant d’identifier d’un simple coup d’œil toute déviation résiduelle par rapport au postulat binomial négatif.

11.2 Implémentation sous Python pour les sciences des données

Dans l’écosystème de programmation Python dédié aux sciences des données et à l’intelligence artificielle, l’analyse probabiliste et la modélisation économétrique s’articulent autour des bibliothèques scientifiques complémentaires scipy.stats et statsmodels. Le sous-module scipy.stats.nbinom met à disposition des data scientists l’ensemble des méthodes formelles de la classe des distributions discrètes, notamment le calcul de l’espérance, de la variance, de la fonction de masse (pmf) et de la fonction de répartition cumulative (cdf).

Pour l’ajustement de régressions sur des jeux de données complexes, la bibliothèque statsmodels s’impose comme la référence méthodologique équivalente aux standards de R. L’interface déclarative par formules permet d’ajuster des modèles selon deux approches distinctes : soit via la classe des Modèles Linéaires Généralisés (statsmodels.api.GLM) en stipulant la famille sm.families.NegativeBinomial(alpha=...), soit de manière plus robuste et flexible au moyen de la classe dédiée statsmodels.discrete.discrete_model.NegativeBinomial.

Dans cette seconde classe dédiée, le paramètre de dispersion α (correspondant scrupuleusement au paramètre de surdispersion quadratique du modèle NB2) est directement intégré au processus global d’optimisation numérique via l’algorithme de quasi-Newton BFGS. La commande summary() extrait un tableau de résultats complet restituant pour chaque covariable le coefficient estimé, son erreur-type asymptotique, la statistique de Wald (valeur z), la probabilité critique associée (p-value) ainsi que l’intervalle de confiance à 95 %. L’extraction des ratios de taux d’incidence s’effectue par l’application vectorielle de la fonction exponentielle numpy.exp(modele.params).

11.3 Bonnes pratiques de rapport scientifique des résultats

La communication scientifique des analyses de comptage exige le respect d’une rigueur méthodologique sans faille, conforme aux directives des manuels de publication reconnus, tels que les normes édictées par l’American Psychological Association (APA, 7e édition). Les chercheurs omettent trop fréquemment de rapporter les paramètres structurants de la dispersion, limitant indûment leur communication aux seuls coefficients de régression et à leur niveau de significativité binaire.

Un rapport d’analyse irréprochable doit impérativement expliciter les éléments quantitatifs suivants au sein du texte principal ou de tableaux synthétiques normalisés :

• Justification formelle du modèle : Exposer explicitement les arguments empiriques ayant guidé le rejet du modèle de Poisson, en rapportant la statistique formelle du test du score de Cameron et Trivedi ou la statistique du test du rapport de vraisemblance composite ajusté pour paramètre à la frontière.
• Métrique de surdispersion : Indiquer la valeur ponctuelle estimée du paramètre de dispersion (α ou θ) accompagnée de son erreur-type ou de son intervalle de confiance à 95 %, en prenant soin de définir avec précision la convention de paramétrisation adoptée (NB1 versus NB2, ou paramétrisation du logiciel utilisé).
• Indices de parcimonie : Fournir les valeurs comparatives de l’AIC et du BIC pour le modèle de Poisson de référence et le modèle binomial négatif final afin de rendre transparente la supériorité de l’ajustement.
• Taille des effets : Privilégier la communication des Ratios de Taux d’Incidence (IRR) plutôt que des coefficients logarithmiques bruts, en explicitant systématiquement leurs intervalles de confiance à 95 % respectifs afin de faciliter l’interprétation clinique directe.
• Reproductibilité : Documenter la version exacte de l’environnement de calcul statistique mobilisé, les bibliothèques tierces exploitées et les paramètres de tolérance de convergence des algorithmes numériques.

12. Extensions complexes : Régression et modèles à inflation de zéros

12.1 La régression binomiale négative multivariée

Dans les contextes d’observation épidémiologique et clinique avancée, l’évaluation de l’effet d’une exposition thérapeutique ou d’un facteur de risque environnemental impose le contrôle simultané d’une multitude de variables confondantes. La régression binomiale négative multivariée s’intègre harmonieusement dans cette dynamique en modélisant l’espérance conditionnelle μ_i en fonction d’un vecteur arbitraire de covariables continues, catégorielles ou binaires via la fonction de lien logarithmique.

Une préoccupation méthodologique fondamentale concerne l’hétérogénéité des temps d’observation ou d’exposition entre les participants. Si certains patients sont suivis sur une période de douze mois tandis que d’autres ne sont observés que durant trois mois, comparer directement la fréquence brute de leurs épisodes de soins induirait un biais d’attrition critique. Pour résoudre cette difficulté de manière formelle, on introduit au sein du prédicteur linéaire un terme déterministe invariable désigné sous le terme d’offset :

ln(μ_i) = ln(T_i) + β₀ + β₁X_i1 + … + β_pX_ip

Où T_i symbolise la durée d’exposition spécifique du sujet i. Sur le plan algébrique, cette formulation revient à modéliser le taux par unité de temps (μ_i / T_i) tout en maintenant la variable dépendante sur son échelle discrète d’origine au sein de la fonction de vraisemblance. L’offset possède un coefficient de régression contraint a priori à la valeur 1, préservant ainsi la nature intrinsèque de l’algorithme d’estimation.

De plus, l’introduction de termes d’interaction au sein du prédicteur linéaire exige une grande prudence conceptuelle. Sur l’échelle multiplicative issue du lien log, l’interaction entre deux variables quantitatives ne représente pas une synergie additive de comptages, mais une modification relative du ratio des taux d’incidence. L’analyste doit impérativement décomposer les pentes simples pour visualiser la modulation du risque à différents niveaux des covariables modératrices.

12.2 Modèles binomiale négative avec excès de zéros (ZINB)

Il est fréquent en recherche comportementale et médicale que les données de comptage présentent une proportion de zéros bien plus massive que ce que même la distribution binomiale négative surdispersée est en mesure d’absorber. Ce phénomène survient lorsqu’une part considérable de la population d’étude n’émet aucun comportement cible au cours de la période d’observation. Pour répondre à cette architecture binaire-discrète, Lambert (1992) puis Greene (1994) ont développé les modèles à excès de zéros, aboutissant au modèle binomial négatif avec excès de zéros (ZINB, pour Zero-Inflated Negative Binomial).

Sur le plan conceptuel, le modèle ZINB postule l’existence d’une hétérogénéité qualitative divisant la population en deux sous-groupes latents distincts générant des zéros par deux processus mécaniques fondamentalement différents :

• Les zéros structurels : Ces observations émanent d’individus qui n’appartiennent tout simplement pas à la sous-population à risque et dont la probabilité d’émettre l’événement est rigoureusement égale à zéro en toute circonstance. C’est le cas par exemple d’individus totalement abstinents face à l’alcool dans une étude sur les épisodes d’ébriété aiguë.
• Les zéros d’échantillonnage (ou accidentels) : Ces zéros proviennent d’individus qui sont de réels consommateurs potentiels (appartenant à la population à risque), mais qui, par pure contingence stochastique durant la fenêtre d’observation limitée de l’étude, n’ont émis aucun épisode d’ivresse.

Mathématiquement, le modèle ZINB combine un modèle logistique binaire (qui prédit la probabilité π d’appartenir au groupe des zéros structurels) et un modèle de comptage binomial négatif traditionnel appliqué à l’ensemble des sujets susceptibles d’émettre l’événement. La fonction de masse composite s’exprime par le système d’équations par morceaux :

P(Y = 0) = π + (1 − π) × P_NB(Y = 0 | μ, α)

P(Y = y) = (1 − π) × P_NB(Y = y | μ, α)   pour tout y > 0

L’arbitrage formel entre un modèle binomial négatif standard et un modèle ZINB ne peut pas s’appuyer sur un simple test du rapport de vraisemblance emboîté classique, car les modèles ne sont pas formellement hiérarchisés. On recourt alors au test non emboîté de Vuong, fondé sur la distribution asymptotique de la différence moyenne des log-vraisemblances individuelles de Kullback-Leibler. Une statistique de Vuong significativement positive (V > 1,96) atteste que l’introduction du processus d’inflation de zéros améliore de façon déterminante la qualité prédictive du modèle par rapport à la binomiale négative usuelle.

12.3 Modèles à troncature et modèles de survie discrets

Une dernière famille d’extensions analytiques concerne les contextes de collecte de données où le zéro est structurellement exclu de l’espace des observations accessibles. Considérons par exemple une enquête épidémiologique menée exclusivement auprès de patients hospitalisés dans laquelle on mesure la durée de séjour exprimée en jours complets révolus. Par définition même du protocole d’échantillonnage, aucun patient admis à l’hôpital ne peut présenter une durée de séjour de zéro jour. Les données sont dites tronquées à gauche en zéro.

L’application d’un modèle binomial négatif classique sur des données tronquées produirait une distorsion sévère de l’estimation de l’espérance, la masse théorique attribuée à zéro tirant faussement l’espérance conditionnelle vers le bas. On recourt alors à la distribution binomiale négative tronquée à zéro (Zero-Truncated Negative Binomial, ZTNB). La fonction de masse s’obtient en renormalisant la distribution originelle par la probabilité complémentaire de non-nullité :

P(Y = y | Y > 0) = P_NB(Y = y) / [ 1 − P_NB(Y = 0) ] = P_NB(Y = y) / [ 1 − (1 + αμ)^−α−1 ], pour tout y ≥ 1

Enfin, la distribution binomiale négative entretient des liens théoriques étroits avec l’analyse de survie en temps discret. Dans de nombreux protocoles longitudinaux, le suivi clinique s’effectue à intervalles de temps réguliers discrets (visites mensuelles ou annuelles). L’apparition d’un événement critique d’intérêt (par exemple le décès, l’abandon thérapeutique ou la rechute clinique) peut être modélisée en envisageant le nombre d’intervalles franchis avec succès avant la survenue de l’événement comme une variable binomiale négative. Cette passerelle méthodologique unifie l’analyse des comptages et l’histoire des événements au sein d’un même corpus statistique rigoureux, ouvrant des horizons prometteurs pour la modélisation des trajectoires cliniques longitudinales.

Références

• American Psychological Association. (2020). Publication manual of the American Psychological Association (7th ed.). American Psychological Association. https://doi.org/10.1037/0000165-000
• Cameron, A. C., & Trivedi, P. K. (2013). Regression analysis of count data (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9781139013567
• Greenwood, M., & Yule, G. U. (1920). An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents. Journal of the Royal Statistical Society, 83(2), 255-279. https://doi.org/10.2307/2341080
• Hilbe, J. M. (2011). Negative binomial regression (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511973420
• Kleiber, C., & Zeileis, A. (2016). Visualizing count data regressions using rootograms. The American Statistician, 70(3), 296-303. https://doi.org/10.1080/00031305.2016.1173590
• Lambert, D. (1992). Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics, 34(1), 1-14. https://doi.org/10.2307/1269547
• McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman & Hall/CRC. https://doi.org/10.1007/978-1-4899-3242-6
• Pólya, G. (1930). Sur quelques points de la théorie des probabilités. Annales de l’Institut Henri Poincaré, 1(2), 117-161. http://www.numdam.org/item/AIHP_1930__1_2_117_0/
• Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S (4th ed.). Springer. https://doi.org/10.1007/978-0-387-21706-2
• Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica, 57(2), 307-333. https://doi.org/10.2307/1912557