Tutoriel : GLM sur données de comptage (régression de Poisson) avec R

glm données de comptage logiciel R

Les GLM (modèles linéaires généralisés) sur données de comptage, ou régression de Poisson, sont des approches statistiques qui doivent être employées lorsque la variable à analyser résulte d’un processus de comptage (comme un nombre d’œufs pondus, un nombre de buts marqués, ou encore un nombre de visites sur un site internet). Ces approches sont indispensables, car dans cette situation les hypothèses des modèles linéaires classiques ne sont plus satisfaites.

De manière plus précise, les modèles de régression de Poisson, sont des GLM, comportant une fonction de lien log et une structure d’erreur de type Poisson.

Dans cet article, je vous propose quelques rappels d’ordre théorique concernant les GLM sur données de comptage, suivis d’un tutoriel pour mettre en oeuvre cette approche, en pas à pas, avec R, dans deux situations courantes : la régression linéaire, et l’ANOVA.

 

1. Rappels

1.1 Pourquoi les modèles linéaires classiques ne sont pas adaptés

Les modèles linéaires classiques ne sont pas adaptés pour analyser des variables à expliquer (ou réponses) de type “comptage”, notamment parce qu’ils supposent que celles-ci sont distribuées selon une loi Normale. Cette hypothèse conduit alors à considérer que la variance des résidus est homogène, autrement dit constante, quelle que soit la valeur des comptages moyens prédits par le modèle. Or, les données de type comptage ne sont pas distribuées selon une loi Normale, mais selon une loi de Poisson. Et compte tenu de cette loi de distribution ,la variance des résidus n’est pas constante mais proportionnelle aux comptages moyens prédits par le modèle.

Cette différence n’est pas une nuance, elle est capitale ! Utiliser un modèle linéaire classique pour traiter des réponses de type “comptage” peut entraîner, entre autres, une estimation biaisée de l’erreur standard des paramètres du modèle, ce qui peut conduire à une estimation erronée de la p-value, et à une conclusion inexacte.

De plus, l’utilisation de modèles linéaires classiques pour analyser des données de type comptage peut conduire à des prédictions négatives. Voir un exemple ici .

 

1.2 La distribution de Poisson

On dit qu’une variable aléatoire Y suit une distribution de Poisson de paramètre Lambda, si elle prend pour valeur y = 0 ,1,2,3,… avec une probabilité P définie par :

$$Pr({Y=y}) = \frac{e^{-\lambda}\; \lambda^y}{y\;!} $$

La distribution de Poisson est ainsi définie par un seul paramètre : Lambda.

Pour fixer les idées, voici quelques exemples de distribution, pour des valeurs de Lambda variant entre 0.1 et 30.

distributions de poisson en fonction de Lambda

 

Plus Lambda augmente, plus la distribution de Poisson se rapproche d’une loi Normale. On dit qu’une loi de Poisson peur être approximée par une loi Normale quand Lambda est égal 20 ou 30 ( ça dépend des auteurs !).

La distribution de Poisson possède deux éléments remarquables :

  • L’espérance ( ou moyenne) d’une variable aléatoire distribuée selon une loi de poisson est égale à Lambda :

$$ E(y) = \lambda $$

  • La variance d’une variable aléatoire distribuée selon une loi de poisson est aussi égale à Lambda :

$$ Var(y) = \lambda $$

1.3 Caractéristiques des GLM pour données de comptage

Comme expliqué dans l’article d’introduction aux GLM, ces modèles sont constitués de trois éléments :

  • un prédicteur linéaire,
  • une fonction de lien Log
  • une structure d’erreur de Poisson

Le prédicteur linéaire suppose simplement que les réponses prédites par le modèle linéaire généralisé vont l’être à partir d’une combinaison linéaire des variables prédictives.

$$ µ_y = \sum_{j=1}^{p}\;\beta_j\;X_{ij} $$

On nomme généralement ce prédicteur linéaire par la lettre ƞ (eta):

$$ \eta \; = \sum_{j=1}^{p} \beta_j\;X_{ij} $$

