Modèles linéaires généralisés

Valérie Monbet

IRMAR, Université de Rennes 1

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 1 / 203

Références, supports
Ce cours est essentiellement une compilation des cours de G. Rodriguez (Princeton
University) et A. Lavenu (Université de Rennes 1).

Hardin, J. and Hilbe, J. (2012). Generalized Linear Models and Extensions, 3rd
Edition. College Station, Texas : Stata Press. Un livre avec des exemples et des
applications incluant des analyses avec Stata.
Notes de cours de G. Rodriguez et exemples de codes R :
http://data.princeton.edu/wws509/
Notes de cours de L. Rouvière
http://perso.univ-rennes2.fr/system/files/users/rouviere_l/poly_logi
Notes de cours de F. Bertrand
http://www-irma.u-strasbg.fr/~fbertran/enseignement/
Ecole_Doctorale_SVS_Automne_2008/ED_RegLog.pdf
Hosmer, D.W. and Lemeshow, S. (2013). Applied Logistic Regression, 3rd Edition.
New York : John Wiley and Sons. Une discussion détaillée sur le modèle logistique
avec des applications.
McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models, 2nd Edition.
London : Chapman and Hall. La "bible" des modèles linéaires généralisés. Très
intéressant, mais plutôt destinée à des étudiants avancés.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 2 / 203

 Introduction

Outline

1 Introduction
Modèles linéaires pour les données continues
Modèles linéaires pour les données discrètes

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 3 / 203

 Introduction Modèles linéaires pour les données continues

Outline

1 Introduction
Modèles linéaires pour les données continues
Modèles linéaires pour les données discrètes

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 3 / 203

 Introduction Modèles linéaires pour les données continues

Modèles linéaires pour les données continues

Les modèles linéaires sont utilisés pour étudier comment une variable continue
dépend d’un ou plusieurs prédicteurs ou variables explicatives. Les prédicteurs
peuvent être quantitatifs ou qualitatifs.
Exemple : données des efforts des plannings familiaux en Amérique du Sud (Mauldin
and Berelson, 1978)
Le niveau social et les efforts des plannings familiaux sont mesurés par une
combinaison d’indices. Plus l’indice est élevé plus le niveau social (resp. l’effort) est
élevé.
Niv. social Effort Déclin du tx nat.
Bolivia 46 0 1
Brazil 74 0 10
Chile 89 16 29
Colombia 77 16 25
CostaRica 84 21 29
Cuba 89 15 40
Dominican Rep 68 14 21
Ecuador 70 6 0
El Salvador 60 13 13
. . . .
. . . .
. . . .
Dans cet exemple on cherche à comprendre comment le niveau social et les efforts
de planification influent sur le taux de natalité.
Notons Y le taux de natalité et X1 et X2 respectivement le niveau social et l’effort de
planification.
On dispose de n = 20 observations {(xi1 , xi2 , yi )}, i = 1, · · · , n.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 3 / 203
 Introduction Modèles linéaires pour les données continues

Exemple (suite)
On observe que le déclin du taux de natalité augmente avec le niveau social et avec
l’effort de planification.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 4 / 203

 Introduction Modèles linéaires pour les données continues

Exemple (suite) et vocabulaire

Si on suppose que Yi est une variable gaussienne de moyenne µi et de variance σ 2 ,

on pourra écrire le modèle

µi = β0 + β1 xi1 + β2 xi2 = β0 + xTi β

Calcul matriciel

La densité de probabilité de Yi est donnée par

!
1 1 (y − µi )2
ϕ(y ) = √ exp −
2πσ 2 2 σ2

Une manière plus standard d’écrire ce modèle linéaire est

Y = β0 + β1 X1 + β2 X2 +  = β0 + XT β + 

avec  une variable aléatoire de loi de Gauss centrée et de variable σ 2 .

X est la matrice de design et XT β est le prédicteur linéaire.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 5 / 203

 Introduction Modèles linéaires pour les données continues

Inférence

L’estimation des paramètres du modèle β0 , β, σ se fait par minimisation d’un critère

des moindres carrées (OLS)
n
X
(β0∗ , β ∗ ) = arg min (yi − β0 − xTi β)2
(β0 ,β)
i=1

σ ∗ = Var (yi − β0∗ − xTi β ∗ )

ou par maximum de vraisemblance (ML).
Dans le cas du modèle linéaire, ces deux estimateurs sont identiques et

b = (XT X)−1 XT y
β

On dispose ensuite de tests (test de Fisher, test de Wald, etc) pour valider et
interpréter le modèle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 6 / 203

 Introduction Modèles linéaires pour les données continues

Régression linéaire simple

> m1<- lm(DecTxNat ~ NivSocial,data=fpe)

> summary(m1)

Call:
lm(formula = DecTxNat ~ NivSocial, data = fpe)

Residuals:
Min 1Q Median 3Q Max
-13.239 -6.260 0.787 6.678 17.162

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -22.1254 9.6416 -2.295 0.03398 *
NivSocial 0.5052 0.1308 3.863 0.00114 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.973 on 18 degrees of freedom

Multiple R-squared: 0.4532,Adjusted R-squared: 0.4228
F-statistic: 14.92 on 1 and 18 DF, p-value: 0.001141

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 7 / 203

 Introduction Modèles linéaires pour les données continues

Régression linéaire simple

Le test t permet de tester l’hypothèse H0 : βj = 0 pour chaque variable j.

Le test de Fisher permet de tester plusieurs paramètres simultanément. C’est
notamment important quand on travaille avec des variables catégorielles recodées en
variables binaires.
Ici les deux tests conduisent à la même conclusion : la variable explicative a un effet
significatif.

> anova(m1)
Analysis of Variance Table

Response: DecTxNat
Df Sum Sq Mean Sq F value Pr(>F)
NivSocial 1 1201.1 1201.08 14.919 0.001141 **
Residuals 18 1449.1 80.51
---

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 8 / 203

 Introduction Modèles linéaires pour les données continues

Régression linéaire simple

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 9 / 203

 Introduction Modèles linéaires pour les données continues

Régression linéaire multiple

> m2 <- lm(DecTxNat ~ NivSocial+Effort,data=fpe)
> summary(m2)

Call:
lm(formula = DecTxNat ~ NivSocial + Effort, data = fpe)

Residuals:
Min 1Q Median 3Q Max
-10.3475 -3.6426 0.6384 3.2250 15.8530

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.4511 7.0938 -2.037 0.057516 .
NivSocial 0.2706 0.1079 2.507 0.022629 *
Effort 0.9677 0.2250 4.301 0.000484 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.389 on 17 degrees of freedom

Multiple R-squared: 0.7381,Adjusted R-squared: 0.7073
F-statistic: 23.96 on 2 and 17 DF, p-value: 1.132e-05

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 10 / 203

 Introduction Modèles linéaires pour les données continues

Régression linéaire simple

> anova(m2)
Analysis of Variance Table

Response: DecTxNat
Df Sum Sq Mean Sq F value Pr(>F)
NivSocial 1 1201.08 1201.08 29.421 4.557e-05 ***
Effort 1 755.12 755.12 18.497 0.0004841 ***
Residuals 17 694.01 40.82
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 11 / 203

 Introduction Modèles linéaires pour les données continues

Cas des prédicteurs qualitatifs : analyse de la variance

Dans certains cas, les prédicteurs sont tous qualitatifs. Par exemple, si le niveau
social est donné par une variable discrète à 3 modalités (inf à 70, entre 70 et 80, sup
à 80).
Les notations sont un peu différentes de celles du modèle de régression.
On note K le nombre de niveaux dans le facteur et nk le nombre d’observations dans
le niveau k . yki est la réponse de l’individu i au niveau k .
Dans l’exemple, yki sera donc le déclin du taux de natalité du pays i pour le niveau
social k . k = 1, · · · , 3, K = 3, ...
On écrit alors le modèle de l’analyse de la variance à un facteur : Yki ∼ N (µki , σ 2 )
avec
µki = µ + αk

où µ est un effet moyen (commun à taux les niveaux du facteur) et αk représente

l’effet spécifique du niveau k .

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 12 / 203

 Introduction Modèles linéaires pour les données continues

ANOVA, un facteur

> m1g <- lm(DecTxNat ~ NivSocial.g)

> summary(m1g)

Call:
lm(formula = DecTxNat ~ NivSocial.g)

Residuals:
Min 1Q Median 3Q Max
-14.750 -6.579 -1.161 5.250 16.400

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.571 3.498 2.164 0.04497 *
NivSocial.gMedium 1.029 5.420 0.190 0.85173
NivSocial.gHigh 16.179 4.790 3.377 0.00358 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.256 on 17 degrees of freedom

Multiple R-squared: 0.4505,Adjusted R-squared: 0.3858
F-statistic: 6.967 on 2 and 17 DF, p-value: 0.006167

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 13 / 203

 Introduction Modèles linéaires pour les données continues

ANOVA, un facteur

Dans ce cas, on voit bien la différence entre le test t et le test de Fisher.

Le test de Fisher a plus de sens car il considère la variable explicative dans son
ensemble et non modalité par modalité.

> anova(m1g)
Analysis of Variance Table

Response: DecTxNat
Df Sum Sq Mean Sq F value Pr(>F)
NivSocial.g 2 1193.8 596.89 6.9672 0.006167 **
Residuals 17 1456.4 85.67
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 14 / 203

 Introduction Modèles linéaires pour les données continues

Ce petit exemple permet de rappeler les différentes étapes de la modélisation linéaire

1. Analyse descriptive, visualisation des données
2. Choix du modèle : transformation de variables, variable à introduire dans le modèle
3. Inférence : estimation des paramètres du modèle, tests, intervalles de confiance
4. Sélection des variables explicatives, choix de leur paramétrisation
(continues/discrétisées)
5. Interprétation du modèle
6. Validation et comparaisons de modèle(s)

On boucle souvent sur les étapes 2 à 6 jusqu’à obtenir un bon modèle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 15 / 203

 Introduction Modèles linéaires pour les données discrètes

Outline

1 Introduction
Modèles linéaires pour les données continues
Modèles linéaires pour les données discrètes

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 16 / 203

 Introduction Modèles linéaires pour les données discrètes

Introduction d’une fonction de lien

Quand la variable à prédire est discrète, elle ne suit plus une loi de Gauss mais une
loi de Bernoulli, une loi binomiale (ou multinomiale), une loi de Poisson, etc.
De façon similaire au modèle linéaire, on va écrire que le paramètre de localisation
de la loi varie avec les prédicteurs.

Exemple des variables dichotomiques : loi de Bernoulli

Yi ∈ {0, 1} alors Yi ∼ B(πi ) avec
 
g(πi ) = β0 + XTi β ou πi = g −1 β0 + XTi β

La fonction g est une fonction de lien définie de [0, 1] sur R. En effet, β0 + XTi β ∈ R
alors que 0 < πi < 1 est une probabilité.

Exemple des variables de comptage : loi de Poisson

Yi ∈ {0, 1, 2, 3, · · · } alors Yi ∼ P(λi ) avec

exp(λi ) = β0 + XTi β

En effet, β0 + XTi β ∈ R alors que λi ∈ R+ .

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 16 / 203

 Introduction Modèles linéaires pour les données discrètes

Introduction d’une fonction de lien

Pour choisir un modèle linéaire généralisé (GLM) il faut

- choisir la loi de Y |X = x dans la famille exponentielle des GLM.
- choisir une fonction de lien inversible g.
Pour utiliser un modèle GLM il faudra donc estimer les paramètres x T β . Une fois
cette estimation réalisée, on disposera d’une estimation de η(x) ainsi que de
E(Y |X = x)) = g −1 (x T β).
A chaque loi (ou type de variable Y ), on associe une fonction de lien canonique

Modèle Logistique Log-linéaire Linéaire

Loi de Y |X = x Bernoulli Poisson Gauss
E(Y |X = x) logit(E(Y |X = x)) = x T β log(E(Y |X = x)) = x T β E(Y |X = x) =

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 17 / 203

 Introduction Modèles linéaires pour les données discrètes

Alternative : introduction d’une variable latente

Une alternative à l’approche précédente consiste à introduire une variable latente qui
suit un modèle linéaire.
Exemple de la loi de Bernoulli
On définit alors
Yi∗ = β0 + XTi β + ui

avec ui une variable à densité symétrique et de fonction de répartition F . On

n’observe pas directement Yi∗ mais on observe Yi et

Yi = 1 si Yi∗ > 0

Yi = 0 si Yi∗ ≤ 0

Les deux formulations sont équivalentes. En effet, ce second modèle permet d’écrire
la probabilité

πi = P(Yi = 1) = P(Yi∗ > 0) = P(ui > −β0 − XTi β) = 1 − F (−β0 − XTi β)

F définit de façon unique la fonction de lien.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 18 / 203

 Introduction Modèles linéaires pour les données discrètes

Données censurées (ou tronquées)

Remarque : dans le cas de données censurées on utilise naturellement la formulation

avec variable latente.
Exemple de données censurées
Un exemple typique de données censurées à droite sont des données de suivi de
cohorte (ex : cancer) avec lesquelles on cherche à estimer une durée de survie. On
n’observe pas la date de décès des patients encore vivant à la fin de l’étude.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 19 / 203

 Introduction Modèles linéaires pour les données discrètes

Inférence

Dans les modèles linéaires généralisés, on n’a pas accès aux moindres carrés (Yi∗
n’est pas observée) et l’estimation des paramètres du modèle se fait par maximum
de vraisemblance.
Comme dans le cas du modèle linéaire, on dispose d’outils statistiques pour valider
et interpréter un modèle ou pour comparer des modèles entre eux.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 20 / 203

 Plan du cours
1 Introduction
Modèles linéaires pour les données continues
Modèles linéaires pour les données discrètes
2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique
3 Inférence pour le modèle logistique
Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs
4 Diagnostiques de régression pour les données binaires
5 Variantes des modèles logistiques
Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale
6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Monbet, 12/2016 (- M2)
Inférence GLM, M2 Pharma. 21 / 203
 Regression logistique

Outline

1 Introduction

2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 22 / 203

 Regression logistique Rappels, vocabulaire

Outline

1 Introduction

2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 22 / 203

 Regression logistique Rappels, vocabulaire

Variables catégorielles

Quand la variable à expliquer est catégorielle, le modèle naturel est la régression

logistique.

Variable binaire ou dichotomique : variable qualitative qui ne peut prendre que deux
valeurs (généralement codée 0 ou 1).
Exemple : malade (oui/non), sexe (H/F), . . .

Variable polytomique : variable qualitative qui peut prendre comme valeur une
modalité parmi K .
Exemple : CSP, partis politiques, symptômes d’une maladie. . .

Variable ordinale : variable qualitative qui peut prendre comme valeur une modalité
parmi K qui sont ordonnées.
Exemple : intensité d’une douleur, classes d’âge, . . .

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 22 / 203

 Regression logistique Distribution de Bernoulli

Outline

1 Introduction

2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 23 / 203

 Regression logistique Distribution de Bernoulli

Exemple : diabète

