Modèle linéaire généralisé avec R

Définition

Il s'agit de modèles linéaires plus souple que le modèle linéaire "traditionnel", qui vont être utiles lorsque la variable à expliquer n’est pas une variable quantitative continue. En effet, le modèle linéaire général ne convient pas pour certaines variables à expliquer, comme par exemple pour les données de comptage ou pour les probabilités.

Par exemple, les valeurs prédites peuvent en effet n’être pas entière (pour les données de comptage) ou hors de l’intervalle [0,1] (pour une probabilité).

En outre, le modèle linéaire général repose sur l’hypothèse de normalité du terme d’erreur, et la variable à expliquer y est ainsi aussi considérée comme distribuée normalement (hypothèse très forte). (voir Regression linéaire avec R). Or, il paraît souvent abusif de considérer la variable à expliquer (et son erreur) comme suivant une loi normale puisque cela supposerait une distribution symétrique (alors que par exemple ce n’est très probablement pas le cas pour des données de comptage).

Enfin, avec l’hypothèse d’espérance nulle du terme d’erreur, on a : <math> E(y/x) = b0+b1 \times x1+b2 \times x2. </math>
Le modèle linéaire général considère ainsique les effets des variables explicatives sur la variable à expliquer sont linéaires, ce qui n’est très souvent pas le cas.

Le modèle linéaire généralisé considère plutôt une relation du type : <math> E(y/x)= g(b0+b1 \times x1+b2 \times x2) </math>

La fonction inverse de g(...), est appelée la fonction de liaison, et prend la forme : <math>f(E(y/x)) = b0 + b1 \times X1 + b2 \times X2 + ... + bk \times Xk </math>

Enfin, le modèle linéaire généralisé va tenir compte de la non-normalité du terme d’erreur, et choisir une autre distribution pour la modéliser.

La fonction GLM()

Généralités

La fonction glm() est utilisée dans R, avec la syntaxe générale : glm(formula, family=familytype(link=linkfunction), data=dataset)

où :

• formula est la formule classique pour écrire le modèle à estimer (de type : y ~ x1 + x2 + ...) où y est la variable à expliquer et x1,x2,... sont les variables explicatives.
• family représente la distribution du terme d’erreur et link la forme de la fonction lien spécifiée.

Quelques exemples :

?family

• binomial(link = "logit") #Modèle LOGIT
• gaussian(link = "identity")
• Gamma(link = "inverse")
• inverse.gaussian(link = "1/mu^2")
• poisson(link = "log")
• quasi(link = "identity", variance = "constant")
• quasibinomial(link = "logit")
• quasipoisson(link = "log")

A noter que family=« gaussian » renvoie à la fonction lm et au modèle linéaire général puisque la fonction de lien est l’identité. Parmi les modèles les plus utilisés on trouve :

• le modèle LOGIT et PROBIT (lorsque la variable à expliquer est de type binaire (typiquement 0-1), et que l’on cherche à expliquer une probabilité : par exemple la probabilité pour un composant d’être défectueux…
• le modèle de Poisson, utilisé pour des données de comptage (ex : nombre de visites chez le médecin, nombre de buts marqués,...).

Modèle LOGIT et PROBIT :

Présentation

Il est sont l’un des usages les plus courants de la fonction glm(). Il sont utilisés lorsque la variable à expliquer est de type binaire, et que les variables explicatives sont continues ou binaires. Tous deux utilisent une loi binomiale (family).

<math> Y=g(b0+b1 \times x1 +b2 \times x2 +....)</math>

L'inverse de la fonction de lien (notée g précédemment) peut être une loi logistique (ce qui donne un modèle LOGIT) ou celle de la loi Normale centrée réduite (modèle PROBIT). A noter que le modèle LOGIT est aujourd’hui le plus utilisé (l’avantage initial du PROBIT était la fonction de lien plus « simple » à calculer à la main, ce qui n’est plus important aujourd’hui avec l’utilisation de logiciels).

Application : modèle

Pour illuster, on utilise le jeu de données contenu dans R : mtcars, qui contient 32 observations et 11 variables. On souhaite déterminer si la probabilité pour une voiture d’être dotée d’une boîte à vitesses manuelle (variable am) dépend où non du nombre de chevaux (variable hp) et de son poids (variable wt). On a un exemple ici où le modèle linéaire général peut produire une valeur estimée de cette probabilité (p=E(am/wt,hp)) ou des effets marginaux hors de l’intervalle [0,1].

##### LOGIT

regLOGIT=glm(am ~ hp + wt, data=mtcars, family=binomial(link=logit))
summary(regLOGIT)

##### PROBIT

regPROBIT=glm(am ~ hp + wt, data=mtcars, family=binomial(link=probit))
summary(regPROBIT)

Fonction (affichage, extraction,..)

Les fonctions présentées pour la fonction lm() (voir : Regression linéaire avec R) sont aussi valables pour la fonction glm(), à savoir pour les plus utilisées :

• summary() pour afficher la sortie complète
• coef() pour extraire les coefficients
• confint() pour extraire leur intervalle de confiance (à 95%)
• fitted() pour extraire les valeurs estimées
• resid() pour extraire les résidus.

Comme le modèle est estimé par maximum de vraisemblance, il peut être utile d’extraire la vraisemblance du modèle, avec la fonction logLik() (NB : Cela est aussi possible pour un objet de classe lm, mais peu utilisé)

#vraissemblance
logLik(regLOGIT)

Interprétation et prédiction

De plus, bien que possible, il n’est pas d’usage d’interpréter les coefficients de tels modèles. En effet, en raison de la fonction qui définit la relation entre la variable à expliquer et les variables explicatives, les coefficients présentés dans le tableau précédent n’ont pas de signification en eux-mêmes (sans transformation), et l'effet marginal (i.e le coefficient associé à un variable) va dépendre des valeurs prises par l'ensemble des variables explicatives, et pas seulement par la variable concernée. (A noter qu'il est d'usage de prendre alors des valeurs moyennes pour les autres variables).

En effet, pour le logit la fonction de lien utilisée étant une fonction logistique, on a : <math> f(y)=log(\frac{p}{1-p})=18.87+0.036 \times hp-8.08 \times wt </math>
où p = probabilité pour la voiture d'être dotée d'un boîte de vitesses manuelles (P(am)=1).

En revanche, on peut utiliser le modèle pour la prédiction :

##Prédiction (modèle logit) :
predict(regLOGIT, newdata=data.frame(hp=100,wt=2.88), type="response")

Une voiture avec une puissance de 100 chevaux et un poids de 2.88 lb/1000 va donc avoir une probabilité estimée d’être équipée d’une boîte manuelle de 31% environ.