Cependant, contrairement aux modèles linéaires classiques, les valeurs prédites par le prédicteur linéaire ne correspondent pas à la prédiction moyenne d’une observation, mais à la transformation (par une fonction mathématique) de celle-ci. Dans le cas de la régression de Poisson il s’agit de la transformation log. Ce lien est défini par l’équation suivante :

$$log(µ_y) =\sum_{j=1}^{p} \beta_j\;X_{ij} $$

En pratique, cela signifie que les valeurs du prédicteur linéaire sont obtenues en transformant préalablement les valeurs observées par la fonction de lien. Autrement dit, les beta sont estimés après transformation des réponses par la fonction Log.

Pour obtenir la prédiction moyenne, il est alors nécessaire d’appliquer la fonction inverse du Log, c’est à dire la fonction exponentielle :

$$ µ_y = e^{\sum_{j=1}^{p} \beta_j\;X_{ij}} $$

Les modèles de régression de Poisson ont une structure d’erreur de Poisson. Cette structure d’erreur, permet notamment de spécifier correctement la relation entre la moyenne et la variance. Cette relation est utilisée par l’approche de maximum de vraisemblance pour estimer les coefficients et les erreurs standard des paramètres du modèle.

1.4 Les conditions de validité du modèle de régression de Poisson

Les résultats issus des modèles de régression de Poisson sont valides
si:

  • les réponses sont indépendantes.
  • les réponses sont distribuées selon une loi de Poisson, de paramètre Lambda.
  • il n’existe pas de surdispersion

 1.4.1 L’indépendance des réponses

L’indépendance des réponses signifie qu’elles ne doivent pas être corrélées entre elles. Par exemple, il ne doit pas avoir de lien entre une réponse et celle de la ligne suivante, ou précédente, dans le jeu de données. On voit cela, lorsque des données sont répétées sur des sujets, ou des unités expérimentales identiques.

L’indépendance des réponses se valide par l’étude du plan expérimental : pour réaliser une régression de Poisson, il ne doit pas avoir de données répétées. Si c’est le cas, il faut utiliser un autre type de modèle, comme un modèle linéaire généralisé à effet mixte (Generalized Linear
Mixed Models, ou GLMM).

1.4.2 La distribution des réponses selon une loi de Poisson

La distribution des réponses (selon une loi de Poisson) est généralement supposée. Néanmoins, cette hypothèse peut être explorée en comparant la distribution des comptages observés à la distribution théorique des comptages sous une loi de poisson de paramètre Lamda, estimé par la moyenne des comptages observés. Si cette hypothèse est rejetée, un autre structure d’erreur devra être utilisée.

1.4.3 Absence de surdispersion

En théorie (c’est à dire selon la loi de Poisson), la variance des réponses est égale au paramètre Lambda. Il y a surdispersion (overdispersion en anglais) lorsque la variance des réponses observée est supérieure à cette variance supposée (ou théorique). Dans cette situation, l’erreur standard des paramètres est sous estimée. Ceci peut conduire à une p-value excessivement faible, et donc aboutir à une conclusion erronée.

La surdispersion est généralement évaluée par le ratio de la déviance résiduelle sur le nombre de degrés de liberté du modèle :

$$ \widehat{\phi} = \frac{{deviance\; residuelle}} {nddl} $$

Si ce ratio est supérieur à 1, alors il y a surdispersion. Cette supériorité à 1 est un peu subjective, dans le sens où il n’existe pas de seuil (1.5, 2 ?) à partir duqel on considère qu’il ya surdispersion.

Les principales causes des surdispersion sont :

  • une corrélation entre les réponses,
  • l’absence d’une variable explicative importante,
  • un sur-représentation des valeurs zéro par rapport à ce qui est attendue selon la distribution de Poissson de paramètre Lambda.

En cas de surdispersion, il est nécessaire d’utiliser d’autres structures d’erreur, telles que les structures “quasi Poisson” ou “négative binomiale”.Celles-ci vont conduire à une augmentation de l’erreur standard des paramètres du modèle, par un facteur :

$$ \sqrt(\widehat{\phi}) $$

 

2.Les Data