Reaven et Miller (1979) ont collecté des données de 145 adultes non obèses
diagnostiqués "chimical diabetic", "overt diabetic" ou "normal". L’étude avait pour
objectif de déterminer des relations entre des données biologiques (prise de sang) et
le statut diabétique du patient.
Les variables mesurées sont
I rw : poids relatif, ie ratio entre le poids et le poids attendu sachant la taille
I fpg : glycémies à jeun
I ga : glycémies test (mesure de l’intolérance au glucose)
I ina : insuline durant le test (mesure de la réponse insuline au glucose administré par voie
orale )
I sspg : concentration de glucose plasmatique à l’équilibre
I cc : groupe clinique (3 =overt diabetic", 2= chemical diabetic, 1=normal)

ident rw fpg ga ina sspg cc

1 0.81 80 356 124 55 Normal
2 0.95 97 289 117 76 Normal
3 0.94 105 319 143 105 Normal
4 1.04 90 356 199 108 Normal
.. .. .. .. .. .. ..
. . . . . . .
Reaven, G. M. and Miller, R. G. (1979). An attempt to define the nature of chemical diabetes using a
multidimensional analysis. Diabetologia 16, 17-24.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 23 / 203

 Regression logistique Distribution de Bernoulli

Exemple : diabète
On observe une dépendance claire entre les glycémies test et le groupe clinique.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 24 / 203

 Regression logistique Distribution de Bernoulli

Exemple : diabète
Ajuster un modèle de régression linéaire (courbe grise) n’a pas de sens !

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 25 / 203

 Regression logistique Distribution de Bernoulli

Exemple : diabète, cas binaire

On recode la variable groupe pour se ramener à une variable binaire.
Sur la figure, on trace en rouge la proportion de patients diabétiques en fonction de la
mesure d’intolérance au glucose.
Pour l’obtenir on calcule la moyenne des yi sur des groupes de d’individus ayant un
test de glucose similaire.
On observe que la courbe rouge croit selon une forme logistique entre 0 et 1.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 26 / 203

 Regression logistique Distribution de Bernoulli

Exemple : diabète, cas binaire

Notons Y la variable de groupe (variable à
prédire).
Pour un patient i, Yi suit une loi de
Zoom sur les glycémies<800
Bernoulli de paramètre πi .
La figure montre bien que πi varie en
fonction des covariables.
En pratique, on va modéliser la relation
entre les covariables Xi et πi = E(Yi ).

g(πi ) = β0 + XTi β

La fonction g(.) = logit(.) est une fonction

de lien définie de [0, 1] sur R.
g est l’inverse de la fonction logistique (en
bleu) . Cette fonction est parfois appelée
fonction sigmoide.

exp(x)
g −1 (x) =
1 + exp(x)

Autrement dit,
Monbet, 12/2016 (- M2) T GLM, M2 Pharma. 27 / 203
 Regression logistique Lien logit

Outline

1 Introduction

2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 28 / 203

 Regression logistique Lien logit

D’où vient le lien logit ?

On a vu qu’on ne peut pas modéliser directement l’évolution de πi comme une
fonction linéaire des covariables Xi .
On va donc transformer ces probabilités pour ne plus avoir la contrainte d’être dans
un intervalle restreint.
On définit alors le ratio des cas favorables sur les cas défavorables (risque relatif).
πi
oddsi =
1 − πi

Il est parfois plus naturel de parler en terme de risque relatif plutôt que de probabilité.
Ensuite, on prend le log avant de s’affranchir de la contrainte de positivité.
πi
ηi = logit(πi ) = log
1 − πi

On note que si πi = 1/2, logit(πi ) = 0.

La transformation logit() est bijective. Son inverse est la fonction sigmoide

eηi
πi = logit −1 (ηi ) =
1 + eηi

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 28 / 203

 Regression logistique Lien logit

Le lien logit
Dans un premier temps, on peut retenir le modèle suivant
T
ex β
P(Y = 1|X = x) =
1+ e xT β

Allure de la courbe pour différentes valeurs de β

Lorsque β augmente,P(Y = 1|X = x) est souvent proche de 0 ou 1.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 29 / 203

 Regression logistique Lien logit

Odds et odds ratio

On peut être tenté de dire : plus β est grand, mieux on discrimine.

Prudence : tout dépend de l’échelle de x (si x change d’échelle, β va également
changer...)
Les coefficients du modèle logistique sont souvent interprétés en terme d’odds ratio.
L’odds (chance) pour un individu x d’obtenir la réponse Y = 1 est défini par :

P(Y = 1|X = x)
odds(x) =
1 − P(Y = 1|X = x)

L’odds ratio (rapport des odds) (rapport des chances) entre deux individus x et x̃ est

odds(x)
OR(x, x̃) =
odds(x̃)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 30 / 203

 Regression logistique Lien logit

Odds et odds ratio : exemple

Supposons qu’à un moment donné un cheval x a une probabilité p(x) = 3/4 de

perdre. Cela signifie que sur 4 courses disputées, il peut espérer en gagner une et en
perdre 3. L’odds vaut 3/1 (3 défaites contre 1 victoire, on dit également que ce cheval
est à 3 contre 1).
Pour la petite histoire, si l’espérance de gain était nulle, cela signifierait que pour 10
euros joués, on peut espérer 30 euros de bénéfice si le cheval gagne. Le rapport
p(x)/(1 − p(x)) varie entre 0 (que des victoires) et l’infini (que des défaites) en
passant par 1 (une victoire pour une défaite).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 31 / 203

 Regression logistique Lien logit

Odds et odds ratio

Il faut être prudent avec l’interprétation des OR : ils sont très souvent utilisés mais
pas toujours bien interprétés.
Comparaison de probabilités de succès entre deux individus :
OR(x, x̃) > 1 ⇔ P(Y = 1|X = x) > P(Y = 1|X = x̃)

OR(x, x̃) = 1 ⇔ P(Y = 1|X = x) = P(Y = 1|X = x̃)

OR(x, x̃) < 1 ⇔ P(Y = 1|X = x) < P(Y = 1|X = x̃)

Interprétation en termes de risque relatif : dans le cas où P(Y = 1|X = x) et

P(Y = 1|X = x̃) sont très petits par rapport à 1, on peut faire l’approximation
P(Y = 1|X = x)
OR(x, x̃) ≈
P(Y = 1|X = x̃)
et interpréter "simplement".
Mesure de l’impact d’une variable : pour le modèle logistique
pβ (x) = β0 + β1 x1 + · · · + βp xp
il est facile de vérifier que
OR(x, x̃) = exp(β1 (x1 − x̃1 )) exp(β2 (x2 − x̃2 )) · · · exp(βp (xp − x̃p ))
Si on considère 2 observation qui diffèrent seulement par la jième variable, alors
OR(x, x̃) = exp(βj (xj − x̃j )) mesure l’influence de cette variable.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 32 / 203
 Regression logistique Distribution binomiale et modèle logistique

Outline

1 Introduction

2 Regression logistique
Rappels, vocabulaire
Distribution de Bernoulli
Lien logit
Distribution binomiale et modèle logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 33 / 203

 Regression logistique Distribution binomiale et modèle logistique

Exemple : données de contraception

Little (1978) étudie des données de sondage aux Fidji portant 1607 femmes. Il
cherche à modéliser l’utilisation (ou non) d’une contraception en fonction de l’âge, du
niveau d’éducation et du désir d’enfant.
Les variables sont toutes discrètes. On observe donc le nombre de femmes utilisant
ou non une contraception dans les différents groupes.

Age Education Desires More Contraceptive use Total

Children ? No Yes
< 25 Lower Yes 53 6 59
No 10 4 14
Upper Yes 212 52 264
No 50 10 60
25-29 Lower Yes 60 14 74
No 19 10 29
Upper Yes 155 54 209
No 65 27 92
30-39 Lower Yes 112 33 145
No 77 80 157
Upper Yes 118 46 164
No 68 78 146
40-49 Lower Yes 35 6 41
No 46 48 94
Upper Yes 8 8 16
No 12 31 43
Total 1100 507 1607

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 33 / 203

 Regression logistique Distribution binomiale et modèle logistique

Loi binomiale
On peut définir

1 si la femme i utilise une contraception
Yi =
0 sinon

Yi suit une loi de Bernoulli de paramètre πi et on a

y
P(Yi = yi ) = πi i (1 − πi )1−yi

Supposons maintenant qu’on s’intéresse à un groupe i, on définit alors

Yi = nombre de femmes utilisant une contraception dans le groupe i.

Yi prend ses valeurs dans 0,1,· · · ,ni . Si les ni observations de chaque groupe sont
indépendantes et ont la même probabilité de succès πi , alors Yi suit une loi
binomiale de paramètres πi et ni : Yi ∼ B(ni , πi )
 
ni y
P(Yi = yi ) = πi i (1 − πi )ni −yi
yi

On remarque que si ni = 1 on retrouve la loi de Bernoulli.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 34 / 203

 Regression logistique Distribution binomiale et modèle logistique

Données individuelles ou agrégées

Données individuelles Données agrégées

X Y X Y n
individu désir d’enfant Contraception désir d’enfant Contraception Nombre de femmes
(oui/non) (oui/non) (oui/non) (oui/non) dans ce groupe

Loi de Bernoulli Loi de Binomiale

P(Yi = yi ) = πi Yi ∼ B(ni , πi )

E(Yi ) = πi E(Yi ) = ni πi

Var (Yi ) = πi (1 − πi ) Var (Yi ) = ni πi (1 − πi )

 
y ni
P(Yi = yi ) = πi i (1 − πi )1−yi P(Yi = yi ) =
y
πi i (1 − πi )ni −yi
yi

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 35 / 203

 Regression logistique Distribution binomiale et modèle logistique

Modèle logistique
Soient y1 , · · · , yk , k observations indépendantes d’une variable aléatoire Yi telle que
Yi ∼ B(ni , πi ), ni ≥ 1.
On suppose de plus que
logit(πi ) = XTi β

avec Xi le vecteur des covariables (incluant la variable constante égale à 1 pour

l’intercept).
Ce modèle est appelé, modèle linéaire généralisé avec réponse binomiale et lien
logit.
Les coefficients β peuvent être interprétés de façon similaire à ceux modèle linéaire,
mais il faut garder en tête qu’ils donnent la variation du logit et non directement de la
moyenne.
On remarque que
πi
= exp(XTi β)
1 − πi
Imaginons que la variable j augmente de 1 : xTi β → xTi β + βj . En prenant
T
β βj
l’exponentiel, on a exi e
T
β βj
exi e
eβj =
xT β
e i

représente un odds ratio ie l’évolution du risque relatif quand la variable augmente

de 1.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 36 / 203
 Regression logistique Distribution binomiale et modèle logistique

Modèle logistique

On peut aussi se rappeler que

T
e xi β
πi = Tβ
1 + exi
Le terme de gauche est une moyenne, comme on en a l’habitude. Mais le terme de
droite est plus difficile à interpréter.
Pour les prédicteurs continus, on peut regarder la dérivée de cette expression

∂πi
= βj πi (1 − πi )
∂xij

On voit que l’effet de la variable j dépend de la valeur du prédicteur et de celle de la

probabilité.
En pratique, dans le cas où on prédit une variable binaire, on interprète les odds ratio
(rapport de risques relatif) et/ou la probabilité associée à la moyenne des variables
explicatives dans chaque groupe.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 37 / 203

 Regression logistique Distribution binomiale et modèle logistique

Exemple des données de diabète

Dans le cas des données diabète,

logit(P(Y = ”Diabete”)) = −90 + 0.21ga

et exp(0.21) = 1.23 signifie que quand le

taux de glycémie test augmente de 1 point,
le rapport des risques relatifs augmente de
23%.
On peut aussi calculer P(Y = ”Diabete”)
pour la moyenne de la glycémie test dans
chaque groupe.

¯ Diabet ) = 1
P(Y = ”Diabete”|ga = ga
La pente de la courbe rouge est d’autant plus
forte que la partie linéaire est grande.
¯ Normal ) = 10−5
P(Y = ”Diabete”|ga = ga

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 38 / 203

 Regression logistique Distribution binomiale et modèle logistique

Identifiabilité : cas des variables qualitatives

Tout comme pour le modèle d’analyse de variance, une variable qualitative est
représentée par les indicatrices associées aux différentes modalités.
Considérons un modèle où la seule variable explicative est le sexe

π(x) = β0 + βF IF (x) + βH IH (x)

ce qui s’écrit aussi

π(x) = β0 + βF + (βH − βF )IH (x)
et il y a une infinité d’écriture possible.
La 1ère écriture correspond à une matrice de design X à 3 colonnes. Les colonnes 2
et 3 sont liées ce qui rend l’estimation impossible.
Une solution pour pallier cette difficulté consiste à mettre une contrainte sur les
coefficients βF et βH .
La solution souvent utilisée par les logiciels est de supprimer une des colonnes de la
matrice X, ce qui revient à considérer que le coefficient de la modalité associée à
cette colonne est nul. Cette modalité est alors prise comme modalité de référence.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 39 / 203

 Regression logistique Distribution binomiale et modèle logistique

Exemple du modèle à un facteur

Revenons à l’exemple des données d’usage de contraception. On considère le
facteur âge.
Chaque groupe a son propre logit

logit(πi ) = η + αi

Pour éviter la redondance, on adopte la méthode qui consiste à supposer que

α1 = 1.
Dans ce cas, η correspond au logit du groupe de référence et αi mesure la différence
de logit entre le groupe de référence et le groupe i.
Paramètre Symbole Estimation Std. err.
Constant η -1.507 0.130
25-29 α2 0.461 0.173
30-39 α3 1.048 0.154
40-49 α4 1.425 0.194
Le logit de la constante, -1.51, pour les femmes de moins de 25 ans correspond à un
odds de 0.22 ie que le modèle prédit que 22% des femmes de moins de 25 ans
utilisent une contraception. La valeur empirique est aussi 22%.
En prenant l’exponentiel des coefficients des groupes d’âge, on obtient des odds ratio
de 1.59, 2.85 et 4.16 ie que le risque d’utiliser une contraception augmente avec
l’âge jusqu’à être multiplié par 4.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 40 / 203

 Inférence pour le modèle logistique

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 41 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 41 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Méthode du maximum de vraisemblance

Soit Yi une variable aléatoire de loi f (.; beta) : P(Yi = y ) = f (y ; beta)
Soient y1 ,· · · ,yn , n réalisations indépendantes de Y1 ,· · · ,Yn
La vraisemblance du modèle f (.; beta) sachant l’échantillon y1 , · · · , yn est la
probabilité d’observer cet échantillon pour un vecteur β donné
Autrement dit,
L(β) = Pβ (Y1 = y1 , · · · , Yn = yn )

Par l’indépendance des variables Y1 ,· · · ,Yn , on a

n
Y n
Y
L(β) = Pβ (Yi = yi ) = f (yi ; βi )
i=1 i=1

Et en prenant le logarithme, on obtient

n
X n
X
log L(β) = log Pβ (Yi = yi ) = log f (yi ; beta)
i=1 i=1

On note que le log est une fonction croissante qui ne change pas la position des
optima locaux d’une fonction.
Le vecteur βb qui réalise le maximum de vraisemblance est celui qui rend l’échantillon
le plus vraisemblable.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 41 / 203
 Inférence pour le modèle logistique Maximum de vraisemblance

Comment trouver le maximum de vraisemblance

L’approche standard pour trouver le maximum d’une fonction de plusieurs variables

consiste à annuler son gradient (dérivée première) et à vérifier que son hessien
(dérivée seconde) est défini négatif.
Pour obtenir l’estimateur du maximum de vraisemblance on résout donc le système
d’équations suivant  ∂ log L(β)

 ∂β1
=0

..

 .
 ∂ log L(β)


∂β
=0
p

Ce système est appelée système d’équations du score.

C’est un système de p équations à p inconnues.
Dans le cas où les équations sont non linéaires (comme pour le modèle logistique),
ce système n’admet pas de solution analytique et on approche sa solution en utilisant
un algorithme itératif.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 42 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Propriétés de l’estimateur du maximum de vraisemblance

Pour les modèles linéaires généralisés, on peut montrer que l’estimateur du

maximum de vraisemblance existe et qu’il est unique.
Il a les propriétés suivantes
- Il est consistent : il tend vers la vraie valeur du paramètre quand le nombre
d’observations tend vers l’infini.
- Il est asymptotiquement efficace : c’est l’estimateur qui a la plus petit variance
possible, sous réserve que le nombre d’observations soit assez grand.
- Il est asymptotiquement distribué suivant une loi de Gauss : quand le nombre
d’observations tend vers l’infini, l’estimateur du maximum de vraisemblance tend, en
loi, vers une variable gaussienne.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 43 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Modèle logistique

Dans le modèle logistique, l’estimation se fait par maximum de vraisemblance.

Soient (y1 , x1 ), · · · ,(yn , xn ), n observations indépendantes du vecteur (Y , X).
La log vraisemblance est donnée par
n
X
log L(β) = yI log(πi ) + (ni − yi ) log(1 − πi )
i=1

où πi dépend des covariables xi et d’un vecteur de paramètres β.

Vraisemblance

Pour maximiser la log-vraisemblance, on peut calculer les dérivées premières et

secondes. Et utiliser l’algorithme de Newton-Raphson pour approcher l’estimateur du
maximum de vraisemblance. On obtient ainsi le score et l’information de Fisher.
Newton

En pratique, on utilise l’algorithme du "iteratively re-weighted least square" (IRLS).

On peut montrer que c’est une façon d’écrire l’algorithme de Newton-Raphson.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 44 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

IRLS pour le modèle logistique

Algorithme
Choisir un β0 initial
Répéter jusqu’à convergence
b = XT β k −1 avec k le numéro d’itération
xxxCalculer η
b = logit −1 (b
xxxCalculer µ η ), probabilité prédite de Yi = 1
dη ni
xxxCalculer zi = η̂i + (yi − µ̂i ) dµi = η̂i + µ̂ (n − (yi − µ̂i )
i i i µ̂i )
xxxRégresser z sur les covariables

β k ← (XT WX)−1 XT Wz

xxx avec une matrice de poids diagonale de terme

wii = µ̂i (ni − µ̂i )/ni

L’estimateur obtenu est consistent (ie qu’il tend en moyenne vers la vraie valeur du
paramètre quand n tend vers l’infini) et sa variance est donnée par

b = (XT WX)−1
Var (β)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 45 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Exemple sur données simulées

On suppose que pour tout Yi soit une loi de Bernoulli de paramètre
πi = 5X1i − 2X2i

avec X1 et X2 des variables aléatoires indépendantes de loi de Gauss centrée et de

variance 0.09.
On simule 500 réalisations indépendantes.
Vraisemblance en fonction de β
Scatter plot

Vraies valeurs de β (noir), estimation

Y=0 (noir), Y=1 (rouge) xxx
(jaune)
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 46 / 203
 Inférence pour le modèle logistique Maximum de vraisemblance

Exemple sur données simulées

Code R
> n = 500
> x=matrix(rnorm(2*n),n,2)
> beta = matrix(c(5,-2),2,1)*.3
> pi = exp(x%*%beta)/(1+exp(x%*%beta))
> u = runif(n)
> Y = u<pi
> mod.glm = glm(Y~x,family=binomial,control=list(trace=1))
Deviance = 534.4288 Iterations - 1
Deviance = 532.3176 Iterations - 2
Deviance = 532.3024 Iterations - 3
Deviance = 532.3024 Iterations - 4
> summary(mod.glm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.00211 0.10581 0.020 0.984
x1 4.50960 0.46758 9.645 < 2e-16 ***
x2 -2.05054 0.36946 -5.550 2.86e-08 ***
---
> pi.hat = exp(mod.glm$coefficients[1]+
+ x%*%mod.glm$coefficients[2:3])/(1+exp(x%*%mod.glm$coefficients[2:3]))
> table(pi>.5,pi.hat>.5)
FALSE TRUE
FALSE 244 5
TRUE 1 250
L’algorithme itératif converge en 4 itérations et on retrouve des valeurs de coefficients proches
des vrais paramètres.
En prédiction, 5 individus sont mal classés.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 47 / 203
 Inférence pour le modèle logistique Maximum de vraisemblance

Cas (pathologique) des nuages de points complètement séparables

(1/2)
Dans le cas où le nuage de points (y1 , x1 ) ,· · · ,(yn , xn ) est complètement séparable
l’algorithme IRLS ne converge pas.
On dit qu’un nuage de points (y1 , x1 ) ,· · · ,(yn , xn ), avec xi ∈ Rp pour tout i, est
complètement séparable s’il existe un β ∈ Rp tel que
pour tout i tel que Yi = 1, xTi β > 0 et
pour tout i tel que Yi = 0, xTi β < 0
On dit qu’un nuage de points (y1 , x1 ) ,· · · ,(yn , xn ), avec xi ∈ Rp pour tout i, est
quasi-complètement séparable s’il existe un β ∈ Rp tel que
pour tout i tel que Yi = 1, xTi β ≥ 0 et
pour tout i tel que Yi = 0, xTi β ≤ 0

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 48 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

Cas (pathologique) des nuages de points complètement séparables

(2/2)

Dans le modèle logistique, la valeur de β

b est liée à la pente de la frontière entre les
groupes.
On voit bien que dans le cas complètement séparable, on peut trouver une infinité de
frontières possibles.
Dans ce cas, l’estimateur du maximum de vraisemblance n’existe pas.
Quand on travaille avec des données réelles, il est rare que les nuages de points
soient complètement séparables.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 49 / 203

 Inférence pour le modèle logistique Maximum de vraisemblance

# On regroupe les cas de diabète

> w = which(chemdiab$cc!="Normal")
> chemdiab$cc <- ordered(chemdiab$cc, levels = c("Normal", "Diabetic"))
> chemdiab$cc[w] = "Diabetic"

# On ajuste un modèle de régression logistique

> mod = glm(cc~ga,data=chemdiab,family=binomial)
Warning message:
glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

> boxplot(chemdiab$ga~chemdiab$cc)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 50 / 203

 Inférence pour le modèle logistique Prédiction et intervalles de confiance

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 51 / 203

 Inférence pour le modèle logistique Prédiction et intervalles de confiance

Prédiction
Sachant la valeur des variables explicatives x, l’estimation de la moyenne de y est µ̂
avec g(µ̂) = xT β
b pour une fonction de lien g.
Par exemple, pour le lien logit
logit(µ̂) = xT β

On peut construire un intervalle de confiance de cet estimateur pour indiquer sa

précision. Pour obtenir cet intervalle de confiance, on a besoin de la loi de µ̂
La variance du prédicteur linéaire est xT βb est
 −1
Var (xT β)
b = xT xT W x x

Un intervalle de confiance [µ− , µ+ ] au risque α est obtenu pour µ̂ via

q
b − Φ−1 (1 − α/2) xT xT W x −1 x
logit(µ− ) = xT β

avec Φ−1 (1 − α/2) le q

quantile de la loi de Gauss centrée réduite. µ+ est obtenu en

−1
ajoutant Φ−1 (1 − α/2) xT xT W x x
Intervalle de confiance

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 51 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 52 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Déviance

Quand on ajuste un modèle, on souhaite savoir s’il décrit bien les données.
La déviance mesure un écart entre les valeurs observées yi et ni − yi et les valeurs
estimées µ̂i et ni − µ̂i .
n     
X yi ni − yi
D=2 yi log + (ni − yi ) log
µ̂i ni − µ̂i
i=1

où yi est la valeur observée et µ̂i la valeur prédite pour l’observation i.

Si l’ajustement est parfait le rapport des valeurs observées sur les valeurs prédites
est égal à 1 et son log est nul. L’ajustement est donc d’autant meilleur que la
déviance est faible.
Dans le cas de données groupées, la distribution de la déviance tend en loi vers une
variable de loi chi2 à n − p degrés de liberté quand n tend vers l’infini, où n est le
nombre de groupes et p le nombre de paramètres.
On peut en déduire un test de qualité d’ajustement.
Pour des données individuelles, la déviance ne tend pas en loi vers un chi2 et elle
sert uniquement d’indice de qualité.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 52 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Test de Pearson

Pour tester la qualité de l’ajustement, on peut aussi utiliser le test de Pearson.

Dans le cas de données binomiales, sa statistique de test s’écrit
n
X ni (yi − µ̂i )2
χ2p =
µ̂i (ni − µ̂i )
i=1

Chaque terme de la somme est un écart au carré entre la valeur observée et la

valeur prédite divisée par la variance de yi .
Pour les données groupées, la statistique χ2p suit approximativement une loi du chi 2
à n − p degrés de liberté.
Pour les données individuelles, l’approximation n’est pas de très bonne qualité et le
test de Pearson n’est donc pas considéré comme une bonne mesure.
La statistique de Pearson est aussi utilisée pour calculer des résidus.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 53 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Test d’Hosmer Lemeshow

Ce test permet de vérifier l’adéquation d’un modèle en présence de données

individuelles.
L’idée est de se rapprocher du cas de données répétées en créant ces répétitions. Le
test s’effectue de la manière suivante (voir Hosmer & Lemeshow (2000), chapitre 5
pour plus de précisions).
1. Les probabilités µi sont ordonnées par ordre croissant
2. Ces probabilités ordonnées sont ensuite séparées en K groupes de taille égale
(on prend souvent K = 10 si n est suffisamment grand). On note
- mk les effectifs du groupe k ;
- ok le nombre de succès (Y = 1) observé dans le groupe k ;
- µk la moyenne des µi dans le groupe k . La statistique de test est alors

K
X (ok − mk µk )
C2 =
mk µk (1 − µk )
k =1

Le test se conduit de manière identique au test de déviance, la statistique C 2 suivant

approximativement une loi du chi 2 à K − 1 degrés de liberté.
Sous R : hoslem.test {ResourceSelection}

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 54 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Autres critères basés sur la vraisemblance

Pour comparer des modèles, il est courant d’utiliser des critères basés sur une
pénalisation de la vraisemblance.
Le critère d’Akaike est défini par

AIC = −2 log L + 2k

avec k le nombre de paramètres à estimer.

Si l’on considère un ensemble de modèles candidats, le modèle choisi est celui qui
aura la plus faible valeur d’AIC. Ce critère repose donc sur un compromis entre la
qualité de l’ajustement et la complexité du modèle.
Le critère "Bayes Information criterion" est défini par

BIC = −2 log L + log(n)k

avec n le nombre d’observations.

L’AIC pénalise le nombre de paramètres moins fortement que le BIC.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 55 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Tests d’hypothèses pour des modèles emboités

Pour tester des modèles emboités, on utilise des tests basés sur des rapports de
vraisemblance.
Pour fixer les idées, on peut partitionner la matrice de design en deux composantes
de dimension p1 et p2  
β1
X = (X1 , X2 ) et β = .
β2

On considère le test d’hypothèse nulle

H0 : β 2 = 0

qui signifie que la composante X2 n’a pas d’effet sur la réponse.

On note D(X1 ) la déviance du modèle incluant uniquement X1 et D(X1 + X2 ) la
déviance du modèle complet. Alors la statistique

χ2 = D(X1 ) − D(X1 + X2 )

tend en loi vers une variable de loi du chi 2 à p2 degrés de liberté quand le nombre
d’observations tend vers l’infini.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 56 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Déviance pour les modèles emboités

H0 : βrw = 0
> logistic.fit1 = glm(cc~fpg,family = binomial,data=chemdiab)
> deviance(logistic.fit1)
[1] 121.944
> aic(logistic.fit1)
lik infl vari aic
-10.740278 7.470741 6.751109 36.422037
> logistic.fit2 = glm(cc~fpg+rw,family = binomial,data=chemdiab)
> deviance(logistic.fit2)
[1] 110.1438
> aic(logistic.fit2)
lik infl vari aic
-9.296681 13.896896 12.820934 46.387154
> 1-pchisq( deviance(logistic.fit1)- deviance(logistic.fit2),df=1)
[1] 0.00059
Selon les déviances, on conclut, au risque 5%, que l’ajout de la variable rw permet
d’améliorer significativement le modèle.
> p = predict(logistic.fit2,chemdiab,typ="respons")
> table(chemdiab$cc=="Diabetic",p>.5)
FALSE TRUE
FALSE 66 10
TRUE 20 49
Pourtant le critère AIC et les erreurs de prédictions sont dégradées : 30/145. ? ? ?
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 57 / 203
 Inférence pour le modèle logistique Qualité d’ajustement

Remarques sur la déviance

Dans le modèle linéaire, la déviance est égale à la somme des carrés des résidus.
Les tests de rapport de vraisemblance sont construit à partir des déviances
normalisées. Dans le cas du modèle linéaire la normalisation est σ 2 .
Pour les données binomiales, la déviance joue encore le rôle d’une somme de carrés
de résidus. Dans ce cas, la normalisation est égale à 1.
Le test de Pearson ne peut pas être généralisé pour les modèles emboités. Ceci
explique qu’on utilise plus facilement les déviances.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 58 / 203

 Inférence pour le modèle logistique Qualité d’ajustement

Tests d’hypothèses pour les paramètres

On considère le test dont l’hypothèse nulle est

H0 : βj = 0

Comme dans le cas de la régression linéaire, on utilise le test de Wald dont la

statistique de test est
β̂j
z= q .
Var
d (β̂j )

Cette statistique tend en loi vers une variable de loi de Gauss centrée et réduite
quand le nombre d’observations tend vers l’infini.
La statistique du test de Wald est aussi utilisée pour calculer des intervalles de
confiance de niveau de confiance 100(1 − α)% pour les paramètres
q
β̂j ± Φ−1 (1 − α/2) Var
d (β̂j )

avec Φ la fonction de répartition de la loi de Gauss centrée réduite 1 .

Intervalle de confiance

Les intervalles de confiance des effets en terme de logit sont obtenus en prenant les
exponentiels des bornes de l’intervalle précédent.
1. Pour α = .05, Φ−1 (1 − α/2) = 1.96
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 59 / 203
 Inférence pour le modèle logistique Qualité d’ajustement

Tests pour les paramètres

H0 : βj = 0

> logistic.fit3 = glm(cc~fpg+rw+sspg,family = binomial,data=chemdiab)

> summary(logistic.fit3)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -18.042959 3.359824 -5.370 7.86e-08 ***
fpg 0.104834 0.024475 4.283 1.84e-05 ***
rw 5.422535 2.195915 2.469 0.01353 *
ina 0.009150 0.003086 2.965 0.00303 **

Null deviance: 200.675 on 144 degrees of freedom

Residual deviance: 96.771 on 141 degrees of freedom

Ici on conclut que chacun des paramètres associé aux variables fpg, rw et ina est
significatif.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 60 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 61 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Données de contraception et test d’homogénéité

Considérons la table à 2 entrées suivante extraite de l’ensemble des données de
contraception
Yi ∼ B(ni , πi )

Desires Using Not Using All

i yi ni − yi ni
Yes 219 753 972
No 288 347 635
All 507 1100 1607
Test d’homogénéité (les 2 groupes ont la même probabilité)
Ce test correspond au modèle nul ie au modèle où les 2 groupes ont le même logit

logit(πi ) = η

# Préparation des données

> cuse <- data.frame(matrix(c(
+ 0, 219, 753,
+ 1, 288, 347), 2, 3, byrow=TRUE))
> names(cuse) <- c("nomore","using","notUsing")
> cuse$n <- cuse$using + cuse$notUsing
> cuse$Y <- cbind(cuse$using, cuse$notUsing)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 61 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Ajustement du modèle nul

Ajustement du modèle nul.

> m0 <- glm( Y ~ 1, data=cuse, family=binomial)
> m0

Call: glm(formula = Y ~ 1, family = binomial, data = cuse)

Coefficients:
(Intercept)
-0.7746

Degrees of Freedom: 1 Total (i.e. Null); 1 Residual

Null Deviance: 91.67
Residual Deviance: 91.67 AIC: 107.5

On vérifie que le paramètre (intercept) correspond à la proportion de femmes utilisant

une contraception.
> p <- sum(cuse$using)/sum(cuse$n) # probabilité empirique
> log(p/(1-p)) # logit de la probabilité empirique
[1] -0.7745545

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 62 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Ecart-types et déviance
Pour obtenir la déviance et les écart-types, on utilise la fonction summary
> summary(m0)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.77455 0.05368 -14.43 <2e-16 ***

Null deviance: 91.674 on 1 degrees of freedom

Residual deviance: 91.674 on 1 degrees of freedom
AIC: 107.54
> 1-pchisq(91.674,df=1)
[1] 0
La déviance du modèle nul est 91.67 pour 1 ddl et la pvalue associée est nulle. On
rejette donc clairement l’hypothèse selon laquelle les 2 groupes sont homogènes.
On vérifie que l’écart type du logit observé est la racine carrée de 1/succès+1/échec
> sqrt( 1/sum(cuse$using) + 1/sum(cuse$notUsing) )
[1] 0.0536794
Il peut être aussi instructif de calculer la déviance "à la main".
> obs <- cuse$Y
> fit <- cbind(p*cuse$n, (1-p)*cuse$n)
> 2*sum(obs*log(obs/fit))
[1] 91.6744

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 63 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Intervalles de confiance

On peut obtenir les intervalles de confiance.

Pour revenir en échelle "probabilité", on utilise la fonction inverse du logit plogis
> confint(m0)
2.5 % 97.5 %
-0.8804716 -0.6700014
> plogis(confint(m0))
2.5 % 97.5 %
0.2930801 0.3384965

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 64 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Odds-ratio (rapport de cotes), modèle à un facteur

On ajuste maintenant le modèle à un facteur "ne veut plus d’enfant".

logit(πi ) = η + βi

Ce modèle est saturé pour ce jeu de données (2 paramètres pour 2 probabilités) il a

donc une déviance nulle.
> m1 <- glm(Y ~ nomore, family=binomial, data=cuse)
> summary(m1)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.23499 0.07677 -16.086 <2e-16 ***
nomore 1.04863 0.11067 9.475 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Null deviance: 9.1674e+01 on 1 degrees of freedom

Residual deviance: 1.7986e-14 on 0 degrees of freedom
AIC: 17.87
La constante correspond au log-odds de "utilise une contraception" parmi les femmes
qui veulent encore des enfants et le coefficient "nomore" est la différence en log-odds
entre les 2 groupes.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 65 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Interprétation de l’odds-ratio

L’estimation de η est le logit de la proportion observée de femme utilisant une

contraception parmi les femmes qui veulent encore des enfants logit(219/972) =
-1.235. L’estimation de β est la différence entre les logits des 2 groupes,
logit(288/635) - logit(219/972) = 1.049.
En prenant l’exponentiel, on trouve un rapport de cote d’environ 3.

> exp(coef(m1)["nomore"])
nomore
2.853737
Ca ne signifie pas que "les femmes qui ne veulent plus d’enfant ont 3 fois plus de
chance d’utiliser une contraception que les autres"
C’est le rapport des cotes qui est égal à 3 et non la probabilité.
En effet, on a par définition
πi
= exp(β)
(1 − πi )

Ici on observe que les proportions sont 0.454 (= 288/635) et 0.225 (=213/972), et le
rapport est 2.01. Ainsi les femmes qui ne veulent plus d’enfant ont 2 fois plus de
chance d’utiliser une contraception que les autres.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 66 / 203

 Inférence pour le modèle logistique Exemple 1 : comparaison de 2 groupes

Test de Wald et intervalle de confiance

Modèle à un facteur
> summary(m1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.23499 0.07677 -16.086 <2e-16 ***
nomore 1.04863 0.11067 9.475 <2e-16 ***
Le test du chi2 associé à la statistique de Wald égale à 89.78 permet de rejeter
l’hypothèse selon laquelle le coefficient est nul.
> b <- coef(m1)
> se <- sqrt(diag(vcov(m1)))
> (b[2]/se[2])^2
nomore
89.77765

On peut aussi calculer un intervalle de confiance (qui ne contient pas 0) au risque

5%.
> exp(confint(m1,"nomore"))
2.5 % 97.5 %
2.298942 3.548111
> exp(confint.default(m1,"nomore"))
2.5 % 97.5 %
nomore 2.297258 3.545015
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 67 / 203
 Inférence pour le modèle logistique Exemple 2 : comparaison de plus de 2 groupes

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 68 / 203

 Inférence pour le modèle logistique Exemple 2 : comparaison de plus de 2 groupes

Données

On considère encore les données de contraception

> cuse <- data.frame(matrix(c(
+ 1, 72, 325,
+ 2, 105, 299,
+ 3, 237, 375,
+ 4, 93, 101), 4, 3, byrow=TRUE))
> names(cuse) <- c("age","using","notUsing")
> cuse$n <- cuse$using + cuse$notUsing
> cuse$age <- factor(cuse$age,
labels=c("< 25","25-29","30-39","40-49"))
> cuse
age using notUsing n
1 < 25 72 325 397
2 25-29 105 299 404
3 30-39 237 375 612
4 40-49 93 101 194

> cuse$Y <- cbind(cuse$using, cuse$notUsing)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 68 / 203

 Inférence pour le modèle logistique Exemple 2 : comparaison de plus de 2 groupes

Modèle à un facteur
On considère le modèle en fonction de l’âge, variable à 4 modalités. Il est bien sûr
saturé (4 observations pour 4 paramètres).

logit(πi ) = η + βi

avec β1 = 0. Chaque i correspond à une classe d’âge.

Ce modèle est analogue à un modèle d’analyse de la variance.
> mag <- glm( Y ~ age, family=binomial, data=cuse)
> summary(mag)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.5072 0.1303 -11.571 < 2e-16 ***
age25-29 0.4607 0.1727 2.667 0.00765 **
age30-39 1.0483 0.1544 6.788 1.14e-11 ***
age40-49 1.4246 0.1940 7.345 2.06e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Null deviance: 7.9192e+01 on 3 degrees of freedom

Residual deviance: 2.1405e-13 on 0 degrees of freedom
AIC: 32.647
La déviance du modèle nul correspond à un test de rapport de vraisemblance et
permet de rejeter l’hypothèse selon laquelle la probabilité d’utiliser une contraception
est la même dans toutes les classes d’âge.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 69 / 203
 Inférence pour le modèle logistique Exemple 2 : comparaison de plus de 2 groupes

Coefficients et test de Wald

Le logit de la modalité de référence -1.51 correspond à un odds de 0.22.
On observe que les chances de succès augmentent avec l’âge de 59% et 185%
quand on passe aux classes d’âge 25-29 et 30-39, et sont quadruplés pour la classe
40-49, en comparaison de la première classe d’âge 25.
> b <- coef(mag)
> exp(b[-1])
age25-29 age30-39 age40-49
1.585145 2.852778 4.156353
Test de Wald
H0 : β2 = β3 = β4 = 0

Statistique de test :
b T var −1 (β)
β b βb

> V <- vcov(mag)

> t(b[-1]) %*% solve(V[-1,-1]) %*% b[-1]
[,1]
[1,] 74.35663
La statistique de test vaut 74.35 pour un test du chi2 à 3 ddl. On rejette l’hypothèse
nulle.
Le test basé sur la déviance (D0 − D1 = 79) et le test de Wald ne sont pas
équivalents mais ils conduisent à la même conclusion.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 70 / 203
 Inférence pour le modèle logistique Exemple 3 : modèle à une variable

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 71 / 203

 Inférence pour le modèle logistique Exemple 3 : modèle à une variable

Données

On considère encore les données de contraception mais on traite l’âge comme une
variable continue : on attribue à chaque femme du groupe i l’âge moyen de ce
groupe.
Le modèle s’écrit
logit(πi ) = α + βxi

> cuse$agem <- c(20, 27.5, 35, 45)[as.numeric(cuse$age)]

> cuse[,c("age","agem")]
age agem
1 < 25 20.0
2 25-29 27.5
3 30-39 35.0
4 40-49 45.0

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 71 / 203

 Inférence pour le modèle logistique Exemple 3 : modèle à une variable

Ajustement du modèle
On estime alors les paramètres
> mam <- glm(Y ~ agem, family=binomial, data=cuse)
> summary(mam)
Deviance Residuals:
1 2 3 4
-0.3697 -0.3739 1.0845 -0.9750

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.672667 0.233249 -11.458 <2e-16 ***
agem 0.060671 0.007103 8.541 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Null deviance: 79.1917 on 3 degrees of freedom
Residual deviance: 2.4034 on 2 degrees of freedom
AIC: 31.05

> b <- coef(mam)

> exp(b["agem"])
agem
1.062549
Le paramètre de pente indique que le logit de la probabilité d’utiliser une
contraception augmente de 0.061 chaque année. Ca correspond à une augmentation
du rapport de cote de exp(0.06) = 1.063 soit 6.3%.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 72 / 203

 Inférence pour le modèle logistique Exemple 3 : modèle à une variable

Test de significativité

Pour tester la significativité de l’âge,

H0 : β = 0

on peut utiliser le test de Wald dont la statistique vaut 8.54 (voir page précédente).
Le carré de cette statistique peut être comparé au quantile du chi2 à un degré de
liberté.
> qchisq(0.95,df=1) # quantile
[1] 3.841459

> 1-pchisq(8.54,df=1) # pvalue

[1] 0.003474256
On peut aussi utiliser un test de rapport de vraisemblance pour comparer ce modèle
au modèle nul. La différence de déviance est égale à 76.8.
Dans les 2 cas, on rejette l’hypothèse nulle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 73 / 203

 Inférence pour le modèle logistique Exemple 3 : modèle à une variable

Test de linéarité
On peut tester formellement l’hypothèse de linéarité (ce qui est équivalent à tester la
significativité de l’âge) : H0 : β = 0
>79.2-2.4 = 76.8
La statistique vaut 76.8 pour 1 degré de liberté, elle est donc très significative. Ceci
indique qu’on ne peut pas rejeter l’hypothèse de linéarité.
Ca implique aussi qu’on peu, si on le souhaite, modéliser l’âge par une variable
quantitative et ainsi gagner 2 degrés de liberté.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 74 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

Maximum de vraisemblance
Prédiction et intervalles de confiance
Qualité d’ajustement
Exemple 1 : comparaison de 2 groupes
Exemple 2 : comparaison de plus de 2 groupes
Exemple 3 : modèle à une variable
Exemple 4 : modèle à deux préditeurs

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 75 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Données

On modélise l’évolution de la probabilité d’utiliser une contraception en fonction de

l’âge et du désir d’enfant.

Yij ∼ B(nij , πij )

i ∈ {1, · · · , 4} et j ∈ {1, 2}

Age Desires Using Not Using All

i j yij nij − yij nij
<25 Yes 58 265 323
No 14 60 74
25-29 Yes 68 215 283
No 37 84 121
30-39 Yes 79 230 309
No 158 145 303
40-49 Yes 14 43 57
No 79 58 137
Total 507 110 1607

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 75 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Table de déviance (1/3)

En construisant successivement plusieurs modèles logistiques, on construit une table
de déviance
Modèle Notation logit(πij ) Déviance dll
Nul Φ η 145.7 7
Age A η + αi 66.5 4
Désir D η + βi 54.0 6
Additif A+D η + αi + βi 16.8 3
Saturated AD η + αi + βi + (αβ)ij 0 0
On note que le modèle nul ne permet pas de décrire les données : la déviance de
145.7 avec 7 ddl est très largement supérieure à 14.1, le quantile de la loi du chi2 à 7
ddl.
L’introduction de l’âge réduit la déviance à 66.5 pour 4 ddl. La différence 145.7-66.5 =
79.2 pour 3 dll est une statistique de test pour l’effet âge. L’effet est très significatif.
On obtient exactement le même résultat que pour les données contraception-âge de
l’exemple 2.
Cette équivalence souligne une propriété des données binomiales : l’effet de l’âge
sur l’utilisation d’une contraception est complètement contenu dans les distributions
marginales contraception-âge.
Les déviances elles-mêmes varient car elles dépendent du contexte. Seule la
différence reste identique.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 76 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Table de déviance (2/3)

Modèle Notation logit(πij ) Déviance dll

Nul Φ η 145.7 7
Age A η + αi 66.5 4
Désir D η + βi 54.0 6
Additif A+D η + αi + βi 16.8 3
Saturated AD η + αi + βi + (αβ)ij 0 0

Le modèle basé sur le désir d’enfant indique de nouveau un effet significatif du désir
d’enfant sur l’utilisation d’une contraception.
Le fait que la différence de déviance pour le désir d’enfant est égale à 91.7 pour 1 ddl
alors qu’elle est seulement de 79.2 pour 3 ddl pour l’âge indique que l’effet du désir
d’enfant est plus fort que celui de l’âge.
Cependant cette comparaison est informelle : les modèles ne sont pas emboités.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 77 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Table de déviance (3/3)

Modèle Notation logit(πij ) Déviance dll
Nul Φ η 145.7 7
Age A η + αi 66.5 4
Désir D η + βi 54.0 6
Additif A+D η + αi + βi 16.8 3
Saturated AD η + αi + βi + (αβ)ij 0 0

Dans le modèle additif, on suppose α1 = β1 = 0 et on interprète les coefficients de la

façon suivante
- η est la probabilité qu’une femme de moins de 25 ans désirant plus d’enfant utilise
une contraception ;
- αi , i = 2, 3, 4 représente l’effet net de la tranche d’âge comparé à la tranche -25
ans dans la même catégorie de désir d’enfant ;
- β2 représente l’effet net du désir d’enfant dans la même tranche d’âge.
On peut comparer les déviances des modèles à un facteur à celle du modèle additif.
Par exemple, quand on passe du modèle D au modèle A + D, la déviance décroit de
37.2 et on perd 3 ddl. C’est la statistique de test observée de l’effet net de l’âge après
avoir pris en compte l’effet du désir d’enfant

H 0 : αi = 0

et le test est très significatif.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 78 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Paramètres du modèle additif

La table de déviance permet de choisir le modèle.

Ensuite, on doit interpréter les paramètres.
Les paramètres du modèle additif sont donnés ci-dessous.
Paramètre Symbole Estimation Std err z-ratio
Constant η -1.694 0.135 -12.53
Age (25-29) α2 0.368 0.175 2.10
Age (30-39) α3 0.808 0.160 5.06
Age (40-49) α4 1.023 0.204 5.01
Désir (No) β2 0.824 0.117 7.04
Les estimations des αi montrent un effet monotone de l’âge. Cet effet pourrait varier
avec la modalité de "Désir".
De même β̂2 montre un fort effet du désir d’enfant. Cet effet pourrait lui aussi varier
avec la tranche d’âge.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 79 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Modèle avec interaction

On considère maintenant le modèle qui inclut une interaction
logit(πij ) = η + αi + βj + (αβ)ij

avec α1 = β1 = (αβ)1j = (αβ)i1 = 0)

On interprète (αβ)i2 pour i = 2, 3, 4 comme l’effet additionnel ne pas désirer
d’enfant, comparé à celui de désirer un autre enfant, pour les femmes du groupe i
comparé aux femmes de moins de 25 ans.
On interprète parfois β2 + (αβ)i2 pour simplifier.
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.5193 0.1450 -10.480 < 2e-16 ***
age25-29 0.3682 0.2009 1.832 0.06691 .
age30-39 0.4507 0.1950 2.311 0.02082 *
age40-49 0.3971 0.3402 1.168 0.24298
nomoremore 0.0640 0.3303 0.194 0.84637
age25-29:nomoremore 0.2672 0.4091 0.653 0.51366
age30-39:nomoremore 1.0905 0.3733 2.921 0.00349 **
age40-49:nomoremore 1.3672 0.4834 2.828 0.00468 **
Les résultats indiquent que l’usage d’une contraception, parmi les femmes qui
désirent des enfants, varie peu avec l’âge. En revanche, l’effet de ne pas vouloir
d’enfant augmente beaucoup avec l’âge.
On peut résumer les résultats : la contraception utilisée pour espacer les naissances
ne varie pas beaucoup avec l’âge au contraire de l’usage de la contraception pour
limiter la fertilité qui augmente très fortement avec l’âge.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 80 / 203
 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Analyse de covariance

Le modèle avec interaction est saturé, ce qui signifie que le modèle reproduit les
données. Si on considère l’âge comme une variable continue, on peut espérer
trouver un modèle plus parcimonieux.
Les modèles proposés sont listés dans le tableau ci-dessous
Modèle Notation logit(πij Déviance ddl
Une ligne X α + βxi 68.88 6
Lignes parallèles X+D αj + βxi 19.99 5
Deux lignes XD αj + βj xi 9.14 4
Le dernier modèle inclut une interaction entre le désir d’enfant et l’âge, il permet donc
à l’effet du désir d’enfant de varier avec l’âge.
La réduction de déviance (de 9.9 pour 1 ddl) correspond au test (même pente)

H0 : β1 = β2

Ce test est rejeté avec une p-value de 0.002.

La déviance du dernier modèle, 9.14 pour 4 ddl, est juste en dessous du seuil critique
9.49 au risque 5%. On ne peut donc pas rejeter ce modèle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 81 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Modèle XD, choix de la paramétrisation

Ici, il est utilise de réfléchir au choix de la paramétrisation du modèle.

L’application directe de la méthode dans laquelle on choisit une modalité de
référence conduit à une paramétrisation à 4 variables : une variable égale à 1, une
variable x pour l’âge, une variable d qui prend la valeur 1 pour les femmes qui ne
veulent plus d’enfant et une variable dx égale au produit de d par x.
On obtient alors les paramètres suivants : constante et pente pour les femmes qui
veulent un autre enfant et la différence en constante et pente pour les femmes qui ne
veulent plus d’enfant. La différence peut-être difficile à interpréter.
L’alternative consiste à définir les variables : d, 1 − d, dx et (1 − d)x.
> cuse$agem <- c(20, 27.5, 35, 45)[as.numeric(cuse$age)]
> cuse$agec <- cuse$agem - 30.6
> ancova = list(
+ one = glm(Y ~ agem, family=binomial, data=cuse),
parallel = glm(Y ~ agem+nomore, family=binomial, data=cuse),
+ two = glm(Y ~ agec*nomore, family=binomial, data=cuse) )
> cuse$more <- 1 - cuse$nomore
> cuse$age.more <- cuse$agec * cuse$more
> cuse$age.nomore <- cuse$agec * cuse$nomore
> alt <- glm(Y ~ more + age.more + nomore + age.nomore - 1,
+ family=binomial, data=cuse)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 82 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Modèle XD, paramètres

On obtient finalement les estimations suivantes

Désir Age Symbole Estimation Std. err. P(>|z|)

More Constante α1 -1.194 0.079 < 2e-16 ***
Pente β1 0.0218 0.0103 0.035 *
Nomore Constante α2 -0.437 0.093 2.70e-06 ***
Pente β2 0.06981 0.01144 1.04e-09 ***

On retrouve que l’effet de l’âge est d’autant

plus fort que la femme ne veut plus
d’enfant.
En comparant les logit théoriques aux
logits empiriques, on observe qu’un
modèle quadratique serait pertinent.
En effet, l’introduction d’un effet
quadratique de l’âge en interaction avec le
désir d’enfant conduit à un très bon
ajustement.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 83 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Modèle multifacteurs
G. Rodriguez présente une analyse complète des données de contraception dans
ses notes de cours http://data.princeton.edu/wws509/notes/c3.pdf,
section 3.6.
Il propose deux approches
1. Une sélection "forward" qui consiste à partir du modèle nul et à le complexifier petit
à petit en ajoutant des effet individuels puis des interactions.
2. Une sélection "backward" qui consiste à partir du modèle contenant toutes les
interactions et à éliminer pas à pas les effets non significatifs.
Dans les deux cas, le modèle retenu est le modèle avec une interaction entre l’âge et
de désir d’enfant + un effet de l’éducation.
Dans le graphique ci dessous L="low education", U="Upper education", Y = desir for
more children, N = no desir for more children.

On observe l’interaction Y/N-age et pas

d’interaction L/U-age (lignes parallèles). L’effet
du non désir d’enfant augmente fortement avec
l’âge (pente forte).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 84 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Autres critères basés sur la vraisemblance

Pour comparer des modèles, il est courant d’utiliser des critères basés sur une
pénalisation de la vraisemblance.
Le critère d’Akaike est défini par

AIC = −2 log L + 2k

avec k le nombre de paramètres à estimer.

Si l’on considère un ensemble de modèles candidats, le modèle choisi est celui qui
aura la plus faible valeur d’AIC. Ce critère repose donc sur un compromis entre la
qualité de l’ajustement et la complexité du modèle.
Le critère "Bayes Information criterion" est défini par

BIC = −2 log L + log(n)k

avec n le nombre d’observations.

L’AIC pénalise le nombre de paramètres moins fortement que le BIC.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 85 / 203

 Inférence pour le modèle logistique Exemple 4 : modèle à deux préditeurs

Comparaison des différents modèles pour les données de contraception.

En noir pour l’âge sous forme de variable discrète, en rouge pour l’âge en variable
continue.
Les résultats des modèles BIC et AIC sont très proches pour ces modèles.
Ils conduisent à choisir le modèle quadratique.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 86 / 203

 Diagnostiques de régression pour les données binaires

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 87 / 203

 Diagnostiques de régression pour les données binaires

Diagnostiques de régression

Les diagnostiques sont aussi importants pour la régression logistique que pour la
régression classique.
Les diagnostiques se basent aussi sur les résidus ie les différences entre les valeurs
observées et les valeurs prédites.
La principale différence par rapport au modèle linéaire est que pour les données
binaires, on ne suppose plus que les données sont identiquement distribuées.
On utilise deux types de résidus : les résidus de Pearson et les résidus basés sur la
déviance.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 87 / 203

 Diagnostiques de régression pour les données binaires

Résidus de Pearson

L’approche la plus simple pour obtenir des résidus est de calculer la différence entre
les valeurs observées et les valeurs prédites et de diviser par l’écart-type des valeurs
observées
yi − µ̂i
pi = p
µ̂i (ni − µ̂i )/ni

où la µi sont les valeurs prédites et le dénominateur est donné par

µ̂ µ̂
var (yi ) = ni πi (1 − πi ) ' ni n i (1 − n i ) .
i i
Ces résidus sont appelés résidus de Pearson parce que la racine carrée de pi est la
contribution de l’individu à la statistique du chi2 de Pearson.
Pour les données groupées, les résidus de Pearson suivent approximativement une
loi de Gauss. Ce n’est pas vrai pour les données individuelles.
Dans les deux cas, un individu qui a un résidu |pi | > 2 doit nécessiter une attention
particulière.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 88 / 203

 Diagnostiques de régression pour les données binaires

Résidus basés sur la déviance

Une alternative consiste à considérer les résidus basés sur la déviance.

Dans ce cas, on s’intéresse au rapport de yi /µ̂i à la place de la différence yi − µ̂i .
Les résidus sont alors définis par
s     
yi ni − yi
di = 2 yi log + (ni − yi ) log
µ̂i n̂i − µi

Si on prend le carré et qu’on somme tout ces résidus, on obtient la déviance.

Les observations telles que di > 2 peuvent indiquer un défaut d’ajustement.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 89 / 203

 Diagnostiques de régression pour les données binaires

Exemple

> additive <- glm(Y ~ age + education + nomore, family=binomial,

+ data=cuse)
> dr <- residuals(additive)
> sum(dr^2)
[1] 29.91722 # Deviance
> pr <- residuals(additive, type="pearson")
> sum(pr^2)
[1] 28.28834
> cbind(
+ cuse[,c("age","education","nomore")],
+ obs = cuse$using/cuse$n, fit = fitted(additive),
+ dr, pr)[abs(dr)>2,]
age educ. nomore obs fit dr pr
4 <25 high 1 0.167 0.308 -2.515 -2.375
8 25-29 high 1 0.293 0.397 -2.065 -2.026
13 40-49 low 0 0.146 0.315 -2.491 -2.325

Le modèle additif est particulièrement mauvais pour les jeunes femmes éduquées qui
veulent plus d’enfants et pour les femmes plus agées qui ont un faible niveau
d’éducation.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 90 / 203

 Diagnostiques de régression pour les données binaires

Résidus studentisés
Les résidus définis ci-dessous ne sont pas standardisés.
Ils prennent en compte le fait que les observations ont des variances différentes,
mais ils ne tiennent pas compte de la variance supplémentaire due à l’estimation des
paramètres ; dans le cas du modèle linéaire les résidus studentisés incluent cette
variance supplémentaire.
On peut approcher la variance des résidus

Var (yi − µ̂i ) ≈ (1 − hii )Var (yi )

avec hii l’effet levier qui est aussi le terme diagonal de la matrice

H = X(XT WX)−1 XT W

et W une matrice diagonale de terme général wii = µi (ni − µi )/ni évaluée à

l’estimateur maximum de vraisemblance.
Ainsi, un résidu studentisé est donné par
pi
si = p
1 − hii

L’ effet levier hii mesurent l’impact de yi dans l’estimation de µ̂i

Un point i a un effet levier si Hii > 2(p + 1)/n
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 91 / 203
 Diagnostiques de régression pour les données binaires

Effet levier et influence

Les éléments diagonaux de la matrice H peuvent être interprétés comme des effets
levier ie qu’ils mesurent l’importance du rôle que joue yi dans l’estimation de µ̂i .
On pourrait calculer les distances de Cook qui comparent l’estimateur du maximum
de vraisemblance β bàβ d (i) , l’estimateur obtenu sans l’observation i. Ce calcul est
coûteux.
Mais on peut approcher ces distances par

1 d b T T b T 2 hii
Di = (β(i) − β) X WX(β
d(i) − β) ' pi
p (1 − hii )2 p

en faisant juste une itération de l’algorithme IRLS sans l’observation i en partant de

β.
b
La distance de Cook mesure l’influence d’une observation sur l’ensemble des
prévisions en prenant en compte l’effet levier et l’importance des résidus.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 92 / 203

 Diagnostiques de régression pour les données binaires

Exemple

> obs <- cuse$Y

> p <- fitted(additive)
> pfit <- p # for clarity
> pobs <- cuse$using / cuse$n
> lev <- hatvalues(additive)
> cd <- cooks.distance(additive)
> i <- order(-lev) # sort descending
> cbind(cuse[,c("age","education","nomore")],
+ pobs,pfit,lev,cd)[i[1:3],]
age educ nomore obs pfit lev cd
3 <25 high 0 0.1969697 0.1623053 0.6696332 2.385867366
7 25-29 high 0 0.2583732 0.2223899 0.5774811 0.843656898
14 40-49 low 1 0.5106383 0.5140024 0.5601446 0.002054933

Les 3 cellules qui ont potentiellement le plus d’influence sur l’estimation sont les
jeunes femmes avec un haut niveau d’éducation et qui veulent des enfants, et les
femmes plus âgées avec un faible niveau d’éducation et qui ne veulent plus d’enfant.
Le groupe le plus jeune a le plus d’influence.
Le groupe des femmes plus âgées n’a pas d’influence : distance de Cook nulle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 93 / 203

 Variantes des modèles logistiques

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 94 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 94 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Fonctions de lien

Tous les modèles étudiés jusqu’ici ont une fonction de lien logit.
En réalité, toutes les transformations qui envoient une probabilité sur la droite réelle
peuvent être utilisées, sous réserve qu’elles soient bijectives, continues et
différentiables.

En particulier, si F est la fonction de répartition d’une variable aléatoire définie sur

tout R, on peut écrire
πi = F (ηi )
ou de façon équivalente
ηi = F −1 (πi )

Les choix les plus courants sont la loi de Gauss, la loi logistique et les lois de valeurs
extrèmes.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 94 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Formulation avec variable latente

On note Yi une variable aléatoire définie sur {0, 1}.
Supposons qu’il existe une variable continue Yi∗ (définie sur R) non observable telle
que
Yi∗ > c

1 si
Yi =
0 sinon

avec c un certain niveau.

On dit que Yi∗ est la réponse latente.
Les biostatisticiens interprètent souvent Yi∗ comme une dose et Yi comme une
réponse. On parle parfois de modèle dose-réponse.
On définit
πi = P(Yi = 1) = P(Yi∗ > c)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 95 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Modèle dose-réponse
Le modèle définit ci-dessus ne change pas si on multiplie Yi∗ et c par une même
constante ou si on leur ajoute une même constante. Ainsi, pour que le modèle soit
identifiable on suppose que Yi∗ est une variable aléatoire centrée et réduite et on
pose c = 0.
Quand la réponse dépend d’un vecteur de covariables x, on peut supposer

Yi∗ = xTi β + Ui

avec Ui un terme d’erreur qui a une fonction de répartition F .

Avec ce modèle, la probabilité d’observer Yi = 1

πi = P(Yi = 1) = P(Ui > −ηi ) = 1 − F (−ηi )

avec ηi = Xi T β el prédicteur linéaire.

Si la distribution de l’erreur est symétrique autour de 0, F (u) = 1 − F (−u) et

πi = F (ηi )

Cette expression définit un modèle linéaire généralisé avec une réponse de Bernoulli
et le lien ηi = F −1 (πi )

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 96 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Modèle Probit

Le modèle probit correspond au cas où Ui ∼ N (0, 1).

On parle alors de lien probit pour ηi = Φ−1 (π)
avec Φ la fonction de répartition de la loi de Gauss standard.
On peut aussi considérer le cas plus général tel que Ui ∼ N (0, σ 2 )
!
Ui xT β
πi = P(Yi∗ > 0) = P(Ui > −xTi β) =P >− i
σ σ
! !
xTi β xTi β
= 1−Φ 1− =Φ
σ σ

On remarque qu’on ne peut pas identifier β et σ séparément. On fixe donc σ = 1.

Il y a un petit inconvénient à utiliser la distribution de Gauss pour la fonction de lien :
elle n’a pas de forme explicite.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 97 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Exemple des données de contraception

Considérons un modèle probit avec comme covariables l’âge (v.a. continue) et le

désir d’enfant et avec une interaction.
> probit.model <- glm(Y ~ agec * nomore,
+ family=binomial(link="probit"), data=cuse)
On obtient les esdtimations suivantes
Paramètre Symbol Estimation Std. Err. P(> |z|)
Constante α1 -0.72 0.046 < 2e-16
Age β1 0.01 0.006 0.033
Désir α2 − α1 0.46 0.073 3.94e-10
Age x Désir β2 − β1 0.03 0.009 0.0009
Pour interpréter ce modèle, on imagine une variable latente continue qui mesure la
motivation de la femme à utiliser une contraception.
Pour un âge moyen (30.6 ans), le fait de ne plus vouloir d’enfant accroit la motivation
de presque 0.5 écart-type.
Chaque année d’âge est associée à un accroissement de la motivation de 0.01
écart-type si la femme veut encore des enfants et à 0.03 si elle n’en veut plus.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 98 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Modèle Logistique

Une alternative au modèle Gaussien est le modèle logistique.

e ηi
πi = F (ηi ) =
1 + eηi

pour −∞ < ηi < +∞

La transfomatiopn inverse est la logit.
πi
ηi = F −1 (πi ) = log
1 − πi

La distribution logistique standard est symétrique,

elle a pour moyenne 0 et pour variance π 2 /3. Sa
forme est très proche de celle de la loi normale
mais elle a des queues plus lourdes.
Les deux sont presque linéaires entre 0.1 et 0.9 et
elles conduisent à des résultats très proches.
La valeurs des coefficients changent à cause de la
variance standard qui change.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 99 / 203

 Variantes des modèles logistiques Autres fonctions de lien

Transformation log-log complémentaire

Un autre choix est la transformation c-log-log

ηI = log (− log(1 − πi )) .

C’est la transformation inverse de la distribution de valeurs extrèmes (Gumbel)

η
F (ηi ) = 1 − e−ei .

Pour les petites valeurs de πi cette transformation est proche du logit.

La transformation log-log correspond au cas où −Ui a une distribution de valeur
extrême reverse standard −ui
F (ui ) = e−e .

Cette distribution est disymétrique avec une queue lourde vers à droite.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 100 / 203

 Variantes des modèles logistiques Loi multinomiale

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 101 / 203

 Variantes des modèles logistiques Loi multinomiale

Quand la variable à prédire est multinomiale

Dans certains problèmes, la variable à prédire est multinomiale.

Dans l’exemple des données de diabète, la variable Y prend 3 modalités "Normal",
"Chemical Diabet" et "Overt Diabet".
Dans l’exemple de données de contraceptions ci dessous, Y prend les valeurs
"Stérilisation", "Autre" et "Aucune".
Age Ster. Autre Aucune Total
15-19 3 61 232 296
20-24 80 137 400 617
25-29 216 131 301 648
30-34 268 76 203 547
40-44 150 24 164 338
45-49 91 10 183 284
Total 1005 489 1671 3165
Données du Salvador (1985) pour des femmes mariées.
On cherche à modéliser les probabilités de Y d’appartenir à chacun des groupes
sachant la classe d’âge.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 101 / 203

 Variantes des modèles logistiques Loi multinomiale

Quand la variable à prédire est multinomiale

La variable à prédire Yi prend alors ses valeurs dans {1, 2, · · · , J} avec J ≥ 2 et on

définit
XJ
πij = P(Yi = j) avec πij = 1
j=1

PJ
La contrainte j=1 πij = 1 implique que le modèle a J − 1 vecteurs de paramètres.
Loi de Y  
ni y y
P(Yij = yij , · · · , YiJ = yiJ ) = πi1i1 · · · πiJiJ
yi1 , · · · , yiJ

Le cas particulier où J = 2 correspond à la loi binomiale.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 102 / 203

 Variantes des modèles logistiques Loi multinomiale

Modèle logistique multinomial

L’approche usuelle pour les données multinomiales consiste à choisir une catégorie
de référence et à modéliser les log-odds pour les autres modalités.
πij
ηij = log = αj + xTi β j
πiJ
avec α une constante et β j un vecteur de paramètres, j = 1, · · · , J − 1.
On obtient donc le modèle suivant
T
eαj +xi βj
P(Yi = j) = PJ−1 αk +xT β
1 + j=1 e i k

1
P(Yi = J) = PJ−1 αk +xT β
1+ j=1 e
i k

On peut vérifier que

J
X
P(Yi = j) = 1
j=1

et on a bien
P(Yi = j)
ηij = log = αj + xTi β j
P(Yi = J)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 103 / 203

 Variantes des modèles logistiques Loi multinomiale

Exemple pour les données de diabète

Pour les données de diabète, on trace en pointillés les probabilités empiriques et en

trait plein les probabilités théoriques associées au modèle multinomial.

P(Yi = ”CD”)
log = −67.25+0.16∗ga
P(Yi = ”No”)

P(Yi = ”OD”)
log = −87.12+0.19∗ga
P(Yi = ”No”)

Le modèle sous estime un peu les probabilités de "Chemical Diabetic" par rapport
aux observations. mais on retrouve bien les tendances.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 104 / 203

 Variantes des modèles logistiques Loi multinomiale

Données de contraception

Pour les modèles à réponse multinomiale, on considère un autre jeu de données de

contraception.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 105 / 203

 Variantes des modèles logistiques Loi multinomiale

Modèle logistique multinomial

On peut aussi définir des modèles non linéaires.

Dans l’exemple des données de contraception, les logits sont des fonctions
quadratiques de l’âge. On pose le modèle

ηij = αj + βj ai + γj ai2

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 106 / 203

 Variantes des modèles logistiques Loi multinomiale

Inférence

Pour estimer les paramètres du modèle par maximum de vraisemblance, on

maximise la vraisemblance
 
ni y y
P(Yij = yij , · · · , YiJ = yiJ ) = πi1i1 · · · πiJiJ
yi1 , · · · , yiJ

avec les probabilités πij qui dépendent des paramètres αj et β j .

L’optimisation requiert d’utiliser des procédures numériques. Les algorithmes de
Fisher scoring et de Newton-Raphson sont généralement les plus performants.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 107 / 203

 Variantes des modèles logistiques Loi multinomiale

Exemple

Pour les données de contraception, l’estimation par maximum de vraisemblance

conduit à une déviance de 20.5 pour 8 ddl. La p-value est égale à 0.009 ce qui
montre un mauvais ajustement.
> cuse$Y <- as.matrix(cuse[,c("none","ster","other")])
> msat <- multinom(Y ~ ageg, data=cuse); msat
L’effet quadratique a un rapport de vraisemblance de 500.6 pour 4 ddl ce qui est très
significatif.
Les paramètres estimés sont
Contraste
Paramètres Ster. vs None Other vs None
Constante -12.62 -4.552
Linéaire 0.7097 0.2641
Quadratique -0.0097 -0.0047
Ils sont utilisés pour tracer la figure vue plus haut. Elle montre que l’ajustement n’est
pas si mauvais à l’exception du groupe d’âge 15-19 pour lequel la probabilité de
stérilisation est sur-estimée.
On pourrait pousser plus loin en ajoutant un effet cubique. Vous pouvez essayer car
ce modèle devrait passer les test !

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 108 / 203

 Variantes des modèles logistiques Modèle logistique conditionnel

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 109 / 203

 Variantes des modèles logistiques Modèle logistique conditionnel

Modèle logistique conditionnel

Le modèle logistique conditionnel est une extension du modèle multinomial

particulièrement bien adapté pour modéliser des comportements face à un choix.
Les variables explicatives peuvent contenir des attributs du choix (prix) et des
informations individuelles (salaire).
Soit Yi une variable discrète qui représente un choix parmi J. Soit Uij la valeur ou
l’utilité du jème choix pour le ième individu.
Les Uij sont des variables indépendantes telles que

Uij = ηij + ij

avec ηij un effet fixe et ij un effet aléatoire.

Les individus sont supposés avoir un comportement rationnel et faire leur choix de
façon à maximiser son utilité ie que l’individu i va choisir j si Uij est le plus grand
parmi Ui1 , · · · , UiJ
Mais Uij comporte un terme aléatoire, on traduit donc l’hypothèse en terme de
probabilité
πij = P(Yi = j) = P(Uij = max(Ui1 , · · · , UiJ ))

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 109 / 203

 Variantes des modèles logistiques Modèle logistique conditionnel

Modèle logistique conditionnel

On peut montrer que si ij suit une loi de valeur extrême de type I (ie Gumbel), de
densité
f () = exp (− − exp(−))

alors
exp(ηij )
πij = PJ
k =1 exp(ηik )

On reconnait l’équation qui définit le modèle multinomial.

Si J = 2, l’individu i choisit Ui1 si Ui1 − Ui2 > 0. Et on peut montrer que la différence
de deux utilités aléatoires qui suivent des lois de Gumbel indépendantes a une
distribution logistique. On retrouve le modèle logistique standard.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 110 / 203

 Variantes des modèles logistiques Modèle logistique conditionnel

Quand le modèle logistique conditionnel ne marche pas bien...

Il existe un cas de figure où le modèle multinomial ne marche pas bien appelé le "bus
rouge ou bleu".
Choix de transport : train/bus rouge/bus bleu.
Hypothèse : la moitié des gens prennent le train et l’autre le bus ; et les voyageurs qui
choisissent le bus sont indifférents à la couleur. Les probabilités de choix sont donc
π = (0.50, 0.25, 0.25) et l’espérance des utilités est η = (log 2, 0, 0).
Hypothèse supplémentaire : le service du bus bleu est interrompu. On s’attend à ce
que les voyageurs qui prennent le bus prennent tous le bus rouge et donc à une
répartition 1 :1. Mais à cause des utilités log 2 et 0 ont obtient une répartition 2 :1.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 111 / 203

 Variantes des modèles logistiques Modèle logistique conditionnel

Logits

Dans le cas des logit multinomials, les utilités moyennes ηij sont modélisées en
fonction des caractéristiques individuelles

ηij = xTi β j

Dans le cas des logit conditionnels, on modélise les utilités moyenne en fonction des
caractéristiques des alternatives.
Si zj représente un vecteur des caractéristiques de la jème alternative, alors

ηij = zTj γ

On peut combiner les deux modèles

ηij = xTi β j + zTij γ

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 112 / 203

 Variantes des modèles logistiques Modèle logistique hiérarchique

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 113 / 203

 Variantes des modèles logistiques Modèle logistique hiérarchique

Modèle logistique hiérarchique

Quand la réponse est multinomiale, on choisit une des modalités qui sert de
référence.

Une alternative consiste à hiérarchiser les modalités (ou les choix) en 2 sous
ensembles et d’ajuster une modèle logistique ordinaire pour chaque comparaison.

Pour les données de contraception, par exemple, on peut considérer

1. le rapport de cote de "utiliser une forme de contraception" contre "ne pas utiliser de
contraception" (1494 :1671)
2. le rapport de cote de stérilisation contre toutes sortes de contraception
(1005 :489).

Cette approche est utile si on considère que les individus font leur choix de façon
séquentielle. Par exemple, pour l’usage d’une contraception, la femme décide
d’abord si elle va utiliser une contraception ou non. Puis, elle doit choisir un moyen de
contraception.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 113 / 203

 Variantes des modèles logistiques Modèle logistique hiérarchique

Exemple
La figure ci-dessous montre les logarithmes des rapports de cote empiriques pour
- l’utilisation d’une contraception contre aucune contraception (u) et
- la stérilisation contre toutes sortes de contraception (s)
en fonction de l’âge.
On observe que l’usage d’une contraception augmente jusqu’à l’âge de 35 ans
environ puis décroit. Alors que l’usage spécifique de la stérilisation continue de
croître jusqu’à 50 ans.

Les données suggèrent donc le modèle

ηij = αj + βj ai + γj ai2

où ai est le centre de la classe d’âge i. j = 1 pour l’équation correspondant à l’usage

d’une contraception et j = 2 pour l’équation correspondant au choix de la méthode.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 114 / 203
 Variantes des modèles logistiques Modèle logistique hiérarchique

Inférence
Une caractéristique importante et intéressante du modèle logistique hiérarchique est
que la vraisemblance s’écrit sous la forme d’un produit de vraisemblances binomiales
qui peuvent être maximisées séparément.
Prenons de nouveau l’exemple des données de contraception. La contribution de
l’individu i à la vraisemblance est
y y y
Li = πi1i1 πi2i2 πi3i3 .

où les πij sont les probabilités et les yij sont les nombres correspondants pour les
femmes stérilisées, utilisant une autre méthode et n’utilisant aucune méthode
respectivement.
Ceci s’écrit aussi
 yi1  yi1
πi1 πi2 y
Li = (πi1 + πi2 )yi1 +yi2 πi3i3 .
πi1 + πi2 πi1 + πi2

En notant ρi1 = πi1 + πi2 la probabilité d’utiliser une contraception dans le groupe
d’âge i et ρi2 = πi1 /(πi1 + πi2 ) la probabilité conditionnelle d’être stérilisée sachant
qu’on utilise une contraception,
y y +y
Li = ρi2i1 (1 − ρi2 )yi2 ρi1i1 i2 (1 − ρi1 )yi3 .
| {z }| {z }
prob. binomiale prob. binomiale

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 115 / 203

 Variantes des modèles logistiques Modèle logistique hiérarchique

Inférence

On ajuste alors les 2 modèles logistiques séparément.

> cuse$A <- cbind(cuse[,"ster"] + cuse[,"other"], cuse[,"none"])
> sla <- glm(A ~ age + agesq, data=cuse, family=binomial)
> cuse$S.A <- cbind(cuse[,"ster"], cuse[,"other"])
> sls <- glm(S.A ~ age + agesq, data=cuse, family=binomial)
> x2 <- deviance(sla) + deviance(sls); x2
[1] 16.89298
> pchisq(x2, 8, lower.tail=FALSE)
[1] 0.03124295
On obtient
Paramètre Contraste
Use vs No use Ster vs Other
Constante -7.180 -8.869
Linéaire 0.4397 0.4942
Quadratique -.0063 -0.0056
Le modèle est assez bien ajusté avec très peu de paramètres (p-value = .03).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 116 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

Autres fonctions de lien
Loi multinomiale
Modèle logistique conditionnel
Modèle logistique hiérarchique
Modèles pour une réponse ordinale

6 Régression de Poisson

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 117 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Réponse ordinale

Exemples de réponse ordinale

- Qualité de vie d’un patient (excellente, bonne, assez bonne, mauvaise)
- Douleur (aucune, faible, forte, sévère)
- Diagnostique (absolument normal, probablement normal, équivoque, probablement
anormal, absolument anormal)
- Tendance politique (très libéral, légèrement libéral, modéré, légèrement conservatif,
très légèrement)

Les modèles multinomiaux pourraient être appliqués pour les réponses ordinales,
mais on perdrait la notion d’ordre.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 117 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèles pour réponse ordinale

Dans le cas d’une réponse ordinale, on va donc construire des modèles spécifiques.
Nous considérons l’exemple des données de diabètes. La réponse est ordonnée :
(1) "normal" < (2) "chemical" < (3) "overt".
Modéliser Yi = j peut être interprété comme la modélisation d’une variable continue
Yi∗ qui prend ses valeurs dans une suite d’intervalles

yi = j si θj−1 ≤ yi∗ < θj , j = 2, · · · , J + 1

En pratique, on va s’intéresser à la probabilité cumulée

γij = P(Yi ≤ j) = P(Yi∗ < θj ).

Et on va relier les probabilités cumulées aux variables explicatives X .

Supposons que Y ∗ = −XT β +  avec E() = 0. La moyenne de Y ∗ sachant
XT = xT est donc xT β et on a

γij = P(i ≤ θj + XTi β)

La distribution de  détermine la forme du modèle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 118 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèle cumulatif et lien logit (1/2)

Si la loi de  est logistique, on a
1
P( ≤ t) =
1 + e−t

d’où
1
γij = P(i ≤ θj + Xi β) = T β)
1 + e−(θj +Xi
et
log(γij ) = θj + XTi β, j = 1, · · · , J

Comme on écrit explicitement la constante, on ne suppose plus que la matrice de

design inclut une colonne de 1.
θj est une constante qui représente la valeur de base de la probabilité cumulée
transformée pour la catégorie j. C’est le seul paramètre qui dépend de j. En pratique
ceci implique que les rapports de cote sont proportionnels (voir exemple ci-dessous).
β représente l’effet des covariables sur la probabilité cumulée transformée.
Le rapport de cote de Y ≤ j pour deux valeurs de x est
T
eθj +x1 β
T
−x2 T )β
= e(x1
θj +x2 T β
e

L’effet est le même pour toutes les modalités j.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 119 / 203
 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèle logistique cumulatif et lien logit (2/2)

Si la fonction de lien est le logit

logit(γij ) = logit(P(Yi ≤ j))

 
P(Yi ≤ j)
= log
1 − P(Yi ≤ j)
!
πi1 + · · · πij
= log
πij+1 + · · · + πiJ

On retrouve donc un modèle de régression logistique avec une variable réponse

binaire pour laquelle une modalité est formée par les modalités 1 à j de Y et l’autre
modalité est formée des modalités j + 1 à J.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 120 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

On considère les données suivantes

http://www.stat.ufl.edu/~aa/ordinal/R_examples.pdf

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 121 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

Il est utile, quand c’est possible, de tracer les données. Par exemple, ici, on aimerait
savoir si les rapports de cotes sont proportionnels ou non.
> trauma = matrix(c(1,2,3,4,59,48,44,43,25,21,14,4,46,44,54,
+ 49,48,47,64,58,32,30,31,41),4,6)
> colnames(trauma) <- c("dose","y1","y2","y3","y4","y5")
> trauma = as.data.frame(trauma)
> trauma
dose y1 y2 y3 y4 y5
1 1 59 25 46 48 32
2 2 48 21 44 47 30
3 3 44 14 54 64 31
4 4 43 4 49 58 41
> plot(1:5,log( trauma[1,2:6]/apply(trauma[2:4,2:6],2,sum)),
+ pch=20,col="blue",ylim=c(-1.3,3),xlab="response",ylab="logit")
> lines(1:5,log( trauma[1,2:6]/apply(trauma[2:4,2:6],2,sum)),
+ col="blue")
> points(1:5, log( apply(trauma[1:2,2:6],2,sum)/apply(trauma[3:4,2:6
+ pch=18,col="green")
> lines(1:5,log(apply(trauma[1:2,2:6],2,sum)/apply(trauma[3:4,2:6],2
> points(1:5,log(apply(trauma[1:3,2:6],2,sum)/trauma[4,2:6]),
+ pch=17,col="red")
> lines(1:5,log(apply(trauma[1:3,2:6],2,sum)/trauma[4,2:6]),
+ col="red")
> grid()
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 122 / 203
 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

On observe que les rapports de cote sont assez proches de rapport proportionnels
(ils évoluent de la même façon). L’hypothèse de linéarité semble raisonnable.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 123 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

On commence par ajuster un modèle logistique à rapports de cote proportionnels
 
P(Yi ≤ j)
log = θj + xij β
1 − P(Yi ≤ j)

>library(VGAM)
> fit = vglm(cbind(y1,y2,y3,y4,y5)~dose,
family = cumulative(parallel=TRUE),data=trauma)
> summary(fit)
Call:
vglm(formula = cbind(y1, y2, y3, y4, y5) ~ dose,
family = cumulative(parallel = TRUE), data = trauma)
Pearson residuals:
logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3]) logit(P[Y<=4])
1 -0.8742 1.4325 -0.17949 -0.9222
2 -0.6845 1.2073 0.19701 -0.2754
3 -0.1411 -0.8381 -0.08163 1.2159
4 1.9782 -1.9899 0.04382 -0.1565
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.71917 0.15881 -4.528 5.94e-06 ***
(Intercept):2 -0.31860 0.15642 -2.037 0.04167 *
(Intercept):3 0.69165 0.15793 4.380 1.19e-05 ***
(Intercept):4 2.05700 0.17369 11.843 < 2e-16 ***
dose -0.17548 0.05632 -3.116 0.00183 **

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 124 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

On obtient l’ajustement suivant

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 125 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

On commence par ajuster un modèle logistique à rapports de cote proportionnels

 
P(Yi ≤ j)
log = θj + xij β
1 − P(Yi ≤ j)

On fait un test de rapport de vraisemblance pour comparer ce modèle au modèle nul.

>library(VGAM)
> fit = vglm(cbind(y1,y2,y3,y4,y5)~dose,
family = cumulative(parallel=TRUE),data=trauma)
> fit = vglm(cbind(y1,y2,y3,y4,y5)~dose,family = cumulative(parallel=TRUE),dat
> fit0 = vglm(cbind(y1,y2,y3,y4,y5)~1,family = cumulative(parallel=TRUE),data=
> statRV = -2*(logLik(fit0)-logLik(fit))
> 1-pchisq(statRV,df=length(coef(fit))-length(coef(fit0)))
[1] 0.001932651
On conclut que ...

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 126 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse, modèle multinomial saturé

Comparaison avec le modèle multinomial saturé (ie à rapports de cotes non

proportionnels)
Dans ce cas, on choisit une modalité de référence. Les rapports de cotes ne sont
donc pas les mêmes et on ne peut pas comparer directement les paramètres estimés
des deux modèles.
> trauma = matrix(c(1,2,3,4,59,48,44,43,25,21,14,4,46,44,54,49,48,
+ 47,64,58,32,30,31,41),4,6)
> count = c(t(trauma[,2:6]))
> dose = c(rep(1,5),rep(2,5),rep(3,5),rep(4,5))
> response = c(matrix(1:5,5,4))
> trauma2 = data.frame(dose=dose,response=response,count=count)
> y = factor(trauma2$response)
> msat = multinom(y~dose,data=trauma2,weight=count)
> summary(msat)
Coefficients:
(Intercept) dose
2 -0.3454744 -0.3544391
3 -0.3664693 0.1470186
4 -0.3718410 0.1945529
5 -0.8458643 0.1914717
> statRV = -2*(logLik(msat)-logLik(fit))
> statRV
’log Lik.’ 2351.399 (df=8)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 127 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse, modèle multinomial saturé

Ajustement du modèle multinomial.
On ne suppose plus que la réponse est ordinale.
Le nombre de paramètres est plus important.
> -2*logLik(msat)+2*8 # modèle multinomial saturé
’log Lik.’ 2465.145 (df=8)
> -2*logLik(fit)+2*5 # modèle à rapports de cotes proportionnels
[1] 107.7456
Avantage : facile à mettre en oeuvre.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 128 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse, rapport de cotes non proportionnels

Si on veut donner de la souplesse au modèle, il est plus naturel de comparer le
modèle à rapports de cotes proportionnels au modèle à rapports de cotes non
proportionnels qu’au modèle mutinomial.
> fit2 = vglm(cbind(y1,y2,y3,y4,y5)~dose,family = cumulative,data=trauma)
> summary(fit2)
...
Pearson residuals:
logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3]) logit(P[Y<=4])
1 0.5289 -0.48790 -0.13879 -0.35261
2 -0.2835 0.64231 0.20238 -0.06988
3 -0.8415 0.03446 -0.09066 1.07587
4 0.5233 -0.11885 0.03690 -0.67979

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.86459 0.19423 -4.451 8.53e-06 ***
(Intercept):2 -0.09374 0.17849 -0.525 0.5994
(Intercept):3 0.70625 0.17558 4.022 5.76e-05 ***
(Intercept):4 1.90867 0.23838 8.007 1.18e-15 ***
dose:1 -0.11291 0.07288 -1.549 0.1213
dose:2 -0.26890 0.06832 -3.936 8.29e-05 ***
dose:3 -0.18234 0.06385 -2.856 0.0043 **
dose:4 -0.11926 0.08470 -1.408 0.1592
> statRV = -2*(logLik(fit)-logLik(fit2))
> 1-pchisq(statRV,df=length(coef(fit2))-length(coef(fit)))
[1] 0.002487748

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 129 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse, rapport de cotes non proportionnels

On teste alors l’hypothèse des rapports de cotes proportionnels

H0 : β j = β, ∀j ∈ {1, · · · , J}

La différence des déviances suit une loi du chi2 à (J − 1)p degrés de liberté.
> pchisq(deviance(fit)-deviance(fit2),
df=df.residual(fit)-df.residual(fit2),lower.tail=FALSE)
[1] 0.002487748
Une alternative consiste à utiliser la statistique de score ie ∇β `(β̃) où β̃ est la valeur
de l’estimateur du maximum de vraisemblance.
Cette statistique suit une loi du chi2 à p(J − 2) degrés de liberté.
L’avantage de cette seconde approche est qu’on n’a pas besoin d’ajuster le modèle
sous l’hypothèse alternative.
AIC
> -2*logLik(msat)+2*8
’log Lik.’ 2465.145 (df=8)
> -2*logLik(fit)+2*5
[1] 107.7456
> -2*logLik(fit2)+2*8
[1] 99.41481

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 130 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse

L’amélioration est significative d’un point de vue statistique. Mais peut-être n’est-elle
pas vraiment utile ? pvalue = 2e-3.
> sqrt(mean(sum((fitted(fit)-trauma)^2)))
[1] 190.9919
> sqrt(mean(sum((fitted(fit2)-trauma)^2)))
[1] 190.9899
On le vérifie que la figure (rapports de cotes proportionnels : traits pleins ; rapports de
cotes non proportionnels : tirets)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 131 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèle cumulatif, autre formulation

On peut utiliser une définition basée sur une variable latente comme pour les
modèles vus précédemment.
On peut alors choisir d’autres fonctions de lien. Ici encore les liens usuels sont le lien
probit et le lien c-log-log.
> fit2 = vglm(cbind(y1,y2,y3,y4,y5)~dose,family = cumulative,data=trauma)
> summary(fit2)
...
Pearson residuals:
logit(P[Y<=1]) logit(P[Y<=2]) logit(P[Y<=3]) logit(P[Y<=4])
1 0.5289 -0.48790 -0.13879 -0.35261
2 -0.2835 0.64231 0.20238 -0.06988
3 -0.8415 0.03446 -0.09066 1.07587
4 0.5233 -0.11885 0.03690 -0.67979

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.86459 0.19423 -4.451 8.53e-06 ***
(Intercept):2 -0.09374 0.17849 -0.525 0.5994
(Intercept):3 0.70625 0.17558 4.022 5.76e-05 ***
(Intercept):4 1.90867 0.23838 8.007 1.18e-15 ***
dose:1 -0.11291 0.07288 -1.549 0.1213
dose:2 -0.26890 0.06832 -3.936 8.29e-05 ***
dose:3 -0.18234 0.06385 -2.856 0.0043 **
dose:4 -0.11926 0.08470 -1.408 0.1592
...
Residual deviance: 3.8516 on 8 degrees of freedom
> pchisq(deviance(fit)-deviance(fit2), df=df.residual(fit)-df.residual(fit2),
+lower.tail=FALSE)
[1] 0.002487748
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 132 / 203
 Variantes des modèles logistiques Modèles pour une réponse ordinale

Exemple : modèle dose-réponse, lien probit

Pour le modèle probit à rapports de cote proportionnels

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 133 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Réponse ordinale, lien c-log-log

Si on suppose que  suit une loi d’extrème, on obtient le modèle

log − log(1 − γij ) = θj + xTi β

Ce modèle s’appelle aussi le modèle de Cox pour données groupées (discrete

proportional hazards model) et il est utilisé pour les données de survie.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 134 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèle à risques adjacents

Dans le cas d’une réponse ordinale, une autre approche usuelle consiste à
considérer les rapports de cotes des modalités adjacentes
!
πij
g = θi + xTi β
πij+1

> fit.adj = vglm(cbind(y1,y2,y3,y4,y5)~dose,

family = acat(parallel=TRUE),data=trauma)
> summary(fit.adj)
Pearson residuals:
logit(P[Y=2]/P[Y=1]) logit(P[Y=3]/P[Y=2]) logit(P[Y=4]/P[Y=3]) logit(P[Y=5]/
1 1.05417 -1.2784 0.11029 0.8
2 0.80369 -1.1253 -0.25909 0.2
3 0.07975 0.8077 0.03749 -1.2
4 -2.19733 1.7747 0.09302 0.2

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -1.27326 0.15264 -8.342 < 2e-16 ***
(Intercept):2 0.93341 0.15406 6.059 1.37e-09 ***
(Intercept):3 -0.05933 0.11471 -0.517 0.60503
(Intercept):4 -0.66471 0.12614 -5.270 1.37e-07 ***
dose 0.06998 0.02263 3.092 0.00198 **

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 135 / 203

 Variantes des modèles logistiques Modèles pour une réponse ordinale

Modèle à risques séquentiels

Dans le cas d’une réponse ordinale, une autre approche usuelle consiste à
considérer les rapports de cotes des modalités "séquentielles"
!
πi1 + · · · + πij
g = θi + xTi β
πij+1
ou !
πij
g = θi + xTi β
πij+1 + · · · + πiJ

> fit.seq = vglm(cbind(y1,y2,y3,y4,y5)~dose,

+family = cratio( parallel=TRUE),data=trauma)
> summary(fit.seq)
Pearson residuals:
logit(P[Y>1|Y>=1]) logit(P[Y>2|Y>=2]) logit(P[Y>3|Y>=3]) logit(P[Y>4|Y>=4])
1 -0.121816 -1.5414 0.8151 1.2524
2 -0.008446 -1.4041 0.1730 0.4764
3 0.564893 0.5629 -0.5312 -1.3353
4 -0.454814 2.6652 -0.4582 -0.1990
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 0.83085 0.13441 6.181 6.35e-10 ***
(Intercept):2 1.82657 0.16939 10.783 < 2e-16 ***
(Intercept):3 0.27102 0.14206 1.908 0.05641 .
(Intercept):4 -0.81712 0.16030 -5.097 3.44e-07 ***
dose 0.12759 0.04418 2.888 0.00388 **
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 136 / 203
 Régression de Poisson

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 137 / 203

 Régression de Poisson Distribution de Poisson

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 137 / 203

 Régression de Poisson Distribution de Poisson

Données de comptage

Dans ce chapitre, nous allons parler des modèles log-linéaires pour les données de
comptage sous l’hypothèse d’une erreur qui suit une loi de Poisson.
Ces modèles ont de nombreuses applications pour analyser des nombres
d’évènements (nombre de cas de cancers) mais aussi pour les tables de contingence
et les données de survie.
Comme exemple, nous considérons des données d’un sondage sur la fertilité dans
les Fidji. On s’intéresse au nombre d’enfants nés de femmes mariées de race
Indienne en fonction du temps écoulé depuis leur 1er mariage (6 catégories), le lieu
de résidence (Suva, Urban, Rural) et le niveau d’éducation (none, lower primary,
upper mrimary, secondary or higher).
Dans chaque cellule du tableau ci dessous on donne la moyenne, la variance et
l’effectif observés.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 137 / 203

 Régression de Poisson Distribution de Poisson

Exemple : nombre de naissances

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 138 / 203

 Régression de Poisson Distribution de Poisson

Exemple : nombre moyen de naissances

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 139 / 203

 Régression de Poisson Distribution de Poisson

Distribution de Poisson

On dit qu’une variable Y a une distribution de Poisson de paramètre µ > 0 si elle

prend des valeurs entières y = 0, 1, 2, · · · avec la probabilité

e−µ µy
P(Y = y ) =
y!

La moyenne et la variance d’une loi de Poisson sont égales à µ

E(Y ) = var (Y ) = µ

Ainsi, dans la loi de Poisson, la variance augmente avec la moyenne. C’est une des
caractéristiques essentielle de la loi de Poisson.
Exemples connus de données qui suivent une loi de Poisson : nombre de bombes qui
ont touché Londres pendant la 2nde guerre mondiale par unité de surface,
désintégrations radioactives, échanges de chromosomes dans une cellule, nombre
de bactéries dans différentes parties d’une boite de Petri.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 140 / 203

 Régression de Poisson Distribution de Poisson

Exemple : variance nombre de naissances

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 141 / 203

 Régression de Poisson Distribution de Poisson

Propriétés de la distribution de Poisson

La distribution de Poisson peut être obtenue comme la limite d’une distribution

binomiale du nombre de succès parmi un très grand nombre d’essais de Bernoulli
avec une faible probabilité de succès.
Si Y ∼ B(n, π) alors la distribution de Y quand n tend vers l’infini et π avec µ = nπ
tend vers une loi de Poisson de paramètre µ.
Autrement dit la loi de Poisson est une loi qui compte les évènements rares.

La somme de variables aléatoires de loi de Poisson suit aussi une loi de Poisson.
Si Y1 et Y2 suivent des lois de Poisson de paramètres µ1 et µ2 alors Y1 + Y2 suit une
loi de Poisson de paramètre µ1 + µ2 .
Ce résultat se généralise à plus de 2 variables.

Une conséquence de ce résultat est qu’il est équivalent d’analyser des données
individuelles et des données groupées.
Sous une hypothèse d’indépendance, si les données individuelles Yij ∼ Pois(µi )
pour j = 1, · · · , ni alors le total du groupe Yi ∼ Pois(ni µi ).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 142 / 203

 Régression de Poisson Modèle log-linéaire

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 143 / 203

 Régression de Poisson Modèle log-linéaire

Modèle log-linéaire
Soient n observations y1 , · · · , yn considérées comme des réalisations indépendantes
de variables aléatoires de Poisson telles que Yi ∼ Pois(µi ).
On suppose de plus que la moyenne µi dépend de covariables xi .
On aurait envie de poser le modèle

µi = xTi β

Mais le terme de droite peut prendre des valeurs négatives... alors que la moyenne
d’une loi de Poisson est forcément positive.
On pose alors
log(µi ) = xTi β

Dans ce modèle βj représente le changement attendu pour le log de la moyenne par

unité de changement du prédicteur xj .
En prenant l’exponentiel, on obtient un modèle multiplicatif pour la moyenne

µi = exp(xTi β)

Quand xj augmente de 1 point, la moyenne est multipliée par exp(βj ).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 143 / 203

 Régression de Poisson Modèle log-linéaire

Effet du log

Un des effets avantageux du log est qu’il ramène les données d’évènements rares à
une échelle plus linéaire.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 144 / 203

 Régression de Poisson Données hétéroscédastiques

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 145 / 203

 Régression de Poisson Données hétéroscédastiques

Données surdispersées
Les données de nombre d’enfants étaient modélisées traditionnellement par un
modèle linéaire gaussien.
Ce choix peut être discuté car le nombre d’enfants varie de 0 à 6... il ne peut pas
suivre une loi de Gauss.
Cependant le point le plus critique, n’est pas la loi des erreurs mais le fait que la
variance n’est pas constante.
On trace la variance de chaque cellule du tableau en fonction de la moyenne en
échelle log-log.

On voit que l’hypothèse de variance constante

n’est pas du tout vérifiée.
La variance est proche de la moyenne
→ plus cohérent avec le modèle de Poisson que
le modèle Gaussien.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 145 / 203

 Régression de Poisson Données hétéroscédastiques

Données groupées et "offset"

A t-on besoin des données individuelles pour la régression de Poisson ?

NotonsPYijk ` le nombre d’enfants nés de la femme ` dans le groupe ijk et
Yijk = ` Yijk ` le total du groupe.
Si chaque observation de ce groupe est issue d’une variable aléatoire de Poisson de
moyenne µijk alors le total suit une loi de Poisson de moyenne nijk µijk
Et dans ce cas

log E(Yijk ) = log(nijk µijk ) + xTijk β

= log(nijk ) + log(µijk ) + xTijk β

Ainsi, le total du groupe suit un modèle log-linéaire avec exactement le même

coefficient β que le modèle pour les données individuelles, à l’exception de l’offset
log(nijk )
L’offset est fréquent dans les modèles log-linéaires. Il représente souvent le log d’une
mesure d’exposition (ici le nombre de femmes).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 146 / 203

 Régression de Poisson Données hétéroscédastiques

Exemple
Les données du nombre d’enfants par femme ne donnent pas d’information
individuelle. Mais on peut travailler avec les données groupées si on introduire un
offset, qui est ici le log nombre de femme dans chaque groupe.
> require(foreign)
> ceb <- read.dta("http://data.princeton.edu/wws509/datasets/ceb.dta")
> head(ceb)
i dur res educ mean var n
1 1 0-4 Suva None 0.50 1.14 8
2 2 0-4 Suva Lower primary 1.14 0.73 21
3 3 0-4 Suva Upper primary 0.90 0.67 42
>ceb$y <- round(ceb$mean*ceb$n, 0)
>ceb$os = log(ceb$n)
>m0 <- glm( y ~ offset(os), data=ceb, family=poisson)
>summary(m0)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.376346 0.009712 141.7 <2e-16 ***

Null deviance: 3731.9 on 69 degrees of freedom

Residual deviance: 3731.9 on 69 degrees of freedom
AIC: 4163.3
>exp(coef(m0))
(Intercept)
3.960403
> sum(ceb$y)/sum(ceb$n))
[1] 3.960403

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 147 / 203

 Régression de Poisson Inférence

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 148 / 203

 Régression de Poisson Inférence

Modèle log-linéaire

Le modèle décrit ci-dessous est un modèle linéaire généralisé avec une erreur
Poisson et un lien log.
La log vraisemblance de n observations de Poisson indépendantes s’écrit
n
X
log L(β) = (yi log µi − µi ) + cte
i=1

où µi dépend des covariables xi et du vecteur de paramètres β et cte est une

constante indépendante des paramètres β
On montre facilement qu’en dérivant la log-vraisemblance en fonction de β et en
annulant ces dérivées, on obtient les équations d’estimation

XT y = XT µ
b

où X est la matrice de design (incluant une colonne de 1), y est la réponse observée
et µ
b est la réponse prédite avec l’estimateur du maximum de vraisemblance.
Les équations d’estimation ont cette forme dans tous les modèles linéaires
généralisés avec la fonction de lien canonique.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 148 / 203

 Régression de Poisson Inférence

Équations d’estimation (1/2)

Les équations d’estimation du modèle de Poisson s’écrivent
XT y = XT µ
b

Considérons un exemple pour mieux comprendre : dans les données des

(U)
naissances, on considère le modèle Yi ∼ Pois(µi ) avec µi = exp(β0 + β1 Xi ) où
X (U) est une variable qui vaut 1 si la mère vit en ville et 0 sinon.
Le modèle inclut un intercept, X contient une colonne de 1 : 1. Or
n
X
1y = yi
i=1

et la première équation s’écrit

n n n
(U)
(U)
X X X
yi = exp(β0 + β1 xi ) = e β0 e β 1 xi = eβ0 (n0 + n1 eβ1 )
i=1 i=1 i=1

avec n1 le nombre de mères qui vivent en ville et n0 le nombre des autres mères.
La seconde colonne de X correspond à la variable dichotomique et l’équation
associée s’écrit donc
(U)
X X
yi = exp(β0 + β1 xi ) = eβ0 n1 eβ1 xU (i)
(U) (U)
{i|xi =1} {i|xi =1}

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 149 / 203

 Régression de Poisson Inférence

Équations d’estimation (2/2)

On obtient ainsi un système d’équations à 2 équations et deux inconnues eβ0 et eβ1

e 0 (n0 + n1 eβ1 ) = ni=1 yi = Ntotal

 β P

eβ0 n1 eβ1 = {i|xU (i)=1} yi = Nné en ville

P

En faisant le rapport des deux équations, on simplifie le terme eβ0 et on obtient

Pn
n0 + n1 eβ1 i=1 yi
Ntotal
β
= P =
n1 e 1 {i|xU (i)=1} yi N né en ville
d’où
n0 Nné en ville n0 N1
e β1 = =
n1 Ntotal − Nné en ville n1 N0

On obtient eβ0 en utilisant la seconde équation

Nné en ville Ntotal − Nné en ville N0

e β0 = = =
n1 eβ1 n0 n0

Ce résultat se généralise à plus de variables et des variables à plus de modalités.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 150 / 203

 Régression de Poisson Inférence

Estimation en pratique
Le résultat précédent se généralise à plus de variables et des variables à plus de
modalités.
Cependant, en général, l’estimation se fait en utilisant l’algorithme numérique Iterated
Re-weighted Least Square (IRLS) comme pour la régression logistique avec ici
yi − µ̂i
zi = ηi +
µ̂i

et les poids
wii = µ̂i

Cet algorithme est initialisé en prenant le log des réponses yi et en le régressant sur
les variables observées par moindres carrés. Si certaines réponses sont nulles, leur
log n’est pas défini alors on leur ajoute une petite constante.
> m1r <- glm( y ~ res + offset(os), data=ceb, family=poisson)
> summary(m1r)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.20460 0.02499 48.199 < 2e-16 ***
resUrban 0.14429 0.03245 4.447 8.72e-06 ***
resRural 0.22806 0.02783 8.194 2.52e-16 ***

Null deviance: 3731.9 on 69 degrees of freedom

Residual deviance: 3659.3 on 67 degrees of freedom
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 151 / 203
 Régression de Poisson Inférence

Goodness-of-fit
La déviance permet de mesurer l’écart entre le modèle et l’observation.
n    
X yi
D=2 yi log − (yi − µ̂i )
µ̂i
i=1

Le premier terme est identique à la déviance binomiale.

Le second terme est égal à zero puisque qu’une des propriétés de la régression de
Poisson est de reproduire les totaux marginaux (voir exemple ci-dessus).
Si n est grand, la déviance suit approximativement une loi du chi2 à n − p degrés de
liberté avec n le nombre d’observations et p le nombre de degrés de liberté.
> m1r <- glm( y ~ res + offset(os), data=ceb, family=poisson)
> summary(m1r)
Null deviance: 3731.9 on 69 degrees of freedom
Residual deviance: 3659.3 on 67 degrees of freedom
AIC: 4094.8
> qchisq(.95,df=67)
[1] 88.25016
> anova(m0,m1r)
Model 1: y ~ offset(os)
Model 2: y ~ res + offset(os)
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 69 3731.9
2 67 3659.3 2 72.572 < 2.2e-16 ***

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 152 / 203

 Régression de Poisson Inférence

Goodness-of-fit, exemple

Modèle
Nb naissance ∼ Pois(os exp(β0 + β1 Résidence))

Le test d’ajustement montre que la variable "Résidence" ne suffit pas à bien prédire
le nombre de naissances.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 153 / 203

 Régression de Poisson Inférence

Goodness-of-fit, exemple

Modèle
Nb naissance ∼ Pois(exp(β0 + β1 Résidence))

Le graphe des résidus montre qu’on explique bien la moyenne mais pas la variance.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 154 / 203

 Régression de Poisson Inférence

Autres tests d’ajustement

Pour tester la qualité d’ajustement,

H0 Le modèle permet de reproduire les observations | H1 Le modèle ne permet pas
de reproduire les observations
on peut aussi utiliser la statistique de Pearson
n
X (yi − µ̂i )2
χ2p =
µ̂i
i=1

qui suit approximativement une loi du chi2 à n − p degrés de libertés si n est grand.
> res.pearson = sum((ceb$y-m1r$fitted.values)^2/m1r$fitted.values)
[1] 3304.612
>

On peut aussi construire un test de rapport de vraisemblance ou un test de Wald.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 155 / 203

 Régression de Poisson Inférence

Déviance
On compare les différents modèles à l’aide des déviances.

On prédit le nombre d’enfants en fonction

Modèle Déviance d.d.l. de la durée du mariage, la résidence et le
Null 3731.52 69 niveau d’éducation.
Modèles à 1 facteur
Durée 165.84 64 Le modèle nul a une déviance de 3731.52
Résidence 3659.23 67 pour 69 ddl ce qui revient à rejeter
Education 2661.00 66 l’hypothèse selon laquelle le nombre
Modèles à 2 facteurs moyen d’enfants est le même pour tous les
D+R 120.68 62 groupes.
D+E 100.01 61
DR 108.84 52 L’introduction de la durée du mariage dans
DE 84.46 46 le modèle permet de faire décroitre
Modèles à 3 facteurs significativement la déviance pour une
D+R+E 70.65 59
D + RE 59.89 53
faible décroissance du nombre de ddl. Ceci
E + DR 57.06 49 reflète le fait que le nombre moyen d’enfant
R+DE 54.91 44 est fortement dépendant du temps pendant
DR+RE 44.27 43 lequel la femme est "exposée" à une
DE+RE 44.60 38 maternité possible.
DR+DE 42.72 34
DR + DE + RE 30.95 28 Il est clair qu’il n’est pas raisonnable de
considérer un modèle ne tenant pas
compte de cette variable.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 156 / 203

 Régression de Poisson Inférence

Interprétation : effet de l’éducation (1/2)

Pour tester l’effet brut de l’éducation, on

Modèle Déviance d.d.l. peut comparer la déviance du modèles
Null 3731.52 69 incluent uniquement l’éducation avec celle
Modèles à 1 facteur du modèle nul
Durée 165.84 64
Résidence 3659.23 67
Education 2661.00 66 3731 − 2661 = 1071/pour 3 ddl
Modèles à 2 facteurs
D+R 120.68 62
D+E 100.01 61 C’est très significatif !
DR 108.84 52 Mais pour ce modèle ça n’a pas de sens
DE 84.46 46 d’exclure la durée du mariage, on préfère
Modèles à 3 facteurs
D+R+E 70.65 59 donc tester l’effet de l’éducation celui de la
D + RE 59.89 53 durée étant pris en compte ie comparer les
E + DR 57.06 49 modèles D+E et D.
R+DE 54.91 44
DR+RE 44.27 43 165 − 100 = 65 pour 3 ddl
DE+RE 44.60 38
DR+DE 42.72 34 soit
DR + DE + RE 30.95 28
pvalue = 2.5e − 14

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 157 / 203

 Régression de Poisson Inférence

Interprétation : effet de l’éducation (2/2)

Modèle Déviance d.d.l.

Null 3731.52 69
Modèles à 1 facteur
Durée 165.84 64
Résidence 3659.23 67 Les femmes ayant un haut niveau
Education 2661.00 66
Modèles à 2 facteurs
d’éducation ont tendance à être plus
D+R 120.68 62 jeunes, l’effet de l’éducation est donc
D+E 100.01 61 amplifié.
DR 108.84 52
En comparant, D+R et D+E+R, la
DE 84.46 46
Modèles à 3 facteurs différence de déviance est un peu plus
D+R+E 70.65 59 faible : 50. Elle reste significative.
D + RE 59.89 53 Les femmes éduquées ont tendance à
E + DR 57.06 49
R+DE 54.91 44 vivre en ville. E et R sont donc en partie
DR+RE 44.27 43 redondantes.
DE+RE 44.60 38
DR+DE 42.72 34
DR + DE + RE 30.95 28

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 158 / 203

 Régression de Poisson Inférence

Interprétation : interactions

Modèle Déviance d.d.l.

Null 3731.52 69
Modèles à 1 facteur
Le niveau d’éducation fait-il plus de
Durée 165.84 64
Résidence 3659.23 67 différence à la campagne qu’en ville ?
Education 2661.00 66 La décroissance de déviance entre D+R+E
Modèles à 2 facteurs et D+RE est seulement de 10 environ pour
D+R 120.68 62 6 ddl ; ce n’est pas significatif (pvalue =
D+E 100.01 61 0.09).
DR 108.84 52
DE 84.46 46 L’effet de du niveau d’éducation croit-il avec
Modèles à 3 facteurs la durée du mariage ?
D+R+E 70.65 59 La décroissance de déviance entre D+R+E
D + RE 59.89 53
E + DR 57.06 49
et D+RE est seulement de 15 environ pour
R+DE 54.91 44 15 ddl ; (pvalue = 0.54).
DR+RE 44.27 43 Mêmes remarques pour les autres
DE+RE 44.60 38
interactions → on garde le modèle D+R+E.
DR+DE 42.72 34
DR + DE + RE 30.95 28

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 159 / 203

 Régression de Poisson Inférence

Modèle additif D+R+E

La constante représente le log du nombre

Param. Estim. Std Err z-ratio d’enfant pour les cellules de références.
Constant -0.1173 0.0549 -2.14 e−0.1173 = 0.89, les femmes mariées
Durée 0-4 - depuis moins de 5 ans vivant à Suva et
5-9 0.9977 0.0528 18.91 n’ayant pas d’éducation ont donc en
10-14 1.3705 0.0511 26.83
moyenne 0.89 enfants.
15-19 1.6142 0.0512 31.52
20-24 1.7855 0.0512 34.86 On observe que la durée du mariage
25-29 1.9768 0.0500 39.50 accroît le nombre moyen d’enfant par
Résidence Suva - femme. En passant de "moins de 5 ans" de
Urbaine 0.1123 0.0325 3.46
Rurale 0.1512 0.0283 5.34 mariage à "enter 5 et 9 ans", le log de la
Education Aucune - moyenne croit de 1. e0.9977 = 2.71.
Faible 0.0231 0.0227 1.02 Tout lieux de résidence et niveaux
Élevée -0.1017 0.0310 -3.28 d’éducation confondus, les femmes
Sec+ -0.3096 0.0552 -5.61 mariées depuis plus de 10 ans ont en
moyenne 7.22 plus d’enfants que les
femmes mariées depuis moins de 5 ans.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 160 / 203

 Régression de Poisson Inférence

Modèle additif D+R+E

L’effet de la résidence montre que les

femmes de Suva ont la fertilité la plus
faible. Quelque soit la durée du mariage et
Param. Estim. Std Err z-ratio le niveau d’éducation, les femmes vivant
Constant -0.1173 0.0549 -2.14
dans un autre zone urbaine ont 12% fois
Durée 0-4 -
5-9 0.9977 0.0528 18.91 plus d’enfants (exp0.1123 = 1.12), et dans
10-14 1.3705 0.0511 26.83 les zones rurales 16%.
15-19 1.6142 0.0512 31.52 Plus le niveau d’éducation est élevé, plus le
20-24 1.7855 0.0512 34.86
25-29 1.9768 0.0500 39.50 nombre moyen d’enfants par femme est
Résidence Suva - faible, quelque soit la durée du mariage et
Urbaine 0.1123 0.0325 3.46 le lieu résidence.
Rurale 0.1512 0.0283 5.34
Education Aucune -
On a vu plus haut que les interactions
Faible 0.0231 0.0227 1.02 n’apportent pas beaucoup au modèle.
Élevée -0.1017 0.0310 -3.28 Cependant, le modèle de Poisson est un
Sec+ -0.3096 0.0552 -5.61 modèle additif en échelle log. Dans
l’échelle d’origine le modèle est
multiplicatif a .
β β
a. β1 log(x1 ) + β2 log(x2 ) = log(x1 1 x2 2 )

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 161 / 203

 Régression de Poisson Inférence

Modèle de Poisson et effet multiplicatif

Un modèle additif en échelle log traduit des interactions dans l’échelle d’origine.
Prédiction du nombre moyen d’enfants
Durée du mariage 0-4 5-9 10-14 15-19 20-24 25-29
Sans éducation 0.89 2.41 3.50 4.47 5.30 6.42
Sec + 0.65 1.77 2 .57 3.28 3.89 4.71
Différence 0.24 0.64 0.93 1.19 1.41 1.71
Les femmes qui ont un niveau d’éducation secondaire ou plus ont en moyenne 27%
(1 − e−0.3096 ) plus d’enfants que les femmes "sans éducation".
Le tableau montre que la différence est croissante avec la durée du mariage. Ainsi
l’effet de l’éducation mesuré dans l’échelle d’origine dépend de la durée du mariage.
On rappelle par ailleurs que le modèle de Poisson permet à la variance d’évoluer
avec la moyenne : on a plus de variabilité dans les cellules de moyenne élevée.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 162 / 203

 Régression de Poisson Sur-dispersion

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson
Distribution de Poisson
Modèle log-linéaire
Données hétéroscédastiques
Inférence
Sur-dispersion

7 Validation, sélection de modèles

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 163 / 203

 Régression de Poisson Sur-dispersion

Sur-dispersion

Une des caractéristiques clé de la distribution de Poisson est l’égalité de la moyenne

et de la variance.
E(Y ) = Var (Y ) = µ

Cependant, les données réelles présentent parfois de la sur-dispersion ie une

variance plus importante que la moyenne.
Il existe plusieurs modèles/solutions permettant de prendre ne compte la
sur-dispersion
- modèle quasi-Poisson,
- modèle de Poisson avec estimation robuste de la variance,
- modèle à réponse binomiale négative.
Supposons que la variance est proportionnelle à la moyenne

Var (Y ) = φE(Y ) = φµ

Si φ > 1 on a une sur-dispersion et si φ < 1 un sous-dispersion (mais ce second cas

est rare en pratique).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 163 / 203

 Régression de Poisson Sur-dispersion

Inférence (1/2)
Dans le cas où la variance est proportionnelle à la moyenne Var (Y ) = φE(Y ) = φµ,
il est facile d’adapter l’algorithme IRLS.
Il suffit en effet de remplacer les poids wii par
µ̂
wii∗ =
φ

Ce sont les poids du modèle de Poisson divisés par φ. La constante φ disparait

quand on calcule
 −1
XT W X XT W z

de telle sorte que l’estimateur du maximum de vraisemblance du modèle de Poisson

est un estimateur du maximum de la quasi-vraisemblance pour le modèle avec
variance proportionnelle à la moyenne.
La variance de β b est alors
 −1
b = φ XT W X
Var (β)

avec W = diag(µ1 , · · · , µn )
Ainsi les écart-types d’estimation du modèle de Poisson sont conservatifs en cas de
sur-dispersion.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 164 / 203
 Régression de Poisson Sur-dispersion

Inférence (2/2)

On a cependant besoin d’estimer la constante φ.

L’approche classique s’appuie sur la statistique du chi2 de de Pearson
n n
X (yi − µi )2 X (yi − µi )2
χ2p = =
var (yi ) φµi
i=1 i=1

Si le modèle est correct l’espérance de cette statistique est n-p.

Par la méthode des moments, on a donc

χ2p
φ̂ =
n−p

Remarque - D’habitude une valeur importante de la statistique de Pearson traduit un

mauvais ajustement. Donc, quand on utilise la statistique il faut être assez sûr que
l’erreur d’ajustement vient bien de la sur-dispersion.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 165 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative (1/3)

Quand la sur-dispersion est due à une hétérogénéité de la moyenne, une alternative
consiste à ajouter un effet multiplicatif aléatoire pour décrire l’hétérogénéité dans le
modèle de Poisson. Ceci conduit au modèle de régression binomiale négative.
On suppose que la loi conditionnelle de Y sachant une variable non observée θ suit
une loi de Poisson de moyenne et de variance µθ :

Y ∼ P(µθ)

θ représente un facteur latent qui a tendance à faire croître Y par rapport à ce qu’on
attend sachant les covariables observées si θ > 1 et décroitre sinon.
En pratique on fixe E(θ) = 1. Si on suppose de plus que θ suit une loi Gamma de
paramètres α et γ, alors Y suit une loi binomiale négative

Γ(α + y ) γ α µy
P(Y = y ) =
y !Γ(α) (µ + β)α+y

C’est la loi qui modélise habituellement le nombre d’échec avant le k ième succès
dans une série de tirage de Bernoulli indépendants de paramètre π.

α = k , π = γ/(µ + γ)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 166 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative (2/3)

La distribution binomial négative avec α = γ = 1/σ 2 est telle que

E(Y ) et Var (Y ) = µ(1 + σ 2 µ)

SI σ = 0 on retrouve la distribution de Poisson. Si σ > 0 la variance est plus grande

que la variance dans une distribution de Poisson. La loi binomiale négative est donc
sur-dispersée par rapport à celle de Poisson.
Avec θ telle que E(Y |θ) = var (Y |θ) = θµ , E(θ) = 1 et Var (θ) = σ 2 , on a


E(Y ) = Eθ EY |θ (Y |θ) = Eθ (θµ) = µ

et


Var (Y ) = Eθ VarY |θ (Y |θ) + Varθ (Eθ (θµ))
= Eθ (θµ) + Varθ (θµ)
= µ + µ2 σ 2 = µ(1 + σ 2 µ)

On retrouve les relations précédentes sans faire d’hypothèse sur la loi de θ et on peu
en déduire des estimateurs des paramètres.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 167 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative (3/3)

La distribution de Poisson est un cas particulier de la distribution binomiale négative

avec σ 2 = 0. On peut donc construire des tests de rapport de vraisemblance pour
comparer les deux modèles.
Mais comme le modèle de Poisson (ie l’hypothèse nulle) est sur la frontière de
l’espace des paramètres, la statistique de test ne converge pas vers une loi du chi2 à
1 ddl comme pourrait s’y attendre.
On peut approcher la loi de la statistique de test par simulation ou utiliser le fait que si
on considère la statistique de test comme un chi2 à 1 ddl le test obtenu est
conservatif.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 168 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative - Exemple (1/3)

Assurance au tiers, Australie - On répertories le nombre de plaintes sur une période

de 12 mois entre 1984 et 1986 dans 176 aires géographiques en New South Wales.
Les autre variables mesurées sont le nombre d’accidents, le nombre de personnes
tuées ou accidentées et la population.

Le log du nombre de plaintes varie linéairement avec le log du nombre d’accidents.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 169 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative - Exemple (2/3)

Pour se faire une idée de la sur-dispersion, on peut estimer la moyenne et la variance
du nombre plaintes en fonction de cinq groupes du nombre d’accidents correspond à
cinq quantiles.
Nombre d’accidents
0-138 139-267 268-596 597-1810 ≥1811 Overall
Mean 30 68 178 497 2176 587
Variance 397 1206 25 622 31 751 1.8×106 1.0×106
Variance/Mean 13 18 144 64 847 1 751
Coefft of var 0.66 0.51 0.90 0.36 0.62 1.73
n 36 35 35 35 35 176

La variance est bien plus importante que la moyenne, un modèle de régression de

Poisson n’est donc pas approprié.
Le modèle de régression binomiale négative avec lien log s’écrit

Y ∼ BN(k , µ), log µ = log(n) + xT β

Dans le cas des plaintes de tiers, on propose

Y ∼ BN(k , µ), log µ = log(n) + β1 + β2 log(z)

avec z le nombre d’accidents et n la population d’une aire géographique.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 170 / 203

 Régression de Poisson Sur-dispersion

Régression binomiale négative - Exemple (3/3)

Pour le modèle

Y ∼ BN(κ, µ), log µ = log(n) + β1 + β2 log(z),

on obtient les résultats suivants
Paramètre ddl β se chi2 pvalue
Intercept 1 -6.954 0.162 1836.69 < 0.0001
log accidents 1 0.254 0.025 100.04 < 0.0001
κ 0.172 0.020

On a donc que le nombre moyen de plaintes varie avec le nombre d’accidents selon
l’équation
µ̂ = ne−6.954+0.254 log(z)

Ainsi, si le nombre d’accident augmente d’un facteur a alors le taux µ/n augmente
d’un facteur e0.254 log(a) = a0.254 .
Par exemple, si le nombre d’accidents augmente de 10%, a = 1.1, on estime un effet
sur le nombre moyen de plaintes de 1.10.254 = 1.02 ie 2%.
Si on fixe κ = 0, on retrouve le modèle de Poisson. Ce modèle a une déviance de
15836.7 pour 174 ddl, ce qui traduit un très mauvais ajustement.
En comparaison, le modèle de binomial négatif a une déviance de 192.3 pour 174
ddl.
La statistique de test du rapport de vraisemblances des deux modèles est égale à
15027. On rejette donc très largement l’hypothèse κ = 0.
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 171 / 203
 Régression de Poisson Sur-dispersion

Régression binomiale négative ou quasi-vraisemblance ?

Dans le cas des plaintes de tiers, on ajuste le modèle basé sur la
quasi-vraisemblance.
On peut alors comparer les différents modèles

Modèle Variance φ β̂1 (se) β̂2 (se)

Poisson µ φ=1 -7.094 (0.027) 0.259 (0.003)
Quasi-vrais. φµ φ̂= 91.02 -7.094 (0.258) 0.259 (0.032)
Bin. nég. µ(1 + κµ) - -6.954 (0.162) 0.254 (0.025)

Le modèle basé sur la

quasi-vraisemblance conduit aux
mêmes estimateurs que le modèle de
Poisson mais avec un écart-type plus
grand qui traduit l’inflation de la
variance.
Le modèle binomial négatif conduit
aussi à des estimations très proches.
Le modèle binomial négatif est
intuitivement plus satisfaisant car il
modélise vraiment le mécanisme de
sur-dispersion. Cependant le modèle
basé sur la quasi-vraisemblance conduit
à des résultats très similaires.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 172 / 203

 Validation, sélection de modèles

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Performance en prédiction
Estimation par validation croisée
Méthodes boostrap, "out of bag"
Sélection de variables

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 173 / 203

 Validation, sélection de modèles

Validation, sélection de modèles

Validation : Est-ce que le modèle sélectionné est "bon" ? En statistique cette question
peut être abordée de différentes façons :
- Est-ce que la qualité d’ajustement globale est satisfaisante : le modèle décrit-il bien
les valeurs observées ? Ce type de question fait l’objet des tests d’ajustement ou
d’adéquation (goodness of fit) → tests basés sur la déviance).
- Est-ce que le modèle est un bon prédicteur ? → performances en prédiction.

Sélection : étant donnés M modèles M1 , · · · , MM , comment choisir le "meilleur" à

partir de l’échantillon dont on dispose.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 173 / 203

 Validation, sélection de modèles

Performance en prédiction
Il convient de définir au préalable une règle de prévision (on se restreint dans un
premier temps au modèle logistique).
Un modèle Mβ fournit une estimation p̂β (x) = pβ̂ (x), il est naturel de définir une
règle de prévision δ̂β à partir de cette estimation :

1 si p̂β (x) ≥ s
δ̂β (x) =
0 sinon

où s ∈ [0, 1] est un seuil fixé par l’utilisateur.

Il existe plusieurs façons de choisir ce seuil, les logiciels statistiques prennent
généralement par défaut la valeur 0.5.
Etant donné δ̂ : Rp+1 → {0, 1} une règle de prévision construite à partir d’un
échantillon Dn = {(X1 , Y1 ), · · · , (Xn , Yn )}, on définit la probabilité d’erreur de δ̂ par
 
L(δ̂) = P δ̂(X ) 6= Y |Dn

Pour comparer K règles (ou modèles), l’approche consiste à estimer les probabilités
d’erreur de toutes les règles candidates à l’aide de l’échantillon puis à choisir la règle
qui possède la plus petite estimation.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 174 / 203

 Validation, sélection de modèles Performance en prédiction

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Performance en prédiction
Estimation par validation croisée
Méthodes boostrap, "out of bag"
Sélection de variables

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 175 / 203

 Validation, sélection de modèles Performance en prédiction

Estimation de l’erreur de prédiction (ou de classement)

La difficulté est de trouver un "bon" estimateur de l’erreur de prédiction L(δ̂).

Une première idée est d’utiliser
n
1X
L̂n = Iδ̂(X )6=Y
n i i
i=1

Comme le modèle saturé ajuste de manière parfait les données, cette procédure
revient à choisir de manière quasi systématique le modèle saturé.
Exemple caricatural : prédiction par le plus proche voisin.
Considérons le modèle suivant : pour prédire la classe d’un individu x, on cherche le
plus proche voisin de x parmi X1 , · · · , Xn (notons le Xi ∗ ) et on pose δ̂(x) = Yi ∗ .
Dans l’approche précédente pour estimer l’erreur de prédiction, on cherche pour
chaque Xi son plus proche voisin parmi X1 , · · · , Xn , on trouve bien sûr Xi et l’erreur
de prédiction est donc nulle.
La faiblesse de cette approche vient du fait qu’on utilise le même échantillon pour
construire le modèle et pour estimer la probabilité d’erreur. Ceci introduit un biais
dans l’estimation de la probabilité d’erreur.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 175 / 203

 Validation, sélection de modèles Performance en prédiction

Estimation de l’erreur de prédiction par apprentissage-validation

La procédure apprentissage-validation s’affranchit de ce problème en séparant de

manière aléatoire les données Dn = {(X1 , Y1 ), · · · , (Xn , Yn )} en deux parties
distinctes :
- D` = {(Xi , Yi ), i ∈ I` }, un échantillon d’apprentissage de taille ` qui est utilisé pour
ajuster les modèles ;
- Dm = {(Xi , Yi ), i ∈ Im }, un échantillon de validation de taille m = n − ` qui est
utilisé pour estimer les erreurs de prédiction

1 X
Ln (δ̂) = Iδ̂(X )6=Y ;
m i i
i∈Im

avec I` ∪ Im = {1, · · · , n} et I` ∩ Im = ∅.
Ln (δ̂) est un estimateur sans biais de L(δ̂).
Exemple : on montre que l’erreur de prédiction est sous estimée par l’estimation en
resubstitution.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 176 / 203

 Validation, sélection de modèles Performance en prédiction

Estimation de l’erreur de prédiction par apprentissage-validation

Exemple
> file = ’http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.csv’
> df=read.table(file, sep=’,’, header=T)
> df$www = as.factor(df$www) ; df$male = as.factor(df$male) ; df$trust = as.fa
> head(df)
trust belief educate income age male www
1 0 2 15 9.0 48 0 0
2 1 3 14 27.5 39 0 1
3 0 0 14 27.5 25 0 1
> appri = sample(nrow(df),nrow(df)*2/3) ; testi = setdiff(1:nrow(df),appri)
> blm<-glm(trust~educate+income+age+male+www,data=df, family=binomial)
> predr = predict(blm,df[testi,],type="response")
> (c.r=table(predr>.5,df$trust[testi]))
0 1
FALSE 168 94
TRUE 50 80
> (err.r = 1-sum(diag(c.r))/length(testi))
[1] 0.3673469
> blm<-glm(trust~educate+income+age+male+www,data=df,
family=binomial,subset=appri)
> predcv = predict(blm,df[testi,],type="response")
> (c.cv = table(predcv>.5,df$trust[testi]))
0 1
FALSE 174 102
TRUE 44 72
> (err.cv = 1-sum(diag(c.cv))/length(testi))
[1] 0.372449

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 177 / 203

 Validation, sélection de modèles Performance en prédiction

Matrice de confusion, sensibilité, spécificité

Pour certaines applications, on évalue d’autres indices que l’erreur de classement.
Pour un seuil s donné, on construit par exemple une matrice de confusion
Prévision Observation Total
Y=1 Y=0
ŷi = 1 n11 (s) n01 (s) n1+ (s)
ŷi = 0 n01 (s) n00 (s) n0+ (s)
Total n+1 n+0 n
- sensibilité (ou taux de vrais positifs) : n11 (s)/n+1 (s)
- spécificité (ou taux de vrais négatifs) : n00 (s)/n+0 (s)
Il est possible de choisir un seuil s qui permette d’atteindre une sensibilité fixé ou une
spécificité fixée.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 178 / 203

 Validation, sélection de modèles Performance en prédiction

Liste d’indices

Observation
Prévision event no event
event A B
no event C D

Sensitivity = A/(A+C)
Specificity = D/(B+D)
Prevalence = (A+C)/(A+B+C+D)
PPV = (sensitivity * Prevalence)/((sensitivity*Prevalence) +
((1-specificity)*(1-Prevalence)))
NPV = (specificity * (1-Prevalence))/(((1-sensitivity)*Prevalence) +
((specificity)*(1-Prevalence)))
Detection Rate = A/(A+B+C+D)
Detection Prevalence = (A+B)/(A+B+C+D)
Balanced Accuracy = (Sensitivity+Specificity)/2

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 179 / 203

 Validation, sélection de modèles Performance en prédiction

Courbe ROC
Tous les indices listés requièrent le choix d’un seuil.
Il existe des indicateurs plus flexibles (toujours basés sur la prévision) qui n’imposent
pas de fixer le seuil.
La courbe ROC fait partie de ces critères.
C’est une courbe paramétrée par le seuil

x(s) = 1 − sp(s)
y (s) = se(s)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 180 / 203

 Validation, sélection de modèles Estimation par validation croisée

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Performance en prédiction
Estimation par validation croisée
Méthodes boostrap, "out of bag"
Sélection de variables

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 181 / 203

 Validation, sélection de modèles Estimation par validation croisée

K fold cross validation

On a vu que a validation croisée est une bonne méthode pour estimer une erreur de
prédiction.
Le plus souvent, on combine plusieurs estimations en apprentissage-validation.
1 2 3 4 5

Train Train Validation Train Train

K
1 X 1 X
CV (f̂ ) = Iδ̂(X )6=Y
K [n/K ] i i
k =1 i∈Ik

Quelle valeur de K choisir ?

I Leave-one-out : K = N → estimation de ErrA ; dans ce cas l’estimateur est
approximativement sans biais, mais peut avoir une grande variance.
I K = 5 ou 10 → estimation de Err ; les ensembles d’apprentissages sont assez différents de
l’ensemble d’origine. La variance est faible, mais le biais peut-être un problème si chaque
ensemble est "petit".

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 181 / 203

 Validation, sélection de modèles Estimation par validation croisée

Validation croisée
P10
Exemple : X suit une loi uniforme sur [0, 1]×20 et Y = 1 si p=1 X > 5 et 0 sinon.
L’ensemble d’apprentissage est composé de 80 individus.
On estime l’erreur de prédiction (en noir) sur un ensemble de test très grand et on la
compare à l’erreur estimée par validation croisée 10 folds (en rouge). Les bâtons
représentent un intervalle de plus ou moins un écart-type de l’erreur moyenne
individuelle de chaque sous ensemble.

Erreur de classement moyenne

0.4
0.3
0.2
0.1
0.0

5 10 15

Nombre de prédicteurs

L’estimation de l’erreur par validation croisée 10-folds sous-estime l’erreur de

prédiction pour les petites valeurs de d mais elle identifie bien le nombre de
composantes.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 182 / 203

 Validation, sélection de modèles Estimation par validation croisée

Estimation de l’erreur de prédiction, K-fold

Exemple
> library(boot)
> blm<-glm(trust~educate+income+age+male+www, data=df, family=binomial(link="l
> predr = predict(blm,df[testi,],type="response")
> c.r=table(predr>.5,df$trust[testi])
> (err.r = 1-sum(diag(c.r))/length(testi))
[1] 0.3673469
> mean(abs((as.numeric(trust)-1)-prev_prob)>.5)
[1] 0.3673469
> cout = function(trust,prev_prob){
c = mean(abs(trust-prev_prob)>.5)
return(c)
}
> cv.glm(df,blm,cout)$delta # leave-one-out
[1] 0.3492334 0.3498624
> cv.glm(df,blm,cout,K=5)$delta
[1] 0.3509370 0.3483838

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 183 / 203

 Validation, sélection de modèles Estimation par validation croisée

Courbe ROC, K-fold

Exemple
K = 10
nk = floor(nrow(df)/K)
pred = NULL
Yobs = NULL
rk=list()
ii = sample(nrow(df),nrow(df)) # permutation aléatoire des indices
for (k in 1:K){
testi = ii[((k-1)*nk+1):(k*nk)]
appri = setdiff(1:nrow(df),testi)
blm<-glm(trust~educate+income+age+male+www, data=df, family=binomial(link="log
pr = predict(blm,df[testi,],type="response")

pred = c(pred,pr)
Yobs = c(Yobs,df$trust[testi])
rk[[k]]=roc(df$trust[testi],pr)
}
add = FALSE
for (k in 1:K){
if (k>1) {add=TRUE}
plot.roc(rk[[k]],col="gray",add=add,lwd=1)
}
r=roc(Yobs,pred)
plot(r,main=paste("AUROC =",round(r$auc*100)/100),add=TRUE)
grid()

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 184 / 203

 Validation, sélection de modèles Estimation par validation croisée

Courbe ROC, K-fold

On peut, par exemple, tracer les K courbe ROC et leur moyenne.

On a ainsi une idée de la variabilité des résultats.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 185 / 203

 Validation, sélection de modèles Estimation par validation croisée

Bien utiliser la validation croisée

Un processus de sélection/validation de modèle comporte généralement plusieurs

étapes :
1. Sélection d’un sous ensemble de "bons" prédicteurs ;
2. Construction d’un modèle basé sur ce sous ensemble et estimation des
paramètres ;
3. Estimation de l’erreur de prédiction.
La validation croisée doit être appliquée à l’ensemble du processus ie qu’on tire
aléatoirement les sous échantillons avant de réaliser les étapes 1 à 3.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 186 / 203

 Validation, sélection de modèles Méthodes boostrap, "out of bag"

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Performance en prédiction
Estimation par validation croisée
Méthodes boostrap, "out of bag"
Sélection de variables

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 187 / 203

 Validation, sélection de modèles Méthodes boostrap, "out of bag"

Méthodes bootstrap

L’objectif du bootstrap est d’estimer Ln (δ̂) mais il estime bien l’erreur de prédiction
attendue L(δ̂).
Algorithme
Pour b = 1,· · · , B
xxxx Tirer aléatoirement avec remise un ensemble Ib de taille napp
xxxxxx dans l’ensemble des données
xxxx Sélectionner les prédicteurs
xxxx Ajuster le modèle
xxxx Estimer l’erreur de prédiction avec les individus qui ne sont pas dans Ib .
fin
Calculer la moyenne des erreurs de prédiction
En pratique on choisi souvent B = 100.
On peut aussi choisir de faire un tirage sans remise d’un ensemble de taille
prédéfinie.
Si la taille n de l’échantillon est assez grande, on peut le décomposer aléatoirement
en 3 sous ensembles
I Ensemble d’apprentissage [50%] → estimation des paramètres
I Ensemble de validation [25%] → sélection de variables
I Ensemble de test [25%] → estimation de Ln (δ̂)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 187 / 203

 Validation, sélection de modèles Méthodes boostrap, "out of bag"

Comparaison de modèles
library(rpart)
data(kyphosis)
summary(kyphosis)
n = nrow(kyphosis)
ntest = round(n/5)
p.add <- p.mul <- Yobs <- NULL
B = 100
for (b in 1:B){
testi = sample(1:n,ntest)
appri = setdiff(1:n,testi)
mod = glm(Kyphosis~.,data=kyphosis,family=binomial,subset=appri)
pmod = predict(mod,kyphosis[testi,])
mod2 = glm(Kyphosis~.^2,data=kyphosis,family=binomial,subset=appri)
pmod2 = predict(mod2,kyphosis[testi,])
Yobs = c(Yobs,kyphosis$Kyphosis[testi])
p.add = c(p.add,pmod)
p.mul = c(p.mul,pmod2)
}
r.add=roc(Yobs,p.add)
plot(r.add,main=paste("AUROC =",round(r.add$auc*100)/100))
r.mul=roc(Yobs,p.mul)
plot(r.mul,add=TRUE,col="blue",lty=2)
legend(.4,.4,legend=c("sans interactions","avec interactions"),lty=1:2,c
grid()
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 188 / 203
 Validation, sélection de modèles Méthodes boostrap, "out of bag"

Comparaison de modèles

L’ajout des interactions permet d’améliorer légèrement le modèle selon la courbe ROC.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 189 / 203

 Validation, sélection de modèles Méthodes boostrap, "out of bag"

Remarques

Les estimations du risque empirique considérées (validation croisée, bootstrap) sont

asymptotiquement équivalentes et il n’est pas possible de savoir laquelle sera la
meilleure à n fini.

La validation par bootstrap est de plus en plus utilisée et remplace petit à petit la
validation croisée car elle permet de mieux tenir compte de la variabilité des données.

L’estimation d’une erreur de prévision est une opération délicate aux conséquences
importantes. Il est donc nécessaire d’utiliser le même estimateur pour comparer
l’efficacité de deux méthodes ou modèle.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 190 / 203

 Validation, sélection de modèles Sélection de variables

Outline

1 Introduction

2 Regression logistique

3 Inférence pour le modèle logistique

4 Diagnostiques de régression pour les données binaires

5 Variantes des modèles logistiques

6 Régression de Poisson

7 Validation, sélection de modèles

Performance en prédiction
Estimation par validation croisée
Méthodes boostrap, "out of bag"
Sélection de variables

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 191 / 203

 Validation, sélection de modèles Sélection de variables

Problème de la sélection de variables

On a vu au travers des différentes partie du cours qu’une des questions importantes

en modélisation est de choisir les variables explicatives.
Cette question devient même cruciale quand on travaille avec des données
comportant un grand nombre de variables (spectrométrie, données omics, ...).
Solutions proposées plus haut
- tests sur des modèles imbriqués (ex : déviance, AIC, ... )
- comparaison de modèles basée sur l’erreur de classification ou la courbe ROC.
Méthodes pas à pas
Méthode descendante (ou backward) : on ajuste le modèle introduisant toutes les
variables disponibles, puis à chaque étape on retire la variable ayant le moins fort
pouvoir explicatif (choix fait d’après un test de Fisher). On s’arrête lorsque toutes les
variables restantes sont significatives (ou que le fait de retirer la variable dégrade trop
le modèle).
Méthode ascendante (ou forward) : on ajuste le modèle introduisant uniquement une
constante, puis à chaque étape on ajoute la variable ayant le plus fort pouvoir
explicatif (choix fait d’après un test de Fisher). On s’arrête lorsque l’ajout d’une
nouvelle variable n’améliore pas le modèle.
Méthode combinée : on peut combiner les méthodes descendantes et ascendantes.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 191 / 203

 Validation, sélection de modèles Sélection de variables

Méthode descendante

> df<-read.table(’http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.csv’, se
> df$www = as.factor(df$www)
> df$male = as.factor(df$male)
> df$trust = as.factor(df$trust)
> blm<-glm(trust~(educate+income+age+male+www)^2, data=df, family=binomial)
> summary(blm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.115e+00 2.283e+00 -2.240 0.0251 *
educate 1.329e-01 1.525e-01 0.871 0.3836
income 9.801e-03 8.048e-02 0.122 0.9031
age 5.025e-02 3.205e-02 1.568 0.1170
male1 7.423e-01 9.783e-01 0.759 0.4480
www1 3.456e-01 1.261e+00 0.274 0.7841
educate:income 1.088e-03 5.010e-03 0.217 0.8280
educate:age -9.966e-04 2.019e-03 -0.494 0.6215
educate:male1 4.877e-02 5.300e-02 0.920 0.3575
educate:www1 9.999e-03 7.495e-02 0.133 0.8939
income:age -4.522e-05 8.101e-04 -0.056 0.9555
income:male1 3.502e-04 2.456e-02 0.014 0.9886
income:www1 1.097e-02 2.722e-02 0.403 0.6870
age:male1 -1.732e-02 1.007e-02 -1.721 0.0853 .
age:www1 1.787e-03 1.241e-02 0.144 0.8855
male1:www1 -5.715e-01 3.375e-01 -1.693 0.0904 .

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 192 / 203

 Validation, sélection de modèles Sélection de variables

Méthode descendante
> blm.st = step(blm,direction="backward")
> anova(blm.st)
Analysis of Deviance Table
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev
NULL 1173 1596.6
educate 1 64.810 1172 1531.8
income 1 17.994 1171 1513.8
age 1 29.989 1170 1483.8
male 1 4.381 1169 1479.5
www 1 11.506 1168 1467.9
age:male 1 2.484 1167 1465.5
male:www 1 2.170 1166 1463.3
>summary(blm.st)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.515585 0.552140 -9.989 < 2e-16 ***
educate 0.150519 0.026237 5.737 9.65e-09 ***
income 0.030862 0.011609 2.658 0.007852 **
age 0.036382 0.006986 5.207 1.91e-07 ***
male1 1.348318 0.525191 2.567 0.010250 *
www1 0.777063 0.231423 3.358 0.000786 ***
age:male1 -0.016686 0.009651 -1.729 0.083802 .
male1:www1 -0.479713 0.326321 -1.470 0.141544
Null deviance: 1596.6 on 1173 degrees of freedom
Residual deviance: 1463.3 on 1166 degrees of freedom
AIC: 1479.3

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 193 / 203

 Validation, sélection de modèles Sélection de variables

Régression Lasso

Les méthodes pas à pas sont optimales si les variables sont orthogonales (aucune
dépendance entre elles).
Une alternative est la régression pénalisée comme la régression Lasso. L’idée est de
pénaliser la vraisemblance par un terme qui devient grand si on met trop de variables
(inutiles) dans le modèle.
En régression Lasso, on minimise alors le critère pénalisé suivant
p
X
− log L(β) + λ |βj |
j=1

La constante λ joue le rôle d’une pondération et donne plus ou moins d’importance à

chacun des deux termes.
En pratique, la pénalisation va forcer les paramètre βj à rester nul si la variable Xj
n’apporte pas d’information utile.
Un des avantages de cette approche est qu’on peut l’utiliser même si le nombre de
variables est plus grand que le nombre d’observation (ex cas des données OMIC).

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 194 / 203

 Validation, sélection de modèles Sélection de variables

Lasso
df<-read.table(’http://www.indiana.edu/~statmath/stat/all/cdvm/gss_cdvm.csv’, sep=
df1 = df[,-2]
df1$EducInc = df$educate*df$income
df1$EducAge = df$educate*df$age
df1$EducMale = df$educate*df$male
df1$EducWww = df$educate*df$www
df1$IncAge = df$income*df$age
df1$IncMale = df$income*df$male
df1$IncWww = df$income*df$www
df1$AgeMale = df$age*df$male
df1$AgeWww = df$age*df$www
df1$MaleWww = df$male*df$www
library(glmnet)
blm.lasso = glmnet(as.matrix(df1[,-1]),df1[,1],family="binomial")
plot(blm.lasso,xvar="lambda") # Valeurs des coefficients en fonction de lambda

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 195 / 203

 Validation, sélection de modèles Sélection de variables

Lasso

Un des points délicats de la régression Lasso est de bien choisir la valeur de λ.

On peut le faire par validation croisée.
blm.cv = cv.glmnet(as.matrix(df1[,-1]),df1[,1],family="binomial")
> blm.lasso = glmnet(as.matrix(df1[,-1]),df1[,1],family="binomial",lambda=blm.
> blm.lasso$beta
15 x 1 sparse Matrix of class "dgCMatrix"
s0
educate 0.025559246
income .
age 0.006306436
male .
www .
EducInc 0.001687337
EducAge 0.001201131
EducMale 0.014927644
EducWww 0.015936374
IncAge .
IncMale .
IncWww 0.006594953
AgeMale .
AgeWww 0.002718073
MaleWww .

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 196 / 203

 Validation, sélection de modèles Sélection de variables

Validation croisée et sélection de variables

Si on utilise la validation croisée (Kfold ou bootstrap) pour comparer des

méthodes/modèles l’étape de la sélection de variables doit être incluse dans la
boucle de validation (ie réalisée pour chaque échantillon d’apprentissage).
n = nrow(df)
ntest = round(n/5)
p.st <- p.lasso <- Yobs <- p.st <- NULL
B = 100
for (b in 1:B){
testi = sample(1:n,ntest)
appri = setdiff(1:n,testi)
blm = glm(trust~(educate+income+age+male+www)^2, data=df, family=binomial,subs
blm.st = step(blm,direction="backward")
pmod = predict(blm.st,df[testi,])
blm.cv = cv.glmnet(as.matrix(df1[appri,-1]),df1[appri,1],family="binomial")
blm.lasso = glmnet(as.matrix(df1[appri,-1]),df1[appri,1],family="binomial",lam
pmod2 = predict(blm.lasso,as.matrix(df1[testi,-1]),type="response")
Yobs = c(Yobs,df$trust[testi])
p.st = c(p.st,pmod)
p.lasso = c(p.lasso,pmod2)
}
r.st=roc(Yobs,p.st)
plot(r.st)
r.lasso=roc(Yobs,p.lasso)
plot(r.lasso,add=TRUE,col="blue",lty=2)

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 197 / 203

 Validation, sélection de modèles Sélection de variables

Validation croisée et sélection de variables

Dans cette exemple, les deux méthodes conduisent au même résultat, avec une aire
sous la courbe ROC d’environ 0.62.

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 198 / 203

 Validation, sélection de modèles Sélection de variables

Conclusion

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 199 / 203

 Validation, sélection de modèles Sélection de variables

Calcul Matriciel

Transposée

x11 x12 ··· x1p x11 x21 ··· xn1

   
 x21 x22 ··· x2p   x12 x22 ··· xn2 
X= .. ..  , XT =  .. ..
   

 . .   . . 
xn1 xn2 ··· xnp x1p x2p ··· xnp

Produit matrice vecteur

x11 x21 · · · xn1 β1 β1 x11 + β2 x12 + · · · + βp x1p
    
 x12 x22 · · · xn2  β2   β1 x21 + β2 x22 + · · · + βp x2p 
XT β =  . .. .. = ..
    
 ..
 
.  .   . 
x1p x2p · · · xnp βp β1 xn1 + β2 xn2 + · · · + βp xnp

Back

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 200 / 203

 Validation, sélection de modèles Sélection de variables

Vraisemblance pour le modèle logistique

On suppose qu’on observe n couples (yi , xi ) avec yi ∈ {0, 1} et xi ∈ Rp

Modèle logistique avec lien binomial
T
y e xi β
P(Yi = yi ) = πi i (1 − πi )1−yi et πi = π(xi ) = Tβ
1 + exi

Log vraisemblance
n
X
log L(β) = log Pβ (Yi = yi )
i=1
n  
y
X
= log πi i (1 − πi )1−yi
i=1
n
X
= yi log πi + (1 − yi ) log(1 − πi )
i=1

Back

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 201 / 203

 Validation, sélection de modèles Sélection de variables

Intervalle de confiance (1/2)

Soit θ̂n un estimateur du paramètre θ tel que θ̂n vérifie un théorème central limit ie
que quand n tend vers l’infini,
θ̂n − θ
q →Z
Var (θ̂n )

avec Z une variable aléatoire de loi de Gauss centrée et réduite.

En général, les estimateurs du maximum de vraisemblance vérifient le théorème

central limit.

On définit un intervalle de confiance au risque α pour θ̂n à partir des bornes

−z1−α/2 et −z1−α/2 telles que
 
θ̂n − θ
P −z1−α/2 < q < z1−α/2  = 1 − α
 
Var (θ̂n )

Monbet, 12/2016 (- M2) GLM, M2 Pharma. 202 / 203

 Validation, sélection de modèles Sélection de variables

Intervalle de confiance (2/2)

Si n est assez grand, on peut supposer que √θ̂n −θ suit approximativement une loi
Var (θ̂n )
de Gauss et on a donc
 
θ̂n − θ
P −z1−α/2 < q < z1−α/2  = Φ(z1−α/2 ) − Φ(−z1−α/2 )
 
Var (θ̂n )
= 2Φ(z1−α/2 ) − 1

avec Φ la fonction de répartition de la loi de Gauss centrée réduite. On a donc

2Φ(z1−α/2 ) − 1 = 1 − α

soit

z1−α/2 = Φ−1 (1 − α/2)

Les bornes de l’intervalle de confiance pour θ̂n s’écrivent alors

q
µ− = θ̂n − Φ−1 (1 − α/2) Var (θ̂n )
q
µ+ = θ̂n + Φ−1 (1 − α/2) Var (θ̂n )

Back
Monbet, 12/2016 (- M2) GLM, M2 Pharma. 203 / 203