Dans la suite de cet article, je vais vous montrer comment réaliser une régression linéaire simple lorsque la variable réponse est de type comptage. Puis dans un second temps, je vous montrerai, dans la même situation, comment réaliser une analyse de type ANOVA à un facteur.

Pour illustrer ces deux types d’analyse, je vais utiliser le jeu de données “Orstein” du package “car”.

 

Ce jeu de données concerne les 248 compagnies les plus importantes du Canada dans les années 70. Il comporte 4 variables :

  • interlocks : nombre de postes d’administrateurs partagés avec d’autres grandes entreprises. Il s’agit de la variable de type comptage à expliquer.
  • assets : actifs de la compagnie en millions de dollars ; il s’agit d’une variable de type numérique continue
  • sector : secteur industriel de la compagnie, il s’agit d’une variable catégorielle à 10 modalités
  • nation : nation contrôlant la compagnie, il s’agit d’une variable catégorielle à 4 modalités.

 

3 Tutoriel régression de Poisson de type régression linéaire simple

Dans cette analyse, nous allons étudier la relation entre la variable réponse “interlocks” et la variable explicative “assets”.

Cette modélisation correspond à l’équation suivante :

$$log(\mu_y) = \beta_0 + \beta_1*log10(assets)_{i} $$

 3.1 Visualisation de la relation

La première étape consiste à visualiser cette relation par un scatterplot :

 

regression de poisson avec R

La variable assets est extrêmement étendue, nous allons alors utiliser une transformation log10.

 

régression de poisson avec R

La transformation log10 permet également de linéariser la relation, elle est donc nécessaire.

 

3.2 Ajustement du modèle de régression linéaire de Poisson

L’ajustement est réalisé à l’aide de la fonction glm() et de l’argument family=”poisson”.

 

3.3 Evaluation des hypothèses de la régression de Poisson

3.3.1 Indépendance des réponses.

Nous allons considérer que cette hypothèse est satisfaite.

 

3.3.2 Distribution des réponses selon une loi de Poisson

Une méthode simple pour évaluer cette hypothèse, consiste à comparer la distribution des comptages réellement observés, avec leur distribution théorique, selon la loi de Poisson de paramètre Lambda (estimé par la moyenne des comptages).

Pour cela, on commence par estimer la moyenne des comptages :

 

La moyenne des comptages est égale à 13.58.

Puis on simule des comptages selon une distribution de Poisson de paramètres Lambda = 13.58.

 

glm poisson count with r

Les comptages observés (variable interlock) sont en bleu et les comptages théoriques en rouge.On peut voir que les deux distributions sont très différentes. La variable interlock ne suit donc pas une distribution de Poisson de paramètre lambda = 13.58. On peut également voir une sur-représentation des valeurs zéro dans les comptages observés ; il est donc très probable qu’une surdispersion soit mise en évidence.

3.3.3 Evaluation de la surdispersion

 

Le ratio residual deviance / ddl est égal à 1904.7 /246, soit 7.74. Ce ratio est très largement supérieur à 1 et permet de mettre en évidence la présence d’une surdispersion. Il est donc nécessaire d’utiliser une autre structure d’erreur dans le modèle de régression.

3.4 Utilisation d’autres structures d’erreurs

3.4.1 Structure d’erreur quasi poisson

Pour cela, il suffit de remplacer family=”poisson” par family=”quasipoisson” :

 

Comme expliqué précédemment, l’utilisation d’une structure d’erreur “quasi poisson” à pour conséquence d’augmenter l’erreur standard des paramètres.

 

Ici l’erreur standard associé à log10(assets) est égale à 0.07, alors qu’elle était estimée à 0.02 avec la structure d’erreur de Poisson.

 

 3.4.2 Structure d’erreur négative binomiale

En cas de surdispersion, il est également possible d’utiliser une structure d’erreur quasi binomiale. Pour cela, on utilise la fonction glm.nb() du package MASS. Par défaut, la fonction de lien utilisée est également la fonction log :

 

Avec cette structure d’erreur, l’erreur standard associée à log10(assets) est estimée à 0.09. Cette correction entraîne une modification très importante de la z value, qui passe de 44.22 avec une structure d’erreur de Poisson à 11.257 ici. Dans d’autres situations, cela pourrait également avoir un impact très important sur la p-value, et donc sur les conclusions.

 

3.5 Interprétation des résultats.

En reprenant les résultats obtenus avec la structure d’erreur “quasi poisson”, les résultats sont obtenus à l’aide de la fonction summary() :

 

Comme rappelé précédemment, dans une régression de poisson, c’est le log de la moyenne des comptages qui est modélisé :

 

$$log(\mu_y) = \beta_0 + \beta_1*log10(assets)_{i} $$

Ici, Beta_0 = -1.04 et Beta_1 = 1.05.

L’interprétation de l’intercept (beta_0) ne présente généralement pas d’intérêt particulier. On s’intéresse alors uniquement au coefficient associé à la variable prédictive.

Le coefficients 1.05 associé à log10(asset) signifie qu’une augmentation d’une unité de log10 des actifs de la compagnie, est associée avec une augmentation de 1.05 du log du nombre moyen d’interlocks (c’est à dire du nombre moyen de postes d’administrateurs partagés).

Pour exprimer le coefficient de régression dans l’échelle originale, c’est à dire en nombre d’interlocks, plutôt qu’en log du nombre d’interlocks, il est nécessaire d’appliquer la transformation inverse du log, c’est à dire la transformation exponentielle.

$$ e^{1.05} = 2.86 $$

Ainsi, le nombre de postes d’administrateurs partagés augmente en moyenne de 2.86 lorsque le log10 des actifs augmente d’une unité.

Au final, en exprimant également les actifs dans leur échelle originale, on peut dire que le nombre de postes d’administrateurs partagés augmente en moyenne de 2,86 lorsque les actifs de la compagnie sont multipliés par 10 (puisqu’un unité en log10, correspond à un facteur 10).

3.6 Prédictions

La fonction predict.glm() permet d’estimer le nombre de postes d’administrateurs partagés prédit, en moyenne, par le modèle, pour un ou plusieurs niveaux d’actifs souhaités.

Pour cela, on crée un data frame nommé “mydf” (par exemple) contenant les niveaux d’actifs pour lesquels on souhaite obtenir une estimation, ici par exemple, 100,1000,5000,10000, 50000,et 100000. Il est très important d’utiliser le même nom de variable que dans le fichier de données utilisé lors de l’ajustement du modèle.

 

Ce dataframe est ensuite passé dans l’argument newdata de la fonction predict.glm(). L’argument type=”response” permet d’obtenir les prédictions dans l’échelle originale (c’est à dire en nombre de postes).

 

Par défaut, la fonction predict.glm() utilise l’argument type=”link”, ce qui renvoi des prédictions dans l’échelle log :

 

3.7 Visualisation du modèle

Le principe est de créer un vecteur de valeurs de assets compris entre les valeurs min et max observées, puis de le passer dans un data frame.

 

Dans un second temps, la fonction predict.glm() est employée pour prédire le nombre moyen d’interlock pour chaque valeurs du vecteur créé précédemment. Ces prédictions sont ajoutées au data frame.

 

Enfin le plot est réalisé :

 

glm données de comptage logiciel r

 

4 Tutoriel régression de Poisson de type ANOVA à un facteur :

Pour illustrer ce cas de figure, nous allons comparer le nombre moyen de postes d’administrateurs, en fonction de la nation contrôlant la compagnie. Cette analyse n’est pas particulièrement judicieuse, car, ici, seule la variable explicative nation est inclue. Cette analyse sert seulement d’illustration.

Pour plus d’information sur l’ANOVA à un facteur, consultez cet article.

4.1 Visualisation

Une première visualisation permet d’appréhender les données :

 

anova sur données de comptages logiciel r

4.2 Calcul des moyennes et de leur intervalle de confiance

Les moyennes peuvent être facilement calculées en couplant les fonctions group_by() et summarise() du package dplyr.

 

Les intervalles de confiance peuvent être obtenus en employant une approche bootstrap par l’intermédiaire de la fonction slipper() du package slipper.

 

C’est un peu fastidieux de le faire pour chaque nation, mais au final on obtient :

 

Certains intervalles de confiance ne se chevauchent pas. Ceci laisse supposer qu’il existe des différences du nombre de postes moyens entre les nations.

4.3 Ajustement du modèle ANOVA de Poisson

Comme précédemment, on utilise simplement la fonction glm() avec l’argument family=”poisson” :

 

4.4 Evaluation des hypothèses de la régression de Poisson

4.4.1 Indépendance des réponses.

Nous allons considérer que cette hypothèse est satisfaite.

4.4.2 Distribution des réponses selon une loi de Poisson

La représentation graphique réalisée précédemment, dans le chapitre dédié à la régression linéaire simple, nous a déjà permis de conclure que la variable interlocks n’est pas distribuée selon une loi de Poisson de paramètre Lambda = 13.58.

4.4.3 Evaluation de la surdispersion

 

Le ratio residual deviance / ddl est égal à 3064.5 /244, soit 12.56. Ce ratio est très largement supérieur à 1 et permet de mettre en évidence la présence d’une surdispersion. Il est donc nécessaire d’utiliser une autre structure d’erreur dans le modèle de régression.

Comme précédemment, nous allons pouvoir utiliser les structures d’erreur de type “quasi Poisson” et “negative binomial”.

4.5 Utilisation d’autres structures d’erreurs

4.5.1 Structure d’erreur quasi poisson

 

4.5.2 Structure d’erreur negative binomiale

 

Comme attendu, dans les deux cas, les erreurs standard des paramètres ont été largement augmentées.

 

 4.6 Interprétation des résultats.

Pour obtenir la table d’analyse de variance, on utilise la fonction Anova() du package “car” :

 

Ici, les résultats montrent que deux moyennes, au moins, sont  différentes.

4.7 Comparaisons multiples

Afin d‘identifier qu’elles sont les différences, nous pouvons réaliser toutes les comparaisons deux à deux à l’aide du package multcomp et de sa fonction glht(). Pour plus de détail, consultez cet article.

 

Les comparaisons peuvent alors être visualisées :

 

comparaisons multiples glm

Le package multcomp, contient également une fonction cld() qui permet, dans le cadre du test de Tukey, d’indiquer par des lettres la significativité des comparaisons. Lorsque deux modalités partagent une même lettre, cela signifie que leurs différences ne sont pas significativement différentes. A l’inverse, lorsque deux modalités ne partagent pas de lettres en commun, alors cela signifie que leurs moyennes sont significativement différentes.

 

Comme expliqué dans un de mes précédents articles, on peut alors utiliser ces lettres pour les ajouter sur un graph réalisé avec ggplot2.

 

 

comparaisons multiples avec R

NB : les losanges noirs représentent les moyennes Remarques :

 

On peut également obtenir les intervalles de confiance des niveaux moyens, comme ceci:

 

Le modèle ayant été ajusté sans intercept (grâce au -1), les coefficients fournis pas la fonction summary() correspondent au log des niveaux moyens de chaque nation. On peut également les obtenir par :

 

Pour les obtenir dans l’échelle originale, c’est à dire en nombre de postes partagés, il est nécessaire d’utiliser la transformation exponentielle :

 

Les intervalles de confiance s’obtiennent alors comme ceci :

 

Ils sont légèrement différents de ceux obtenus par bootstrap, cela n’est pas étonnant, car l’approche par bootstrap est basée sur des tirages aléatoires.

 

J’espère que cet article aura répondu à un grand nombre de questions que vous vous posiez au sujet des GLM pour données de comptage, qu’elles soient d’ordre théorique ou pratique. Si certaines interrogations subsistent, indiquez-les-moi en commentaire, et j’essaierai de vous apporter une réponse.

En attendant, si cet article vous a plus, n’oubliez pas de la partager 😉

Crédit photo : Clker

Prolongez votre lecture :

 

Partager l'article
  •  
  •  
  •  
  •  
  •  
    17
    Partages
  • 17
  •  
  •  
  •  
  •  

5 commentaires

  1. YAO ANICET G KOUAME Répondre

    Bonjour Claire. Félicitations et merci pour cet article très bien élaboré qui, à coup sûr, améliorera la compréhension de beaucoup sur les G LM.

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *