Université Rennes 2, UFR Sciences

Sociales

Régression logistique avec R

Laurent Rouvière

Université Rennes 2
Place du Recteur H. le Moal
CS 24307 - 35043 Rennes
Tel : 02 99 14 18 06
Mel : [Link]@[Link]
Table des matières

1 Introduction 5
1.1 Rappels sur le modèle linéaire..................................................................................................5
1.2 Le modèle linéaire généralisé : GLM........................................................................................6
1.2.1 La régression logistique.................................................................................................6
1.2.2 La régression log-linéaire...............................................................................................9
1.2.3 Généralisation : GLM..................................................................................................10
1.3 Exemples de fonctions de liens pour la régression d’une variable binaire...........................11

2 Analyse discriminante logistique 13

2.1 Le modèle statistique................................................................................................................13
2.1.1 L’échantillonnage..........................................................................................................13
2.1.2 Identifiabilité du modèle.............................................................................................14
2.2 L’estimateur du maximum de vraisemblance.........................................................................15
2.2.1 La vraisemblance..........................................................................................................15
2.2.2 Comportement asymptotique de l’estimateur du maximum de vraisemblance 16
2.3 L’algorithme IRLS....................................................................................................................18
2.3.1 Rappel sur l’algorithme de Newton-Raphson............................................................18
2.3.2 Calcul des estimateurs..................................................................................................18
2.4 Dimensions explicatives, variables explicatives......................................................................19
2.4.1 Variable explicative quantitative.................................................................................19
2.4.2 Variable explicative qualitative...................................................................................20
2.4.3 Interactions....................................................................................................................21
2.5 Interprétation des coefficients β..............................................................................................22
2.6 Précision des estimateurs et tests...........................................................................................24
2.6.1 Loi asymptotique...........................................................................................................24
2.6.2 Intervalles de confiance...............................................................................................24
2.6.3 Tests de nullité de q coefficients libres.......................................................................26
2.7 Le schéma d’échantillonnage rétrospectif................................................................................28
2.8 Un exemple avec R..................................................................................................................30
2.8.1 Modèles “simples”........................................................................................................30
2.8.2 Encore d’autres modèles.............................................................................................32

3 Sélection et validation de modèles 35

3.1 Sélection ou choix de modèle..................................................................................................35
3.1.1 Un outil spécifique : la déviance.................................................................................35
3.1.2 Test de déviance entre 2 modèles emboîtés...............................................................38
3.1.3 Critère de choix de modèles........................................................................................38
3.1.4 Apprentissage/validation...............................................................................................39

Modèle logistique Laurent Rouvière

 3.1.5 Validation croisée.........................................................................................................41
3.1.6 Sélection automatique...................................................................................................43
3.2 Validation du modèle...............................................................................................................46
3.2.1 Test d’adéquation de la déviance...............................................................................46
3.2.2 Test d’Hosmer Lemeshow...........................................................................................47
3.2.3 Analyse des résidus......................................................................................................47
3.2.4 Points leviers et points influents..............................................................................52

4 Modèle logistique multi-classes 55

4.1 Le modèle saturé ou modèle multinomial...............................................................................55
4.2 Modèle polytomique ordonné..................................................................................................56
4.2.1 Cas binaire....................................................................................................................56
4.2.2 Généralisation...............................................................................................................57
4.2.3 L’égalité des pentes......................................................................................................59
4.2.4 Le test d’égalité des pentes.........................................................................................60
4.3 Le modèle polytomique nominal.............................................................................................61
4.3.1 Le modèle.....................................................................................................................61
4.3.2 Estimation et interprétation des paramètres...........................................................62

Annexes 65
A.1 Rappels sur la méthode du maximum de vraisemblance......................................................65
A.2 Echantillonnage Rétrospectif....................................................................................................67
A.3 Exercices....................................................................................................................................69
A.4 Correction..................................................................................................................................73

Bibliographie 77

Annales 78

Laurent Rouvière Modèle logistique

Chapitre 1

Introduction

Notations :
— X = (1, X1, . . . , Xp)′ : vecteur aléatoire de dimension p + 1, les marginales Xj sont
les variables explicatives. On note également x = (1, x1, . . . , xp)′ une réalisation de X ;
— Y variable (univariée) à expliquer.
— (X1, Y1), . . . , (Xn, Yn) : un n-échantillon aléatoire (i.i.d. et de même loi que le couple (X, Y
)), tel que Xi = (1, Xi1, . . . , Xip)′ ;
— (x1, y1), . . . , (xn, yn) une réalisation de (X1, Y1), . . . , (Xn, Yn).
— X : la matrice des observations :
 
1 x11 . . . x1p
. . . . 
X = . . . .  .
1 xn1 . . . xnp

1.1 Rappels sur le modèle linéaire

Le contexte
Nous cherchons à expliquer une variable Y par p variables X = (1, X1, . . . , Xp)′. Pour ce faire,
on dispose de n réalisations (x1, y1), . . . , (xn, yn) du couple (X, Y ). Le but est de modéliser la
dépendance de la variable réponse Y sur les variables explicatives X1, . . . , Xp. Plusieurs
raisons peuvent motiver cette modélisation :
— la description : on veut un modèle qui permette de décrire la relation entre Y et X ;
— l’évaluation des contributions relatives de chaque prédicteur pour expliquer Y ;
— la prédiction : prévoir la valeur de Y pour des nouvelles valeurs des variables explicatives.
On rappelle que le modèle linéaire s’écrit :

Y = X ′ β + ǫ = β0 + β1X1 + . . . + βpXp + ǫ,

avec β = (β0, β1, . . . , βp)′ ∈ Rp+1 et ǫ ∼ (0, σ2). On distingue alors deux cas :
— Les variables Xi sont N déterministes (non aléatoires) :

Y ∼ N (X′β, σ2), E[Y ] = X ′ β ;

— Les variables Xi sont aléatoires :

(Y |X) ∼ N (X′β, σ2), E[Y |X] = X′β.

Modèle logistique Laurent Rouvière

6 Introduction

Plaçons nous maintenant dans le cas où la variable à expliquer Y est qualitative (sexe, couleur,
présence ou absence d’une maladie...) et possède un nombre fini de modalités g1, . . . , gm. Le pro-
blème consiste alors à expliquer l’appartenance d’un individu à un groupe à partir des p variables
explicatives, on parlera de discrimination au lieu de régression.

Il est bien entendu impossible de modéliser directement la variable Y par une relation linéaire
(imaginons que Y soit le sexe d’une personne ou la couleur de ses cheveux). Afin de pallier
| = gk X = x). Supposons pour simplifier
cette difficulté, on va s’intéresser aux probabilités P(Y
que la variable Y prenne uniquement deux valeurs : 0 (“groupe 0”) ou 1 (“groupe 1”). La
connaissance de P(Y = 1|X = x) implique celle de P(Y = 0|X = x) : il suffit par conséquent de
modéliser la probabilité p(x) = P(Y = 1|X = x). On peut par exemple envisager une relation de
la forme

pβ(x) = β0 + β1x1 + . . . + βpxp + ε = x′β + ε.

Cette approche possède plusieurs inconvénients :

— Remarquons tout d’abord que la variance de |Y X = x vaut pβ(x)(1 pβ(x)). Contraire-
ment au modèle linéaire traditionnel, cette variance n’est pas constante et par
conséquent l’hypothèse d’homoscédasticité des résidus ne sera pas vérifiée.
— Le fait qu’aucune restriction ne soit effectuée sur les β implique que x′β peut prendre n’im-
porte quelle valeur dans R. Ce fait est gênant pour l’estimation d’une probabilité (imaginer
une estimation du genre Pβˆ(Y = 1|X = x) = —1297.56 ! ! !).
Pour ces raisons, nous devons étendre le modèle linéaire classique aux cas où :
— Y peut être une variable qualitative (présence ou absence d’une maladie, appartenance à
une catégorie...) ;
— les erreurs peuvent ne pas avoir la même variance (s’affranchir de l’hypothèse d’homoscé-
dasticité).

1.2 Le modèle linéaire généralisé : GLM

1.2.1 La régression logistique

Nous nous plaçons tout d’abord dans un contexte de classification binaire, c’est-à-dire que nous
supposons qu’il existe seulement deux groupes à discriminer. Nous verrons dans le chapitre 4
comment étendre les techniques à des modèles multiclasses (plus de deux groupes).

Exemple 1.1 Nous souhaitons expliquer la variable Y présence (1)/ absence (0) d’une
maladie cardio-vasculaire (Chd) par l’âge des patients. Les données sont représentées sur la
figure 1.1.

Laurent Rouvière Modèle logistique

1.2 Le modèle linéaire généralisé : GLM 7

* * * * * * * * * * * * * * * * * * * * * * * * * *

1.
0
0.
8
0.
6
ch
d 0.
4
0.
2

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
0.
0

20 30 40 50 60 70
age

FIGURE 1.1 – Représentation directe de Chd (variable à expliquer Y ) en fonction de

l’âge (variable explicative X).

Cette figure montre qu’il est difficile de modéliser les données brutes, la variabilité de la variable
CHD est élevée pour tout âge. Une méthode permettant de réduire cette variabilité consiste à
regrouper les patients par classe d’âge. Nous obtenons le tableau suivant :

CHD
Age Effectif Absent Présent % Présent
]19 ;29] 10 9 1 0.1
]29 ;34] 15 13 2 0.133333
]34 ;39] 12 9 3 0.25
]39 ;44] 15 10 5 0.333333
]44 ;49] 13 7 6 0.461538
]49 ;54] 8 3 5 0.625
]54 ;59] 17 4 13 0.764706
]59 ;69] 10 2 8 0.8
TaBLE 1.1 – Données regroupées en classe d’âge.

La liaison entre l’âge et la présence de la maladie devient beaucoup plus claire. Il apparaît en effet
que lorsque l’âge augmente, la proportion d’individus atteint par la maladie augmente. La figure
1.2 permet d’évaluer cette liaison : elle apparaît nettement sous la forme d’une courbe sigmoïde
(en forme de “S”). Il semblerait donc “naturel” de modéliser cette proportion de malades par classe
d’âge en fonction de l’âge par une courbe sigmoïde.

Modèle logistique Laurent Rouvière

8 Introduction

1.
0
o

0.
8
o

0.
6
ch
d
o
0.
4

o
0.
2

o
o
0.
0

20 40 60 80
age

FIGURE 1.2 – Fréquence de Chd par classe d’âge en fonction de l’âge.

La colonne “moyenne” du tableau 1.1 fournit une estimation de E[Y |X = x] pour chaque classe
d’âge. Nous pouvons donc proposer une modélisation de l’espérance conditionnelle de E[Y |X =
x] :

E[Y |X = x] = hβ(x)

où l’allure de la courbe représentative de hβ est une sigmoïde.

Plusieurs fonctions hβ peuvent naturellement être utilisées. Pour le modèle logistique on considère
la fonction h(x) = exp(x)/(1 + exp(x)), ce qui donne le modèle
exp(β0 + β1x)
E[Y |X = x] = p(x) = ,
1 + exp(β0 + β1 x)
où encore
p(x)
logit p(x) = log = β + β x,
0 1
1 — p(x)
logit désignant la fonction bijective et dérivable de ]0, 1[ dans R : p'→ log(p/(1 p)) (voir figures
1.3 et 1.4). Nous verrons qu’une telle modélisation permettra de retrouver un grand nombre des
“bonnes” propriétés du modèle linéaire.
3

1.0
2

0.8
1

0.6
0

0.4
−1

0.2
−2
−3

0.0

0.0 0.2 0.4 0.6 0.8 −3 −2 −1 0 1 2 3

1.0

FIGURE 1.3 – Fonction logit . FIGURE 1.4 – Fonction inverse de logit .

Laurent Rouvière Modèle logistique

1.2 Le modèle linéaire généralisé : GLM 9

La loi conditionnelle de la variable d’intérêt diffère entre le modèle logistique et le modèle linéaire.
Dans le modèle de régression linéaire Y = β0 + β1X + ε, on fait l’hypothèse que les résidus ε
suivent une loiN (0, σ2). Il vient Y|X = x ∼(β0 + β1x, σ2). Pour le modèle logistique, pour une
observation x deN la variable explicative, on peut exprimer la variable d’intérêt comme suit :

Y = p(x) + ε.

La variable aléatoire ε peut prendre simplement deux valeurs : si y = 1 alors ε = 1 — p(x) et si

y = 0 alors ε = —p(x). Par conséquent ε prend pour valeur 1 — p(x) avec probabilité p(x) et —
p(x) avec probabilité 1 — p(x) : Y |X = x suit une loi de Bernoulli de paramètre p(x).
Définition 1.1 (Régression logistique) Soit Y une variable à valeurs { dans
} 0, 1 à
expliquer par p variables explicatives X = (1, X1, . . . , Xp) . Le modèle logistique propose une
′

modélisation de la loi de Y |X = x par une loi de Bernoulli de paramètre pβ(x) = Pβ(Y = 1|X =
x) telle que :

log p β( x ) =
β +β +...+β = x′β, (1.1)
x x
0 1 1 p p
1 — pβ (x)
ou encore
logit pβ(x) = x′β,
logit désignant la fonction bijective et dérivable de ]0, 1[ dans R : p '→ log(p/(1 — p)).

L’égalité (1.1) peut également s’écrire

exp(x′β)
pβ(x) = Pβ(Y = 1|X = x) = .
1 + exp(x′β)
Remarque 1.1 Dans un modèle logistique, nous effectuons deux choix pour définir le modèle :
1. le choix d’une loi pour Y |X = x, ici la loi de Bernoulli ;
2. le choix de la modélisation de P(Y = 1|X = x) par

logit Pβ(Y = 1|X = x) = x′β.

La fonction logit est bijective et dérivable. Elle est appelée fonction de lien.
Remarquons également que
,
Eβ[Y |X = x] = Pβ(Y = 1|X = x)


,
Vβ(Y |X = x) = Pβ(Y = 1|X = x) 1 — Pβ(Y = 1|X = x)

ce qui implique que la variance n’est pas constante et varie selon x.

1.2.2 La régression log-linéaire

Dans le modèle logistique la variable à expliquer est une variable binaire. Le modèle log-linéaire
traite le cas d’une variable de comptage. Voici quelques exemples :
— nombre de catastrophes aériennes sur une période donnée ;
— nombre de voitures à un feu rouge ;
— nombre d’accidents par jour sur une autoroute...
Modèle logistique Laurent Rouvière
1 Introduction
0
Définition 1.2 (Régression log-linéaire) Soit Y une variable de comptage à expliquer
par le vecteur X = (1, X1, . . . , Xp)′. Le modèle log-linéaire propose une modélisation de la loi
| de Y
X = x par une loi de poisson de paramètre λ = λ(x) telle que :

log E[Y |X = x] = x′β.

Pour une nouvelle mesure x effectuée, le modèle log-linéaire va donc prédire la valeur exp(x′β).

Remarque 1.2 Ici encore, deux choix sont effectués pour définir le modèle :
1. le choix d’une loi pour Y |X = x, ici la loi de Poisson ;
2. le choix de la modélisation de E(Y |X = x) par

log E[Y |X = x] = x′β.

La fonction log est bijective et dérivable.

1.2.3 Généralisation : GLM

On peut résumer les remarques précédentes par le tableau :

Choix logistique log-linéaire linéaire

loi de Y |X = x Bernoulli Poisson Normale

modélisation
de logit E[Y |X = x] = x′β log E[Y |X = x] = x′β E[Y |X = x] = x′β
E[Y |X = x]

Une généralisation de ces méthodes est appelée GLM (Generalized Linear Model). L’approche
GLM consiste à :
1. choisir une loi pour Y |X = x parmi un ensemble restreint de lois (les lois exponentielles
GLM) ;
2. choisir une fonction de lien g(.) parmi une ensemble réduit de fonctions bijectives et déri-
vables.
| X = x] par la fonction g est ensuite
3. la transformation de l’espérance conditionnelle E[Y
modélisée par une fonction η qui n’est autre qu’une combinaison linéaire des variables
explicatives :
g (E[Y |X = x]) = η(x) = x′β.

On peut résumer un modèle GLM par le schéma suivant :

Laurent Rouvière Modèle logistique

1.3 Exemples de fonctions de liens pour la régression d’une variable binaire 11

A expliquer
composante aléatoire Explicatif
Lien Composante systématique
Y |X = x suit une loi fixée.
E[Y X| = x] dépend de η(x) On modélise η par une com- binaison
au travers de la fonction g linéaire des Xj
appelée fonction de lien p
Σ
g(E[Y |X]) = η(X) η(x) = xjβj
j=0

g est une fonction inversible.

Remarque 1.3 1. Pour choisir un modèle GLM il faut

— choisir la loi de Y | X = x dans la famille exponentielle des GLM.
— choisir une fonction de lien inversible g.
2. Pour utiliser un modèle GLM il faudra donc estimer les paramètres β = (β0, β1, . . . , βp)
′
. Une fois cette estimation réalisée, on disposera d’une estimation de η(x) ainsi que de
E[Y |X = x] = g−1(η(x)).

Le tableau 1.2 donne quelques exemples de GLM.

Loi Nom du lien Fonction de lien

Bernoulli/Binomiale lien logit g(µ) = logit (µ) = log(µ/(1 — µ))
Poisson lien log g(µ) = log(µ)
Normale lien identité g(µ) = µ
Gamma lien réciproque g(µ) = —1/µ

TaBLE 1.2 – Exemples de GLM.

1.3 Exemples de fonctions de liens pour la

régression d’une variable binaire
D’autres fonctions de lien que logit peuvent être utilisées dans le cas où la variable à expliquer Y
est binaire. On trouve notamment dans la littérature les transformations :
— probit, qui n’est autre que l’inverse de la fonction de répartition de la loi normale
centrée réduite : ∫ ǫ
1 1
∀p ∈ [0, 1], probit(p) = ǫavec √ exp 2
du = p.
2π −∞ — u
2
— log-log définie par :
∀p ∈ [0, 1], log-log(p) = log(— log(1 — p)).

Ces transformations sont représentées sur la figure 1.5.

Modèle logistique Laurent Rouvière

1 Introduction
2

4
2
0
−
2
−
4

0.0 0.2 0.4 0.6 0.8 1.0

FIGURE 1.5 – Fonctions de liens : probit (trait plein), logit (tirets), log-log (pointillés).
Des trois fonctions de lien présentées, la transformation log-log est bien appropriée aux cas où
l’on souhaite modéliser les probabilités de succès de manière asymétrique. Les transformations
logit et probit possèdent des propriétés identiques. Dans de nombreux cas, on préfère utiliser la
transformation logistique. Plusieurs raisons motivent ce choix :
— d’un point de vue numérique, la transformation logistique est plus simple à manipuler
(notamment pour l’écriture des estimateurs du maximum de vraisemblance, voir section
2.2) ;
— on a une interprétation claire des coefficients en terme d’odds ratio pour la transformation
logistique (voir section 2.5).
— le modèle logistique est particulièrement bien adapté à un schéma d’échantillonnage rétros-
pectif (voir section 2.7).
Nous nous focaliserons dans la suite sur le modèle logistique. Les différents résultats obtenus
pour- ront s’étendre aux autres modèles GLM. Il est important de connaître les notations des GLM
présentées dans cette partie. C’est en effet sous cette forme qu’elles sont présentées dans la litté-
rature ainsi que dans la plupart des logiciels statistiques (par exemple sur R).

Laurent Rouvière Modèle logistique

Chapitre 2

Analyse discriminante logistique

Nous rappelons que Y désigne une variable à expliquer binaire (qui prend 2 valeurs 0 ou 1 pour
simplifier) ou un label qui dénote l’appartenance à un groupe et X1, . . . , Xp désignent p variables
explicatives. On souhaite :
— expliquer la variable Y à l’aide des p variables explicatives X = (1, X1, . . . , Xp)′ ;
— étant donnée une nouvelle mesure x des p variables explicatives, prédire le label y associé à
cette nouvelle observation.
Nous avons vu dans le chapitre précédent que le modèle logistique est défini par

logit pβ(x) = β0 + β1x1 + . . . + βpxp = x′β (2.1)

où β = (β0, . . . , βp)′ et x = (1, x1, . . . , xp)′. Nous nous posons maintenant le problème de
l’estima- tion des paramètres β à partir d’un échantillon (x1, y1), . . . , (xn, yn).

2.1 Le modèle statistique

Les paramètres du modèle logistique sont généralement estimés par la méthode du maximum de
vraisemblance (voir annexe A.1). Afin d’écrire “proprement” la vraisemblance, il convient de
définir avec précision le modèle statistique dans lequel nous allons nous placer. On rappelle qu’un
modèle statistique est un couple (Hn, P) où H est l’espace de chaque observation et P est une
famille de lois de probabilité sur Hn muni de sa tribu borélienne.

2.1.1 L’échantillonnage
On note (x1, y1), . . . , (xn, yn) une réalisation d’un n-échantillon (X1, Y1), . . . , (Xn, Yn) i.i.d. et de
même loi que (X, Y ). Sans perte de généralité, nous supposons que l’échantillonnage a été réalisé
en deux temps :
— x1, . . . , xn est la réalisation d’un n-échantillon i.i.d. X1, . . . , Xn. On note PX la loi de X1 ;
— pour i = 1, . . . , n, yi est généré selon une loi de Bernoulli de paramètre pβ(xi) avec
logit pβ(xi) = x′β (β∈ Rp+1 étant le paramètre à estimer).
Nous distinguerons deux structures de données. On parlera de
— données individuelles lorsque tous les xi sont différents. On appellera alors design l’ensemble
{x1, . . . , xn}. Le modèle est alors défini par
}
Mn = {0, 1}n, {B(pβ(x1)) ⊗ . . . ⊗ B(pβ(xn)), β ∈ Rp+1} .

Modèle logistique Laurent Rouvière

1 Analyse discriminante logistique
4
— données répétées lorsque certaines valeurs de xi sont répétées plusieurs fois (ce qui peut se
produire si par exemple PX est discrète à support borné). On note
— x1, . . . , xT les différentes valeurs observées ;
ΣT
— nt, t = 1, . . . , T le nombre de fois où xt a été observé (on a donc t=1 nt = n) ;
— yt, t = 1, . . . , T le nombre de succés (Y = 1) obervé au point xt.
On appellera design l’ensemble {(x1, n1), . . . , (xT , nT ) }. Pour cette structure de données, le
modèle sera défini par
Mn = }{0, . . . n1} × . . . × {0, . . . nT }, {Bin(n1, pβ(x1)) ⊗ . . . ⊗ Bin(nT , pβ(xT )), β ∈ Rp+1}
.
Exemple 2.1 Dans le cadre d’un essai clinique, on teste 2 traitements (A et B) sur n = 220
patients atteints d’une pathologie. On reporte dans le tableau 2.1 le nombre de patients qui ont
guéri au bout de 3 mois de traitements.

Guérisons Non Guérisons

A 40 60
B 90 30
TaBLE 2.1 – Exemple de données répétées.
Nous sommes ici en présence de données répétées. Un individu (x, y) est un patient (x représente
le traitement subi et y vaut 1 si le patient a guéri au bout de 3 mois, 0 sinon). Le design est
{(A, 100), (B, 120)} et on a y1 = 40 et y2 = 90.
Les propriétés des estimateurs sont très proches pour ces deux types de données. Néanmoins,
certains concepts tels que la forme de la vraisemblance où les tests d’adéquation par la déviance
peuvent légèrement différer. Dans ce chapitre, nous nous focalisons sur le cas de données indi-
viduelles (qui est le cas le plus fréquent). Pour une étude plus approfondie du cas des données
répétées, nous renvoyons le lecteur à l’annexe A.3 (pour l’écriture de la vraisemblance) ou aux
ouvrages de Hosmer & Lemeshow (2000) et Collet (2003).

2.1.2 Identifiabilité du modèle

On rappelle qu’un modèle {H, {Pθ, θ ∈ Θ}} est identifiable si θ '→ Pθ est injective. Dit brutale-
ment, un modèle est dit identifiable si deux paramètres différents définissent deux lois différentes.

Exemple 2.2 Plaçons nous dans le cadre de la régression d’une variable binaire par une variable
explicative à deux modalités A et B. On considère le modèle logistique
M }
3
1 = {0, 1}, {B(pβ(x)), β ∈ R }

où
logit pβ(x) = β0 + β11x=A + β21x=B
= β0 + β1(1 — 1x=B) + β21x=B
= (β0 + β1) + 01x=A + (β2 — β1)1x=B
Si on pose β˜ = (β0—β1, 0, β2 —β1) on a β =/ β˜ et pβ(x) = pβ˜(x). M 1 n’est donc pas identifiable.
Une solution classique consiste à mettre une contrainte sur les coefficients. Si par exemple, on pose
β2 = 0 alors le modèle ′ }
M1 = {0, 1}, {B(pβ(x)), β ∈ R2}
où logit pβ(x) = β0 + β11x=A est identifiable.

Laurent Rouvière Modèle logistique

2.2 L’estimateur du maximum de vraisemblance 1
5
Propriété 2.1 Le modèle Mn est identifiable si et seulement si rang(X) = p + 1.

Preuve. Soit β /= β˜. Par définition Mn est identifiable si et seulement si il existe i ≤ n tel que
pβ(xi) /= p ˜(xi). Supposons que : ∀1 ≤ i ≤ n, x′ β = x ′ β˜ . On a alors
β i i

     
 1 x
 .   x1p
α0  .. + .
11
1  .  + . . . + αp  .  = 0
α 1 xn1 xnp

avec α = β — β˜ /= 0. Une colonne de X est donc combinaison linéaire des autres, d’où rang(X)
< p +. 1

2.2 L’estimateur du maximum de vraisemblance

2.2.1 La vraisemblance
On se place dans le cas de données individuelles. La vraisemblance du modèle (identifiable) Mn
est définie par :

Ln : {0, 1}n × Rp+1 → R+

(y1, . . . , yn, β) '→ B(pβ(x1)) ⊗ . . . ⊗ B(pβ(xn))({y1, . . . , yn}).

Remarque 2.1 Si on désigne par Yi une variable aléatoire de loi B(pβ(xi)) alors les variables
Y1, . . . , Yn sont indépendantes mais pas de même loi.

Lorsqu’il n’y aura pas de confusion possible, on commetra l’abus de notation Ln(y1, . . . , yn, β) =
Ln(β). Calculons la vraisemblance. On a
n n
Y Y
Ln(β) = Pβ(Y = yi|X = xi) = pβ(xi)yi (1 — pβ(xi))1−yi .
i=1 i=1

En passant au log nous obtenons

Ln(β) = Σ
n
{yi log(pβ(xi)) + (1 — yi) log(1 — pβ(xi))}
i=1
Σ pβ(xi)
n
= y log + log(1 — p (x )) ,
i β i
i=1
1 — pβ (xi)

ou encore, d’après (2.1),

Ln(β) = Σn i i
{yix′β — log(1 + exp(x′β))}. (2.2)
i=1

Modèle logistique Laurent Rouvière

1 Analyse discriminante logistique
6 h i′
∂L
Le vecteur gradient au point β défini par ∇L (β) = (β), . . . , (β) s’obtient par dérivation
∂L

n n
n
∂β0 ∂βp
n
i
∂L Σ xij
y — exp(x β)
′i
n
(β) = xi ij
∂βj 1 + exp(x′β)
i=1
n
Σ
= [xij(yi — pβ(xi))] .
i=1

Ce qui donne en écriture matricielle

n
Σ
∇Ln(β) = [xi(yi — pβ(xi))] = X′(Y — Pβ)
i=1

où Y = (y1, . . . , yn)′ et Pβ = (pβ(x1), . . . , pβ(xn))′. L’estimateur du maximum de vraisemblance

(si il existe) est solution de l’équation (appelée équation du score) :

S(β) = ∇Ln(β) = X′(Y — Pβ) = 0. (2.3)

On rappelle que si cette équation admet une solution en β notée g(y1, . . . , yn) (et que cette solution
est un maximum de Ln(β)) alors l’estimateur du maximum de vraisemblance est βˆ = g(Y1, . . . , Yn).

Trouver explicitement βˆ n’est pas possible. En effet, l’équation (2.3) se réécrit :

,
exp(β1x11 + . . . + βpx1p) exp(β1xn1 + . . . + βpxnp)
 x11 y1 + . . . + n1
x ny = 11
x 1 + exp(β x + . . . + β x + . . . + xn1 1 + exp(β +...+β x )
) x
1 11 p 1p 1 n1 p np

.
 .
, xp y + . . . + x y = xp exp(β1x11 + . . . + βpx1p) + . . . + x exp(β1xn1 + . . . + βpxnp)
.
1 1 np n 1 + exp(β . . . + β p 1 + exp(β 1 . . . + β
1 +
1 11 x p 1 ) np 1 n x+ p np )
Ce système (qui n’est pas xlinéaire en β) n’admet pas x
de solution analytique. On a donc recours à
des méthodes numériques qui nécessitent de connaître d’éventuelles propriétés sur la régularité de
la fonction à optimiser (en terme de convexité par exemple).

2.2.2 Comportement asymptotique de l’estimateur du maximum

de vrai- semblance
Définition 2.1 Le nuage de points est dit :
— complètement séparable si Eβ ∈ Rp+1 : ∀i tel que Yi = 1 on a xi′β > 0 et ∀i tel que Yi = 0
on a xi′β < 0 ;
— quasi-complètement séparable si Eβ ∈ Rp+1 : ∀i tel que Yi = 1 on a xi′β ≥ 0, ∀i tel que
Yi = 0 on a xi′ β ≤ 0 et {i : xi′ β = 0} =/ ∅ ;
— en recouvrement s’il n’est ni complètement séparable ni quasi-complètement séparable
(voir figure 2.1).

Laurent Rouvière Modèle logistique

 2.2 L’estimateur du maximum de vraisemblance 1
7

FIGURE 2.1 – Exemple de séparabilité complète (gauche), quasi-complète (milieu) et de

recouvrement (droite).
Le théorème suivant nous assure la convergence d’un algorithme itératif vers βˆ et nous donne le
comportement asymptotique de l’estimateur du maximum de vraisemblance. Considérons le jeu
d’hypothèses suivant :
— H1 : rang(X) = p + 1.
— H2 : on est en situation de recouvrement.
— H3 : la matrice E[XX′] existe et est définie positive.

Théorème 2.1 1. Sous les hypothèses H1 et H2 la log-vraisemblance β '→ L n(β) est stricte-
ment concave : βˆ existe et est unique.
2. Sous H3 on a
(a) la convergence en probabilité de βˆ vers β quand n → ∞.
(b) la loi asymptotique de l’estimateur du maximum de vraisemblance
√
n(βˆ — β) → N (0, I(β)−1),

où I(β) est la matrice d’information de Fisher (de dimension (p + 1) (p + 1)) au point

β :
h ∂2 i
I(β)kℓ = — L (Y, β) , 0 ≤ k, ℓ ≤
∂βk ∂βℓ 1
E p,
où avec un léger abus de notation I(β)kℓ désigne le terme de la (k + 1)ème ligne et
(ℓ + 1)ème colonne de I(β).

Pour la preuve de la concavité, on pourra se référer au polycopié de Guyon (2005) ou à l’article de

Albert & Anderson (1984). La loi asymptotique découle de la théorie du maximum de vraisem-
blance (voir annexe A.1 et Antoniadis et al (1992) pour une preuve détaillée). La concavité a une
conséquence numérique importante puisqu’elle justifie qu’un algorithme itératif convergera bien
vers la valeur de βˆ. Il n’y a donc pas de risque de converger vers un maximum local non global et
la convergence de l’algorithme ne dépend pas du point d’initialisation de l’algorithme. Pour la
conver-
gence en probabilité et la normalité asymptotique, on pourra regarder Fahrmeir & Kaufmann (1985).

Modèle logistique Laurent Rouvière

 1 Analyse discriminante logistique
8
2.3 L’algorithme IRLS
Deux algorithmes sont généralement implémentés sur les logiciels de statistique pour calculer les
estimateurs du maximum de vraisemblance : l’agorithme du score de Ficher et l’algorithme IRLS
(Iterative Reweighted Least Square). Nous présentons ici uniquement le second algorithme.

2.3.1 Rappel sur l’algorithme de Newton-Raphson

La méthode de Newton-Raphson permet une résolution numérique des équations du score.
Pour simplifier les notations, nous supposons que β est univarié. On part tout d’abord d’une
valeur initiale arbitraire de β, notée β0 et on désigne par β1 = β0 + h une valeur candidate pour
être solution de S(β) = 0, c’est-à-dire S(β0 + h) = 0. Par un développement limité à l’ordre un
de la fonction S, on obtient l’approximation suivante :

S(β0 + h) ≃ S(β0) + hS′(β0).

Comme S(β0 + h) = 0, on obtient pour h la valeur suivante :

h = —[S′(β0)]−1 S(β0)

et donc
β1 = β0 — [S′(β0)]−1 S(β0).
Dans le cas qui nous concerne β ∈ Rp+1 et S(β) = ∇Ln(β). La formule de récurrence se traduit
par
2
β1 = β0 — [∇ Ln(β0)]−1∇Ln(β0)
où ∇2Ln(β0) désigne la matrice hessienne de la log-vraisemblance au point β0 :

∇2 )
0 kℓ
= ∂2L (β 0 ) , 0 ≤ k, ℓ ≤
Ln(β ∂βk ∂βℓ p

où nous comettons toujours l’abus de désigner par ∇ 2 n(β0)kℓ le terme de la (k + 1)ème ligne et
L
(ℓ + 1)ème colonne de ∇2 n(β0). Le processus est ensuite itéré jusqu’à convergence. Il se résume
L
de la manière suivante :
1. choix d’un point de départ β0 ;
2. On construit βk+1 à partir de βk

βk+1 = βk + Ak∇Ln(βk),

où ∇Ln(βk) est le gradient au point βk et Ak = —(∇2Ln(βk))−1 est la matrice de “pas” de

l’algorithme (l’inverse du hessien de Ln au point βk).

2.3.2 Calcul des estimateurs

Calculons la matrice hessienne 2
∂ 2L n(β)
∇ Ln(β) = :
n
∂β
′
k ∂βℓ 0≤k,ℓ≤p
n
2 Σ Σ
∂ Ln exp(xiβ)
(β) = — i i =— i i β (x i)(1 — β (x i)),
∂βk ∂βℓ (1 + exp(xi′
p
xkxℓ β))2 xkxℓp
i=1 i=1
Laurent Rouvière Modèle logistique
2.4 Dimensions explicatives, variables explicatives 1
9
Algorithme 1 maximisation de la vraisemblance
Require: β0
k→1
repeat
βk+1 → βk + Ak∇Ln(βk)
k→k+1
until βk+1 ≈ βk et/ou Ln(βk+1) ≈ Ln(βk)

en écriture matricielle nous obtenons

Σn
2
∇ Ln(β) = — xixi′ pβ(xi)(1 — pβ(xi)) = X′WβX,
i=1

où Wβ est la matrice diagonale diag (pβ(xi)(1— pβ(xi)), i = 1, . . . , n. Nous pouvons maintenant

exprimer βk+1 en fonction de βk :

βk+1 = βk + (X′Wβk X)−1X′(Y — Pβk )

= (X′Wβk X)−1X′Wβk (Xβk + Wβk −1(Y — Pβk ))
= (X′Wβk X)−1X′Wβk Zk,
où Zk = Xβk +Wβk −1(Y— Pβk ). Cette équation est simplement une régression pondérée du
vecteur Zk où les poids Wβk dépendent de X et βk ; d’où le nom de IRLS pour cette méthode.
Les poids sont réévalués à chaque étape de l’algorithme, une étape étant une simple
régression pondérée.

2.4 Dimensions explicatives, variables explicatives

Les remarques formulées dans cette partie s’appliquent pour la plupart des modèles de régres-
sion (modèles linéaires et d’analyse de variance par exemple). Pour plus de détails, on pourra se
rapporter aux ouvrages de Droesbeke et al (2007) et Cornillon & Matzner-Løber (2007).
On rappelle que la dimension d’un modèle paramétrique (identifiable){H ,{ Pθ, θ∈ Θ}} est la
dimension de l’espace des paramètres Θ. Pour le modèle logistique, les lois sont déterminées par
logit pβ(x) = β0 + β1x1 + . . . + βpxp.
Tout comme pour le modèle linéaire, la dimension d’un modèle logistique s’obtient en sommant
les dimensions explicatives associées aux différentes variables explicatives du modèle, lesquelles
varient suivant la nature de la variable explicative. Nous étudions dans cette partie les dimensions
explicatives pour des variables explicatives quantitatives et qualitatives. Nous étudierons
également les cas d’interactions entre variables.

2.4.1 Variable explicative quantitative

Si on dispose d’une seule variable explicative X quantitative (non regroupée en classe) le
modèle s’écrit
logit pβ(x) = β0 + β1x.
Un seul coefficient (β1) est alloué à X, cette variable est représentée par une seule colonne dans la
matrice du design X. Sa dimension est donc égale à 1.

Modèle logistique Laurent Rouvière

 20 Analyse discriminante logistique

2.4.2 Variable explicative qualitative

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 :
logit pα(x) = α0 + αF 1F (x) + αH1H(x), (2.4)
mais aussi
logit pα(x) = (α0 + αF ) + (αH — αF )1H(x).
Il y a une infinité d’écritures possibles... Le modèle (2.4) correspond à une matrice du design X
à trois colonnes où la première colonne est une colonne de 1 et les deux dernières sont obtenues
en effectuant un codage disjonctif complet pour chaque individu (le ième terme de la 2ème (resp.
3ème) colonne vaut 1 si le ième individu de l’échantillon est une femme (resp. un homme)). Par
conséquent, la somme des deuxième et troisième colonnes vaut 1 ce qui rend l’estimation
impossible puisque la matrice X n’est pas de plein rang (X′WβX n’est pas inversible et le
modèle n’est pas identifiable). Une solution pour pallier cette difficulté consiste à mettre une
contrainte sur les coefficients αH et αF . 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 prise comme modalité de référence
par rapport à laquelle on mesure des déviations. Le choix de cette modalité n’a bien entendu pas
d’influence sur les lois Pβ. Il en a cependant une sur la valeur des coefficients estimés ainsi que
sur leurs écarts types. Ainsi le nombre de coefficients significativement différents de 0 peut
changer suivant le choix de la modalité de référence. Ceci montre clairement que, pour juger
l’apport d’une variable qualitative, il n’est pas pertinent d’utiliser les tests de significativité des
coefficients. Il sera préférable de réaliser un test entre modèles emboîtés (voir page 38).

Exemple 2.3 Considérons le cas d’une variable explicative à trois niveaux g1, g2, g3. Les
observa- tions sont récoltées dans les tableaux suivants (équivalents)

observation X Y
1 g1 1 X #{Y = 1} #{Y = 0}
2 g2 1
g1 2 1
3 g3 1
g2 1 1
4 g1 1
g3 1 0
5 g2 0
6 g1 0

On effectue une régression logistique sur R :

> X <- factor(c("g1","g2","g3","g1","g2","g1"))
> Y <- factor(c(1,1,1,1,0,0))
> model <- glm(Y~X,family=binomial)
> model

Call: glm(formula = Y ~ X, family = binomial)

Coefficients:
(Intercept) Xg2 Xg3
0.6931 -0.6931 17.8729

Laurent Rouvière Modèle logistique

2.4 Dimensions explicatives, variables explicatives 2
1
Degrees of Freedom: 5 Total (i.e. Null); 3 Residual
Null Deviance: 7.638
Residual Deviance: 6.592 AIC: 12.59
La modalité g1 est ici prise comme modalité de référence. Le modèle estimé s’écrit donc :
0.6931 si j =
, 0 si j = 2
logit pˆ(g j ) = logit pβˆ(gj) 1
,
= 0.6931 + 17.8729 = 18.566 si j = 3.
ou encore
, exp(0.6931) = 2/3 si j = 1
 1+exp(0.6391)
pˆ(gj ) = 1/2 si j = (2.5)
, 2
exp(18.566)
1+exp(18.566)
= 1.0000 si j = 3.
Bien évidemment, changer la contrainte modifie l’écriture du modèle mais ne modifie la loi de
probabilité sous-jacente. Prenons par exemple, la modalité g3 comme modalité de référence :
> model1 <- glm(Y~C(X,base=3),family=binomial)
> model1

Call: glm(formula = Y ~ C(X, base = 3), family = binomial)

Coefficients:
(Intercept) C(X, base = 3)1 C(X, base = 3)2
18.57 -17.87 -18.57

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

Null Deviance: 7.638
Residual Deviance: 6.592 AIC: 12.59
Il est par exemple facile de vérifier que les probabilités pˆ(g j ), j = 1, 2, 3 sont identiques à celles de
(2.5).

2.4.3 Interactions
Tout comme en analyse de la variance, on ne peut se contenter de modèles purement additifs.
Reprenons l’exemple développé dans Droesbeke et al (2007) (page 122). Nous considérons le cas
où la variable Y représente le fait de faire (codé 1) ou non (codé 0) de la couture. On dispose de
deux variables explicatives : l’age et le sexe. Le modèle “purement” additif s’écrit :
logit pβ(x) = β0 + β1age + β21femme,
la modalité homme a été choisie comme modalité de référence. Une telle écriture revient à
supposer que les pentes sont identiques pour les hommes et les femmes (voir Figure 2.2).
logit logit
pβ(x) Femmes Femmes
pβ(x)

Hommes

Hommes

age age
FIGURE 2.2 – Modèle additif. FIGURE 2.3 – Modèle avec interaction.
Modèle logistique Laurent Rouvière
22 Analyse discriminante logistique

A priori, il semble que les hommes ne font que très rarement de la couture quel que soit leur
âge, il parait préférable de pouvoir utiliser un modèle du genre (voir Figure 2.3) :
logit p(x) = β0 + β1age + β21femme + β3age1femme.
Ce modèle revient à considérer l’interaction entre les variables age et sexe. On rappelle que
deux variables interagissent si l’effet de l’une sur Y diffère suivant les valeurs de l’autre. Bien
entendu, l’ajout d’une interaction augmente la dimension explicative du modèle. Le nombre de
composantes supplémentaires s’obtient en faisant le produit du nombre de dimensions des
variables qui interagissent (ici les variables sexe et age sont de dimension 1, on rajoute donc une
dimension).

2.5 Interprétation des coefficients β

0.7

0.8
0.2
0.3

β=0 β = 0.5
1.0

1.0
0.0

0.0

β=2 β = 10
FIGURE 2.4 – Pβ(Y = 1|X = x) pour différentes valeurs de β.
Nous avons représenté sur la Figure 2.4 l’allure de la courbe représentative de la fonction x '→
exp(xβ)
1+exp(xβ) pour différentes valeurs du paramètre β. On remarque que :
— pour de faibles valeurs de β on a une large plage de valeurs de x pour lesquelles la
fonction se situe aux alentours de 0.5 (la fonction est même constante (0.5) dans le cas
extrême β = 0). Pour ces valeurs pβ(x) = Pβ|(Y = 1 X = x) sera proche de 0.5 et on peut
donc penser qu’il sera difficile de discriminer ;
— lorsque β augmente, la zone où la fonction est proche de 0.5 diminue et la fonction est
proche de 0 ou 1 pour un grand nombre de valeurs de x. Par conséquent, Pβ(Y = 1| X = x)
sera souvent proche de 1 ou 0, ce qui risque de minimiser d’éventuelles erreurs de
prévisions.
On peut interpréter ainsi : plus β est grand, mieux on discrimine. Cependant une telle interpré-
tation dépend des valeurs que x prend (de son échelle). C’est pourquoi en général l’interprétation
des coefficients β s’effectue en termes d’odds ratio. Les odds ratio sont des outils souvent
appréciés dans le domaine de l’épidémiologie (mais pas toujours bien utilisés !).
Les odds ratio servent à mesurer l’effet d’une variable quantitative ou le contraste entre les effets
d’une variable qualitative. L’idée générale est de raisonner en termes de probabilités ou de rapport

Laurent Rouvière Modèle logistique

2.5 Interprétation des coefficients β 23

de cotes (odds). Si on a, par exemple, une probabilité p = 1/4 de gagner à un jeu, cela signifie
que sur 4 personnes une gagne et les trois autres perdent, soit un rapport de 1 gagnant sur trois
perdants, c’est-à-dire p/(1 — p) = 1/3. Ce rapport p/(1 p) varie entre 0 (0 gagnant) et l’infini
(que des gagnants) en passant par 1 (un gagnant pour un perdant).
Définition 2.2 L’odds (chance) pour un individu x d’obtenir la réponse Y = 1 est défini par :
p(x)
odds(x) = , où p(x) = P(Y = 1| X = x).
1 — p(x)
L’odds ratio (rapport des chances) entre deux individus x et x˜ est
odds(x) p(x)
1−p(x)
OR(x, x˜ ) = = .
odds(x˜) p(x˜)
1−p(x˜
)

Exemple 2.4 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).
Les odds ratio peuvent être utilisés de plusieurs manières :
1. Comparaison de probabilités de succès entre deux individus (voir Tableau 2.2) ;

OR(x, x˜ ) > 1 ⇐⇒ p(x) > p(x˜)

OR(x, x˜ ) = 1 ⇐⇒ p(x) = p(x˜)
OR(x, x˜ ) < 1 ⇐⇒ p(x) < p ( x˜ )
TaBLE 2.2 – Règles d’interprétation des odds ratio.
2. Interprétation en termes de risque relatif : dans le cas où p(x) et p(x˜ ) sont
très petits par rapport à 1, comme dans le cas d’une maladie très rare, on peut faire
~
l’approximation OR(x, x˜ ) p(x)/p(x˜) et interpréter simplement. Par exemple si OR(x,
x˜ ) = 4 alors la réponse (maladie) est 4 fois plus probable dans le cas où X = x que dans
le cas où X = x˜ .
3. Mesure de l’impact d’une variable : pour le modèle logistique
logit pβ(x) = β0 + β1x1 + . . . + βpxp,
il est facile de vérifier que
OR(x, x˜ ) = exp(β1(x1 — x˜ 1 )) . . . exp(βp(xp — x˜ p )).
Pour mesurer l’influence d’une variable sur l’odds ratio, il suffit de considérer deux obser-
vations x et x˜ qui diffèrent uniquement par la jème variable. On obtient alors
OR(x, x˜ ) = exp(βj(xj — x˜ j )).
Ainsi une variation de la jème variable d’une unité (sur l’échelle de cette variable)
corres- pond à un odds ratio exp(βj) qui est uniquement fonction du coefficient βj. Le
—
coefficient βj permet de mesurer l’influence de la jème variable sur le rapport pβ(x)/(1
pβ(x)) lorsque xj varie d’une unité, et ceux indépendamment de la valeur de xj. Une telle
analyse peut se révéler intéressante pour étudier l’influence d’un changement d’état d’une
variable qualita- tive.

Modèle logistique Laurent Rouvière

24 Analyse discriminante logistique

Exemple 2.5 Reprenons l’exemple des cotes pour les courses de chevaux. On cherche a
expliquer la performance d’un cheval en fonction du jockey qui le monte. Pour simplifier on
suppose que l’on a que deux jockeys A et B. On désigne par Y la variable aléatoire qui prend pour
valeurs 0 si le cheval remporte la course, 1 sinon. On considère le modèle logistique

logit pβ(x) = β0 + β11x=B.

On a OR(B, A) = exp(β2). Imaginons que pour un échantillon de taille n on obtienne une estima-
tion βˆ2 = log(2). On a alors OˆR(B, A) = 2, ce qui signifie que la cote du cheval est multipliée par
2 lorsqu’il est monté par B par rapport à A.

2.6 Précision des estimateurs et tests

2.6.1 Loi asymptotique
Nous avons obtenu dans le théorème 2.1 le comportement asymptotique de l’estimateur du maxi-
mum de vraisemblance βˆ :
√
n(βˆ — β) →L N (0, I(β)
−1 ),

où I(β) est la matrice d’information de Fisher au point β. On déduit

(βˆ — β)′nI(β)(βˆ — β) →L χ2 p+1 .

Un tel résultat n’est pas utilisable tel quel puisque la matrice I (β) est inconnue. On remarque que
d’après la loi des grands nombres
n
1 Σ Σ
n 2 2 2
∂ 1 ∂ 1 ∂
Iˆ (β kℓ
) = — L 1 (Y i , β) = — L 1 (Y i; β) = — L n(Y 1, . . . , Yn , β)
n i=1 ∂βk ∂βℓ n ∂βk ∂βℓ i=1 n ∂βk ∂βℓ
1
= (X′W X) ,
β kℓ
n
converge presque surement vers I (β)kℓ. Comme βˆ converge faiblement vers β, on obtient grace au
théorème de Slutsky et aux opérations classiques sur la convergence en loi

(βˆ — β)′ Σˆ (βˆ —

2
χp+1 (2.6)
β) →L
avec Σˆ = (X
′
WβˆX).

2.6.2 Intervalles de confiance

On déduit du paragraphe précédent qu’un estimateur de la variance de βˆ j est donné par le jème
terme de la diagonale de Σˆ − 1 . Notons jσˆ 2 cet estimateur. Il vient
(β ˆ ˆ — βj L
— βj) L
β
2 2
j
2 → χ1 ou encore σˆ j → N (0, 1). (2.7)
σˆ j
j

Un intervalle de confiance (asymptotique) de niveau 1 — α pour βj est donc donné par

h i
ˆ ˆ
IC1−α(βj) = β j — u1−α/2 σˆj ; β j + u1−α/2 σˆj ,
Laurent Rouvière Modèle logistique
 2.6 Précision des estimateurs et tests 25

où u1−α/2 représente le quantile de niveau (1 — α/2) de la loi normale N (0, 1).

La validité de ces intervalles est relative puisqu’il s’agit d’une approximation valable asymptoti-
quement. Il est toujours possible de compléter cette étude par un bootstrap afin d’obtenir d’autres
intervalles de confiance dans le cas ou ceux-ci sont particulièrement importants. Cela dit, en pra-
tique, on se contente de l’intervalle de confiance bâti grâce à la matrice d’information de Fisher.

On déduit également de (2.7) les tests de nullité des coefficients du modèle. On note H0 : βj = 0
et H 1: β j/= 0, alors sous H 0, βˆj/ σˆj → N (0, 1). On rejettera H0 si la valeur observée de βj ˆ /jσˆ
L

dépasse en valeur absolue le quantile d’ordre 1 — α/2 de la loi N (0, 1).

Exemple 2.6 Reprenons l’exemple du TP1. Un chef d’entreprise souhaite vérifier la qualité
de ces machines en fonction de l’âge et de la marque des moteurs. Il dispose
— d’une variable binaire Y (1 si le moteur a déjà connu une panne, 0 sinon) ;
— d’une variable quantitative age repésentant l’âge du moteur ;
— d’une variable qualitative a 3 modalités marque représentant la marque du moteur.
On souhaite expliquer la variable Y à partir des deux autres variables. On construit un modèle
logistique permettant d’expliquer Y par l’age du moteur et sa marque.
> model <- glm(panne~.,data=donnees,family=binomial)

On obtient les variances estimées des estimateurs via

> P <- predict(model,type="response")
> W <- diag(P*(1-P))
> X1 <- rep(1,n)
> X2 <- donnees$age
> X3 <- [Link](donnees$marque==1)
> X4 <- [Link](donnees$marque==3)
> X <- matrix(c(X1,X2,X3,X4),ncol=4)
> #ecarts types estimes
> sqrt(diag(solve(Sigma)))
[1] 0.83301450 0.09398045 0.81427979 1.05357830

Bien évidemment, on peut retrouver ces écarts types ainsi que les probabilités critique des tests de
nullité des coefficients avec
> summary(model)

Call:
glm(formula = panne ~ ., family = binomial, data = donnees)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.4232 -1.2263 0.9082 1.1062 1.5982

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.47808 0.83301 0.574 0.566
age 0.01388 0.09398 0.148 0.883
marque1 -0.41941 0.81428 -0.515 0.607
marque3 -1.45608 1.05358 -1.382 0.167

(Dispersion parameter for binomial family taken to be 1)

Modèle logistique Laurent Rouvière

26 Analyse discriminante logistique

Null deviance: 45.717 on 32 degrees of freedom

Residual deviance: 43.502 on 29 degrees of freedom
AIC: 51.502

2.6.3 Tests de nullité de q coefficients libres

La théorie du maximum de vraisemblance nous donnant la loi (asymptotique) des estimateurs,
il est possible de tester la significativité des variables explicatives. Pour cela, trois tests sont
généralement utilisés :
— Le test de Wald ;
— Le test du rapport de vraisemblance ou de la déviance.
— Le test du score ;
Les hypothèses s’écrivent :

H0 : βj1 = βj2 = . . . = βjq = 0 contre H1 : Ek ∈ {1, . . . , q} : βjk /= 0.

Pour alléger les notations, nous supposons sans perte de généralité que nous testons la nullité des
q premiers coefficients du modèle

H0 : β0 = β1 = . . . = βq−1 = 0 contre H1 : Ek ∈ {0, . . . , q — 1} : βk /= 0.

Test de Wald Il est basé sur (2.6). On note β0,...,q−1 le vecteur composé des q premières com-
posantes de β et Σˆ0,...,q
−1
la matrice bloc composée des q premières lignes et colonnes de Σˆ − 1 . Il est
facile de voir que sous H0
βˆ0,...,q−1Σ
ˆ
′ ˆ
0,...,q β 0,...,q−1
L
χ2q.
→

Test du rapport de vraisemblance ou de la déviance La statistique de test est

basée sur la différence des rapports de vraisemblance entre le modèle complet et le modèle sous
H0. On note
βˆH0 l’estimateur du maximum de vraisemblance contraint par H0 (il s’obtient en supprimant les q
premières variables du modèle). On a alors sous H0
L 2
2(Ln (βˆ) — Ln (βˆH0 )) → χq .

Test du score On cherche ici à vérifier si la fonction de score (gradient de la log-

vraisemblance) est proche de 0 sous H0. Sous H0 on a

S (βˆ )′ Σˆ −1 S (βˆ ) L→ χ2,

H0 H0 H0 q

où
= X′W βˆ X.
Σˆ H 0 H0

Pour ces 3 tests, on rejette l’hypothèse nulle si la valeur observée de la statistique de test dépasse
le quantile d’ordre 1 α de la loi χ2. Pour la preuve des convergences des statistiques du
—
maximum de vraisemblance et du
q score, on pourra se référer à l’annexe D de Antoniadis et al
(1992). La figure 2.5 permet de visualiser les trois tests. Le test du score revient à tester la nullité
de la pente
Laurent Rouvière Modèle logistique
2.6 Précision des estimateurs et tests 27
ˆ ˆ ˆ ˆ
en βH 0 (β sous H0 ), le test de Wald la nullité de la distance entre β et βH et le test du rapport
0

de vraisemblance la nullité de la différence entre les vraisemblances en ces deux points.

Modèle logistique Laurent Rouvière

 28 Analyse discriminante logistique

Ln(β)

Test du rapport des vraisemblances

Test de Wald
Test du score
ℓmax

ℓ0
Log-vraisemblance

βˆH βˆ0 β
FIGURE 2.5 – Rapport de vraisemblance, score, test de Wald.
Remarque 2.2 — La PROC LOGISTIC sous SAS réalise les trois tests pour H0 : β1 = β2 =
. . . = βp = 0.
— Pour les tests “variable par variable” ou paramètre par paramètre

H0 : β j = 0 contre H1 : βj /= 0,

la PROC LOGISTIC utilise le test de Wald.

Exemple 2.7 Reprenons l’exemple précédent. Le modèle s’écrit

logit pβ(x) = β0 + β1age + β21marque=1 + β31marque=3,

la modalité 0 de la variable marque est prise comme modalité de référence. On éffectue les 3 tests
présentés ci-dessus pour

H0 : β 1 = β 2 = β 3 = 0 H1 : Ej ∈ {1, 2, 3} βj /= 0.
Le calcul des statistiques de test et des probabilités critiques s’obtient avec les commandes
> #test de Wald
> Sigma <- t(X)%*%W%*%X
> inv_Sigma <- solve(Sigma)
> inv_SigmaH0 <- inv_Sigma[2:4,2:4]
> betaH0 <- coef(model)[2:4]
> statWald <- t(betaH0)%*%solve(inv_SigmaH0)%*%betaH0
> pcWald <- 1-pchisq(statWald,df=3)
> pcWald
[,1]
[1,] 0.5655233
> #test du rapport de vraisemblance
> modelH0 <- glm(panne~1,data=donnees,family=binomial)
> statRappvrais <- 2*(logLik(model)-logLik(modelH0))
> pcRappvrais <- 1-pchisq(statRappvrais,df=3)

Laurent Rouvière Modèle logistique

2.6 Précision des estimateurs et tests 29

> pcRappvrais[1]
[1] 0.528955
> #test du score
> prevH0 <- predict(modelH0,type="response")
> scoreH0 <- t(X)%*%([Link](donnees$panne)-1-prevH0)
> WH0 <- diag(prevH0*(1-prevH0))
> SigmaH0 <- t(X)%*%WH0%*%X
> statscore <- t(scoreH0)%*%solve(SigmaH0)%*%scoreH0
> pcscore <- 1-pchisq(statscore,df=3)
> pcscore
[,1]
[1,] 0.5392691

Ces 3 tests acceptent l’hypothèse nulle.

2.7 Le schéma d’échantillonnage rétrospectif

Jusqu’à présent nous avons considéré un échantillon (X1, Y1), . . . , (Xn, Yn) i.i.d. de même loi que
(X, Y ). Cette phase d’échantillonnage n’est pas nécessairement toujours la mieux adaptée. Consi-
dérons l’exemple suivant.

Exemple 2.8 Une clinique cherche à mesurer l’effet du tabac sur le cancer du poumon. Elle
prélève parmi ses patients un échantillon composé de n1 = 250 personnes atteintes par le cancer et
n0 = 250 personnes ne présentant pas la maladie. Les résultats (simulés !)de l’étude sont présentés
dans le tableau suivant :

Fumeur Non fumeur

Non malade 48 202
Malade 208 42
TaBLE 2.3 – Résultats de l’enquête.
Le statisticien responsable de l’étude réalise un modèle logistique. Les sorties sur R sont :

Call: glm(formula = Y ~ X, family = binomial)

Coefficients:
(Intercept) Xnon_fumeur
1.466 -2.773

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

Null Deviance: 692.3
Residual Deviance: 499.9 AIC: 503.9

On obtient pβˆ(fumeur) = 0.812, ce qui peut paraître un peu élevé. La valeur surprenante d’une
telle estimation vient du fait que l’échantillonnage n’est pas fait selon protocole précédent : il est
fait conditionnellement à Y . Il est facile de voir que les répartitions d’individus selon la variable
Y ne sont pas les mêmes dans la population et dans l’échantillon. Ceci va entrainer un biais au
niveau des estimateurs. On peut modéliser le schéma d’échantillonnage rétrospectif de la manière
suivante.

Modèle logistique Laurent Rouvière

2.730
Le schéma d’échantillonnage rétrospectif 29
Analyse discriminante logistique

Le schéma d’échantillonnage On s’intéresse toujours à la loi conditionnelle |de Y X qui

est une bernoulli de paramètre pβ(x) telle que

logit pβ(x) = β0 + β1x1 +.....+ βpxp.

On désigne toujours par PX la loi de X. La différence ici est que l’échantillon n’est pas i.i.d. de
même loi que (X, Y ). On désigne par S une variable aléatoire qui prend pour valeur 0 et 1 et
par τ1 et τ0 les taux de sondage P(S| = 1 Y = 1) et P(S = 1 Y = 0). Pour un individu (x, y)
|
généré selon PX,Y , on tire une réalisation s de S selon une loi de Bernoulli de paramètre τy, si s
= 1 on garde (x, y) dans l’échantillon, sinon on le jette. On répète le protocole jusqu’à obtenir n
individus.

Le modèle logistique Ici un individu est représenté par un triplet (X, Y, S) et l’échantillon
constitué s’écrit (x1, y1, 1), . . . , (xn, yn, 1). Le modèle étudié pour ce schéma d’échantillonnage est
donc Mn = {{0, 1}n, B(pγ(x1) ⊗ × ⊗B(pγ(xn)), γ ∈ Rp+1} avec

logit pγ(x) = logit Pγ(Y = 1|X = x, S = 1) = γ0 + γ1x1 +. + γpxp.

Ainsi, il suffit d’appliquer tout ce qui a été vu dans ce chapitre pour obtenir les estimateurs
des γj (on déduit également les intervalles de confiance, test...). La question qui se pose est :
comment retrouver pβ(x) en connaissant pγ(x) ?
Propriété 2.2 Avec les notations ci dessus, on a

par conséquent
τ1
logit pγ(x) = logit (x) + log ,
pβ τ0

logit p
β (x) = log τ1
γ — +γ x +...+γ x .
0 1 1 p p
τ0
Preuve. Il suffit de remarquer que grâce au théorème de Bayes on a

pβ(x)P(S = 1|Y = 1, X = x) pβ(x)P(S = 1|Y = 1)

pγ (x) = P(S = 1) = P(S = 1) .

□ Ce résultat est intéressant puisqu’il implique que

le biais dû au mode d’échantillonnage est exclusivement concentré sur la constante du modèle. Si
de plus on connait les taux de sondage, alors on peut corriger ce biais. Il s’agit là d’une propiété
spécifique au modèle logistique.

Ce schéma d’échantillonnage est souvent utilisé lorsque les probabilités πi = P(Y = i) (probabilité
d’appartenance au groupe i) sont très différentes les unes des autres. Dans de tels cas, le schéma
classique conduit à travailler avec des effectifs trop petits dans certains groupes, alors que le
schéma rétrospectif permet de travailler avec des effectifs comparables. Par exemple, en
diagnostic médical, on a souvent des problèmes de discrimination entre deux groupes où l’un des
groupes (1 par exemple) est associé à une maladie et l’autre est caractéristique de l’absence de
Laurent Rouvière Modèle logistique
~
cette maladie. Dans de telles situations, on a bien sûr π1 beaucoup plus petit que π0. L’usage
consiste alors à étudier deux échantillons de taille à peu près équivalente (n1 n0), le premier étant
celui des malades, le second celui des individus sains.

Modèle logistique Laurent Rouvière

30 Analyse discriminante logistique

Exemple 2.9 Reprenons l’exemple 2.8. Des études préalables ont montré que les
probabilités P(Y = 1) = π1 et P(Y = 0) = π0 pouvaient être estimées par 0.005 et 0.995.
L’échantillon- nage a été effectué de manière à travailler avec des échantillons équilibrés,
nous pouvons donc estimer les probabilités P(Y = 1|S = 1) et P(Y = 0|S = 1) par 1/2. Ainsi en
remarquant que

τ1 P(Y = 1|S = 1) P(Y = 0)

τ0 = P(Y = 0|S = 1) P(Y = 1)
on peut estimer τ1
τ0
par 0.995/0.005 ≈ —5.293. On a donc

logit pβˆ(x) = (1.466 — 5.293) — 2.7731non fumeur(x).

D’où pβˆ(fumeur) = 0.0213.

2.8 Un exemple avec R

Le traitement du cancer de la prostate change si le cancer a atteint ou non les nœuds lymphatiques
entourant la prostate. Pour éviter une investigation lourde (ouverture de la cavité abdominale) un
certain nombre de variables sont considérées comme explicatives de la variable Y binaire : Y = 0
le cancer n’a pas atteint le réseau lymphatique, Y = 1 le cancer a atteint le réseau lymphatique.
Le but est d’expliquer Y par les variables suivantes :
— âge du patient au moment du diagnostic : age ;
— le niveau d’acide phosphatase sérique : acide ;
— Le résultat d’une analyse par rayon X, 0= négatif, 1=positif : rayonx ;
— La taille de la tumeur, 0=petite, 1=grande : taille ;
— L’état pathologique de la tumeur déterminé par biopsie (0=moyen, 1=grave) : grade ;
— Le logarithme népérien du niveau d’acidité : [Link].
age acide rayonx taille grade [Link].
1 66 0.48 0 0 0 -0.73396918
2 68 0.56 0 0 0 -0.57981850
3 66 0.50 0 0 0 -0.69314718
4 56 0.52 0 0 0 -0.65392647
5 58 0.50 0 0 0 -0.69314718
6 60 0.49 0 0 0 -0.71334989
7 65 0.46 1 0 0 -0.77652879
8 60 0.62 1 0 0 -0.47803580
9 50 0.56 0 0 1 -0.57981850
10 49 0.55 1 0 0 -0.59783700
TaBLE 2.4 – Représentation des dix premiers individus.

2.8.1 Modèles “simples”

Nous sommes en présence de 6 variables explicatives X1, . . . , X6 avec :
— X1, X2 et X6 quantitatives ;
— X3, X4 et X5 qualitatives (2 niveaux pour chacune).

Laurent Rouvière Modèle logistique

 2.8 Un exemple avec R 31

Premier modèle
Considérons tout d’abord les trois variables explicatives qualitatives X = (X3, X4, X5) :

logit pβ(x) = β0 + β31{x3=1} + β41{x4 =1} + β51{x5 =1}.

Ce modèle possède 4 paramètres. Les sorties du logiciel R sont :

> model_quali<-glm(Y~rayonx+taille+grade,data=donnees,family=binomial)
> model_quali

Call: glm(formula = Y ~ rayonx + taille + grade, family = binomial, data = donnees)

Coefficients:
(Intercept) rayonx1 taille1 grade1
-2.1455 2.0731 1.4097 0.5499

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

Null Deviance: 70.25
Residual Deviance: 52.78 AIC: 60.78

Si par exemple (x3, x4, x5) = (1, 0, 1), on aura alors :

logit pβˆ(x) = βˆ0 + βˆ3 + βˆ5 = —2.1455 + 2.0731 + 0.5499 = 0.4785

et
exp(0.4785)
pβˆ(x) = = 0.6174.
1 + exp(0.4785)
Ainsi, dans un contexte de prévision, nous serons tentés d’assigner le label 1 à la nouvelle obser-
vation x.

Deuxième modèle
Considérons maintenant le modèle uniquement composé de variables quantitatives,

logit pβ(x) = β0 + β1x1 + β2x2 + β6x6.

> model_quanti<-glm(Y~age+acide+[Link].,data=donnees,family=binomial)
> model_quanti

Call: glm(formula = Y ~ age + acide + [Link]., family = binomial, data = donnees)

Coefficients:
(Intercept) age acide [Link].
12.34700 -0.02805 -9.96499 10.54332

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

Null Deviance: 70.25
Residual Deviance: 59.95 AIC: 67.95

Troisième modèle
Le modèle “complet” à 6 variables s’écrit

logit pβ(x) = β0 + β1x1 + β2x2 + β31{x3 =1} + β41{x4 =1} + β51{x5 =1} + β6x6.

Modèle logistique Laurent Rouvière

 32 Analyse discriminante logistique

> model_complet<-glm(Y~.,data=donnees,family=binomial)
> model_complet

Call: glm(formula = Y ~ ., family = binomial, data = donnees)

Coefficients:
(Intercept) age acide rayonx1 taille1 grade1
10.08672 -0.04289 -8.48006 2.06673 1.38415 0.85376
[Link].
9.60912

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

Null Deviance: 70.25
Residual Deviance: 44.77 AIC: 58.77

2.8.2 Encore d’autres modèles...

Comme dans le cas du modèle “linéaire” on peut également considérer des interactions entre
les variables explicatives. On dit qu’il y a interaction entre deux variables F 1 et F 2 sur une
variable Y si l’effet de l’une des variables diffère selon la modalité de l’autre. Remarquons que
cette notion n’a rien à voir avec celle de corrélation qui ne concerne que deux variables alors
que l’interaction met en jeu une troisième variable Y .
Exemple 2.10 (Construction d’interaction) On s’intéresse à l’effet de deux traitements
X1 et X2 sur le rhume. Le traitement X1 consiste à prendre à intervalle de temps réguliers deux
verres de cognac et X2 représente un traitement aux antibiotiques (il n’est pas difficile de
comprendre l’intérêt d’envisager une interaction). La variable réponse Y correspond à l’état du
patient (1 si malade, 0 si bonne santé). N’ayant pas encore trouvé suffisamment de volontaires
pour réaliser l’étude, on simule un échantillon suivant le modèle
1. deux facteurs X1 et X2 à deux niveaux équiprobables ;
2. la loi de Y conditionnellement à X1 et X2 est donnée dans le tableau 2.5.

X1
0 1
X2
0 B(0.95) B(0.05)
1 B(0.05) B(0.95)
TaBLE 2.5 – Lois conditionnelles de Y (B désigne la loi de Bernoulli).
On estime les pourcentages de mal classés sur un échantillon indépendant (voir section 3.1.4) et
on reporte dans le tableau suivant les pourcentages de mal classés pour les modèles sans et avec
interaction. Nous voyons l’intérêt d’inclure une interaction pour cet exemple.

Sans 0.54
Avec 0.065
TaBLE 2.6 – Pourcentages de mal classés.

Pour l’exemple du cancer de la prostate, le modèle avec toutes les interactions d’ordre 2 s’écrit :
> model_inter<-glm(Y~.^2,data=donnees,family=binomial)
Warning message:
des probabilités ont été ajustées numériquement à 0 ou 1 in: [Link](x = X, y = Y,

Laurent Rouvière Modèle logistique

2.8 Un exemple avec R 33

weights = weights, start = start, etastart = etastart,

> model_inter

Call: glm(formula = Y ~ .^2, family = binomial, data = donnees)

Coefficients:
(Intercept) age acide rayonx1
2.843e+17 -4.229e+15 -3.117e+17 -5.453e+16
taille1 grade1 [Link]. age:acide
2.516e+16 -5.778e+15 2.026e+17 4.665e+15
age:rayonx1 age:taille1 age:grade1 age:[Link].
2.077e+13 -5.245e+13 -1.670e+14 -2.869e+15
acide:rayonx1 acide:taille1 acide:grade1 acide:[Link].
5.572e+16 -2.420e+16 2.336e+16 -5.687e+15
rayonx1:taille1 rayonx1:grade1 rayonx1:[Link]. taille1:grade1
1.129e+15 -1.176e+15 -4.004e+16 -5.496e+15
taille1:[Link]. grade1:[Link].
8.625e+15 -1.228e+16
Degrees of Freedom: 52 Total (i.e. Null); 31 Residual
Null Deviance: 70.25
Residual Deviance: 504.6 AIC: 548.6

On peut vérifier que ce modèle nécessite l’estimation de 22 paramètres (1 + 6 + 26 ). Bien entendu,

d’autres sous-modèles avec interactions peuvent être utilisés. De plus, nous pouvons nous
demander si toutes les variables sont bien explicatives ? Dés lors, des méthodes sélection et
validation de modèles doivent être envisagées.

Modèle logistique Laurent Rouvière

Chapitre 3

Sélection et validation de modèles

Ce chapitre se divise en deux parties :

1. Sélection : Etant donnés M modèlesM1, . . . MM , comment choisir le “meilleur” à partir
de l’échantillon dont on dispose.
2. 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).
— L’ajustement peut être aussi regardé observation par observation (individus aber-
rants) par des méthodes graphiques (analyse des résidus) ou analytiques.
• Est-ce que les hypothèses sont vérifiées ? Les méthodes sont essentiellement graphiques
(analyse des résidus).
• L’influence des observations sur l’estimation des paramètres peut être aussi envisagée
(distance de Cook, robustesse).
Dans ce chapitre nous allons traiter ces questions à travers l’exemple du modèle logistique. Nous
noterons Mβ le modèle logistique défini par le vecteur de paramètre β. L’ensemble des méthodes
présentées peut s’étendre à d’autres problématiques de sélection-validation de modèles.

3.1 Sélection ou choix de modèle

Si on se restreint à des modèles logistiques, sélectionner un modèle revient à choisir les variables
(interactions inclues) qui vont constituer le modèle.

3.1.1 Un outil spécifique : la déviance

Il est difficile de se faire une idée sur l’ajustement en se basant sur la valeur de la vraisemblance
puisqu’elle dépend (entre autres) de la taille de l’échantillon. Pour la régression logistique, un outil
spécifique est introduit : la déviance. Elle compare la vraisemblance obtenue à celle d’un modèle
de référence : le modèle complet (ou modèle saturé). Ce modèle ne place pas de contrainte sur la
forme du paramètre p(x) de la loi de Y |X = x (cette loi est une Bernoulli de paramètre p(x)).

Modèle logistique Laurent Rouvière

 36 Sélection et validation de modèles

Modèle saturé en présence de données individuelles

En présence de données individuelles (x1, y1), . . . , (xn, yn), il est facile de voir que, sous le
modèleMsaturé {{sat = } 0,{1 n, B(p(x⊗1)) .⊗. . B(psat(xn)), psat(xi) ∈R }}
, l’estimateur du maximum de
vraisemblance pˆsat (x i ) de psat(xi) est donné par pˆsat (x i ) = yi la log-vraisemblance maximisée
L
sat vaut ainsi 0. Ce modèle contient autant de paramètres que de données et reconstitue “parfaitement”
l’échantillon.

Modèle saturé en présence de données répétées

En présence de données répétées de design {(x1, n1), . . . , (xT , nT )} il est facile de voir que les
estimateurs du maximum de vraisemblance du modèle saturé sont donnés par pˆsat (x t ) = y¯ t =
yt/nt, t = 1, . . . , T . Dans ce cas, la vraisemblance maximisée est donnée par :
T
Y n t nt−yt
Lsat = psat (x( t))yt (1 — sat (xt )) .
t=1
yt ˆ pˆ
On obtient la log-vraisemblance
maximisée
T
Σ nt T
Σ
Lsat = log + yt log pˆsat (x t ) + (nt — yt) log(1 — pˆsat (xt )).
yt t=1
t=1

Que ce soit en présence de données individuelles ou répétées, le modèle saturé possède autant
de paramètres que de points du design. Ce modèle est le modèle le plus complexe (en terme de
nombre de paramètres) puisqu’il propose un coefficient différent pour chaque point du design. Tous
les autres modèles logistiques sont emboîtés dans celui-ci.
Définition 3.1 La déviance d’un modèle Mβ est définie par

DMβ = 2(Lsat — L n (βˆ)).

En présence de données individuelles on a DMβ = — 2L ˆ

n (β ) et en présence de données répétées, la
déviance s’écrit "
=2 T #
nt y¯t y¯ t 1 — y¯
p )ˆ(x + (1 — y¯t ) 1 — p ˆ(x ) )
t
D Mβ Σ
log
log
β β t
t=1
" t
# t t
n —
t y— t,
Σ n t — t t
T
y¯ y) β
= yt log
t
+ log n p ˆ(x
2 t=1 pβˆ (x t (n )
) une différence de log-vraisemblance. Elle constitue un
où y¯t = yt/nt. La déviance est égale à 2 fois
écart en terme de log-vraisemblance entre le modèle saturé d’ajustement maximum et le modèle
considéré.

Ajustement
parfait
bon moyen mauvais Qualité d’ajustement
✲

Laurent Rouvière Modèle logistique

3.1 Sélection ou choix de modèle 37
0 Déviance

Modèle logistique Laurent Rouvière

 38 Sélection et validation de modèles

Exemple 3.1 (calcul de déviance) Considérons l’exemple du cancer de la prostate et

calculons d’abord la déviance pour le modèle Y~age+acide. Nous somme ici en présence de données
indivi- duelles, on obtient la déviance via les commandes :
> mod1 <- glm(Y~age+acide,data=donnees,family=binomial) #construction du modele
> #calcul de la vraisemblance
> prev <- mod1$[Link] #on obtient les pi
> vrais <- rep(0,nrow(donnees))
> vrais[donnees$Y==1] <- prev[donnees$Y==1]
> vrais[donnees$Y==0] <- 1-prev[donnees$Y==0]
> vrais <- prod(vrais) #vrais est la vraisemblance du modele
> dev <- -2*log(vrais) #pour ce modele, il n’y a pas de repetitions aux points du
design, donc Lsat=0
> dev
[1] 65.72393

Bien entendu, le logiciel peut retourner directement la valeur de la déviance

> mod1$deviance
[1] 65.72393

Si maintenant on considère le modèle Y~age+taille, nous somme en présence de données répétées.

Les données se trouvent dans le fichier “donnees_bin_age_taille.txt” dont voici les
premières lignes :
"age" "taille" "Y1" "Y0"
49 "0" 0 1
50 "0" 1 0
51 "0" 0 2
52 "0" 0 1
56 "0" 1 3
58 "0" 0 2
Les deux premières colonnes représentent les valeurs des variables explicatives. On retrouve
ensuite (colonne Y1) le nombre de réponses Y=1 et (colonne Y0) le nombre de réponses Y=0. Le
modèle est construit via la commande :
> donnees1 <- [Link]("donnees_bin_age_taille.txt",header=T)
> model1<-glm(cbind(Y1,Y0)~age+taille,data=donnees1,family=binomial)

La déviance se calcule comme suit

> prev <- model1$fitted #calcul des pi

> ni <- apply(donnees1[,3:4],1,sum)
> ti <- donnees1$Y1
> ybi <- donnees1$Y1/ni
> #calcul des termes combinatoires (facultatif)
> vect_comb <- rep(0,nrow(donnees1))
> for (i in 1:nrow(donnees1)){
+ vect_comb[i]<- log(prod(1:ni[i])/(max(c(prod(1:ti[i]),1))*
max(c(prod(1:(ni[i]-ti[i])),1))))}
> vect <- ni*(ybi*log(prev)+(1-ybi)*log(1-prev))
> vrais_model1 <- sum(vect_comb+vect) #log_vraisemblance du modele
> #modele sature
> vect_sat <- ni*(ybi*log(ybi)+(1-ybi)*log(1-ybi))
> vect_sat[[Link](vect_sat)] <- 0
> vrais_modelsat <- sum(vect_comb+vect_sat)

Laurent Rouvière Modèle logistique

3.1 Sélection ou choix de modèle 39

> #on deduit la deviance

> 2*(vrais_modelsat-vrais_model1)
[1] 37.15260

On retrouve cette valeur directement

> model1$deviance
[1] 37.15260

3.1.2 Test de déviance entre 2 modèles emboîtés

Rappelons que par définition un modèle est emboîté dans un autre plus général (ou plus grand)
lorsqu’il est un cas particulier de ce modèle plus général.

Exemple 3.2 Les modèles logistiques définis par

logit pβ(x) = β0 + β1x1 + β2x2

et

logit pγ(x) = γ0 + γ1x1 + γ2x2 + γ3x3

sont emboîtés l’un dans l’autre.

Il est facile de voir que faire un test entre modèles emboités est équivalent à tester la nullité de
certains coefficients du grand modèle. On peut ainsi utiliser les tests de Wald, du maximum de
vraisemblance ou du score pour tester deux modèles emboités. Si par exemple, M 1 ⊂ M2, alors
on a
—2(Ln ( 1) — Ln ( 2 )) 2
χp −p 2 1

M →
L
M
ˆ ˆ

où pj désigne la dimension du modèle Mj.

— DM1 DM2 , c’est pourquoi ce test est

Remarque 3.1 La statistique de test peut s’écrire
également appelé test de déviance.

3.1.3 Critère de choix de modèles

Le test que nous venons d’étudier permet de sélectionner un modèle parmi deux modèles
emboîtés. Or, à partir de p variables explicatives, il est possible de définir un grand nombre de
modèles logistiques qui ne sont pas forcément emboîtés. L’utilisation d’un simple test de déviance
se révèle alors insuffisante. On a recours à des critères de choix de modèles qui permettent de
comparer des modèles qui ne sont pas forcément emboîtés les uns dans les autres.

Les critères AIC et BIC sont les plus utilisés. Ces critères sont basés sur la philosophie suivante :
plus la vraisemblance est grande, plus grande est donc la log-vraisemblance et meilleur est le
modèle (en terme d’ajustement). Cependant la vraisemblance augmente avec la complexité du
modèle, et choisir le modèle qui maximise la vraisemblance revient à choisir le modèle saturé. Ce
modèle est clairement sur-paramétré, il “sur-ajuste” les données (overfitting).
Modèle logistique Laurent Rouvière
 40 Sélection et validation de modèles

Exemple 3.3 On considère un échantillon de taille n = 100 simulé suivant le modèle :

Xi ∼ N (0, 1), 1U ≤0.25 si Xi ≤ 0

~ U[0, 1], et Yi
i

1Ui≥0.25 si Xi ≥ 0.
Ui =

Les données sont représentées sur la figure 3.1 : environ 3/4 des labels valent 0 pour les valeurs
de Xi négatives et 1 pour les valeurs positives. Le modèle saturé ajuste parfaitement les
observations. Nous voyons cependant qu’il est difficile, pour ne pas dire impossible à utiliser dans
un contexte de prévision. De plus le modèle saturé possède ici n = 100 paramètres tandis que le
modèle logis- tique n’en possède que 2. Ceci est nettement plus avantageux pour expliquer Y d’un
point de vue descriptif.

Pour choisir des modèles plus parcimonieux, une stratégie consiste à pénaliser la vraisemblance
par une fonction du nombre de paramètres.
— Par définition l’AIC (Akaike Informative Criterion) pour un modèle M de dimension p est
défini par

AIC(M) = —2 L n ( Mˆ ) + 2p.

— Le critère de choix de modèle le BIC (Bayesian Informative Criterion) pour un modèle M

de dimension p est défini par

BIC(M) = —2 L n ( Mˆ ) + p log(n).

On choisira le modèle qui possède le plus petit AIC ou BIC. L’utilisation de ces critères est
simple. Pour chaque modèle concurrent le critère de choix de modèles est calculé et le modèle qui
présente le plus faible est sélectionné.

Remarque 3.2 — Remarquons que certains logiciels utilisent— AIC et BIC il est donc
prudent de bien vérifier dans quel sens doivent être optimisés ces critères (maximisation ou
minimisation). Ceci peut être fait aisément en comparant un modèle très mauvais (composé
uniquement de la constante par exemple) à un bon modèle et de vérifier dans quel sens
varient les critères de choix.
— Les pénalités ne sont pas choisies au hasard... Le critère BIC est issue de la théorie Bayé-
sienne. En deux mots, les modèles candidats sont vus comme des variables aléatoires sur
lesquels on met une loi de probabilité a priori. Dans ce contexte, choisir le modèle qui
possède le plus petit BIC revient à choisir le modèle qui maximise la probabilité a
posteriori.

3.1.4 Apprentissage/validation
Ce critère mesure la performance d’un modèle en terme de prévision. Il convient donc de définir
au préalable une règle de prévision pour un modèle logistique. Un modèle Mβ fournit une
estimation
pˆβ (x) = pβˆ(x), il est naturel de définir de définir une règle de gˆ β à partir de cette
prévision estimation :
( ≥
gˆ β x) = 1 si pˆβ (x) s
0 sinon,

Laurent Rouvière Modèle logistique

3.1 Sélection ou choix de modèle 41
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.

Modèle logistique Laurent Rouvière

42 Sélection et validation de modèles

+ + + + + ++ + ++ +++ ++ +++ + ++ + + + + + ++ + ++ +++ ++ +++ + ++

1.

1.
0

0
+ + + + + +

Y
0.

0.
5

5
+ + + + ++ +++++ +++ + ++ + + + + + + +++++++ +++ + ++ + +
0.

0.
0

0
+ +
−2 0 2 −2 0 2

X X

FIGURE 3.1 – Gauche : Représentation des observations (gauche). Droite : Représentation des
modèles saturés (pointillés) et logistique (trait plein).
Définition 3.2 Etant donné gˆ : Rp+1 → {0, 1} une règle de prévision construite à partir d’un
échantillon Ðn = (X1, Y1), . . . , (Xn, Yn), on définit la probabilité d’erreur de gˆ par
L(gˆ) = P(gˆ(X) /= Y |Ðn).
Etant donné K règles gˆ k , k = 1, . . . , K, l’approche consiste à
1. estimer les probabilités d’erreur de toutes les règles candidates à l’aide de l’échantillon ;
2. choisir la règle qui possède la plus petite estimation.
Toute la difficulté est de trouver un “bon” estimateur de L(gˆ). Une première idée serait d’utiliser
n
1Σ
ˆ
L (gˆ) = 1 {gˆ(X )/=Y }.
n i i

i=1

Comme le modèle saturé ajuste de manière “parfaite” les données, on a Lˆ ( gˆ s a t ) = 0. Cette

procédure conduirait à sélectionner de manière quasi-systématique le modèle saturé. La faiblesse
de cette approche tient du fait que le même échantillon (X1, Y1), . . . , (Xn, Yn) est utilisé pour :
— construire le modèle (ou la règle) ;
— estimer la probabilité d’erreur.
Ceci introduit un biais dans l’estimation de la probabilité d’erreur. La procédure apprentissage-
validation s’affranchit de ce problème en séparant de manière aléatoire les données (X1, Y1), . . . ,
(Xn, Yn) en deux parties distinctes :
— Ðℓ = {(Xi, Yi), i ∈ℓ Iun échantillon d’apprentissage de taille ℓ ;
— }Ðm = {(Xi, Yi), i ∈m Iun échantillon de validation de taille m tel que ℓ + m = n,
où Iℓ ∪}Im = 1,
{ . . . , n} et Iℓ ∩ Im = . L’échantillon d’apprentissage est utilisé pour construire les
∅
modèles concurrents (pour estimer les paramètres des différents modèles logistiques envisagés) à
partir desquels on définit les règles de prévision candidates. L’échantillon de validation est ensuite
utilisé pour estimer les probabilités d’erreur de chaque règle (voir figure 3.2).
Plus précisément, une fois les paramètres des modèles candidats estimés sur l’échantillon
d’appren- tissage, on construit les règles de prévision gˆ k , k = 1, . . . , K associées à ces modèles.
La probabilité d’erreur d’une règle gˆ k
L(gˆk ) = P(gˆk (X) /= Y |Ðℓ)

Laurent Rouvière Modèle logistique

 3.1 Sélection ou choix de modèle 43

est ensuite estimée à l’aide de l’échantillon de validation :

1 Σ
ˆ
L (gˆ )=
k (X )/ .
{gˆ 1
m i∈£ =Y } i
k i

m
On choisira la règle qui minimise Lˆ ( gˆ k ) . Il est facile de voir que Lˆ ( gˆ k ) est un estimateur sans
biais de L(gˆk ) puisque, conditionnellement à Ðℓ, m Lˆ ( gˆ k ) suit une loi Binomiale Bin(m, L(gˆk )).

Apprentissage

Séparation
Toutes les variables
Estimations des modèles
Y
concurrents
Y X

Y Validation
Uniquement les X
Données de départ

Valeurs observées Valeurs prédites (pour tous les modèles concurrents)

FIGURE 3.2 – Procédure d’apprentissage/validation.

Le tableau 3.1 compare les probabilités d’erreur estimées sur les données de l’exemple de la figure
3.1. La procédure qui utilise un seul échantillon pour estimer ces probabilités va sélectionner
le modèle saturé, ce n’est pas le cas de la procédure Apprentissage-Validation qui fournit des
estimations des taux d’erreurs plus précises et qui sélectionnera le modèle logistique.
Saturé Logistique
Sans AV 0 0.146
Avec AV 0.244 0.160
TaBLE 3.1 – Pourcentages de mal classés des modèles saturés et logistique de l’exemple de la
Figure
3.1 avec et sans la procédure apprentissage-validation (les deux échantillons sont de
même taille). Cette procédure semble la plus indiquée pour choisir un modèle. Il faut néanmoins
la nuancer car elle requiert beaucoup de données
— dans l’échantillon d’apprentissage pour estimer convenablement les paramètres du modèle
et ainsi ne pas trop pénaliser les modèles avec beaucoup de variables dont les coefficients
seront moins bien estimés ;
— dans l’échantillon de validation pour bien évaluer la probabilité d’erreur des règles.
De plus il n’existe pas de règle pour choisir les tailles des deux échantillons.

3.1.5 Validation croisée

Lorsque l’on n’a pas assez de données pour l’apprentissage/validation, on peut avoir recours à une
procédure de validation croisée. Le principe est de “moyenner” le pourcentage de mal classés à

Modèle logistique Laurent Rouvière

44 Sélection et validation de modèles

l’aide de plusieurs découpages de l’échantillon. Plus précisément, on divise l’échantillon initial en

K sous échantillons Ek de même taille et on effectue K procédures apprentissage-validation
pour lesquelles :
— l’échantillon test sera constitué d’une division Ek ;
— l’échantillon d’apprentissage sera constitué de l’ensemble des autres divisions — E Ek
(voir figure 3.3).

Y X

E1
E2

Ek

EK
FIGURE 3.3 – Découpage de l’échantillon pour la validation croisée. L’échantillon
d’apprentissage correspond à la partie hachurée.
On obtient ainsi une prévision pour chaque individu de la division Ek et une fois les K procédures
apprentissage-validation effectuées, on a une prévision pour tous les individus de l’échantillon.
Il suffit alors de comparer ces prévisions aux valeurs observées pour obtenir une estimation de
la probabilité d’erreur de la règle. Le modèle retenu sera le modèle qui conduit à l’estimation
minimale.
Bien entendu le choix du nombre K de parties n’est pas anodin.
— Plus K est faible, plus la capacité de prévision sera évaluée dans de nombreux cas puisque
le nombre d’observations dans la validation sera élevé, mais moins l’estimation sera
précise puisque moins de données seront utilisées pour estimer les paramètres du modèle ;
— Au contraire, un K élevé conduit à peu d’observations dans la validation et donc à une plus
grande variance dans l’estimation de la probabilité d’erreur.

Remarque 3.3 Sous R, la librairie boot permet d’estimer le pourcentage de mal classées par
validation croisée. Si, par exemple, on considère le modèle composé des 6 variables explicatives sur
les données du cancer de la prostate, on obtient :

> modele <- glm(Y~.,data=donnees,family=binomial)

> library(boot)
> cout <- function(Y_obs,prevision_prob){
+ return(mean(abs(Y_obs-prevision_prob)>0.5))}
> [Link](donnees,modele,cout)$delta[1]
1
0.3396226

Laurent Rouvière Modèle logistique

3.1 Sélection ou choix de modèle 45

Modèle de départ

Modèle en cours = M0

AIC M0 moins bon Ajout d’un coefficient

M1 devient M0
Choix parmi tous les modèles (+ petit AIC)

Modèle sélectionné =M1

Comparaison AIC modele M0 et modele M1

AIC M0 meilleur

Modèle courant M0 retenu

FIGURE 3.4 – Technique ascendante utilisant l’AIC.

3.1.6 Sélection automatique

Les procédures que nous venons d’étudier permettent de sélectionner un modèle à partir d’une
famille de modèles donnée. Une autre approche de la sélection de modèle consiste à chercher
parmi les variables X1, . . . , Xp, celles qui “expliquent le mieux” Y . Par exemple, pour la régression
logistique, nous pourrions nous poser le problème de chercher le meilleur sous-ensemble des
p variables explicatives pour un critère C donné (AIC, BIC...). Le nombre de sous ensembles
de p variables étant 2p, nous serions en présence de 2p modèles logistiques possibles, c’est-à-
dire 2p modèles différents. Bien entendu, nous sélectionnerions le modèle qui optimiserait le
critère C. Cependant, dans de nombreuses situations, p est grand et par conséquent le nombre
de modèles considérés est “très grand”. Les algorithmes d’optimisation du critère C
deviennent très coûteux en temps de calcul. On préfère alors souvent utiliser des méthodes de
recherche pas à pas.

Recherche pas à pas, méthode ascendante (forward selection)

A chaque pas, une variable est ajoutée au modèle.

— Si la méthode ascendante utilise un test de déviance, nous rajoutons la variable Xj dont la
valeur p (probabilité critique) associée à la statistique de test de déviance qui compare les
2 modèles est minimale. Nous nous arrêtons lorsque toutes les variables sont intégrées ou
lorsque la valeur p est plus grande qu’une valeur seuil.
— Si la méthode ascendante utilise un critère de choix, nous ajoutons la variable Xj dont
l’ajout au modèle conduit à l’optimisation la plus grande du critère de choix. Nous nous
arrêtons lorsque toutes les variables sont intégrées ou lorsque qu’aucune variable ne permet
l’optimisation du critère de choix (voir aussi Figure 3.4).

Modèle logistique Laurent Rouvière

46 Sélection et validation de modèles

Recherche pas à pas, méthode descendante (backward selection)

A la première étape toutes les variables sont intégrées au modèle.
— Si la méthode descendante utilise un test de déviance, nous éliminons ensuite la variable
Xj dont la valeur p associée à la statistique de test de déviance est la plus grande. Nous
nous arrêtons lorsque toutes les variables sont retirées du modèle ou lorsque la valeur p est
plus petite qu’une valeur seuil.
— Si la méthode descendante utilise un critère de choix, nous retirons la variable Xj dont le
retrait du modèle conduit à l’augmentation la plus grande du critère de choix. Nous nous
arrêtons lorsque toutes les variables sont retirées ou lorsque qu’aucune variable ne permet
l’augmentation du critère de choix.

Recherche pas à pas, méthode progressive (stepwise selection)

Idem que l’ascendante, sauf que l’on peut éliminer des variables déjà introduites. En effet, il peut
arriver que des variables introduites au début de l’algorithme ne soient plus significatives après
introduction de nouvelles variables. Remarquons qu’en général la variable “constante” est toujours
présente dans le modèle.

Exemple 3.4 Reprenons l’exemple des données du cancer de la prostate. Nous allons
sélectionner des modèles par les différentes approches pas à pas.
1. Méthode ascendante : le modèle initial est constitué uniquement de la variable age.
> model_age<-glm(Y~age,data=donnees,family=binomial)
> model_asc<-step(model_age,direction="forward",scope=list(upper=
formula("Y~(age+acide+rayonx+taille+grade+[Link].)")))
> model_asc

Call: glm(formula = Y ~ age + rayonx + taille + [Link]., family = binomial,

data = donnees)

Coefficients:
(Intercept) age rayonx1 taille1 [Link].
2.65636 -0.06523 2.08995 1.75652 2.34941

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

Null Deviance: 70.25
Residual Deviance: 47.68 AIC: 57.68

2. Méthode descendante : le modèle initial est ici constitué de toutes les variables
(sans interactions).
> modelcomplet<-glm(Y~.,data=donnees,family=binomial)
> model_des<-step(modelcomplet,direction="backward")
> model_des

Call: glm(formula = Y ~ acide + rayonx + taille + [Link]., family = binomial,

data = donnees)

Coefficients:
(Intercept) acide rayonx1 taille1 [Link].
9.067 -9.862 2.093 1.591 10.410

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

Laurent Rouvière Modèle logistique

3.1 Sélection ou choix de modèle 47

Null Deviance: 70.25

Residual Deviance: 46.43 AIC: 56.43

3. Méthode progressive : le modèle initial est ici constitué de toutes les variables
(sans interactions).

> model_pro<-step(modelcomplet,direction="both")
> model_pro

Call: glm(formula = Y ~ acide + rayonx + taille + [Link]., family = binomial,

data = donnees)

Coefficients:
(Intercept) acide rayonx1 taille1 [Link].
9.067 -9.862 2.093 1.591 10.410

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

Null Deviance: 70.25
Residual Deviance: 46.43 AIC: 56.43

On peut également mettre des variables d’interactions parmi les variables candidates.

> model_pro1<-step(modelcomplet,direction="both",scope=list(upper=formula("Y~(age+acide+
rayonx+taille+grade+[Link].)^2"),lower=formula("Y~1")))
> model_pro1

Call: glm(formula = Y ~ acide + rayonx + taille + grade + [Link]. + taille:grade +

taille:[Link]. + acide:grade, family = binomial,data = donnees)

Coefficients:
(Intercept) acide rayonx1 taille1
49.385 -49.186 3.135 -2.635
grade1 [Link]. taille1:grade1 taille1:[Link].
1.227 53.329 -14.264 -21.719
acide:grade1
17.629

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

Null Deviance: 70.25
Residual Deviance: 26.47 AIC: 44.47

Nous voyons sur cet exemple que suivant le choix de la méthode pas à pas et du modèle initial, les
modèles sélectionnés diffèrent. La sélection d’un seul modèle peut s’effectuer en deux temps :

1. On sélectionne un nombre faible (entre 5 et 10 par exemple) de modèles candidats via des
algorithmes pas à pas ;

2. On choisit le modèle qui minimise un critère de choix (AIC, BIC, ou minimisation de la

probabilité d’erreur).

Une fois le modèle choisi, il est nécessaire de mener une étude plus approfondie de ce dernier qui
permettra de le “valider” ou de l’affiner (suppression de points aberrants, analyse des résidus....).

Modèle logistique Laurent Rouvière

48 Sélection et validation de modèles

3.2 Validation du modèle

3.2.1 Test d’adéquation de la déviance
Ce test permet de mesurer l’ajustement d’un modèle. L’idée est simple. Nous avons vu que le
modèle saturé était le meilleur modèle en terme de qualité d’ajustement. Pour mesurer
l’ajustement d’un modèle Mβ, nous allons le comparer au modèle saturé en effectuant un test de
rapport de vraisemblance (appelé déviance ici). Ce test n’est valable qu’en présence de données
répétées.

Généralement, on écrit les hypothèses d’un test entre modèles emboités en terme de nullité de
certains coefficients (ceux du grand modèle qui ne sont pas dans le petit). Il est difficile de
présenter dans un cadre général une telle écriture puisque l’on a pas d’écriture générale de modèle
saturé. C’est pourquoi, on commettra l’abus d’écrire les hypothèses sous cette forme
— H0 : le modèle considéré à Mβ de dimension p est adéquat ;
— H1 : le modèle saturé adéquat.
Cependant, dans la plupart des cas, on pourra écrire les hypothèses de ce test en terme de nullité
de certains coefficients du modèle saturé.

Exemple 3.5 Considérons l’exemple où on dispose de deux variables explicatives X1 et X2 à

deux modalités 0 et 1. On souhaite tester le modèle Mβ

logit pβ(x) = β0 + β11x1=1 + β21x2 =1

contre le modèle saturé qui (sous réserve que tous les croisements entre X1 et X2 soient observés)
va ici s’écrire
logit psat(x) = γ0 + γ11x1=1 + γ21x2=1 + γ31x1 =11x2=1.

Sur cet exemple, les hypothèses du test de deviance vont s’écrire H0 : γ3 = 0 contre H1 : γ3 /= 0.

Le test de la déviance compare le modèle saturé au modèle considéré au moyen de la déviance.

Nous savons que
— si la déviance est grande, alors le modèle considéré est loin du modèle saturé et que par
conséquent il n’ajuste pas bien les données ;
— Par contre si la déviance est proche de 0, le modèle considéré sera adéquat.
Pour quantifier cette notion de “proche de 0” et de “grande déviance”, la loi de la déviance sous
H0 (le modèle considéré est le vrai modèle) va nous être utile. En effet si H0 est vraie, le modèle
considéré est vrai par définition. La déviance sera répartie sur R + mais avec plus de chance d’être
proche de 0. Par contre si H0 n’est pas vraie la déviance sera répartie sur R + mais avec plus de
chance d’être éloignée de 0. Il nous faut donc connaître la loi de la déviance sous H0.

En présence de données répétées, si le nombre de répétitions nt de chaque point du design tend

vers ∞, alors sous H0 la déviance DMβ = 2(Lsat — L n (βˆ )) converge en loi vers un χT 2−p où p est la
dimension de Mβ. Ainsi, on niveau α, on rejettera H0 si la valeur observée de DMβ est supérieure
au quantile d’ordre 1 — α de la loi χ2T −p (voir Figure 3.5).

Laurent Rouvière Modèle logistique

3.2 Validation du modèle 47

H0 conservé H0 repoussé

densité
[Link].150.200.25
0 2 4 6 8 10
D
FIGURE 3.5 – Test de déviance, la droite horizontale représente le seuil de rejet Dc = q1−α(T —
p).

Remarque 3.4 La validité de la loi et donc du test n’est qu’asymptotique, il est nécessaire
d’avoir un peu de recul quant aux conclusions. Ce test ne peut être utilisé qu’en présence de
données répétées. En effet, l’approximation de la loi de la déviance par une loi du χ2 est d’autant
plus valable lorsque le nombre de répétitions aux points du design est grand. En présence de
données
individuelles (aucune répétition sur les points du design), DMβ ne suit pas une loi du χ2 : le test
d’adéquation d’Hosmer Lemeshow est alors conseillé.

3.2.2 Test d’Hosmer Lemeshow

Ce test permet de vérifier l’adéquation d’un modèle Mβ en présence de données individuelles.
L’approche consiste à 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 pˆβ (x i ) sont ordonnées par ordre croissant (pˆβ (x i ) est la probabilité
Pβ(Y = 1|X = xi) estimée par le modèle) ;
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 pˆβ (x i ) dans le groupe k.
La statistique de test est alors
K 2
2 Σ
C = (ok — mkµk)
mkµk(1 — µk) .
k=1

Le test se conduit de manière identique au test de déviance, la statistique C2 suivant approxima-

tivement sous H0 un χ2 à K — 1 degrés de liberté.

3.2.3 Analyse des résidus

On se donne Mβ un modèle logistique. Pour simplifier les notations on écrira pi = pβ(xi) et

Modèle logistique Laurent Rouvière

 48 Sélection et validation de modèles
pˆi = pβˆ(xi).

Laurent Rouvière Modèle logistique

3.2 Validation du modèle 49

Les différents types de résidus

A l’image de la régression plusieurs types de résidus sont proposés par les logiciels. Le premier,
le plus simple à calculer est tout simplement Yi —pˆ i . Ces résidus sont appelés résidus bruts. Ils
permettent de mesurer l’ajustement du modèle sur chaque observation. Ces résidus n’ayant pas la
même variance, ils sont difficiles à comparer. En effet, on rappelle que Vβ(Y| X = xi) = pi(1 pi).
Par conséquent, la variance de tels résidus risquent d’être élevées pour des valeurs de pi proches
de 1/2. Un moyen de pallier à cette difficulté est de considérer les résidus de Pearson
Y pˆ
RPi = √ ii
— i i .
pˆ (1 — pˆ ) (3.1)

Par définition on standardise les résidus par la variance théorique de Yi.

Remarque 3.5 En présence de données répétées ces résidus s’écrivent
t
—
RPt = √ t t nt pt ˆt
Y
nˆ pˆ (1 —
.
pˆ )
ΣT
Pour nt grand, RPt ∼ N (0, 1) si le modèle est valide et la statistique t=1 RP t2 ∼ χT2 −p où p est
la dimension du modèle. Cette statistique est appelée Statistique de Pearson, d’où le nom de résidu
de Pearson. Cette statistique peut être utilisée pour construire un test d’adequation du modèle en
question en présence de données répétées.
Cependant, comme pˆi est aléatoire, il est évident que Vβ(Yi — pi) /= Vβ(Yi — pˆi ). En effet, en notant
,
 εi = Yi — pi
,
εˆi = Yi — pˆi
on a
Hypothèses Réalité
E(εi) = 0 E(εˆi ) ≃ 0

V(εi) = pi(1 — pi) V ( εˆi ) ≃ pi(1 — pi)(1 — hii)

où hii est l’élément de la ième ligne et de la ième colonne de la matrice H = X(X′WβˆX)−1X′Wβˆ.
Il est par conséquent intéressant de considérer la version standardisée des résidus de Pearson
Yi — p ˆ i
√ i i ii .
pˆ (1 — pˆ )(1 — h )
Ces résidus seront en effet plus facile à analyser (leur distribution étant “presque” centrée réduite).

Les résidus de déviance sont définis par

q √
ˆ
rdi = signe(Yi — pˆi ) 2L1(Yi, β ) = signe(Yi — pˆi ) 2(Yi log pˆi + (1 — Yi) log(1 — pˆi )).
Là encore pour tenir compte de la variabilité ces résidus sont standardisés :
rdi
√ .
1 — hii
Ces deux types de résidus de déviance sont ceux qui sont en général conseillés.

Modèle logistique Laurent Rouvière

 50 Sélection et validation de modèles

Remarque 3.6 1. En présence de données individuelles, nous n’avons pas d’informations

pré- cises sur la loi de ces résidus. On sait juste qu’ils sont approximativement d’espérance
nulle et de variance 1.
2. En présence de données répétées,
s les résidus de déviance s’écrivent
y¯ t nt — y t
rd = 2 y log + (n — y ) log ;
n t — nt pˆt
t t t t
pˆt
ΣT
on retrouve bien que la déviance est égale à rd2. Par conséquent, lorsque les nombres
t=1 t
de répétitions nt sont grands, les résidus de déviance suivent approximativement une loi
N (0, 1). On peut ainsi les analyser de la même manière que dans le modèle linéaire.
3. Tout comme dans le modèle linéaire, on peut également étudier une version studentisée de
ces résidus en estimant les quantités pˆi par validation croisée.

Examen des résidus

Index plot Pour le modèle logistique les résidus de déviance sont souvent préférés. De nom-
breuses études expérimentales ont montré qu’ils approchent mieux la loi normale que les résidus
de Pearson. Pour cette raison ces résidus prennent généralement des valeurs qui varient entre -2
et 2. Nous pourrons construire un index plot pour détecter des valeurs aberrantes. Ce graphique
ordonne les résidus en fonction du numéro de leur observation. Les points pour lesquels on
observe on résidu élevé (hors de [—2, 2] par exemple) devront faire l’objet d’une étude
approfondie.

Exemple 3.6 Voici un exemple de calcul des résidus avec R sur les données du cancer de la
prostate.

> model<-glm(Y~rayonx+taille+grade+[Link].+taille:grade+grade:[Link].,data=donnees,
family=binomial)
> res_dev <- residuals(model) #residus de deviance
> res_pear <- residuals(model,type="pearson") #residus de Pearson
> res_dev_stand <- rstandard(model) #residu de deviance standardises
> H <- influence(model)$hat #diagonale de la hat matrix
> res_pear_stand <- res_pear/sqrt(1-H) #residu de Pearson standardises

> plot(rstudent(model),type="p",cex=0.5,ylab="Résidus studentisés par VC")

> abline(h=c(-2,2))

Graphique prédiction linéaire/résidus Ce graphique qui représente X ′ βˆ en abscisse et

εˆ en ordonné permet de détecter les valeurs aberrantes mais aussi les structurations suspectes. Si
une structuration suspecte apparaît, il sera peut être adéquat d’ajouter une nouvelle variable afin
de prendre en compte cette structuration. Dans le cas des données individuelles ce type de
graphique donne toujours des structurations (Figure 3.7) et n’est donc pas à conseiller.

> plot(predict(model),rstudent(model),type="p",cex=0.5,xlab="prévision linéaire",

ylab="Résidus studentisés par VC")

Laurent Rouvière Modèle logistique

3.2 Validation du modèle 51

26
34

2 1
Résidus studentisés par

0
VC
−
1
−
2

0 10 20 30 40 50

Index

FIGURE 3.6 – Index plot des résidus de déviance studentisés.

2 1
Résidus studentisés par

0
VC
−
1
−
2

−5 0 5 10

prévision linéaire

FIGURE 3.7 – Graphique prédiction/résidus pour un modèle logistique.

Résidus partiels Les résidus partiels sont définis par
εˆP Yi — pˆi
.j = +jβˆ.j X
pˆ (1i — pˆi)
L’analyse consiste à tracer pour toutes les variables j les points avec en abscisse la variable j et
en ordonnée les résidus partiels. Si le tracé est linéaire alors tout est normal. Si par contre une
tendance non linéaire se dégage, il faut remplacer la variable j par une fonction de celle ci donnant
la même tendance que celle observée.
Exemple 3.7 Nous reprenons l’exemple du modèle précédent et traçons sur la figure 3.8 les
résidus partiels pour la variable [Link]..
> residpartiels<-resid(model,type="partial")
> prov<-loess(residpartiels[,"[Link]."]~donnees$[Link])
> ordre<-order(donnees$[Link].)
> plot(donnees$[Link].,residpartiels[,"[Link]."],type="p",cex=0.5,xlab="",ylab="")
> matlines(donnees$[Link].[ordre],predict(prov)[ordre])
> abline(lsfit(donnees$[Link].,residpartiels[,"[Link]."]),lty=2)

Modèle logistique Laurent Rouvière

 52 Sélection et validation de modèles

2
26

15
10
5
0
−
5

−0.5 0.0 0.5

FIGURE 3.8 – Résidus partiels pour la variable [Link]., le trait continu représente le résumé
lissé des données par l’estimateur loess, le trait discontinu représente l’estimateur linéaire
par moindre carré.
Même si l’ajustement n’est pas “parfaitement linéaire”, la figure 3.8 montre qu’aucune transforma-
tion n’est nécessaire, les résidus partiels étant répartis le long de la droite ajustée.

Exemple 3.8 Etudions ici les résidus partiels de la variable age pour l’exemple sur les
données de panne du TP 1 (voir exemple 2.6).

> model <- glm(panne~.,data=donnees,family=binomial)

> residpartiel <- residuals(model,type="partial")
> plot(donnees$age,residpartiel[,"age"],cex=0.5)
> est <- loess(residpartiel[,"age"]~donnees$age)
> ordre <- order(donnees$age)
> matlines(donnees$age[ordre],predict(est)[ordre])
> abline(lsfit(donnees$age,residpartiel[,"age"]),lty=2)
3
2
1
residpartiel[,

0
"age"]
−
1
−
2

5 10 15

donnees$age

FIGURE 3.9 – Résidus partiels pour la variable age.

La figure 3.9 suggère de prendre en compte la variable age2 dans le modèle.

Laurent Rouvière Modèle logistique

3.2 Validation du modèle 53

Mallows (1986) propose d’utiliser les résidus partiels augmentés qui dans certaines situations per-
mettent de mieux dégager cette tendance. Les résidus partiels augmentés pour la jème variable
nécessitent un nouveau modèle logistique identique mis à part le fait qu’une variable explicative
supplémentaire est ajoutée : Xp+1 = X2 jla jème variable élevée au carré. Le nouveau vecteur de
coefficient β du modèle est estimé et les résidus partiels sont alors définis comme
Yi — pˆi
εˆP A = + βˆ X + βˆ X2 .
j .j p+1 .j
.j
pˆi(1 — pˆi )
L’analyse des diagrammes est identique à ceux des résidus partiels. Pour une analyse plus
complète sur l’utilisation des résidus, on pourra se reporter au chapitre 5 de l’ouvrage de Collet
(2003).

3.2.4 Points leviers et points influents

Ces notions sont analogues à celles du modèle linéaire (voir Cornillon & Matzner-Løber (2007),
chapitre 4).

Points leviers
Par définition les points leviers sont les points du design qui déterminent très fortement leur
propre estimation. Nous avons vu que l’algorithme d’estimation des paramètres effectue à chaque
étape une régression linéaire et s’arrête lorsque le processus devient stationnaire :
βˆ = (X′W ˆX)−1X′W ˆZ,
β β

et la prédiction linéaire est

X βˆ = X(X′W ˆX)−1X′W ˆZ = HZ,
β β

où H est une matrice de projection selon la métrique Wβˆ. Comme nous transformons X βˆ par
une fonction monotone, des X βˆ extrêmes entraînent des valeurs de pˆ extrêmes. Nous allons
donc utiliser la même méthode de diagnostic que celle de la régression simple avec une nouvelle
matrice de projection H. Pour la ième prédiction linéaire nous avons
Σ
[Xβˆ] i = HiiZi + HijZj.
j/=i

Si Hii est grand relativement aux Hij, j /= i alors la ième observation contribue fortement à
la construction de [Xβˆ] i . On dira que le “poids” de l’observation i sur sa propre estimation
vaut Hii.

Comme H est un projecteur nous savons que 0 ≤ Hii ≤1. Nous avons les cas extrêmes suivants :
— si Hii = 1, pˆi est entièrement déterminé par Yi car Hij = 0 pour tout j.
— si Hii = 0, Yi n’a pas d’influence sur pˆi .
La trace d’un projecteur étant égale à la dimension du sous espace dans lequel on projette, on
a tr(H) =Σi Hii = p. Donc en moyenne Hii vaut p/n. Pour dire que la valeur de Hii contribue trop
fortement à la construction de pˆ i , il faut un seuil au delà duquel le point est un point levier.
Par
habitude, si Hii > 2p/n ou si Hii > 3p/n alors le ième point est déclaré comme un point levier.
En pratique un tracé de Hii est effectué et l’on cherche les points dont le Hii est supérieur à 3p/n
ou 2p/n. Ces points sont leviers et leur valeur influe fortement sur leur propre prévision.

Modèle logistique Laurent Rouvière

54 Sélection et validation de modèles
Pour le modèle de l’exemple 3.6, on a

Laurent Rouvière Modèle logistique

3.2 Validation du modèle 55

> p<-length(model$coefficients)
> n<-nrow(donnees)
> plot(influence(model)$hat,type="h",ylab="hii")
> seuil1<-3*p/n
> abline(h=seuil1,col=1,lty=2)
> seuil2<-2*p/n
> abline(h=seuil2,col=1,lty=3)

34
0.5

9
0.4
0.3
hii

0.2
0.1
0.0

0 10 20 30 40 50

Index

FIGURE 3.10 – Points leviers.

Points influents
Les points influents sont des points qui influent sur le modèle de telle sorte que si on les enlève,
alors l’estimation des coefficients sera fortement changée. La mesure la plus classique d’influence
est la distance de Cook. Il s’agit d’une distance entre le coefficient estimé avec toutes les
observations et celui estimé avec toutes les observations sauf une. La distance de Cook pour
l’individu i est définie par
1 ˆ r2P iHii
Di = (β — βˆ)′ X′ W ˆ X(βˆ — βˆ) ≈ ,
(i) β (i) 2
p p(1 — Hii)

où rPi est le résidu de Pearson pour le ième individu.

Les distances de Cook sont généralement représentées comme sur la figure 3.11. Si une distance
se révèle grande par rapport aux autres, alors ce point sera considéré comme influent. Il convient
alors de comprendre pourquoi il est influent, soit
— il est levier ;
— il est aberrant ;
— (les deux !)
Dans tous les cas il convient de comprendre si une erreur de mesure, une différence dans la popu-
lation des individus est à l’origine de ce phénomène. Eventuellement pour obtenir des conclusions
robustes il sera bon de refaire l’analyse sans ce(s) point(s).
Toujours pour le modèle de l’exemple 3.6, on représente les distances de Cook avec la commande

> plot([Link](model),type="h",ylab="Distance de Cook")

Modèle logistique Laurent Rouvière

56 Sélection et validation de modèles

34

0.8
0.6
Distance de
Cook
0.4 0.2
0.0

0 10 20 30 40 50

Index

FIGURE 3.11 – Distances de Cook.

Laurent Rouvière Modèle logistique

Chapitre 4

Modèle logistique multi-classes

Nous traitons dans ce chapitre le cas où la variable à expliquer Y prend plus de deux
modalités. Pour simplifier les notations, nous supposerons que Y peut prendre K valeurs 1, . . .
, K. Nous cherchons donc à expliquer Y par le vecteur X = (1, X1, . . . , Xp)′ des variables
explicatives (Xj étant qualitatives, quantitatives ou pouvant représenter des interactions). Nous
distinguons deux cas :
— les modalités de Y sont ordonnées : il existe une hiérarchie naturelle entre elles. Par
exemple le degré de satisfaction relativement à un produit, le degré d’adhésion à une
opinion... En biostatistique, il peut s’agir d’un diagnostic sur l’état de santé (très bonne,
bonne, moyenne, mauvais santé), sur le stade d’évolution d’une maladie, ou encore sur la
taille ou la nature d’une tumeur (tumeur absente, bénigne, ou maligne). On parle dans ce
cas de modèle polytomique ordonné ;
— il n’existe pas de relation d’ordre sur les modalités de Y , la variable à expliquer est
purement nominale : accord pour un prêt (oui, non, examen du dossier). On parle dans ce
cas de modèle polytomique nominal où de modèle polytomique multinomial.
Tout comme dans le cas binaire, nous cherchons ici à modéliser la loi | de Y X = x. Si on
{ }
suppose que Y prend ses valeurs dans 1, . . . , K , le problème va consister à mettre une forme
paramé- trique sur les probabilités πj =| P(Y = j X = x), j = 1, . . . , K. Que les données soient
{
individuelles (yi, xi),} i = 1, . . . , {n ou répétées (yit, xt), i = 1, . .}. , nt, t = 1, . . . , T , on
supposera que les ob- servations sont indépendantes. Dans ce chapitre, nous n’aborderons pas
en détails les méthodes d’estimations ainsi que les comportements asymptotiques des estimateurs.
Une fois les conditions d’identifiabilité des modèles vérifiées, l’estimation des paramètres se fera par
maximum de vraisem- blance. Sous des conditions du même type que celles présentées dans le
cas binaire, on admettra les résultats classiques sur le comportement asymptotique des
estimateurs : variance asymptotique se déduisant de la matrice d’information de Fisher,
normalité asymptotique de l’estimateur du maximum de vraisemblance, tests de Wald, du
rapport de vraisemblance et du score...

4.1 Le modèle saturé ou modèle multinomial

Tout comme dans le cas binaire, ce modèle présente un intérêt essentiellement dans le cas de
données répétées. Ici encore, le modèle saturé ne met pas de contrainte sur la forme des
probabilités πj. Il est donc défini par l’ensemble des probabilités

πjt = P(Yit = j|X = xt), i = 1, . . . , nt; t = 1, . . . , T, j = 1, . . . , K — 1

Modèle logistique Laurent Rouvière

56 Modèle logistique multi-classes
ΣK
(j varie de 1 à K — 1 car ∀t, j=1 πjt = 1). Le modèle saturé est donc de dimension T (K — 1). Les
estimateurs du maximum de vraisemblance sont donnés par

1 n
πˆ j t = Σt

nt 1Y it =j ,
i=1

c’est-à-dire par la proportion des Yit prenant la valeur j au point xt. On note π(t) = (π1t, . . . , π(K−1)t)′.

Propriété 4.1 1. πˆ (t) est un estimateur sans biais de π(t).

1
2. V(πˆ(t)) =
n Σt où Σt(j, j) = πjt(1 — πjt) et Σt(j, ℓ) = —πjtπℓt si j /= ℓ.
t
3. nt → ∞ alors
Si √ nt(πˆ(t) — π(t)) →L N (0, Σt).

Les points 1 et 2 se déduisent directement du fait que ∀t = 1, . . . , T le vecteur aléatoire

nt !
nt ′
Σ
Σ 1Yit=1, . . . 1Y it =K−1
i=1
, i=1

suit une loi multinomiale (nt, π(t)). On comprend bien pourquoi ce modèle présente un intérêt
uniquement en présence deMdonnées répétées. En effet, en présence de données individuelles (nt =
1), chaque paramètre est estimé à partir d’une seule observation et on obtient de très faibles
performances pour les estimateurs.

4.2 Modèle polytomique ordonné

4.2.1 Cas binaire
Nous revenons d’abord au cas où Y est binaire (0 ou 1) et supposons, sans perte de généralité, que
nous sommes en présence d’une seule variable explicative X. On introduit ǫ une variable aléatoire
centrée et une variable latente (non observée) Y ⋆ = β˜0 + β1x + ǫ telle| que Y X = x vaut 1
lorsque
la variable latente Y ⋆ est grande (supérieure à un seuil s) et 0 sinon. Nous obtenons :

P(Y = 1|X = x) = P(β˜0 + β1x + ǫ > s) = P(—ǫ < —s + β˜0 + β1x) = F (β0 + β1x)

où F est la fonction de répartition de la variable — ǫ et β0 = s + β˜0 . Pour finir de spécifier le

modèle, il reste à choisir la fonction de répartition F . Si on choisit

F (x) = 1 exp(x)
= , (4.1)
1 + exp(—x) 1+
exp(x)

on obtient le modèle logistique étudié dans les chapitres précédents. Si F est la fonction de
répar- tition associée à la loi normale centrée réduite, nous obtenons alors le modèle probit (voir
Laurent Rouvière Modèle logistique
4.2 Modèle polytomique ordonné 57
section 1.3 et figure 4.1).

Modèle logistique Laurent Rouvière

 58 Modèle logistique multi-classes

1.0
0.8
0.6
0.4
0.2
0.0

−4 −2 0 2 4

FIGURE 4.1 – Fonctions de répartition des lois normale (trait plein) et logistique (tirets).

4.2.2 Généralisation
Le modèle polytomique ordonné peut être présenté comme une généralisation du modèle dichoto-
mique présenté dans la partie précédente, avec cette fois Y prenant K modalités ordonnées. On
se place toujours dans le cas d’une seule variable explicative X. On introduit non plus un seul
seuil, mais plusieurs seuils α1, . . . , αK−1 déterministes tels que :

1 si Y ⋆ <
(Y |X = x) j si αj−1 ≤ Y ⋆ < αj, j = 2, . . . , K —
α avec Y ⋆ = β1x + ǫ.
, 1
= , 1
k si Y ≥ αK−1
⋆
Le choix de la fonction de répartition logistique (4.1) conduit au modèle :

Pβ(Y ≤ j|X = x) = F (αj — β1x), j = 1, . . . , K — 1

ou encore
logit βpj (x) = logit Pβ(Y ≤ j|X = x) = αj — β1x, j = 1, . . . , K — 1. (4.2)
Si on est en présence de p variables explicatives, le modèle devient

logit βpj (x) = αj — β1x1 — . . . — βpxp, j = 1, . . . , K — 1, (4.3)

ou encore
j exp(αj — β1x1 — . . . — βpxp)
p (x) = 1 + exp(αj — β1x1 — . . . — βpxp).
β
Dans ce modèle, seule la constante diffère suivant les différents niveaux de Y . Ce modèle nécessite
ΣK
l’estimation de p+K —1 coefficients (p pentes et K —1 constantes car j=1 Pβ(Y = j|X = x) = 1).

Remarque 4.1 Suivant le logiciel les coefficients estimés peuvent différer. La procédure
LOGISTIC de SAS estime par exemple les pentes — bj = βj. Sous R les fonctions polr, lmr et
vglm des librairies MASS, Design et VGAM permettent de construire des modèles logistiques pour
expliquer une variable qualitative ordinale. Il est important de consulter l’aide de la fonction afin de
connaître la signification des coefficients estimés.

Laurent Rouvière Modèle logistique

4.2 Modèle polytomique ordonné 59

Exemple 4.1 La fonction polr de la librairie MASS utilise un modèle de la forme (4.2) et (4.3).
Reprenons le jeu de données du TP2.
On pose à 195 étudiants la question : si vous trouvez un portefeuille dans la rue contenant de
l’argent et des papiers :
— vous gardez tout (réponse 1) ;
— vous gardez l’argent et rendez le portefeuille (réponse 2) ;
— vous rendez tout (réponse 3).
On construit alors la variable WALLET telle que
— WALLET=1 si l’étudiant répond 1 ;
— WALLET=2 si l’étudiant répond 2 ;
— WALLET=3 si l’étudiant répond 3
; Pour chaque étudiant, on note :
— Le sexe (variable MALE=1 si homme, 0 si femme) ;
— la nature des études suivies (variable BUSINESS= 1 pour les écoles de commerce, 0 pour les
autres écoles) ;
— l’existence de punitions passées (variable PUNISH=1 si puni seulement à l’école primaire, 2
si puni seulement à l’école primaire et secondaire et 3 si puni seulement à l’école primaire,
secondaire et supérieur) ;
— l’explication ou pas par les parents des punitions reçues dans l’enfance (variable
EXPLAIN=1 si les parents expliquaient, 0 sinon).
On cherche à expliquer la variable WALLET par les autres variables. On note Y = WALLET, X1 =
MALE, X2 = BUSINESS, X3 = PUNISH = 1 et X4 = EXPLAIN. Le modèle s’écrit

logit pj = α — β 1 — β21x — β31x — β41x — β51x

β j 1 x =1 =1 =2 =3 =1.
1 2 3 3 4

On obtient les estimations avec polr.

> library(MASS)
> model <- polr(wallet~.,data=donnees)
> model
Call:
polr(formula = wallet ~ ., data = donnees)

Coefficients:
male1 business1 punish2 punish3 explain1
-1.0598227 -0.7388746 -0.6276423 -1.4030892 1.0518775

Intercepts:
1|2 2|3
-2.5678520 -0.7890143

Residual Deviance: 307.3349

AIC: 321.3349

Tout comme dans le cas binaire, on peut tester la nullité d’un sous ensemble de coefficients à
l’aide des statistiques de Wald, du rapport de vraisemblance et du score. Par exemple, on obtient la
probabilité critique du test du rapport de vraisemblance pour le test H0 : β1 = . . . = β5 = 0
contre H1 : Ej ∈ {1, . . . , 5}, βj /= 0 avec les commandes
> model0 <- polr(wallet~1,data=donnees)
> statRV <- -2*(logLik(model0)-logLik(model))
> 1-pchisq(statRV,df=length(coef(model))-length(coef(model0)))

Modèle logistique Laurent Rouvière

 60 Modèle logistique multi-classes

[1] 1.589732e-08
attr(,"df")
[1] 2
attr(,"class")
[1] "logLik"

4.2.3 L’égalité des pentes

Dans le modèle (4.3), les coefficients des variables explicatives sont identiques quel que soit le
niveau de Y (on dit qu’il y a égalité des pentes) tandis que les constantes diffèrent. Ainsi, dans
ce modèle, quelle que soit la modalité j considérée, une variable explicative donnée a la
même influence sur ≤ Pβ(Y j X = x). Dans le cas où l’on dispose d’une seule variable
|
explicative X, le modèle peut être représenté par la figure 4.2.
ηj(x)
=j 3

=j 2

j=1

x
β p (x) = αj — β1x = ηj(x).
FIGURE 4.2 – Représentation du modèle logit j

Si on envisage des valeurs différentes pour les paramètres de pentes βm (m = 1, . . . , p) alors les
droites de la figure 4.2 vont se couper. Il est facile de voir que ceci remet en cause le caractère
ordonné des modalités de la variable expliquée. En effet, supposons qu’au delà d’une valeur x0, la
première droite (j = 1) se situe au-dessus de la seconde (j = 2), on a alors :

∀x > x0, Pβ(Y ≤ 1|X = x) > Pβ(Y ≤ 2|X = x).

Un tel résultat est évidemment gênant : il remet en cause la modélisation retenue, qui distribue
les différentes modalités de Y sur un même axe ordonné. Il faut dans ce cas se tourner vers un
modèle plus général (voir section 4.3).
Dans le cas où la variable Y est ordonnée, on définit l’odds d’un individu xi relativement à l’évè-
nement {Y ≤ j} par
Pβ(Y ≤ j|X = x)
odds(x ) = = exp(α — x′ β).
i j
1 — Pβ(Y ≤ j|X = x)
L’odds ratio entre deux individus xi et xi′ relativement à l’évènement {Y ≤ j} s’écrit donc

OR(xi, xi′ ) = exp (xi′ — xi)′β . (4.4)

Cet odds ratio ne dépend pas de la modalité j, on dit que l’hypothèse des seuils aléatoires se
traduit par une “hypothèse de proportionnalité des odds ratio”.

Laurent Rouvière Modèle logistique

4.2 Modèle polytomique ordonné 61

4.2.4 Le test d’égalité des pentes

On se pose la question de vérifier si la modélisation en terme de seuil aléatoire est “raisonnable”
vis à vis de nos données. Nous avons vu dans la partie précédente que cette hypothèse se traduit
par une hypothèse d’égalité des pentes ou de proportionnalité des odds ratio. Un moyen de vérifier
si cette hypothèse est raisonnable consiste à tester l’hypothèse d’égalité des pentes. Cette approche
consiste tout simplement à comparer le modèle en question (avec égalité des pentes) à un modèle
où on lève l’égalité des pentes. Il s’agit de considérer les hypothèses
H0 : β1m = . . . = βm
K−1
, ∀m = 1, . . . , p
contre
H1 : Em ∈ {1, . . . , p}, E(j, ℓ) ∈ ({1, . . . , K — 1}) 2tels que β j /= β
ℓ
m m
pour le modèle
pΣ
logit βpj (x) = αj — m
βj xm, j = 1, . . . , k — 1 (4.5)
m=1
Sous H0 (égalité des pentes), nous retrouvons le modèle polytomique ordinal. Ce test peut être
mené à l’aide des statistiques de Wald, du rapport de vraisemblance et du score. Par exemple, la
PROC LOGISTIC de SAS effectue le test du score (“Score Test for Proportional Odds
Assumption”). Si H0 est vérifiée, la pente de la log-vraisemblance pour l’estimateur du maximum
de vraisemblance
contraint γˆ = (αˆ 1 , . . . , αˆ k , βˆ 1 , . . . , βˆ p ) doit être proche de 0. On utilise le fait que sous H0
la
statistique du score
∇L n (γˆ)IˆH−1 ∇L n (γˆ)
0

converge en loi vers un χ p(K−2) (Iˆ H 0 est un estimateur de la matrice d’information de Fisher du
2

modèle). Le nombre de degrés de liberté s’obtient en faisant la différence entre la dimension du

modèle (4.5) (K — 1 + p(K — 1)) et celle du modèle sous H0 (K — 1 + p)
La probabilité critique du test d’égalité des pentes pour la statistique de rapport de vraisemblance
s’obtient facilement sur R. Pour le modèle prenant en compte uniquement les variables Male et
EXPLAIN de l’exemple 4.1, on obtient cette probabilité critique avec les commandes
> library(VGAM)
> model1 <- vglm(wallet~male+explain,data=donnees,cumulative(parallel=TRUE)) #modele sous H0
[1] "head(extra$orig.w)"
NULL
> model2 <- vglm(wallet~male+explain,data=donnees,cumulative(parallel=FALSE))
#modele sans egalite des pentes
[1] "head(extra$orig.w)"
NULL
> statRV <- -2*(logLik(model1)-logLik(model2))
> 1-pchisq(statRV,df=length(coef(model2))-length(coef(model1)))
[1] 0.3702765

Un exemple d’interprétation des coefficients

Tout comme pour le modèle dichotomique, les coefficients s’interprètent en terme d’odds ratio.
Considérons le cas où l’on cherche à expliquer la variable Y qui correspond au type de mention
obtenue au BAC (3 modalités : “AB”, “B”, “TB”) par la variable X qui correspond à la moyenne
en maths au cours des deux premiers trimestres. Supposons que le coefficient β estimé pour le
modèle

Modèle logistique Laurent Rouvière

62 Modèle logistique multi-classes
logit p = αj — βx
j

Laurent Rouvière Modèle logistique

 4.3 Le modèle polytomique nominal 61

soit égal à — log 2 et que pour une note fixée x, le modèle fournisse les estimations :
1 1 1
Pβˆ(Y = ”AB”|X = x) = , Pβˆ(Y = ”B”|X = x) = et Pβˆ(Y = ”TB”|X = x) = .
4 2 4
Dans ce cas, on déduit de (4.4) les tableaux suivants :

Proba odds{Y = j} odds{Y ≤ j} odds{Y ≤ j} odds{Y = j} Proba

AB 1/4 1/3 1/3 AB 2/3 2/3 2/5
B 1/2 1 3 B 6 16/19 16/35
TB 1/4 1/3 +∞ TB +∞ 1/6 1/7
TaBLE 4.1 – Odds Ratio pour l’individu x. TaBLE 4.2 – Odds Ratio pour l’individu x + 1.
{ ≤ }
En effet, la formule (4.4) entraîne que les odds relativement aux évènements Y j sont
multipliés par exp β lorsque la variable x est augmentée d’une unité. Ainsi, la première
colonne du tableau
4.2 s’obtient en multipliant la dernière colonne du tableau 4.1 par 2. On pourra par exemple
dire que si pour la note x, j’ai 1 mention AB pour trois autres mentions, alors pour la note x +
1 j’ai 2 mentions AB pour trois autres mentions.

4.3 Le modèle polytomique nominal

Dans le cas du modèle polytomique ordonné, le caractère ordinal de la variable expliquée
permettait d’introduire une variable latente, des seuils, et d’estimer des probabilités cumulées.
Lorsque la variable à expliquer est purement nominale, cette démarche n’est plus possible.

4.3.1 Le modèle
Désignons par {1, . . . , K} les modalités (non ordonnées) de Y et par X = (1, X1, . . . , Xp) les p
variables explicatives. Tout comme pour le modèle dichotomique, nous allons chercher à modéliser
ΣK
les probabilités P(Y = j|X = x) pour j = 1, . . . , K. Comme j=1 P(Y = j|X = x) = 1, il suffit
de modéliser les probabilités P(Y = j| X = x) pour j = 1, . . . , K 1 par exemple. Nous prenons
ici le groupe K comme groupe témoin et définissons le modèle multinomial par
pj P (Y = j|X = x)
(x)
log β = log β = βj + β j x + . . . + β j x = x′βj, j = 1, . . . , K — 1.
0 1 1 p p
pβK(x) Pβ(Y = K|X = x)
On a alors
j
exp(x′βj)
pβ(x) = Σ K−1 (4.6)
1 + k=1 exp(x′βk)
avec pour convention βK = 0.
Remarque 4.2 — Si K = 2, on retombe sur le modèle logistique dans le cas binaire.
— L’inconvénient majeur de ce modèle est la multiplication des paramètres à estimer (un
vecteur de paramètres βj pour K—1 modalités de Y ) qui rend les résultats plus délicats
à interpréter.
— Bien entendu, le choix du groupe de référence à une importance sur la valeur des paramètres
estimés mais pas sur le modèle : quel que soit le groupe de référence choisi, les probabilités
estimés Pβ(Y = j |X = x) seront bien entendu identiques.
— Sous R, les fonctions multinom et vglm des librairies nnet et VGAM permettent d’ajuster
un modèle polytomique nominal. Sous SAS, on utilise la procédure CATMOD ou la procédure
LOGISTIC avec l’option link=glogit.

Modèle logistique Laurent Rouvière

62 Modèle logistique multi-classes

4.3.2 Estimation et interprétation des paramètres

Les paramètres β j , j = 1, . . . —
, k 1 sont estimés par maximum de vraisemblance. Pour une
obser- vation (x, y), on désigne par y1, . . . , yK un codage disjonctif complet de y, i.e., yj = 1
si y = j, 0 sinon. La vraisemblance s’écrit
L(y, β) = p1 (x)y1 . . . pK(x)yK
β β

où pj(x)
| = P(Y = j X = x) est défini par (4.6). La vraisemblance suit donc une loi
M multinomiale (1, p1(x), . . . , pk(x)). C’est pour cela que ce modèle est également appelé
“modèle polytomique multinomial”. Les estimateurs du maximum de vraisemblance s’obtiennent
une fois de plus en an- nulant les dérivées partielles par rapport aux différents paramètres de la
vraisemblance de l’échan- tillon. Comme pour le cas dichotomique, il n’y a pas de solutions
explicites pour les estimateurs et on a recours à des méthodes numériques pour les calculer. Il
n’y a pas de réelles nouveautés par rapport au cas binaire, l’algorithme est simplement plus
délicat à écrire à cause de la multi- plication du nombre de paramètres. Le lecteur pourra
consulter l’ouvrage d’Hosmer & Lemeshow
(2000) pour plus de détails.
Les odds ratio n’apparaissent généralement pas dans les sorties logiciels pour le modèle
multi- nomial : il faut donc les calculer à la main en prenant garde au codage particulier des
variables explicatives qualitatives. On rappelle que pour un individu x, l’odds d’un
évènement Y = j est égal | au rapport P(Y
/ = j X = x)/P(Y = j X = x). Dans le cas du modèle
|
multinomial, on définit l’odds d’un évènement j1 contre un évènement j2 par
Pβ(Y = j1|X = x)
odds(x, Y = j vs Y = j ) = = exp((βj1 — βj2 )′x).
1 2
Pβ (Y = j2|X = x)
Et pour deux individus xi et xi′ , on définit alors l’odds ratio par

OR(x , x ′, Y = j vs Y = j ) = Pβ(Y = j1|X = xi)/Pβ(Y = j2|X = xi)

i i 1 2
Pβ(Y = j1|X = xi′ )/Pβ(Y = j2|X = xi′ )
= exp((βj1 — βj2 )′(xi — xi′ )).

Ainsi si les deux individus xi et xi′ ne diffèrent que d’une unité pour la variable ℓ, on a
j1 j2
OR(xi, xi′ , Y = j1 vs Y = j2) = exp(βℓ — βℓ ).

Exemple 4.2 Nous reprenons le problème de l’exemple 4.1 et considérons le modèle polytomique
nominal
pβj (x)
log = β j + β j1 x =1 + β j1 x =1 + β j1 x =2 + β j1 x =3 + β j1 x =1 , j = 1, 2.
3 0 1 2 3 4 5
pβ (x) 1
2 3 3 4

On obtient les estimations avec la fonction multinom

> library(nnet)
> model3 <- multinom(wallet~.,data=donnees)
# weights: 21 (12 variable)
initial value 214.229396
iter 10 value 151.045327

Laurent Rouvière Modèle logistique

4.3 Le modèle polytomique nominal 63

final value 150.780162

converged
> model3
Call:
multinom(formula = wallet ~ ., data = donnees)

Coefficients:
(Intercept) male1 business1 punish2 punish3 explain1
2 1.299394 -0.09558412 -0.7635377 -0.8959272 -1.788003 0.7957086
3 2.406209 -1.26720191 -1.1791095 -1.1450957 -2.141172 1.5935325

Residual Deviance: 301.5603

AIC: 325.5603
On remarque que la fonction multinom n’ajuste pas directement le modèle écrit précédemment
puisqu’elle prend comme modalité de référence la première modalité de Y . Il faut transformer
la valeur des estimations pour obtenir les βˆj . La fonction vglm prend quant à elle la troisième
modalité de Y comme modalité de référence :
> model4 <- vglm(wallet~.,data=donnees,multinomial)
[1] "head(extra$orig.w)"
NULL
> model4
Call:
vglm(formula = wallet ~ ., family = multinomial, data = donnees)

Coefficients:
(Intercept):1 (Intercept):2 male1:1 male1:2 business1:1
-2.4062095 -1.1068174 1.2672025 1.1716184 1.1791118
business1:2 punish2:1 punish2:2 punish3:1 punish3:2
0.4155709 1.1450946 0.2491692 2.1411709 0.3531734
explain1:1 explain1:2
-1.5935357 -0.7978215

Degrees of Freedom: 390 Total; 378 Residual

Residual Deviance: 301.5603
Log-likelihood: -150.7802

Tout comme pour les autres modèles, on peut tester la nullité d’un sous ensemble de coefficients à
l’aide des statistiques de Wald, du rapport de vraisemblance et du score. Par exemple, on obtient la
probabilité critique du test du rapport de vraisemblance pour le test H0 : βj = . . . = βj = 0 contre
1 5
H1 : E(j, k) ∈ {1, 2} × {1, . . . , 5}, βkj /= 0 avec les commandes
> model5 <- vglm(wallet~1,data=donnees,multinomial)
[1] "head(extra$orig.w)"
NULL
> statRV <- -2*(logLik(model5)-logLik(model4))
> 1-pchisq(statRV,df=length(coef(model4))-length(coef(model5)))
[1] 2.087963e-07

Exemple 4.3 Nous terminons cette partie par un exemple de calcul d’odds ratio sur R. On
consi- dère le problème d’expliquer Y à trois niveaux 1,2 et 3 par X une variable qualitative à
trois modalités a, b et c. Les données sont obtenues par les commandes :
> Y <- factor(c(rep(3,5),rep(2,5),rep(1,5)))
> [Link](189)
> X <- factor(sample(c(rep("a",5),rep("b",5),rep("c",5))))
> donnees <- [Link](X,Y)

Modèle logistique Laurent Rouvière

64 Modèle logistique multi-classes

On utilise la fonction multinom pour ajuster un modèle logistique multinomial :

> library(nnet)
> model <- multinom(Y~X,data=donnees)

Les coefficients du modèle sont obtenus par :

> beta <- coef(model)
> beta
(Intercept) Xb Xc
2 1.0979772850 -10.5259443 -0.4045297
3 -0.0005551724 -0.4047684 0.6939763

Nous remarquons que le niveau de référence de Y est ici 1. On veut calculer l’odds ratio OR(b, c, Y =
2 vs Y = 3). On calcule d’abord odds(b, Y = 2 vs Y = 3) et odds(c, Y = 2 vs Y = 3)
> oddb23 <- exp(beta[1,1]+beta[1,2]-beta[2,1]-beta[2,2])
> oddc23 <- exp(beta[1,1]+beta[1,3]-beta[2,1]-beta[2,3])

On déduit l’odds ratio cherché

> oddb23/oddc23
[1] 0.0001206436

On aurait pu le calculer directement

> exp(beta[1,2]-beta[2,2]-beta[1,3]+beta[2,3])
[1] 0.0001206436

Laurent Rouvière Modèle logistique

Annexes

A.1 Rappels sur la méthode du maximum de

vraisemblance
Pour plus de précisions, le lecteur pourra se reporter au chapitre 6 de l’ouvrage de Lejeune (2004)
ainsi qu’au polycopié de Cadre (2010).
On considère X1, . . . , Xn un n-échantillon i.i.d. dont la loi mère admet pour densité fθ(x), θ
étant le paramètre à estimer. On désigne par θˆ = θˆ(X1, . . . , Xn) un estimateur de θ.

Théorème A.1 (Inégalité de Cramer-Rao) On suppose que ∈ θ R

θˆ est sans biais.
et que Sous certaines conditions de régularité, on a
1
V (θˆ) ≥ ,
θ
nI(θ)

où I(θ) est l’information de Fisher :

∂
I(θ) = E θ 2
∂θ ln f (X) .

Si un estimateur sans biais atteint la borne de Cramer-Rao, on dit qu’il est efficace.

Supposons maintenant que le paramètre θ à estimer est de dimension supérieure à 1. On introduit

la matrice d’information de Fisher I(θ) symétrique d’ordre k dont l’élément en position (i, j) est :

∂ ∂ i
Eh ∂θ ln f (X, θ)∂θ ln f (X, h ∂2
i =— ln f (X, θ) .
∂θ
θ) E ∂θ
i j i j

On montre alors que pour tout estimateur sans biais θˆ et pour tout u ∈ Rk

V(u′θˆ) ≥ u ′ [I (θ)]
−1
u,
n

où V(u′θˆ) représente la variance de la combinaison linéaire u′θˆ.

Définition A.1 On appelle fonction de vraisemblance de θ pour une réalisation donnée x1, . . . , xn
de l’échantillon, la fonction de θ :
n
Y
L(θ; x1, . . . , xn) = fθ(xi).
i=1
66 Modèle logistique multi-classes

Définition A.2 On appelle estimation du maximum de vraisemblance une valeur θˆ,

s’il en existe une, telle que :
L(θˆ) = sup L(θ).
θ∈Θ

Une telle solution dépend de x1, . . . , xn (θˆ = h(x1, . . . , xn)). La statistique θˆ = h(X1, . . . , Xn) est
appelée estimateur du maximum de vraisemblance (EMV).

Théorème A.2 Soit θˆ l’estimateur du maximum de vraisemblance défini ci dessus. Sous certaines
conditions de régularité, on a :
— θˆ converge presque sûrement vers θ (il est donc asymptotiquement sans biais) ;
— θˆ est asymptotiquement normal :
√
n(θˆ — θ) →L N (0, [I(θ)]
−1
).

La matrice de variance-covariance de θˆ se “rapproche” de

n I [ (θ)]
1 −1
lorsque→
n . On dit que
∞
l’estimateur du maximum de vraisemblance est asymptotiquement efficace.

Laurent Rouvière Modèle logistique

A.2 Echantillonnage Rétrospectif 67

A.2 Echantillonnage Rétrospectif

Une clinique cherche à mesurer l’effet du tabac sur le cancer du poumon. Elle prélève parmi
ses patients un échantillon composé de n1 = 250 personnes atteintes par le cancer et n0 = 250
personnes ne présentant pas la maladie. Les résultats (simulés !)de l’étude sont présentés dans le
tableau suivant :

Fumeur Non fumeur

Non malade 48 202
Malade 208 42
TaBLE A.3 – Résultats de l’enquête.
Le statisticien responsable de l’étude réalise un modèle logistique. Les sorties sur R sont :
Call: glm(formula = Y ~ X, family = binomial)

Coefficients:
(Intercept) Xnon_fumeur
1.466 -2.773

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

Null Deviance: 692.3
Residual Deviance: 499.9 AIC: 503.9

On note Y la variable qui prend pour valeur 1 si l’individu est malade, 0 sinon et X la variable
explicative (fumeur ou non fumeur).
1. Ecrire le modèle logistique. Quelle est la probabilité P(Y = 1| X = fumeur) estimée par
ce modèle ?
2. Ce résultat vous paraît-il surprenant ?
Pour un individu (X, Y ), on définit S la variable indicatrice d’appartenance à l’échantillon de
l’individu, i.e.,
1 si l’individu est dans l’échantillon
S
0 sinon.
=
On définit également :
— τ1 = P(S = 1 Y | = 1, X = x) = P(S = 1 Y = 1) : taux de sondage dans le groupe 1
(malade). |
— τ0 = P(S = 1 Y | = 0, X = x) = P(S = 1 Y = 0) : taux de sondage dans le groupe 0
(malade). |
— p0(x) = P(Y = 1|X = x, S = 1).
— p(x) = P(Y = 1|X = x).
3. Montrer à l’aide du théorème de Bayes que :
τ1
logit p (x) = log + logit p(x).
0
τ0

4. On définit
— π1 = P(Y = 1) la probabilité à priori d’appartenance au groupe 1 ;
— π0 = P(Y = 0) la probabilité à priori d’appartenance au groupe 0 ;

Modèle logistique Laurent Rouvière

68 Modèle logistique multi-classes

Exprimer τ1 en fonction de n1/n (la proportion d’individus du groupe 1 dans l’échantillon),

π1 et P(S = 1) (la probabilité qu’un individu soit dans l’échantillon).
5. Des études préalables ont montré que les probabilités π1 et π0 pouvaient être estimées par
0.005 et 0.995. En déduire une estimation de P(Y = 1|X = fumeur).

Commentaires
Cette exercice nous montre que la manière de constituer l’échantillon d’apprentissage (X1, Y1), . . . ,
(Xn, Yn) est un point important qui ne doit pas être ignoré. On distingue essentiellement deux
schémas d’échantillonnage :
Le schéma de mélange : tous les individus ont la même probabilité d’être sélectionnés
dans l’échantillon. Les observations sont tirées aléatoirement sans distinction des groupes
0 et
1. Dans ce cas, les proportions d’individus par groupe sont sensiblement identiques dans
l’échantillon et dans la population. On peut montrer que les estimateurs du maximum de
vraisemblance des πi sont les quantités ni/n.
Le schéma rétrospectif : les proportions d’individus par groupe ne sont pas les mêmes
dans l’échantillon et dans la population. Ce schéma d’échantillonnage est souvent utilisé
lorsque les probabilités πi (probabilité d’appartenance au groupe i) sont très différentes les
unes des autres. Dans de tels cas, le schéma de mélange conduit à travailler avec des
effectifs trop petits dans certains groupes, alors que le schéma rétrospectif permet de
travailler avec des effectifs comparables. Par exemple, en diagnostic médical, on a souvent
des problèmes de discrimination entre deux groupes où l’un des groupes (1 par
exemple) est associé à une maladie et l’autre est caractéristique de l’absence de cette
maladie. Dans de telles situations, on a bien sûr π1 beaucoup plus petit que π0. L’usage
consiste alors à étudier deux échantillons de taille à peu~ près équivalente (n1 n0), le
premier étant celui des malades, le second celui des individus sains (voir exercice
précédent).

Pour un tel schéma, le modèle logistique appliqué directement sur les données ne nous fournit pas
une estimation de P(Y = 1 X| = x) mais de P(Y = 1 X = x, S = 1). Une propriété
remarquable du modèle logistique | est que l’effet de ce changement de proportion entre
la population et l’échantillon n’intervient dans l’expression de la probabilité P(Y | =1X
= x) qu’au niveau de la constante β0 :

τ1
logit P(Y = 1|X = x) = logit P(Y = 1|X = x, S = 1) — log τ0 .

En pratique : lorsque l’échantillon (X1, Y1), . . . , (Xn, Yn) est constitué selon un schéma rétros-
pectif il faut :
1. Construire le modèle logistique sur cet échantillon. On obtient alors une estimation de
logit Pˆ ( Y = 1|X = x, S = 1) = β0 + β1x1 + . . . + βpxp.
2. On en déduit alors l’estimation de la probabilité recherchée

0 τ1 1 1 p p
logit Pˆ ( Y = 1|X = x) = β — log
τ0 +β x +...+β x .

Laurent Rouvière Modèle logistique

A.3 Exercices 69
Un tel schéma nécessite bien entendu une connaissance a priori des probabilités π1 et π0.

Modèle logistique Laurent Rouvière

 70 Modèle logistique multi-classes

A.3 Exercices
Exercice 1 On se place dans le cas où les données sont présentées sous forme binomiale
(répéti- tions au point de design xt). On note :
— T le nombre de points de design différents {x1, . . . , xT } ;
— nt le nombre de répétitions au point xt ;
Σnt
— st le nombre de succès au point xt : st = i=1 yit ;
— y¯t = n le nombre moyen de succès au point xt.
st

xi nb succès nb échecs moyenne succès

x1 s1 n1 — s1 y¯1
. . . .
xt st nt — st y¯t
. . . .
xT sT nT — sT y¯ T
TaBLE A.4 – Données répétées.
Montrer que la log-vraisemblance s’écrit :
Σ
T
nt ΣT
L(β) = log + nt {y¯t log(p(xt)) + (1 — y¯t ) log(1 — p(xt))}
st t=T

Σ
t=1

T nt n
p(xt)
= log + + log(1 — p(xt))
t=1
st i=1
1 — p(xt)
Σ nt y¯t
log

où p(xt) = P(Y = 1|X = xt).

Exercice 2 (Nombre de paramètres identifiables) On se place dans le cas de deux

variables explicatives A et B de type facteur admettant respectivement deux (a1, a2) et trois (b1,
b2, b3) niveaux. On dispose de 10 individus, les données sont :
   
a1 b3 1
a2 b2 0
 
 a2 b1
 
1 
a b2 1
 a11 b1   
, Y = 0 .

a2 b3   0
 a1 b3
 
0 
a2 b1 1
 a2 b2   1 
a1 b3 0

1. Rappeler brièvement l’algorithme d’estimation des paramètres par la méthode du

maximum de vraisemblance pour le modèle logistique.
2. Ecrire la matrice de codage disjonctif complet des variables explicatives.

Laurent Rouvière Modèle logistique

A.3 Exercices 71
3. Discuter du rang de cette matrice et en déduire le nombre de paramètres identifiables de
manière unique du modèle logistique pour ce jeu de données.

Modèle logistique Laurent Rouvière

 72 Modèle logistique multi-classes

4. Identifier ces paramètres sur les sorties R suivantes :

A<-factor(c("a1","a2","a2","a1","a1","a2","a1","a2","a2","a1"))
B<-factor(c("b3","b2","b1","b2","b1","b3","b3","b1","b2","b3"))
donnees<-[Link](A,B)
Y<-factor(c(1,0,1,1,0,0,0,1,1,0))
model<-glm(Y~A+B,data=donnees,family=binomial)
model

Call: glm(formula = Y ~ A + B, family = binomial, data = donnees)

Coefficients:
(Intercept) Aa2 Bb2 Bb3
5.690e-01 1.882e-01 -6.310e-16 -1.716e+00

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

Null Deviance: 13.86
Residual Deviance: 12.12 AIC: 20.12
5. Refaire les questions 2, 3 et 4 en considérant les variables d’interactions.
model1<-glm(Y~A+B+A:B,data=donnees,family=binomial)
model1

Call: glm(formula = Y ~ A + B + A:B, family = binomial, data = donnees)

Coefficients:
(Intercept) Aa2 Bb2 Bb3 Aa2:Bb2 Aa2:Bb3
-19.57 39.13 39.13 18.87 -58.70 -58.01

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

Null Deviance: 13.86
Residual Deviance: 6.592 AIC: 18.59

6. On considère le modèle model2 dont les sorties R sont

model2<-glm(Y~A:B-1,data=donnees,family=binomial)
model2

Call: glm(formula = Y ~ A:B - 1, family = binomial, data = donnees)

Coefficients:
Aa1:Bb1 Aa2:Bb1 Aa1:Bb2 Aa2:Bb2 Aa1:Bb3 Aa2:Bb3
-1.957e+01 1.957e+01 1.957e+01 -2.748e-16 -6.931e-01 -1.957e+01

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

Null Deviance: 13.86
Residual Deviance: 6.592 AIC: 18.59

Identifier les paramètres estimés à l’aide de la matrice de codage disjonctif complet étudiée à
la question précédente. Montrer que les modèles model1 et model2 sont équivalents
(on pourra montrer que les probabilités P(Y | = 1 X = x) estimées par les deux modèles
sont les mêmes pour tout x).

Exercice 3 Le traitement du cancer de la prostate change si le cancer a atteint ou non les nøeuds
lymphatiques entourant la prostate. Pour éviter une investigation lourde (ouverture de la cavité
abdominale) un certain nombre de variables sont considérées comme explicative de la variable Y

Laurent Rouvière Modèle logistique

A.3 Exercices 73

binaire : Y = 0 le cancer n’a pas atteint le réseau lymphatique, Y = 1 le cancer a atteint le réseau
lymphatique. Le but de cette étude est donc d’essayer d’expliquer Y par les variables suivantes.
— âge du patient au moment du diagnostic age
— le niveau d’acide phosphatase sérique acide, que l’on appellera par la suite niveau d’acidité
— Le résultat d’une analyse par rayon X, 0= négatif, 1=positif rayonx
— La taille de la tumeur, 0=petite, 1=grande, taille
— L’état pathologique de la tumeur déterminée par biopsie, 0=moyen, 1=grave, grade
— Le logarithme népérien du niveau d’acidité [Link]

age acide rayonx taille grade [Link].

1 66 0.48 0 0 0 -0.73396918
2 68 0.56 0 0 0 -0.57981850
3 66 0.50 0 0 0 -0.69314718
4 56 0.52 0 0 0 -0.65392647
5 58 0.50 0 0 0 -0.69314718
6 60 0.49 0 0 0 -0.71334989
7 65 0.46 1 0 0 -0.77652879
8 60 0.62 1 0 0 -0.47803580
9 50 0.56 0 0 1 -0.57981850
10 49 0.55 1 0 0 -0.59783700
TaBLE A.5 – Les 10 premiers individus.

On a construit 10 modèles logistiques dont les sorties R sont présentées dans le tableau A.7.
On souhaite donner une prévision pour la probabilité P(Y| = 1 X = x) pour cinq nouveaux
individus. Reporter dans le tableau suivant les probabilités estimées par les différents
modèles.

Modèle
M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
Individu
(61,0.60,1,0,1,-0.51)
(49,0.86,0,0,1,-0.15)
(67,0.72,1,0,1,-0.33)
(51,0.95,1,1,1,-0.05)
TaBLE A.6 – Probabilités estimés par les différents modèles.

Modèle logistique Laurent Rouvière

 74 Modèle logistique multi-classes

MODELE 1 MODELE 2
glm(formula = Y ~ age + acide, ...) glm(formula = Y ~ [Link]. + rayonx, ...)

Coefficients: Coefficients:
(Intercept) age acide (Intercept) [Link]. rayonx1
1.32811 -0.05639 2.14004 -0.3308 2.0363 2.1020

MODELE 3 MODELE 4
glm(formula = Y ~ taille + grade +
glm(formula = Y ~ taille + grade, ...)
taille:grade, ...)
Coefficients:
Coefficients:
(Intercept) taille1 grade1
(Intercept) taille1 grade1 taille1:grade1
-1.6008 1.4253 0.7293
-2.251 2.588 2.657 -2.860

MODELE 5 MODELE 6
glm(formula = Y ~ taille:grade - 1, ...)
glm(formula = Y ~ age + acide + age:acide, ...)
Coefficients:
taille0:grade0 taille1:grade0
Coefficients:
-2.2513 0.3365
(Intercept) age acide age:acide
-2.61203 0.00675 7.93956 -0.09282
taille0:grade1 taille1:grade1
0.4055 0.1335

MODELE 7 MODELE 8
glm(formula = Y ~ taille:grade - 1, ...)
glm(formula = Y ~ taille:age, ...)
Coefficients:
taille0:grade0 taille1:grade0
Coefficients:
-2.2513 0.3365
(Intercept) taille0:age taille1:age
2.79408 -0.07157 -0.04351
taille0:grade1 taille1:grade1
0.4055 0.1335

MODELE 9 MODELE 10
glm(formula = Y ~ taille:[Link].
+ rayonx:grade - 1, ...)
glm(formula = Y ~ taille:(age + grade), ...)
Coefficients:
taille0:[Link]. taille1:[Link]. Coefficients:
3.3406 0.9356 (Intercept) taille0:age taille1:age
3.17016 -0.09228 -0.04758
rayonx0:grade0 rayonx1:grade0
-0.6776 1.2771 taille0:grade1 taille1:grade1
2.79235 -0.23813
rayonx0:grade1 rayonx1:grade1
0.1129 2.3963

TaBLE A.7 – Les 10 modèles.

Laurent Rouvière Modèle logistique

 A.4 Correction 73

A.4 Correction
Exercice 1 La variable| S X = xt suit une loi binomiale Bin(nt, p(xt)). Dès lors la vraisemblance
s’écrit :
n n
Y Y t
n p(x t)nt (1 — p(xt ))nt−st.
P (S = st |X = xt )
i=1 = t=1
st
On obtient donc pour log-vraisemblance :

Σ
T
nt T
Σ
L(β) = log + [nt log p(xt) + (nt — st) log(1 — p(xt))] .
st t=1
t=1

On pourra consulter le livre de Collet (2003) (page 59) pour plus de détails.

Exercice 2 1. La procédure d’estimation par MV consiste à “mettre à jour” les coefficients

via une régression pondérée : βˆ k + 1 = (X′W kX)−1X′W k Z k . La première difficulté consiste
à définir X dans cet exemple. X est une matrice qualitative, il est impossible de faire des
opérations dessus, on a donc recours à un codage disjonctif complet.
2. La matrice de codage s’écrit
 0 0 0 1 
0 1 0 1 0
 1
0 1 1 0 0 
1 0 0 1 0
 1 0 1 0 0 
0 1 0 0 1 .

 1 0 0 0 1 
 0 1 
1 0 0 
0 1 0 1 0
La matrice X utilisée pour la régression
1 vaut
0 :0 0 1

1 0 0 0 1 
 1 0 1 0 1 0
X  1
1 0 1 1 0 0
1 1 0 0 1 0
= 1 1 0 1 0 0 
.
 1 0 1 0 0 1 
1 1 0 0 0 1 
 
 1 0 1 1 0 0
1 0 1 0 1 0
1 1 0 0 0 1
3. Le nombre de paramètres estimés par le modèle est égal à la dimension de X′X, c’est-à-dire
au nombre de colonnes de X. Pour effectuer la régression, nous devons inverser X′X (on
considère que la matrice des poids W vaut l’identité). Les vecteurs colonnes de X ne sont
clairement pas linéairement indépendants (voir deuxième et troisième colonne), la matrice
X n’est donc pas de plein rang, le modèle est surparamétré. Il est évident que :

Modèle logistique Laurent Rouvière

74 Modèle logistique multi-classes

— la connaissance de la 2ème colonne implique celle de la 3ème ;

— la connaissance des 4ème et 5ème colonne implique celle de la 6ème.
On peut donc “supprimer” des colonnes à X sans “perdre” d’information. Le logiciel R
supprime la première colonne correspondant chaque modalité, ce qui revient à considérer la
matrice :

1 1 0 0 0 1
 1 0 1 0 1 0
X  1 0 1 1 0 0 
1 1 0 0 1 0 
= 1 1 0 1 0 0
 1 0 1 0 0 1 
1 1 0 0 0 1

 1 0 1 1 0 0 
1 0 1 0 1 0 

 1 1 0 0 0 1

Cette nouvelle matrice est de plein rang (égal à 4), le modèle est constitué de 4 paramètres
identifiables de manière unique.

4. Le modèle est ainsi défini par :

x b1 b2 b3
a1 β0 β0 + Bb2 β0 + Bb3
a2 β0 + Aa2 β0 + +Aa2 + Bb2 β0 + Aa2 + Bb3
TaBLE A.8 – Valeur de logit p(x) pour les différentes classes.

5. Dans le cas d’intéraction la matrice de codage “totale” ou “surparamétrée” va s’écrire


1 1 0 0 0 1 0 0 1 0 0 0
 1 0 1 0 1 0 0 0 0 0 1 0 
X  1 0 1 1 0 0 0 0 0 1 0 0
1 1 0 0 1 0 0 1 0 0 0 0
= 1 1 0 1 0 0 1 0 0 0 0 0
 1 0 1 0 0 1 0 0 0 0 0 1
1 1 0 0 0 1 0 0 1 0 0 0 
 
 1 0 1 1 0 0 0 0 0 1 0 0 
 1 0 1 0 1 0 0 0 0 0 1 0
1 1 0 0 0 1 0 0 1 0 0 0

On peut montrer que le rang de cette matrice est m1m2 = 6. Dans sa procédure d’estimation,

Laurent Rouvière Modèle logistique

A.4 Correction 75

R considère la matrice :

1 1 0 0 0 1 0 0 1 0 0 0
1 0 1 0 1 0 0 0 0 0 1 0 
 1 0 1 1 0 0 0 0 0 1 0 0
1 1 0 0 1 0 0 1 0 0 0 0 
1 1 0 1 0 0 1 0 0 0 0
X
0 = 
1 0 1 0 0 1 0 0 0 0 0
1 
1 1 0 0 0 1 0 0 1 0 0 

 0
 1 0 1 1 0 0 0 0 0 1 0

0
1 0 1 0 1 0 0 0 0 0 1
0
1 1 0 0 0 1 0 0 1 0 0
0

Le modèle est défini par :

x b1 b2 b3
a1 β0 β0 + Bb2 β0 + Bb3
a2 β0 + Aa2 β0 + Aa2 + Bb2 + Aa2 : Bb2 β0 + Aa2 + Bb3 + Aa2 : Bb3
TaBLE A.9 – Valeur de logit p(x) pour les différentes classes.
6. En terme de matrice de codage, le modèle model2 ne considère que les variables
d’intérac- tions :
11 00 0 10 0 1 0 0 0 
10 10 1 00 0 0 0 1 0
1 0 1 1 0 0 0 0 0 1 0 0 

1 1 0 0 1 0 0 1 0 0 0 0
1 1 0 1 0 0 1 0 0 0 0 
X 0= 
10 10 0 10 0 0 0 0 1 
11 00 0 10 0 1 0 0 0 

1 0 1 1 0 0 0 0 0 1 0 0
1 0 1 0 1 0 0 0 0 0 1 0

11 00 0 10 0 1 0 0 0
Il est défini par :
x b1 b2 b3
a1 Aa1 : Bb1 Aa1 : Bb2 Aa3 : Bb3
a2 Aa2 : Bb1 Aa2 : Bb2 Aa2 : Bb3
TaBLE A.10 – Valeur de logit p(x) pour les différentes classes.
Pour montrer que ces modèles sont équivalents, il suffit de montrer que les logit sont
identitiques pour chaque classe dans les modèles model2 et model3.

Exercice 3 Correction sous R. Les commandes sont :

data<-[Link]("~/COURS/GLM/TP/[Link]",sep=";")
donnees<-data[,names(data)!="Y"]
Modèle logistique Laurent Rouvière
 76 Modèle logistique multi-classes
donnees[,"rayonx"]<-factor(donnees[,"rayonx"])
donnees[,"taille"]<-factor(donnees[,"taille"])
donnees[,"grade"]<-factor(donnees[,"grade"])

Laurent Rouvière Modèle logistique

A.4 Correction 77

Y<-factor(data[,names(data)=="Y"])

model1<-glm(Y~age+acide,data=donnees,family=binomial) #2 continues

model2<-glm(Y~[Link].+rayonx,data=donnees,family=binomial) #1 cont, 1 quali

model3<-glm(Y~taille+grade,data=donnees,family=binomial) #2 quali

model4<-glm(Y~taille+grade+taille:grade,data=donnees,family=binomial) #2 quali+inter

model5<-glm(Y~taille:grade-1,data=donnees,family=binomial) #= modele 4

model6<-glm(Y~age+acide+age:acide,data=donnees,family=binomial) # 2 conti+inter

model7<-glm(Y~taille:grade-1,data=donnees,family=binomial) #1 inter continue

model8<-glm(Y~taille:age,data=donnees,family=binomial) # 1 inter qualit/conti

model9<-glm(Y~taille:[Link].+rayonx:grade-1,data=donnees,family=binomial)

model10<-glm(Y~taille:(age+grade),data=donnees,family=binomial)

X_test<-[Link](matrix(c(61,0.60,1,0,1,-0.51,49,0.86,0,0,1,-0.15,67,0.72,1,0,1,
-0.33,51,0.95,1,1,1,-0.05),ncol=6,byrow=T))
names(X_test)<-names(donnees)

X_test[,"rayonx"]<-factor(X_test[,"rayonx"])
X_test[,"taille"]<-factor(X_test[,"taille"])
X_test[,"grade"]<-factor(X_test[,"grade"])

PRED<-matrix(0,ncol=10,nrow=4)
for (i in 1:10){
PRED[,i]<-eval(parse(text=paste("predict(model",i,",newdata=X_test,type=’response’)",sep="")))}

Le tabeau cherché est la matrice PRED.

round(PRED,digits=2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.30 0.68 0.29 0.60 0.60 0.30 0.60 0.17 0.67 0.58
[2,] 0.60 0.35 0.29 0.60 0.60 0.65 0.60 0.33 0.40 0.81
[3,] 0.29 0.75 0.29 0.60 0.60 0.28 0.60 0.12 0.78 0.45
[4,] 0.62 0.84 0.64 0.53 0.53 0.69 0.53 0.64 0.91 0.62

Modèle logistique Laurent Rouvière

 78 Modèle logistique multi-classes

Bibliographie

Albert A. & Anderson D. (1984). On the existence of maximum likelihood estimates in logistic
regression models. Biometrika, 71, 1–10.

Antoniadis A., Berruyer J. & Carmona R. (1992). Régression non linéaire et applications. Econo-
mica.

Cadre B. (2010). Statistique, cours de m1. Poly. de cours, 40 pages, http ://[Link]-
[Link]/math/people/[Link]/.

Collet D. (2003). Modelling Binary Data. Chapman & Hall.

Cornillon P. & Matzner-Løber E. (2007). Régression : Théorie et Application. Springer.

Droesbeke J., Lejeune M. & Saporta J. (2007). Modèles Statistiques pour Données Qualitatives.
Technip.

Fahrmeir L. & Kaufmann H. (1985). Consistency and asymptotic normality of the maximum
likelihood estimator in generalized linear models. The Annals of Statistics, 13, 342–368.

Guyon X. (2005). Le modèle linéaire et ses généralisations. Poly. de cours, 83 pages,

http ://[Link]/-Xavier-Guyon-.

Hosmer D. & Lemeshow S. (2000). Applied Logistic Regression. Wiley.

Lejeune M. (2004). Statistique, la Théorie et ses Applications. Springer.

Perreti-Watel P. (2002). Régression sur variables catégorielles. Poly. de cours, 72 pages.

Pommeret D. (2008). Régression sur variables catégorielles et sur variables de comptage. Poly. de
cours, 100 pages.

Laurent Rouvière Modèle logistique

Modèle logistique Laurent Rouvière
 Examen de mod`ele lin´eaire g´en

´eralis´e ENSAI - mars 2015

Sans document, calculatrice fournie par l’´ecole, dur´ee : 1h45

Pr´eambule. Le sujet est compos´e de trois exercices ind´ependants. La qualit´e de la r´edaction sera
prise en compte. Toutes les r´eponses seront donn´ees sur la copie (ne pas rendre le sujet). Tous les
mod`eles ´ecrits devront ˆetre identifiables. Bien que cela puisse ˆetre aberrant dans certaines parties du
devoir, on supposera dans tout le devoir que les hypoth`eses sur la convergence de l’estimateur du
maximum de vraisemblance dans le mod`ele logistique sont v´erifi´ees et que le nombre d’observations est
suffisamment grand pour approcher la loi des estimateurs par la loi limite vue en cours.
Rappels. Lorsque cela est n´ecessaire, on pourra utiliser que le quantile d’ordre 0.975 de la loi normale
centr´ee r´eduite vaut 1.96.
Exercice 1
On souhaite expliquer une variable binaire Y par p variables explicatives X1, . . . , Xp. On supposera que
les variables explicatives sont d´eterministes et r´eelles. On dispose de n observations du couple (X, Y
). On se placera dans le cas de donn´ees r ´e p ´e t ´e e s .
1. Qu’est-ce qu’un mod`ele lin´eaire g´en´eralis´e ? Vous donnerez une d´efinition math´ematique pr
´ecise et concise.
2. Ecrire le mod`ele logistique permettant d’expliquer Y par X. On prendra soin d’ˆetre pr´ecis dans les
notations, notamment de faire la distinction entre les quantit´es al´eatoires et d´eterministes.
3. Ecrire la vraisemblance du mod`ele en fonction des observations et des param`etres du mod`ele.
Re- trouver la vraisemblance du mod`ele logistique vue en cours dans le cas de donn´ees
individuelles.
4. Ecrire le mod`ele satur´e. Ecrire sa vraisemblance et en d´eduire les estimateurs du maximum de
vraisem- blance (on se contentera de trouver les param`etres qui annulent le gradient de la log-
vraisemblance).
5. Ecrire le test d’ad´equation de la d´eviance (on pr´ecisera notamment les hypoth`eses, la
statistique de test ainsi que sa loi (asymptotique) sous H0).
On se place maintenant dans le cas de donn´ees individuelles. On consid`ere le mod`ele de r´egression
logistique binaire.
6. Ecrire la matrice d’information de Fisher du mod`ele en fonction des observations et des
param`etres du mod`ele.
7. Ecrire la loi asymptotique des estimateurs du maximum de vraisemblance (on pr´ecisera les
hypoth`eses sous lesquelles cette propri´et´e asymptotique est vraie).
8. On souhaite tester la nullit´e des deux premiers param`etres du mod`ele (on pourra les noter β1 et
β2). Ecrire la statistique de test du rapport de vraisemblance et donner sa loi sous H0.
9. En pr´esence d’un grand nombre de variables explicatives p, expliquer en quoi il peut-ˆetre int
´eressant de s´electionner des variables. La r´eponse devra se baser sur des arguments pr´ecis et ne
devra pas exc´eder 5 lignes de texte.
Exercice 2
On cherche `a expliquer une variable Y `a 3 modalit´es (non ordonn´ees) not´ees 1,2,3 par 2 variables
explicatives qualitatives X1 (qui prend pour valeurs A ou B) et X2 (qui prend pour valeurs C ou D). Les
donn´ees sont r´esum´ees dans le tableau suivant.
X1 X2 Y =1 Y =2 Y =3
A C 3 89 28
A D 89 23 18
B C 27 18 65
B D 47 65 28
On lira : ”on a observ´e 18 fois la valeur Y = 2 lorsque (X1, X2) = (B, C)”. On cherche `a expliquer la variable
2
Y par X.

1
 1. Donner une ´ecriture math´ematique (ou plutˆot probabiliste) du probl`eme p os´e. On prendra
soin de donner des notations pr´ecises aux donn´ees observ´ees, de pr´eciser les variables al´eatoires
associ´ees `a ces donn´ees, la ou les lois des variables al´eatoires que l’on cherche `a reconstruire
ainsi que les param`etres ou quantit´es que l’on cherche `a estimer.
2. On consid`ere tout d’abord le mod`ele satur´e.
(a) Ecrire ce mod`ele. Combien poss`ede t-il de param`etres identifiables ?
(b) Donner les estimations par maximum de vraisemblance de deux de ses param`etres au choix.
(c) Calculer un intervalle de confiance asymptotique de niveau 95% pour la probabilit´e de l’´ev`enement
{Y = 2} lorsque (X1, X2) = (A, D) (on est toujours dans le mod`ele satur´e).
3. On ajuste ensuite un mod`ele logistique multinomial. Ce mod`ele est a just´e `a l’aide de la proc
CAT- MOD. Les sorties se trouvent dans l’annexe 1.
(a) Ecrire le mod`ele. Combien poss`ede t-il de param`etres identifiables ? Donner les estimations
des param`etres.
(b) Ecrire les hypoth`eses concernant le test associ´e `a la ligne X2 dans le tableau ”Analyse de
variance par maximum de vraisemblance”. Quel est le nombre de degr´es de libert´e
manquant dans cette ligne. Quelle est la conclusion du test.
(c) Rappeler ce que signifie une int´eraction entre deux variables explicatives. Le mod`ele
logistique multinomial premettant d’expliquer Y par X1, X2 ainsi que l’interaction entre X1 et X2
est-il satur´e ? Justifier.
(d) Les log-vraisemblances maximis´ees sont de
— -487.1461 pour le mod`ele logistique multinomial sans interaction ;
— -438.5001 pour le mod`ele logistique multinomial avec interaction.
Effectuer le test d’ad´equation de la d´eviance pour le mod`ele logistique multinomial sans
interaction (on donnera les hypoth`eses, la statistique de test, sa loi sous H0 et la conclusion du
test).

Exercice 3
On dispose de n = 5 observations (x1, y1), . . . , (xn, yn) telles que xi ∈ R et yi ∈ {0, 1}. Les xi correspondent
`a des observations d’une variable X suppos´ee d´eterministe et les yi correspondent `a des observations d’une
variable Y . On cherche `a expliquer Y par X. Les donn´ees sont dans le tableau suivant.

xi 0.47 -0.55 -0.01 1.07 -0.71

yi 1 0 0 1 1

1. Ecrire le mod`ele de r´egression logistique permettant d’expliquer Y par X (on incluera la

constante dans le mod`ele).
2. Les estimation pβˆ(xi) des probabilit´es P(Yi = 1) par ce mod`ele sont

pβˆ(x1) = 0.76, pβˆ(x2) = 0.40, pβˆ(x3) = 0.60, pβˆ(x4) = 0.89, pβˆ(x5) = 0.35.

Calculer la log-vraisemblance maximis´ee L n (βˆ) du mod`ele.

3. On trouve ci-dessous un extrait du tableau des coefficients de ce mod`ele a just´e sur R.
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4383 ??? ??? ???
X 1.5063 ??? ??? ???
(a) Compl´eter les colonnes Std. Error et z value . Si vous ne disposez pas de suffisamment de
temps pour faire les applications num´eriques (elles ne sont pas difficiles mais un peu p´enible),
on pourra se contenter de les ´ecrire sans les effectuer.
(b) Proposer deux proc´edures de test permettant de tester la nullit´e du coefficient associ´e `a la
variable explicative X. On pourra prendre un niveau de test de α = 5%.
(c) Effectuer les deux tests.
2
Annexe 1 : exercice 2
proc catmod data=exercice2 ;
model Y=X1 X2 ;
run ;

Le Système SAS

Procédure CATMOD
18:22 Wednesday, March 11, 2015 1

Synthèse des données

Réponse Y Niveaux de réponse 3

Variable de pondération Néant Populations 4

Table EXERCICE2 Total des valeurs 500

Valeur(s) manquante(s) 0 Observations 500

Profils de population
Echantillon X1 X2 Taille échantillon
1A C 120

2A D 130

3B C 110

4B D 140

Profils de
réponse
Réponse Y

1 1

2 2

3 3

3
 Le Système SAS

Procédure CATMOD

18:22 Wednesday, March 11, 2015 2

Analyse du max. de vraisemblance

Convergence des calculs du maximum de vraisemblance.

Analyse de variance par maximum de vraisemblance

Source DDL Khi-2 Pr > Khi-2
Intercept ? 19.30 <.0001

X1 ? 25.12 <.0001

X2 ? 75.46 <.0001

Rapport de vraisemblance ? 97.29 <.0001

Analyse des valeurs estimées du maximum de vraisemblance

Nombre Erreur Khi-
Paramètre de fonctions Valeur estimée type 2 Pr > Khi-2
Intercept 1 0.0349 0.1401 0.06 0.8032

2 0.4598 0.1208 14.49 0.0001

X1 A 1 0.5781 0.1318 19.23 <.0001

A 2 0.5354 0.1176 20.73 <.0001

X2 C 1 -1.1713 0.1400 69.98 <.0001

C 2 -0.3130 0.1192 6.89 0.0087

4
 Examen de mod`ele lin´eaire g´en

´eralis´e ENSAI - f´evrier 2014

Sans document

Calculatrice fournie par l’´ecole, dur´ee : 1h45

Pr´eambule : Le sujet est compos´e de trois exercices ind´ependants. La qualit´e de la r´edaction sera
prise en compte. Toutes les r´eponses seront donn´ees sur la copie (ne pas rendre le sujet). Tous les
mod`eles ´ecrits devront ˆetre identifiables. Bien que cela puisse ˆetre aberrant dans certaines parties du
devoir, on supposera dans tout le devoir que les hypoth`eses sur la convergence de l’estimateur du
maximum de vraisemblance dans le mod`ele logistique sont v´erifi´ees et que le nombre d’observations est
suffisamment grand pour approcher la loi des estimateurs par la loi limite vue en cours.

Exercice 1 (Questions de cours)

On souhaite expliquer une variable Y par p variables explicatives X1, . . . , Xp. On dispose de n
observations (x1, y1), . . . , (xn, yn) telles que xi = (xi1, . . . , xip) ∈ Rp et yi ∈ {0, 1}. Les xi sont suppos´es
tous diff´erents et d´eterministes.
1. Qu’est ce qu’un mod`ele lin´eaire g´en´eralis´e ? Vous donnerez une d´efinition pr´ecise et concise.
2. Ecrire le mod`ele logistique permettant d’expliquer Y par X.
3. Ce mod`ele est-il un mod`ele lin´eaire g´e n´e r a l i s´e ? Si oui, justifier en donnant ses composantes, si
non, justifier ´egalement.
4. Ecrire la vraisemblance du mod`ele en fonction des xi, yi et des param`etres du mod`ele.
5. Ecrire la matrice d’information de Fisher du mod`ele en fonction des xi et des param`etres du mod`ele.
6. Ecrire la loi asymptotique des estimateurs du maximum de vraisemblance.
7. On souhaite tester la nullit´e des deux premiers param`etres du mod`ele (on pourra les noter β1 et
β2). Ecrire la statistique de test du rapport de vraisemblance et donner sa loi sous H0.
Exercice 2
On cherche `a expliquer une variable Y `a 3 modalit´es ordonn´ees 1 < 2 < 3 par 3 variables
explicatives qualitatives X1 (qui prend pour valeurs A ou B), X2 (qui prend pour valeurs C ou D) et X3
(qui prend pour valeurs E ou F ). Les donn´ees sont r´esum´ees dans le tableau suivant.
X1 X2 X3 Y =1 Y =2 Y =3
A C E 1 35 3
A C F 3 31 7
A D E 39 7 1
A D F 31 6 6
B C E 15 4 35
B C F 18 5 33
B D E 22 27 11
B D F 18 33 9
On lira : ”on a observ´e 31 fois la valeur Y = 2 lorsque (X1, X2, X3) = (A, C, F )”.
1. On consid`ere tout d’abord le mod`ele satur´e.
(a) Ecrire ce mod`ele. Combien poss`ede t-il de param`etres identifiables ?
(b) Donner les estimations par maximum de vraisemblance de deux de ses param`etres au choix.
(c) Calculer un intervalle de confiance asymptotique de niveau 95% pour la probabilit´e de l’´ev`enement
Y = 2 lorsque (X1, X2, X3) = (B, D, E) (on est toujours dans le mod`ele satur´e).
2. On ajuste ensuite un mod`ele logistique multinomial `a l’aide de la proc CATMOD. Les sorties se
trouvent dans l’annexe 1.

1
 (a) Ecrire le mod`ele a just´e et donner les estimations des param`etres.
(b) Ecrire les hypoth`eses concernant le test associ´e `a la ligne X3 dans le tableau ”Analyse de
variance par maximum de vraisemblance”. Quels est le nombre de degr´es de libert´e
manquant dans cette ligne. Quelle est la conclusion du test.
(c) Quel test est r´e a l i s´e `a la ligne ”Rapport de vraisemblance” du tableau ”Analyse de
variance par maximum de vraisemblance” ? On ´ecriera les hypoth`eses, la statistique de test,
sa loi sous H0 ainsi que la conclusion du test. On expliquera ´egalement d ’ o u` viennent les 8 degr
´es de libert´e.
(d) Quelle est la probabilit´e estim´ee par le mod`ele que {Y = 2} lorsque (X1, X2, X3) = (B, C, E) ?
(e) On ajoute maintenant `a ce mod`ele l’interaction entre X1 et X2. Les sorties se trouvent
dans l’annexe 2. Rappeler ce que signifie l’interaction entre deux variables explicatives.
(f) Le mod`ele avec interaction est-il satur´e ? Justifier.
(g) Privil´egiez-vous le mod`ele avec interaction au mod`ele sans interaction ? Justifier en utilisant
au moins deux crit`eres.
3. On consid`ere maintenant le mod`ele logistique ordinal.
(a) Ecrire les diff´erentes ´etapes permettant de d´efinir le mod`ele logistique ordinal (on partira
de la mod´elisation par une variable latente pour arriver `a l’´ecriture des mod`eles).
(b) L’annexe 3 contient les sorties de la proc LOGISTIC. Ecrire le mod`ele a just´e par cette proc
´edure et donner les estimations des param`etres.
(c) Quelle est la probabilit´e estim´ee par ce mod`ele que {Y = 2} lorsque (X1, X2, X3) = (B, C, E) ?
(d) Expliquer en quelques mots le test effectu´e dans le tableau ”Test de score pour l’hypoth`ese
des cotes proportionnelles”. Ecrire en particulier les hypoth`eses math´ematiques du test et pr
´eciser le nombre de degr´es de libert´e manquant dans le tableau ainsi que la conclusion du
test.
(e) Ecrire les hypoth`eses des tests r´ealis´es dans le tableau ”Test de l’hypoth`ese nulle
globale BETA=0”. Donner les degr´es de libert´e manquants ainsi que la conclusion des tests.

Exercice 3
Dans cet exercice, on pourra utiliser que le quantile d’ordre 0.975 de la loi normale centr´ee r´eduite vaut 1.96.

On dispose de n = 5 observations (x1, y1), . . . , (xn, yn) telles que xi ∈ R et yi ∈ {0, 1}. Les xi correspondent
`a des observations d’une variable X suppos´ee d´eterministe et les yi correspondent `a des observations d’une
variable Y . On cherche `a expliquer Y par X. Les donn´ees sont dans le tableau suivant.

xi 0.72 0.88 0.76 0.89 0.46

yi 0 0 1 1 1

1. Ecrire le mod`ele de r´egression logistique permettant d’expliquer Y par X (on incluera la

constante dans le mod`ele).
2. Les estimation pβˆ(xi) des probabilit´es P(Yi = 1) par ce mod`ele sont

pβˆ(x1) = 0.64, pβˆ(x2) = 0.46, pβˆ(x3) = 0.60, pβˆ(x4) = 0.45, pβˆ(x5) = 0.86.

Calculer la log-vraisemblance maximis´ee L n (βˆ) du mod`ele.

3. On consid`ere maintenant le mod`ele satur´e permettant d’expliquer Y par X.
(a) Ecrire ce mod`ele et donner les estimations de ses param`etres.
(b) Combien vaut la log-vraisemblance maximis´ee de ce mod`ele ?
(c) En d´eduire la d´eviance du mod`ele de la question 1.
4. On cherche maintenant `a expliquer Y par X `a l’aide d’un mod`ele logistique comprenant
uniquement le coefficient associ´e `a X (sans constante). On notera γ le param`etre de ce mod`ele
et γˆ son estimateur du maximum de vraisemblance.

2
(a) Ecrire le mod`ele.

1
 (b) Exprimer σ^
2 l’estimateur de la variance de γˆ en fonction des x et de γˆ .
i

(c) Les estimation pγˆ(xi) des probabilit´es P(Yi = 1) par ce mod`ele sont

pγˆ(x1) = 0.56, pγˆ(x2) = 0.58, pγˆ(x3) = 0.57, pγˆ(x4) = 0.58, pγˆ(x5) = 0.54.

^
Calculer σ 2.

(d) On a de plus γˆ = 0.3574. D´eduire de la question pr´ec´edente une proc´edure asymptotique

de test de niveau α = 5% pour les hypoth`eses H0 : γ = 1 contre H1 : γ /= 1. Donner la
conclusion du test.
5. Tester au niveau α = 5% la nullit´e du param`etre correspondant `a la constante dans le mod`ele
de la question 1.

4
3
 05:33 mardi, février 25, 2014 2

Le Système SAS

Procédure CATMOD

Annexe 1 : exercice 2
proc catmod data=exercice2;
model Y=X1 X2 X3;
run;

Synthèse des données

Réponse Y Niveaux de réponse 3

Variable de pondération Néant Populations 8

Table EXERCICE2 Total des valeurs 400

Valeur(s) manquante(s) 0 Observations 400

Profils de population
Echantillon X1 X2 X3 Taille échantillon
1 A C E 39

2 A C F 41

3 A D E 47

4 A D F 43

5 B C E 54

6 B C F 56

7 B D E 60

8 B D F 60

Profils de
réponse
Réponse Y
1 1

2 2

3 3

Analyse du max. de vraisemblance

Convergence des calculs du maximum de vraisemblance.

6
 05:33 mardi, février 25, 2014 1

Le Système SAS

Procédure CATMOD

Analyse de variance par maximum de vraisemblance

Source DDL Khi-2 Pr > Khi-2
Intercept ? 23.94 <.0001

X1 ? 38.68 <.0001

X2 ? 57.03 <.0001

X3 ? 0.44 0.8019

Rapport de vraisemblance 8 122.39 <.0001

Analyse des valeurs estimées du maximum de vraisemblance

Nombre Erreur Khi-
Paramètre de fonctions Valeur estimée type 2 Pr > Khi-2
Intercept 1 0.6744 0.1753 14.81 0.0001

2 0.8234 0.1689 23.77 <.0001

X1 A 1 0.9519 0.1697 31.46 <.0001

A 2 0.9536 0.1606 35.26 <.0001

X2 C 1 -1.1659 0.1557 56.05 <.0001

C 2 -0.6083 0.1471 17.11 <.0001

X3 E 1 0.0925 0.1459 0.40 0.5260

E 2 0.0359 0.1388 0.07 0.7960

5
 05:36 mardi, février 25, 2014 2

Le Système SAS

Procédure CATMOD

ANNEXE 2 : exercice 2
proc catmod data=exercice2;
model Y=X1 X2 X1*X2 X3;
run;

Synthèse des données

Réponse Y Niveaux de réponse 3

Variable de pondération Néant Populations 8

Table EXERCICE2 Total des valeurs 400

Valeur(s) manquante(s) 0 Observations 400

Profils de population
Echantillon X1 X2 X3 Taille échantillon
1 A C E 39

2 A C F 41

3 A D E 47

4 A D F 43

5 B C E 54

6 B C F 56

7 B D E 60

8 B D F 60

Profils de
réponse
Réponse Y
1 1

2 2

3 3

Analyse du max. de vraisemblance

Convergence des calculs du maximum de vraisemblance.

8
 05:36 mardi, février 25, 2014 1

Le Système SAS

Procédure CATMOD

Analyse de variance par maximum de vraisemblance

Source DDL Khi-2 Pr > Khi-2
Intercept 2 5.18 0.0748

X1 2 23.06 <.0001

X2 2 34.67 <.0001

X1*X2 2 78.99 <.0001

X3 2 0.37 0.8320

Rapport de vraisemblance 6 8.90 0.1791

Analyse des valeurs estimées du maximum de vraisemblance

Nombre Erreur Khi-
Paramètre de fonctions Valeur estimée type 2 Pr > Khi-2
Intercept 1 0.3394 0.1981 2.94 0.0867

2 0.3967 0.1816 4.77 0.0290

X1 A 1 0.3539 0.1980 3.19 0.0739

A 2 0.8574 0.1816 22.30 <.0001

X2 C 1 -1.1581 0.1980 34.20 <.0001

C 2 -0.4629 0.1816 6.50 0.0108

X1*X2 AC 1 -0.4500 0.1981 5.16 0.0231

AC 2 1.0976 0.1816 36.53 <.0001

X3 E 1 0.0867 0.1449 0.36 0.5494

E 2 0.0389 0.1543 0.06 0.8012

7
 05:39 mardi, février 25, 2014 1

Le Système SAS

Procédure LOGISTIC

ANNEXE 3 : exercice 2
proc logistic data=exercice2;
class X1 X2 X3;
model Y=X1 X2 X3;
run;

Informations sur le modèle

Table WORK.EXERCICE2

Variable de réponse Y

Nombre de niveaux de réponse 3

Modèle logit cumulé

Technique d'optimisation Score de Fisher

Nombre d'observations lues 400

Nombre d'observations utili 400

Profil de réponse
Valeur Fréquence
ordonnée Y totale
1 1 147

2 2 148

3 3 105

Les probabilités modélisées sont cumulées sur les valeurs ordonnées inférieures.

Informations sur le niveau de

classe
Variables
Classe Valeur d'expérience
X1 A 1

B -1

X2 C 1
D -1

X3 E 1
F -1

9
 05:39 mardi, février 25, 2014 2

Le Système SAS

Procédure LOGISTIC

Etat de convergence du modèle

Critère de convergence (GCONV=1E-8) respecté.

Test de score pour

l'hypothèse des cotes
proportionnelles
Khi-2 DDL Pr > Khi-2
34.6463 ? <.0001

Statistiques d'ajustement du
modèle
Constante
Constante et
Critère uniquement covariables
AIC 873.478 773.330

SC 881.461 793.287

-2 Log 869.478 763.330

Test de l'hypothèse nulle globale : BETA=0

Test Khi-2 DDL Pr > Khi-2
Rapport de vrais 106.1478 ? <.0001

Score 83.4044 ? <.0001

Wald 101.3350 ? <.0001

Analyse des effets Type 3

Khi-2
Effet DDL de Wald Pr > Khi-2
X1 1 40.4077 <.0001

X2 1 75.9504 <.0001

X3 1 0.3251 0.5686

10
 05:39 mardi, février 25, 2014 3

Le Système SAS

Procédure LOGISTIC

Estimations par l'analyse du maximum de vraisemblance

Valeur Erreur Khi-2
Paramètre DDL estimée type de Wald Pr > Khi-2
Intercept 1 1 -0.6796 0.1181 33.1357 <.0001

Intercept 2 1 1.2620 0.1304 93.7324 <.0001

X1 A 1 0.6525 0.1026 40.4077 <.0001

X2 C 1 -0.9128 0.1047 75.9504 <.0001

X3 E 1 0.0558 0.0978 0.3251 0.5686

Estimations des rapports de cotes

Valeur
estimée Intervalle de confiance
Effet du point de Wald à 95 %
X1 A vs B 3.687 2.466 5.514

X2 C vs D 0.161 0.107 0.243

X3 E vs F 1.118 0.762 1.640

Association des probabilités prédites et des réponses

observées
Pourcentage concordant 68.0 D de Somers 0.454

Pourcentage discordant 22.6 Gamma 0.501

Pourcentage lié 9.4 Tau-a 0.300

Paires 52731 c 0.727

11
 Examen de R´egression sur variables cat´egorielles (Rattrapage)

ENSAI - 12 avril 2011

Sans document

Calculatrice fournie par l’´ecole, dur´ee : 1h30

Pr´eambule : Le sujet est compos´e de trois exercices ind´ependants. La qualit´e de la r´edaction sera
prise en compte. Toutes les r´eponses seront donn´ees sur la copie (ne pas rendre le sujet). Tous les
mod`eles ´ecrits devront ˆetre identifiables. Bien que cela puisse ˆetre aberrant dans certaines parties du
devoir, on supposera dans tout le devoir que les hypoth`eses sur la convergence de l’estimateur du
maximum de vraisemblance dans le mod`ele logistique sont v´erifi´ees et que le nombre d’observations est
suffisamment grand pour approcher la loi des estimateurs par la loi limite vue en cours.
Le barˆeme est donn´e `a titre indicatif.
Exercice 1 (4.5 points)
On dispose de n observations (xi, yi), i = 1, . . . , n d’une variable (qualitative) `a expliquer Y `a valeurs dans
{0, 1} par une variable explicative X `a valeurs dans R. On consid`ere le mod`ele logistique
P(Y = 1|X = x)
logit p(x) = log = βx.
1 − P(Y = 1|X = x)
1. Exprimer la log-vraisemblance de ce mod`ele en fonction des xi, yi et du param`etre β.
2. Les hypoth`eses n´ecessaires `a la convergence de l’estimateur du maximum de vraisemblance βˆ
´etant v´erifi´ees, on a
ˆ 2
(β − β) L 2
→ χ? .
σˆ 2
(a) Exprimer σˆ en fonction des xi et de βˆ.
2

(b) Quelle est le nombre de degr´es de libert´e de la loi du χ2 ?

(c) En d´eduire une proc´edure de test (asymptotique) de niveau α pour les hypoth`eses H0 : β = 1
contre H1 : β1 > 0 (on donnera la statistique de test, sa loi sous H0 ainsi que la zone de rejet que
l’on pourra exprimer en fonction des quantiles de la loi χ2).?

Exercice 2 (7 points)
On cherche `a expliquer la concentration en ozone d’une ville en fonction de la temp´erature `a 12h et de la
quantit´e de pluie. On note :
– Y la variable prenant pour valeur 1 si la concentration en ozone est forte, 0 si elle est faible.
– Pluie la variable qui prend pour valeur oui si on a enregistr´e plus de 5 millim`etres de pluie dans
la journ´ee, non sinon.
– T12 la temp´erature en degr´es Celcius relev´ee `a 12h.
On dispose des valeurs de ces 3 variables mesur´ees sur une p´eriode de n = 112 jours. On consid`ere les
mod`eles suivants :
logit P(Y = 1|X = (Pluie,T12)) = β0 + β11{Pluie=Oui} + β21{Pluie=Non} + β3T12 (1)
et
logit P(Y = 1|X = (Pluie,T12)) = α0 + α11{Pluie=Oui} + α21{Pluie=Non} + α3T12 (2)
o u` (1) est le mod`ele utilisant la contrainte β1 + β2 = 0 et (2) celui ayant pour contrainte α1 = 0.
1. L’annexe 1 contient les sorties de la proc´edure SAS
proc logistic data=exercice2 descending;
class Pluie;
model Y=Pluie T12;
run;

1
 (a) Donner les estimations βˆj des param`etres βj, j = 0, . . . , 3.
(b) Donner les estimations αˆ j des param`etres αj, j = 0, . . . , 3.
2. Quelle est la probabilit´e estim´ee par le mod`ele (1) que la concentration en ozone soit faible pour
une journ´ee o u` on observe pluie=Non et T12=22.2 ?
3. Ecrire les hypoth`eses des tests r´ealis´es dans le tableau Test de l’hypoth`ese nulle globale :
BETA=O. Quels sont les nombres de degr´es de libert´e manquants ? Quelle(s) est (sont) la (les)
conclusion(s) de ces tests ?
4. Ecrire les hypoth`eses du test r´e a l i s´e sur la ligne Pluie du tableau Analyse des effets de Type
3. Quel est le nombre de degr´es de libert´e manquant ? Quelles sont les conclusions des tests de ce
tableau.
5. On souhaite effectuer un test du rapport de vraisemblance permettant de comparer le model (1) au
mod`ele logistique ayant pour variable explicative uniquement la variable T12 (ainsi que la constante).
(a) Ecrire les hypoth`eses de ce test.
(b) Donner la statistique de test ainsi que sa loi sous H0.
(c) La log-vraisemblance du mod`ele logistique ayant pour variable explicative uniquement la
variable T12 (ainsi que la constante) est ´egale `a -61.240. Calculer la valeur observ´ee de la
statistique de test.
(d) Ecrire la probabilit´e critique de ce test en fonction de la fonction de r´epartition de la loi du χ2
`a p degr´es de libert´e (on notera Fp cette fonction de r´epartition et on pr´ecisera la valeur
de p). Conclure.

Exercice 3 (8.5 points)

Partie 1 : question de cours

On souhaite ici expliquer Y une variable qualitative ordinale `a 3 modalit´es 1 < 2 < 3 par une variable
explicative X `a valeurs dans R.
1. Ecrire les diff´erentes ´etapes qui permettent de d´efinir le mod`ele polytomique ordonn´e : on
partira de la mod´elisation par une variable latente pour arriver `a l’´ecriture du mod`ele. On pr
´ecisera notamment les probabilit´es P(Y = j|X = x), j = 1, 2, 3 en fonction des coefficients du
mod`ele.
2. Expliquer pourquoi cette mod´elisation engendre une hypoth`ese de proportionnalit´e des odds ratio.

Partie 2 : application

On souhaite `a nouveau expliquer la concentration en ozone d’une ville. On consid`ere ici comme variables
explicatives la direction du vent et la quantit´e de pluie. On note :
– Y la variable prenant pour valeur 1 si la concentration en ozone est faible, 2 si elle est moyenne, 3 si
elle est forte.
– Vent la variable `a 4 modalit´es correspondant `a la direction du vent (Est, Nord, Ouest, Sud).
– Pluie la variable qui prend pour valeur oui si on a enregistr´e plus de 5 millim`etres de pluie dans
la journ´ee, non sinon.
On dispose des valeurs de ces 3 variables mesur´ees sur une p´eriode de n = 112 jours. On effectue sous SAS
la proc´edure
proc logistic data=exercice3;
class Vent Pluie;
model Y=Vent Pluie;
run;
1. Ecrire le mod`ele ajust´e par cette proc´edure. Vous pr´eciserez les ´eventuelles contraintes
d’identifiabilit´e des param`etres.
2. Les sorties de la proc´edure se trouvent en annexe 2. Donner les estimations des param`etres du
mod`ele que vous avez ´ecrit `a la question pr´ec´edente.
3. Quelle est la probabilit´e estim´ee par ce mod`ele que la concentration en ozone soit faible pour
2
une journ´ee o u` on observe pluie=Oui et Vent=Est ?

3
4. Expliquer en quelques mots le test effectu´e dans le tableau Score Test for the Proportional
Odds Assumption. Donner le nombre de degr´es de libert´e manquant dans ce tableau (justifier)
ainsi que la conclusion de ce test.
5. On consid`ere maintenant le mod`ele polytomique nominal permettant d’expliquer la variable Y par
la variable Vent uniquement. Le tableau 1 contient les valeurs observ´ees de ces deux variables.

1 2 3 Total
Est 1 1 8 10
Nord 10 13 8 31
Ouest 22 17 11 50
Sud 4 6 11 21
Table 1 – Mesures des variables Y et Vent.

(a) Ecrire ce mod`ele et donner sa dimension.

(b) Ce mod`ele est-il satur´e ? Justifier.
(c) Calculer la valeur de l’estimateur du maximum de vraisemblance de la probabilit´e πNord(j) =
P(Y = j|X = Nord).
(d) Calculer un intervalle de confiance (asymptotique) de niveau 95% de πNord(1) (on pourra ap-
procher le quantile d’ordre 0.975 de la loi normale centr´ee r´eduite par 1.96).

2
5
ANNEXE 1 : exercice 2

proc logistic data=exercice2 descending;

class Pluie;
model Y=Pluie T12;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE2

Response Variable Y
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 112

Number of Observations Used 112

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 56
2 0 56

Probability modeled is Y='1'.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

Pluie Non 1
Oui -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 157.265 119.597

SC 159.983 127.753
-2 Log L 155.265 113.597

4
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 41.6677 ? <.0001

Score 35.8430 ? <.0001
Wald 25.9206 ? <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

Pluie ? 8.4176 0.0037

T12 ? 12.5706 0.0004

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -6.2380 1.6816 13.7611 0.0002

Pluie Non 1 0.7161 0.2468 8.4176 0.0037
T12 1 0.2848 0.0803 12.5706 0.0004

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

Pluie Non vs Oui 4.188 1.592 11.021

T12 1.330 1.136 1.556

Association des probabilités prédites et des réponses observées

Percent Concordant 82.0 Somers' D 0.643

Percent Discordant 17.6 Gamma 0.646
Percent Tied 0.4 Tau-a 0.325
Pairs 3136 c 0.822

7
ANNEXE 2 : exercice 3

proc logistic data=exercice3;

class Vent Pluie;
model Y=Vent Pluie;
run;

The LOGISTIC Procedure

Informations sur le

modèle

Data Set WORK.EXERCICE3

Response Variable Y
Number of Response Levels 3
Model cumulative logit
Optimization Technique Fisher's scoring

Number of Observations Read 112

Number of Observations Used 112

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 37
2 2 37
3 3 38

Probabilities modeled are cumulated over the lower Ordered Values.

Informations sur le niveau de classe

Variables de
Classe Valeur création

Vent Est 1 0 0
Nord 0 1 0
Ouest 0 0 1
Sud -1 -1 -1

Pluie Non 1
Oui -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Score Test for the Proportional Odds Assumption

Khi 2 DF Pr > Khi 2

2.3730 ? 0.6675

6
 Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 250.071 221.320

SC 255.508 237.631
-2 Log L 246.071 209.320

Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 36.7512 4 <.0001

Score 31.4803 4 <.0001
Wald 30.0695 4 <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

Vent 3 10.3163 0.0161

Pluie 1 18.9092 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 1 -1.1417 0.2949 14.9904 0.0001

Intercept 2 1 0.6468 0.2846 5.1666 0.0230
Vent Est 1 -1.4080 0.6103 5.3232 0.0210
Vent Nord 1 0.8324 0.3516 5.6036 0.0179
Vent Ouest 1 0.7883 0.3231 5.9532 0.0147
Pluie Non 1 -0.9034 0.2077 18.9092 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

Vent Est vs Sud 0.303 0.052 1.746

Vent Nord vs Sud 2.844 0.959 8.433
Vent Ouest vs Sud 2.721 0.980 7.556
Pluie Non vs Oui 0.164 0.073 0.371

Association des probabilités prédites et des réponses observées

Percent Concordant 66.6 Somers' D 0.474

Percent Discordant 19.2 Gamma 0.553
Percent Tied 14.2 Tau-a 0.319
Pairs 4181 c 0.737

1
 Examen de R´egression sur variables cat

´egorielles ENSAI - 19 janvier 2011

Document autoris´e : une feuille A4

manuscrite Calculatrice fournie par l’´ecole,

dur´ee : 1h45

Pr´eambule : Le sujet est compos´e de quatre exercices ind´ependants. La qualit´e de la r´edaction sera
prise en compte. Toutes les r´eponses seront donn´ees sur la copie (ne pas rendre le sujet). Tous les
mod`eles ´ecrits devront ˆetre identifiables. Bien que cela puisse ˆetre aberrant dans certaines parties du
devoir, on supposera dans tout le devoir que les hypoth`eses sur la convergence de l’estimateur du
maximum de vraisemblance dans le mod`ele logistique sont v´erifi´ees et que le nombre d’observations est
suffisamment grand pour approcher la loi des estimateurs par la loi limite vue en cours.
Le barˆeme est donn´e `a titre indicatif.
Exercice 1 (5 points)
On dispose de 5 observations (xi, yi), i = 1, . . . , 5 d’une variable (qualitative) `a expliquer Y `a valeurs dans
{0, 1} par une variable explicative X `a valeurs dans R. Les donn´ees se trouvent dans le tableau suivant :

xi 3 8 14 11 5
yi 1 0 0 1 0

1. Ecrire le mod`ele logistique permettant d’expliquer Y par X (on pr´ecisera notamment la loi
condition- nelle de Y |X = x et la mod´elisation du param`etre de cette loi conditionnelle).
2. Exprimer la log-vraisemblance de ce mod`ele en fonction des xi, yi et des param`etres du mod`ele.
3. L’annexe 1 contient certaines sorties de la proc´edure logistique effectu´ee sous SAS. Calculer la
log- vraisemblance du mod`ele. En d´eduire la d´eviance.
4. On effectue un test du rapport de vraisemblance (ou de d´eviance) permettant de comparer le
mod`ele construit au mod`ele compos´e uniquement de la constante.
(a) Ecrire les hypoth`eses (math´ematiques) de ce test.
(b) Donner la statistique de test ainsi que sa loi sous H0. Calculer la valeur observ´ee de cette statis-
tique de test.
(c) Ecrire la probabilit´e critique de ce test en fonction de la fonction de r´epartition de la loi du χ2 `a
p degr´es de libert´e (on notera Fp cette fonction de r´epartition et on pr´ecisera la valeur de p).
5. Calculer la valeur de l’odds ratio manquant dans le tableau Estimation des rapports de cotes de
l’annexe 1.
Exercice 2 (8 points)
On fait d´eguster un vin `a n = 430 personnes `a qui on demande de mettre une note A, B ou C (A si
on a aim´e, B si on a moyennement appr´eci´e, C si on n’a pas aim´e). Le but de l’´etude est de voir si
l’appr´eciation de ce vin est diff´erente selon que l’on a mang´e du fromage avant et/ou fum´e une
cigarette. On d´esigne ainsi par :
– Y la variable prenant pour valeur A, B ou C selon que l’appr´eciation du vin.
– X1 la variable qui prend pour valeur F si on a mang´e du fromage avant de goutter le vin, NF sinon.
– X2 la variable qui prend pour valeur C si on a fu m´e une cigarette, NC
sinon. Les r´esultats de l’´etude se trouvent dans le tableau 1.

8
 A B C Total
FC 5 80 15 100
NF C 23 12 85 120
F NC 59 20 11 90
NF NC 48 53 19 120

TABLE 1 – R´esultats de la d´egustation.

On peut par exemple lire dans ce tableau que 100 personnes ont mang´e du fromage et fu m´e une cigarette
avant de d´eguster le vin. Parmi ces 100 personnes, 5 ont mis la note A, 80 la note B et 15 la note C.
Dans cet exercice, on cherche `a expliquer la variable Y par les variables X1 et X2.
1. On consid`ere tout d’abord le mod`ele satur´e.
(a) Quels sont les param`etres de ce mod`ele ? Donner sa dimension.
(b) Pour deux param`etres au choix, calculer les estimateurs du maximum de vraisemblance.
(c) On consid`ere uniquement les personnes qui ont fum´e et mang´e du fromage avant la d
´egustation (groupe (X1, X2) = (F, C)). Donner la loi asymptotique de l’estimateurs du
maximum de vraisem- blance du vecteur
 
P(Y = A|X1 = F, X2 = C)
πF C = P(Y = B|X1 = F, X2 = C)
P(Y = C|X1 = F, X2 = C)
(d) D´eduire de la question pr´ec´edente un intervalle de confiance (asymptotique) de niveau 95% pour
P(Y = A|X1 = F, X2 = C) (on pourra approcher le quantile d’ordre 0.975 de la loi normale centr
´ee r´eduite par 1.96).
2. On consid`ere maintenant le mod`ele polytomique nominal ayant pour variables explicatives X1 et X2
(sans interaction). Les sorties SAS se trouvent en annexe 2.
(a) Ecrire ce mod`ele et donner sa dimension. Vous pr´eciserez les ´eventuelles contraintes
d’identifia- bilit´e des param`etres.
(b) Donner les estimations des param`etres du mod`ele que vous avez ´ecrit `a la question pr´ec´edente.
(c) Ecrire les hypoth`eses du test de l’effet de X1 dans le tableau "Analyse des effets de
type 3". Quels sont les nombres de degr´es de libert´e manquants dans ce tableau ?
(d) Calculer la probabilit´e estim´ee par ce mod`ele qu’un individu mette la note B en ayant fum´e
une cigarette mais sans avoir mang´e de fromage avant la d´egustation.
3. On consid`ere maintenant le mod`ele polytomique nominal ayant pour variables explicatives X1 et X2
ainsi que l’interaction entre X1 et X2. Les sorties SAS se trouvent en annexe 3.
(a) Donner la dimension de ce mod`ele.
(b) Rappeler dans le contexte du probl`eme ce que signifie l’interaction entre les variables X1 et X2.
(c) Ce mod`ele est-il satur´e ? Justifier.
(d) Proposer un test permettant de tester le mod`ele satur´e contre le mod`ele de la question 2 (le
mod`ele n’ayant pas ´e t ´e ´ecris, on se contentera de d´ecrire les hypoth`eses du test et de
nommer la statistique de test (Wald, rapport de vraisemblance ou score), on donnera la loi de
cette statistique sous l’hypoth`ese nulle, la valeur observ´ee de cette statistique ainsi que la
probabilit´e critique du test... et tant qu’ `a faire la conclusion du test).

Exercice 3 (6 points)
On souhaite expliquer expliquer un probl`eme d’insuffisance coronarienne (variable chd qui prend pour
valeur 1 pour insuffisance 0 sinon) par 5 variables :
– tobaco : variable quantitative qui repr´esente la quantit´e moyenne de cigarettes fum´ees par jour ;
– ldl : variable quantitative qui repr´esente un indicateur l i´e au taux de cholest´erol ;
– famhist : variable qualitative qui prend pour valeurs present si l’insuffisance coronarienne a ´e t ´e observ
´ee dans la famille de l’individu, absent sinon ;
– alcohol : variable quantitative qui repr´esente un indicateur l i´e `a la consommation d’alcool de
l’individu sur un mois ;

3
– age : variable quantitative qui repr´esente l’age du patient.
On souhaite expliquer la variable chd par les 5 variables ci dessus `a l’aide d’un mod`ele logistique. On dispose
de n = 462 individus. Les sorties de la proc´edure

proc logistic data=exercice3 descending;

class famhist;
model chd=tobacco ldl famhist alcohol age;
run;

se trouvent en annexe 4.
1. Ecrire le mod`ele a just´e par cette proc´edure et donner les estimations des param`etres.
2. Ecrire les hypoth`eses des tests r´ealis´es dans le tableau Test de l’hypoth`ese nulle globale :
BETA=O. Quels sont les nombres de degr´es de libert´e manquants ? Quelle(s) est (sont) la (les)
conclusion(s) de ces tests ?
3. Ecrire les hypoth`eses du test r´e a l i s´e sur la ligne famhist du tableau Analyse des effets de Type
3. Quel est le nombre de degr´es de libert´e manquant ? Quelles sont les conclusions des tests de ce
tableau.
4. Quelle serait la premi`ere variable supprim´ee `a ce mod`ele par un algorithme pas `a pas
descendant (obtenu en ajoutant uniquement l’argument selection=backward dans la proc
logistic).
On consid`ere maintenant la mod`ele logistique construit `a partir des 4 variables tobaco, ldl, famhist et
age. La proc´edure SAS ainsi que les sorties se trouvent en annexe 5.
5. On souhaite r´ealiser un test du rapport de vraisemblance permettant de comparer ce mod`ele au
mod`ele de la question 1.
(a) Ecrire les hypoth`eses de ce test.
(b) Donner la statistique de test ainsi que sa loi sous H0 et calculer la valeur observ´ee de
cette statistique.
(c) Ecrire la probabilit´e critique de ce test en fonction de la fonction de r´epartition de la loi du χ2
`a p degr´es de libert´e (on notera Fp cette fonction de r´epartition et on pr´ecisera la valeur de p).
Conclure.
6. A partir des sorties propos´ees en annexe 4 et 5, comparer le mod`ele que l’on vient de construire
(sans la variable alcohol) au mod`ele de la question 1 (vous utiliserez tous les crit`eres possibles
permettant de comparer ces deux mod`eles).

Exercice 4 (3 points)
On reprend les donn´ees de la question 1. Les estimations pˆ(x i ) des probabilit´es P(Y = 1|X = xi) sont
donn´ees par

pˆ(x 1 ) = pˆ(x 2 ) = pˆ(x 3 ) = pˆ(x 4 ) = pˆ(x 5 ) = 0.498.

0.564, 0.399, 0.231, 0.309,

des probabilit´es. Calculer un intervalle de confiance de niveau 95% pour les param`etres du mod`ele
logistique (on pourra approcher le quantile d’ordre 0.975 de la loi N (0.1) par 1.96).

2
5
proc logistic data=exercice1 descending;
model Y=X;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE1

Response Variable Y
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 5

Number of Observations Used 5

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 2
2 0 3

Probability modeled is Y='1'.

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 8.730 10.418

SC 8.340 9.637
-2 Log L 6.730 6.418

Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio ?????? ? ??????

Score 0.3046 ? 0.5810
Wald 0.2931 ? 0.5882

4
 Analyse des estimations de la vraisemblance maximum

ANNEXE 1: exercice Erreur Khi 2

1 Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 0.6538 2.1219 0.0949 0.7580

X 1 -0.1327 0.2451 0.2931 0.5882

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

X ????? 0.542 1.416

Association des probabilités prédites et des réponses observées

Percent Concordant 66.7 Somers' D 0.333

Percent Discordant 33.3 Gamma 0.333
Percent Tied 0.0 Tau-a 0.200
Pairs 6 c 0.667

7
 Test de l'hypothèse nulle globale : BETA=0

proc logistic data=exercice2 descending;

class X1 X2 ;
model Y=X1 X2/ link=glogit;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE2

Response Variable Y
Number of Response Levels 3
Model generalized logit
Optimization Technique Fisher's scoring

Number of Observations Read 430

Number of Observations Used 430

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 C 130
2 B 165
3 A 135

Logits modeled use Y='A' as the reference category.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

X1 F 1
NF -1

X2 C 1
NC -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 943.910 802.598

SC 952.038 826.980
-2 Log L 939.910 790.598

6
ANNEXE 2 : exercice
2
Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 149.3125 4 <.0001

Score 137.1818 4 <.0001
Wald 116.9070 4 <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

X1 ? 49.3245 <.0001
X2 ? 77.8695 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre Y DF Estimation std de Wald Pr > Khi 2

Intercept C 1 -0.2730 0.1661 2.7037 0.1001

Intercept B 1 0.3842 0.1329 8.3550 0.0038
X1 F C 1 -0.8169 0.1585 26.5560 <.0001
X1 F B 1 0.1679 0.1244 1.8208 0.1772
X2 C C 1 1.3738 0.1572 76.3797 <.0001
X2 C B 1 0.7617 0.1331 32.7509 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Y Estimate de Wald

X1 F vs NF C 0.195 0.105 0.363

X1 F vs NF B 1.399 0.859 2.279
X2 C vs NC C 15.607 8.427 28.902
X2 C vs NC B 4.588 2.723 7.730

9
 Test de l'hypothèse nulle globale : BETA=0

proc logistic data=exercice2 descending;

class X1 X2;
model Y=X1 X2 X1*X2/ link=glogit;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE2

Response Variable Y
Number of Response Levels 3
Model generalized logit
Optimization Technique Fisher's scoring

Number of Observations Read 430

Number of Observations Used 430

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 C 130
2 B 165
3 A 135

Logits modeled use Y='A' as the reference category.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

X1 F 1
NF -1

X2 C 1
NC -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 943.910 729.307

SC 952.038 761.817
-2 Log L 939.910 713.307

8
ANNEXE 3 : exercice
2 Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 226.6030 6 <.0001

Score 232.6809 6 <.0001
Wald 176.4404 6 <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

X1 2 28.5933 <.0001
X2 2 50.2397 <.0001
X1*X2 2 64.6155 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre Y DF Estimation std de Wald Pr > Khi 2

Intercept C 1 -0.0502 0.1773 0.0800 0.7773

Intercept B 1 0.2848 0.1669 2.9103 0.0880
X1 F C 1 -0.2404 0.1773 1.8370 0.1753
X1 F B 1 0.5606 0.1669 11.2739 0.0008
X2 C C 1 1.2530 0.1773 49.9209 <.0001
X2 C B 1 0.7762 0.1669 21.6143 <.0001
X1*X2 F C C 1 0.1361 0.1773 0.5887 0.4429
X1*X2 F C B 1 1.1510 0.1669 47.5320 <.0001

11
 Test de l'hypothèse nulle globale : BETA=0

proc logistic data=exercice3 descending;

class famhist;
model chd=tobacco ldl famhist alcohol age;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE3

Response Variable chd
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 462

Number of Observations Used 462

Profil de réponse

Valeur Fréquence
ordonnée chd totale

1 1 160
2 0 302

Probability modeled is chd='1'.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

famhist Absent 1
Present -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 598.108 497.379

SC 602.244 522.192
-2 Log L 596.108 485.379

10
ANNEXE 4 : exercice
3 Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 110.7298 ? <.0001

Score 99.5958 ? <.0001
Wald 80.2518 ? <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

tobacco ? 9.3601 0.0022

ldl ? 9.6321 0.0019
famhist ? 16.9328 <.0001
alcohol ? 0.0654 0.7981
age ? 20.4477 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -3.7699 0.5073 55.2213 <.0001

tobacco 1 0.0794 0.0260 9.3601 0.0022
ldl 1 0.1687 0.0543 9.6321 0.0019
famhist Absent 1 -0.4602 0.1118 16.9328 <.0001
alcohol 1 0.00112 0.00439 0.0654 0.7981
age 1 0.0442 0.00977 20.4477 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

tobacco 1.083 1.029 1.139

ldl 1.184 1.064 1.317
famhist Absent vs Present 0.398 0.257 0.618
alcohol 1.001 0.993 1.010
age 1.045 1.025 1.065

Association des probabilités prédites et des réponses observées

Percent Concordant 78.0 Somers' D 0.562

Percent Discordant 21.8 Gamma 0.563
Percent Tied 0.2 Tau-a 0.255
Pairs 48320 c 0.781

13
 Test de l'hypothèse nulle globale : BETA=0

proc logistic data=exercice3 descending;

class famhist;
model chd=tobacco ldl famhist age;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE3

Response Variable chd
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 462

Number of Observations Used 462

Profil de réponse

Valeur Fréquence
ordonnée chd totale

1 1 160
2 0 302

Probability modeled is chd='1'.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

famhist Absent 1
Present -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 598.108 495.444

SC 602.244 516.122
-2 Log L 596.108 485.444

12
ANNEXE 5 : exercice
3 Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 110.6646 4 <.0001

Score 99.5934 4 <.0001
Wald 80.3225 4 <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

tobacco 1 10.0039 0.0016

ldl 1 9.5638 0.0020
famhist 1 17.1448 <.0001
age 1 20.4333 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -3.7422 0.4945 57.2732 <.0001

tobacco 1 0.0807 0.0255 10.0039 0.0016
ldl 1 0.1676 0.0542 9.5638 0.0020
famhist Absent 1 -0.4621 0.1116 17.1448 <.0001
age 1 0.0440 0.00974 20.4333 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

tobacco 1.084 1.031 1.140

ldl 1.182 1.063 1.315
famhist Absent vs Present 0.397 0.256 0.615
age 1.045 1.025 1.065

Association des probabilités prédites et des réponses observées

Percent Concordant 78.0 Somers' D 0.563

Percent Discordant 21.7 Gamma 0.564
Percent Tied 0.2 Tau-a 0.255
Pairs 48320 c 0.781

1
 Test de l'hypothèse nulle globale : BETA=0

Examen de R´egression sur variables cat

´egorielles ENSAI - 19 janvier 2010

Document autoris´e : une feuille A4

manuscrite Calculatrice fournie par l’´ecole,

dur´ee : 1h45

Pr´eambule : Le sujet est compos´e de quatre exercices ind´ependants. La qualit´e de la r´edaction

sera prise en compte. Toutes les r´eponses seront donn´ees sur la copie (ne pas rendre le sujet). Tous les
mod`eles ´ecrits devront ˆetre identifiables.

Exercice 1 (Question de cours)

On dispose de n + m observations i.i.d. (Xi, Yi), i = 1, . . . , n + m o u` Xi est une variable al´eatoire `a
valeurs dans Rd et Yi est une variable al´eatoire binaire `a valeurs dans {0, 1}. On construit un mod`ele
logistique `a l’aide des n premi`eres observations seulement et on d´eduit de ce mod`ele la r`egle de
discrimination

1 si pˆ(x) ≥ 0.5
gˆ(x) = 0 sinon.
o u` pˆ(x) d´esigne la probabilit´e P(Y = 1|X = x) estim´ee par le mod`ele. On note Ðn l’´echantillon
compos´e des n premi`eres observations (Ðn = (X1, Y1), . . . , (Xn, Yn)). On souhaite estimer la
probabilit´e d’erreur de gˆ d´efinie par
L(gˆ) = P(gˆ(X) /= Y |Ðn).
1. Proposer un estimateur L m (gˆ) sans biais de L(gˆ) (on pourra construire cet estimateur `a partir des m
derni`eres observations seulement).
2. Montrer que, conditionnellement `a Ðn, cet estimateur est sans biais (c’est-`a - dire E(Lm (gˆ)|Ðn ) = L(gˆ)).
3. Exprimer la variance conditionnelle V(Lm (gˆ)|Ðn ) de cet estimateur en fonction de m et de L(gˆ).
Exercice 2
On dispose de 4 observations (xi, yi), i = 1, . . . , 4 d’une variable `a expliquer Y `a valeurs dans {0, 1} par une
variable explicative X `a valeurs dans R. Les donn´ees se trouvent dans le tableau suivant :

xi -1.1 0.8 1.4 3.2

yi 1 1 0 1

1. Ecrire le mod`ele logistique permettant d’expliquer Y par X (on pr´ecisera notamment la loi
condition- nelle de Y |X = x et la mod´elisation du param`etre de cette loi conditionnelle).
2. On effectue la r´egression logistique permettant d’expliquer Y par X sous SAS. Les probabilit´es pˆ(x i ) =
Pˆ ( Y = 1|X = xi) estim´ees par le mod`ele sont
pˆ(x 1 ) = 0.82, pˆ(x 2 ) = pˆ(x 3 ) = pˆ(x 4 ) = 0.67.
0.76, 0.74,
(a) Calculer la log-vraisemblance du mod`ele construit.
(b) Quelle est la vraisemblance du mod`ele satur´e ?
(c) En d´eduire la d´eviance.
Exercice 3
Une banque cherche `a mod´eliser la qualit´e de ses clients `a partir de son historique. La variable `a expliquer
Y est binaire, elle est d´efinie par

1 si le client a eu 0 impay´e au cours de ses remboursements mensuels

Y =
14
 0 sinon (si le client a connu au moins un impay´e).
Pour simplifier, on suppose que l’on dispose de deux variables explicatives uniquement : l’ˆage et le sexe des
clients. On dispose d’un ´echantillon d’apprentissage de taille 100 (fichier exercice3) dans lequel

1
– la premi`ere colonne Y contient les observations de la variable `a expliquer ;
– la deuxi`eme colonne X1 contient les observations de la variable binaire sexe : H “homme”, F ”femme” ;
– la troisi`eme colonne X2 contient les observations de la variable continue age.
On consid`ere les mod`eles suivants :

logit P(Y = 1|X = x) = β0 + β11{x1=F} + β21{x1=H} + β3x2 (1)

et
logit P(Y = 1|X = x) = α0 + α11{x1=F} + α21{x1=H} + α3x2 (2)
o u` x = (x1, x2), (1) est le mod`ele utilisant la contrainte β1 +β2 = 0 et (2) celui ayant pour contrainte α1 = 0.
1. L’annexe 1 contient les sorties de la proc´edure SAS
proc logistic data=exercice3 descending;
class X1;
model Y=X1 X2;
run;
(a) Donner les estimations βˆj des param`etres βj, j = 0, . . . , 3.
(b) Donner les estimations αˆ j des param`etres αj, j = 0, . . . , 3.
2. Donner la probabilit´e critique du test de Wald pour les hypoth`eses H0 : α1 = α2 = 0 contre H1 :
α1 /= 0 ou α2 /= 0.
3. Quelle est la probabilit´e estim´ee par le mod`ele (2) d’avoir au moins un impay´e pour un homme
de 40 ans ?
4. Comparer le mod`ele (2) au mod`ele logistique compos´e uniquement de la constante en utilisant
tous les crit`eres possibles (voir annexe 1).
5. Quelle est la loi de la statistique de test du test de Wald pour l’hypoth`ese H0 : β1 = 0 contre
H1 : β1 /= 0.
6. On souhaite maintenant comparer le mod`ele (1) au mod`ele qui contient uniquement la variable sexe.
On lance donc la proc´edure
proc logistic data=exercice3 descending;
class X1;
model Y=X1;
run;
Parmi les sorties on obtient :
Statistiques d’ajustement du mod`ele
Coordonn´ee Coordonn´ee `a

l’origine Crit`ere `a l’origine uniquement et covariables

AIC 138.663 140.622

SC 141.268 145.833
-2 Log L 136.663 136.622

On effectue un test de d´eviance (ou test du rapport de vraisemblance) permettant de comparer le

mod`ele (1) au mod`ele qui contient uniquement la variable sexe.
(a) Ecrire les hypoth`eses de ce test.
(b) Donner la statistique de test ainsi que sa loi sous H0. Calculer la valeur observ´ee de cette statis-
tique de test.
(c) Ecrire la probabilit´e critique de ce test en fonction de la fonction de r´epartition de la loi du χ2 `a
p degr´es de libert´e (on notera Fp cette fonction de r´epartition, pr´eciser la valeur de p).

2
 7. On compl`ete le mod`ele (1) en ajoutant l’interaction entre les variables X1 et X2. L’annexe 2
contient les sorties de la proc´edure SAS
proc logistic data=exercice3 descending;
class X1;
model Y=X1 X2 X1*X2;
run;
En vous inspirant de (1) et (2), ´ecrire ce nouveau mod`ele. Vous pr´eciserez les contraintes d’identifiabilit
´e ainsi que les valeurs estim´ees des param`etres.
8. Quelle est la probabilit´e estim´ee (par ce nouveau mod`ele) d’avoir au moins un impay´e pour un
homme de 40 ans ?
9. Comparer les mod`eles avec et sans interaction (mod`ele (1)) en utilisant au moins 3 crit`eres.
10. On consid`ere l’´ecriture suivante du mod`ele avec interaction
0 1
logit P(Y = 1|X = x)
= γ0f + γ1f x2 si x1 = F.

γh + γhx2 si x1 = H
En utilisant les sorties de l’annexe 2, calculer les estimations γˆ h , γˆ h , γˆ f et γˆ f des param`etres γh, γh,
f 0 1 0 1 0 1
et γf .
γ0 1
Exercice 4
Nous consid´erons le mˆe me probl`eme que dans l’exercice pr´ec´edent. La banque dispose d’une autre fichier de
donn´ees o u` la variable Y `a expliquer est d´efinie par :
1 si le client a eu 0 impay
Y =  2´e si le client a eu entre 1 et 3 impay´es (3 est inclu)

3 si le client a eu strictement plus de 3 impay´es.
Nous cherchons ici `a expliquer la variable Y par la variable sexe uniquement. On dispose d’un ´echantillon
de taille n = 200. On effectue sous SAS la proc´edure
proc catmod data=exercice4;
class sexe;
model Y=sexe;
run;
1. Ecrire le mod`ele ajust´e par cette proc´edure. Vous pr´eciserez les contraintes d’identifiabilit´e
des pa- ram`etres.
2. Les sorties de la proc´edure se trouvent en annexe 3. Donner les estimations des param`etres du
mod`ele que vous avez ´ecrit `a la question pr´ec´edente.
3. Quelle est la probabilit´e estim´ee (par ce mod`ele) d’avoir au moins un impay´e pour un homme ?
4. D ’ o u` viennent les 0 degr´es de libert´e sur la ligne Likelihood Ratio du tableau Maximum
Likelihood Analysis of Variance de l’annexe 3 ? Conclure.
On effectue maintenant la proc´edure
proc logistic data=exercice4;
class sexe;
model Y=sexe;
run;
5. Ecrire le mod`ele a just´e par cette proc´edure. Vous pr´eciserez les contraintes d’identifiabilit´e
des pa- ram`etres.
6. Les sorties de la proc´edure se trouvent en annexe 4. Donner les estimations des param`etres du
mod`ele que vous avez ´ecrit `a la question pr´ec´edente.
7. Quel test est r´e a l i s´e dans la partie Score test for the Proportional Odds Assumption. Vous pr
´eciserez l’hypoth`ese nulle, donnerez le nombre de degr´es de libert´e de la statistique de test (il
est inutile de donner la statistique de test) ainsi que la conclusion du test.
8. Calculer la valeur de l’odds ratio manquant dans le tableau Estimations des rapports de cotes
de l’annexe 4.
3
4
ANNEXE 1 : exercice
3
proc logistic data=exercice3 descending;
class X1;
model Y=X1 X2;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE3

Response Variable Y
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 100

Number of Observations Used 100

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 57
2 0 43

Probability modeled is Y='1'.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

X1 F 1
H -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 138.663 142.543

SC 141.268 150.358
-2 Log L 136.663 136.543

5
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 0.1200 2 0.9418

Score 0.1199 2 0.9418
Wald 0.1198 2 0.9419

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

X1 1 0.0400 0.8415
X2 1 0.0791 0.7785

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 0.0341 0.9031 0.0014 0.9699

X1 F ? 0.0404 0.2021 0.0400 0.8415
X2 1 0.00508 0.0181 0.0791 0.7785

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

X1 F vs H 1.084 0.491 2.394

X2 1.005 0.970 1.041

Association des probabilités prédites et des réponses observées

Percent Concordant 48.4 Somers' D 0.000

Percent Discordant 48.4 Gamma 0.000
Percent Tied 3.2 Tau-a 0.000
Pairs 2451 c 0.500

6
ANNEXE 2 : exercice
3
proc logistic data=exercice3 descending;
class X1;
model Y=X1 X2 X1*X2;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE3

Response Variable Y
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 100

Number of Observations Used 100

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 57
2 0 43

Probability modeled is Y='1'.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

X1 F 1
H -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 138.663 125.472

SC 141.268 135.893
-2 Log L 136.663 117.472

7
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 19.1906 3 0.0002

Score 18.0623 3 0.0004
Wald 15.8971 3 0.0012

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

X1 1 15.0068 0.0001
X2 1 0.0419 0.8377
X2*X1 1 15.8327 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 0.5560 1.0830 0.2636 0.6076

X1 F 1 -4.1954 1.0830 15.0068 0.0001
X2 1 -0.00449 0.0219 0.0419 0.8377
X2*X1 F 1 0.0872 0.0219 15.8327 <.0001

Association des probabilités prédites et des réponses observées

Percent Concordant 73.6 Somers' D 0.482

Percent Discordant 25.4 Gamma 0.487
Percent Tied 1.0 Tau-a 0.239
Pairs 2451 c 0.741

8
ANNEXE 3 : exercice
4
proc catmod data=exercice4;
model Y=sexe;
run;
The CATMOD Procedure

Récapitulatif sur les données

Response Y Response Levels 3

Weight Variable None Populations 2
Data Set EXERCICE4 Total Frequency 200
Frequency Missing 0 Observations 200

Profils de population

Échantillon sexe Sample Size

ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
1 F 111
2 H 89

Profils de réponse

Réponse Y
ƒƒƒƒƒƒƒƒƒƒƒƒ
1 1
2 2
3 3

Analyse de vraisemblance max.

Maximum likelihood computations converged.

Maximum Likelihood Analysis of Variance

Source DF Khi 2 Pr > Khi 2

ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Intercept 2 11.94 0.0026
sexe 2 49.46 <.0001

Likelihood Ratio 0 . .

Analysis of Maximum Likelihood Estimates

Nombre Erreur
Parameter de fonctions Estimation std
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Intercept 1 0.6096 0.4057
2 1.1179 0.3784
sexe F 1 -2.7046 0.4057
F 2 -1.6547 0.3784

Analysis of Maximum Likelihood Estimates

Nombre Khi-
Parameter de fonctions 2 Pr > Khi 2
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Intercept 1 2.26 0.1330
2 8.73 0.0031
sexe F 1 44.43 <.0001
F 2 19.12 <.0001

9
10
ANNEXE 4 : exercice
4
proc logistic data=exercice4;
class sexe;
model Y=sexe;
run;

The LOGISTIC Procedure

Informations sur le modèle

Data Set WORK.EXERCICE4

Response Variable Y
Number of Response Levels 3
Model cumulative logit
Optimization Technique Fisher's scoring

Number of Observations Read 200

Number of Observations Used 200

Profil de réponse

Valeur Fréquence
ordonnée Y totale

1 1 63
2 2 70
3 3 67

Probabilities modeled are cumulated over the lower Ordered Values.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

sexe F 1
H -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Score Test for the Proportional Odds Assumption

Khi 2 DF Pr > Khi 2

2.0695 ? 0.1503

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 443.074 334.977

SC 449.670 344.872
-2 Log L 439.074 328.977

11
 The LOGISTIC Procedure

Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 110.0972 1 <.0001

Score 93.5282 1 <.0001
Wald 73.8009 1 <.0001

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

sexe 1 73.8009 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 1 -1.1689 0.2113 30.6088 <.0001

Intercept 2 1 1.3737 0.2158 40.5178 <.0001
sexe F 1 -1.6869 0.1964 73.8009 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

sexe F vs H ??? 0.016 0.074

Association des probabilités prédites et des réponses observées

Percent Concordant 58.1 Somers' D 0.555

Percent Discordant 2.6 Gamma 0.914
Percent Tied 39.2 Tau-a 0.372
Pairs 13321 c 0.778

12
 Examen de R´egression sur variables cat

´egorielles ENSAI - 21 janvier 2009

Document autoris´e : une feuille A4

manuscrite Calculatrice fournie par l’´ecole,

dur´ee : 2h

Pr´eambule : Le sujet est compos´e de trois exercices ind´ependants. La qualit´e de la r´edaction sera prise
en compte. Il n’est pas n´ecessaire de traˆıter la totalit´e du sujet pour obtenir la note maximale.

Exercice 1
On dispose de 5 observations (xi, yi) d’une variable `a expliquer Y `a valeurs dans {0, 1} par une variable
explicative X `a valeur dans R. Les donn´ees sont r´esum´ees dans le tableau suivant :

xi -0.16 -0.41 0.33 -0.18 1.45

yi 1 0 0 1 1

1. Que vaut la log-vraisemblance du mod`ele satur´e ?

2. Ecrire le mod`ele logistique permettant d’expliquer Y par X.
3. On note pˆ(x i ), i = 1, . . . , 5 les estimations des probabilit´es P(Y = 1|X = xi) par le mod`ele
logistique. Ecrire la vraisemblance du mod`ele en fonction des pˆ(x i ), i = 1, . . . , 5.
4. Les estimations pˆ(x i ), i = 1, . . . , 5 sont donn´ees dans le tableau suivant

pˆ(x 1 ) pˆ(x 2 ) pˆ(x 3 ) pˆ(x 4 ) pˆ(x 5 )

0.519 0.450 0.651 0.514 0.866

Calculer la log-vraisemblance de ce mod`ele (penser `a prendre le logarithme n´ep´erien sur la calculatrice

!).
5. En d´eduire la d´eviance.
Exercice 2
Dans le cadre d’une surveillance de la population angevine, le CHU d’Angers a men´e une ´etude sur 100 individus
afin de d´eterminer leur aptitude `a ronfler. Les variables consid´er´ees sont :
– AGE : en ann´ee ;
– POIDS : en kg ;
– TAILLE : en cm ;
– ALCOOL : nombre de verres bus par jour (en ´equivalent verre de vin rouge) ;
– SEXE : sexe de la personne (0=homme, 1=femme) ;
– RONFLE : Diagnostic de ronflement (1=ronfle, 0=ne ronfle pas) ;
– TABA : Comportement au niveau du tabac (1=fumeur, 0=non fumeur).

Le but de cette ´etude est d’essayer d’expliquer le ronflement (variable RONFLE) par les six variables explicatives
pr´esent´ees ci-dessus (la variable ALCOOL sera consid´er´ee comme une variable continue).

On consid`ere d’abord le mod`ele logistique construit `a partir des variables X1 =SEXE et X2 =TABA :
logit P(Y = 1|X = x) = β0 + β11{x1=0} + β21{x1=1} + β31{x2=0} + β41{x2=1}, x = (x1, x2). (1)

1. 0n effectue la r´egression logistique sous SAS :

proc logistic data=ronfle descending;
class sexe taba;
model ronfle=sexe taba;
run; 1

12
 Les sorties se trouvent en annexe 1. Quelles sont les contraintes sur les coefficients utilis´ees dans cette
proc´edure ?
2. Donner les estimations βˆj des param`etres βj, j = 0, . . . , 4.
3. Quelle est la probabilit´e de ne pas ronfler estim´ee par ce mod`ele pour une femme qui ne fume pas ?
4. On change la contrainte : on choisit pour les deux variables explicatives 0 comme modalit´e de r´ef
´erence.
On consid`ere donc le mod`ele

logit P(Y = 1|X = x) = α0 + α11{x1=0} + α21{x1=1} + α31{x2=0} + α41{x2=1}

muni des contraintes α1 = 0 et α3 = 0. Donner les estimations αˆ j des param`etres αj, j = 0, . . . , 4
pour ce nouveau jeu de contraintes.
5. Parmi les sorties propos´ees, quelles sont celles qui permettent de comparer le mod`ele d´efini par (1)
au mod`ele logistique d´efini `a partir de la constante uniquement ?
6. Pourquoi les probabilit´es critiques des tests associ´es `a la variable SEXE sont identiques dans les tableaux
Analyse des effets Type 3 et Analyse des estimations de la vraisemblance maximum page 6 ?
7. On compl`ete le mod`ele (1) en ajoutant l’interaction entre les variables SEXE et TABA. Rappeler en
quelques mots la signification de cette interaction.
8. Ecrire la proc´edure SAS qui permet de construire ce nouveau mod`ele.
9. Les sorties se trouvent en annexe 2. Comparer les mod`eles avec et sans interaction en utilisant au moins
trois crit`eres.
On consid`ere maintenant le mod`ele logistique “complet” : celui construit `a partir des 6 variables explicatives
(sans interaction) :
proc logistic data=ronfle descending;
class sexe taba;
model ronfle=age poids taille alcool sexe taba ;
run;
Les sorties se trouvent en annexe 3.
10. En vous inspirant de (1), ´ecrire ce nouveau mod`ele. Vous donnerez ´egalement les contraintes
d’identifia- bilit´e.
11. Quels tests sont r´ealis´es dans le tableau Test de l’hypoth`ese nulle globale : BETA=0 page 10 ?
(on
´ecrira simplement les hypoth`eses ainsi que les lois des statistiques de test : les statistiques de test ne
sont pas demand´ees).
12. Quelle serait la premi`ere variable supprim´ee `a ce mod`ele par un algorithme pas `a pas descendant
(obtenu en ajoutant uniquement l’argument selection=backward dans la proc´edure pr´ec´edente) ?

Exercice 3
On s’int´eresse ici `a un ´echantillon de 371 mineurs travaillant dans des mines de charbon. Ceux-ci sont fr
´equemment atteints d’une maladie des poumons : la pneumoconoisis. Une radio permet de juger si leurs
poumons sont s´ev`erement atteints, l´eg`erement atteints ou sains (variable groupe dont les modalit´es
sont respectivement cod´ees par 3, 2 et 1). On tente d’expliquer le r´esultat de la radio par le nombre d’ann
´ees de travail dans les mines (variable period).
On souhaite discriminer la variable Y (groupe) `a 3 modalit´es par la variable X (period) qui sera consid´er
´ee ici comme continue.
1. Ecrire le mo d`ele polytomique multinomial (on prendra comme modalit´e de r´ef´erence Y = 3,
comme dans la proc´edure catmod). On pr´ecisera ´egalement les probabilit´es P(Y = j|X = x), j =
1, 2, 3 en fonction des coefficients du mod`ele.
2. Les sorties de la proc´edure catmod sont donn´ees en annexe 4. Le tableau Profils de population
page 11 r´esume les observations de la variable explicative. Par exemple, la premi`ere ligne signifie que
l’on a observ´e 98 fois une valeur de 5.8 pour la variable X. Donner les valeurs des estimations des
param`etres du mod`ele que vous avez ´ecrit `a la question pr´ec´edente.
3. Ce mod`ele est-il satur´e ? Justifier.
4. u` proviennent les 12 degr´es de libert´e dans la derni`ere ligne du tableau Maximum Likelihood
D’o
Analysis of Variance page 11 ?
2
5. On souhaite maintenant ajuster un mod`ele polytomique ordonn´e. Ecrire les diff´erentes ´etapes
qui permettent de d´efinir un tel mod`ele : on partira de la mod´elisation par une variable latente pour
arriver
`a l’´ecriture du mod`ele. On pr´ecisera notamment les probabilit´es P(Y ≤ j|X = x), j = 1, 2, 3 en
fonction des coefficients du mod`ele.
6. Expliquer pourquoi cette mod´elisation engendre une hypoth`ese de proportionnalit´e des odds ratio.
7. Les sorties de la proc´edure
proc logistic data=mineurs;
model groupe=period;
run;
se trouvent en annexe 5. La mod´elisation par un mod`ele polytomique ordonn´e est-elle acceptable ?
8. Quelles sont les probabilit´es d’ˆetre l´eg`erement atteints par la maladie pour une personne ayant
travaill´e 30 ans dans la mine estim´ees par les deux mod`eles construits (polytomique multinomial et
ordonn´e).

3
4
Annexe 1 : exercice 2

proc logistic data=ronfle descending;

class sexe taba;
model ronfle=sexe taba;
run;

Le Système SAS

The LOGISTIC

Procedure

Informations sur le modèle

Data Set [Link]

Response Variable ronfle
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 100

Number of Observations Used 100

Profil de réponse

Valeur Fréquence
ordonnée ronfle totale

1 1 35
2 0 65

Probability modeled is ronfle=1.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

sexe 0 1
1 -1

taba 0 1
1 -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 131.489 129.151

SC 134.094 136.966
-2 Log L 129.489 123.151

5
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 6.3388 2 0.0420

Score 5.9940 2 0.0499
Wald 5.6112 2 0.0605

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

sexe 1 4.5833 0.0323

taba 1 2.7662 0.0963

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -0.8754 0.2796 9.8058 0.0017

sexe 0 1 0.6339 0.2961 4.5833 0.0323
taba 0 1 0.3937 0.2367 2.7662 0.0963

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

sexe 0 vs 1 3.553 1.113 11.343

taba 0 vs 1 2.198 0.869 5.558

Association des probabilités prédites et des réponses observées

Percent Concordant 44.9 Somers' D 0.262

Percent Discordant 18.7 Gamma 0.411
Percent Tied 36.4 Tau-a 0.120
Pairs 2275 c 0.631

6
Annexe 2 : exercice 2
The LOGISTIC Procedure

Informations sur le

modèle

Data Set [Link]

Response Variable ronfle
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 100

Number of Observations Used 100

Profil de réponse

Valeur Fréquence
ordonnée ronfle totale

1 1 35
2 0 65

Probability modeled is ronfle=1.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

sexe 0 1
1 -1

taba 0 1
1 -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 131.489 131.011

SC 134.094 141.432
-2 Log L 129.489 123.011

Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi

2 Likelihood Ratio 6.4784 3 0.0905

Score 5.9945 3 0.1119
Wald 5.2920 3 0.1516

7
 Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

sexe 1 4.2070 0.0403

taba 1 2.0759 0.1496
sexe*taba 1 0.1333 0.7150

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -0.9311 0.3283 8.0454 0.0046

sexe 0 1 0.6733 0.3283 4.2070 0.0403
taba 0 1 0.4730 0.3283 2.0759 0.1496
sexe*taba 0 0 1 -0.1199 0.3283 0.1333 0.7150

Association des probabilités prédites et des réponses observées

Percent Concordant 44.9 Somers' D 0.262

Percent Discordant 18.7 Gamma 0.411
Percent Tied 36.4 Tau-a 0.120
Pairs 2275 c 0.631

8
Annexe 3 : exercice 2

proc logistic data=ronfle descending;

class sexe taba;
model ronfle=age poids taille alcool sexe taba ;
run;

Le Système SAS

The LOGISTIC

Procedure

Informations sur le modèle

Data Set [Link]

Response Variable ronfle
Number of Response Levels 2
Model binary logit
Optimization Technique Fisher's scoring

Number of Observations Read 100

Number of Observations Used 100

Profil de réponse

Valeur Fréquence
ordonnée ronfle totale

1 1 35
2 0 65

Probability modeled is ronfle=1.

Informations sur le niveau de classe

Variables
de
Classe Valeur création

sexe 0 1
1 -1

taba 0 1
1 -1

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 131.489 123.418

SC 134.094 141.654
-2 Log L 129.489 109.418

9
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 20.0716 6 0.0027

Score 17.7452 6 0.0069
Wald 14.8584 6 0.0214

Analyse des effets Type 3

Khi 2
Effet DF de Wald Pr > Khi 2

age 1 7.1091 0.0077

poids 1 0.2162 0.6420
taille 1 0.1009 0.7508
alcool 1 7.5462 0.0060
sexe 1 0.9371 0.3330
taba 1 4.6294 0.0314

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 -6.0760 5.9893 1.0292 0.3104

age 1 0.0621 0.0233 7.1091 0.0077
poids 1 -0.0154 0.0332 0.2162 0.6420
taille 1 0.0151 0.0475 0.1009 0.7508
alcool 1 0.2365 0.0861 7.5462 0.0060
sexe 0 1 0.3261 0.3369 0.9371 0.3330
taba 0 1 0.6003 0.2790 4.6294 0.0314

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

age 1.064 1.017 1.114

poids 0.985 0.923 1.051
taille 1.015 0.925 1.114
alcool 1.267 1.070 1.500
sexe 0 vs 1 1.920 0.513 7.189
taba 0 vs 1 3.322 1.113 9.917

Association des probabilités prédites et des réponses observées

Percent Concordant 75.6 Somers' D 0.516

Percent Discordant 24.0 Gamma 0.518
Percent Tied 0.4 Tau-a 0.237
Pairs 2275 c 0.758

10
Annexe 4 : exercice 3

proc catmod data=mineurs;

direct period;
model groupe=period;
run;
Le Système SAS

The CATMOD Procedure

Récapitulatif sur les données

Response groupe Response Levels 3

Weight Variable poids Populations 8
Data Set MINEURS Total Frequency 371
Frequency Missing 0 Observations 22

Profils de population

Échantillon period Sample

Size
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
1 5.8 98
2 15 54
3 21.5 43
4 27.5 48
5 33.5 51
6 39.5 38
7 46 28
8 51.5 11

Profils de réponse

Réponse groupe
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
1 1
2 2
3 3

Analyse de vraisemblance max.

Maximum likelihood computations converged.

Maximum Likelihood Analysis of Variance

Source DF Khi 2 Pr > Khi

2 ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Intercept 2 123.22 <.0001
period 2 60.61 <.0001

Likelihood Ratio 12 13.93 0.3053

Analysis of Maximum Likelihood Estimates

Nombre Erreur Khi-

Parameter de fonctions Estimation std 2 Pr > Khi 2
ƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒƒ
Intercept 1 5.0598 0.5964 71.98 <.0001
2 0.7682 0.7377 1.08 0.2977
period 1 -0.1093 0.0165 44.04 <.0001
2 -0.0257 0.0198 1.69 0.1932

11
Annexe 5 : exercice 3

proc logistic data=mineurs;

model groupe=period;
run;

Le Système SAS

The LOGISTIC

Procedure

Informations sur le modèle

Data Set [Link]

Response Variable groupe
Number of Response Levels 3
Weight Variable poids
Model cumulative logit
Optimization Technique Fisher's scoring

Number of Observations Read 22

Number of Observations Used 22
Sum of Weights Read 371
Sum of Weights Used 371

Profil de réponse

Valeur Fréquence Pondération

ordonnée groupe totale totale

1 1 8 289.00000
2 2 7 38.00000
3 3 7 44.00000

Probabilities modeled are cumulated over the lower Ordered Values.

État de convergence du modèle

Convergence criterion (GCONV=1E-8) satisfied.

Score Test for the Proportional Odds Assumption

Khi 2 DF Pr > Khi 2

0.0002 1 0.9881

Statistiques d'ajustement du modèle

Coordonnée à l'origine
Coordonnée à l'origine et
Critère uniquement covariables

AIC 509.162 422.919

SC 511.344 426.192
-2 Log L 505.162 416.919

12
 Test de l'hypothèse nulle globale : BETA=0

Test Khi 2 DF Pr > Khi 2

Likelihood Ratio 88.2432 1 <.0001

Score 80.7246 1 <.0001
Wald 64.5206 1 <.0001

Analyse des estimations de la vraisemblance maximum

Erreur Khi 2
Paramètre DF Estimation std de Wald Pr > Khi 2

Intercept 1 1 3.9559 0.4096 93.2527 <.0001

Intercept 2 1 4.8691 0.4437 120.4349 <.0001
period 1 -0.0959 0.0119 64.5206 <.0001

Estimations des rapports de cotes

Point 95% Limites de confiance

Effet Estimate de Wald

period 0.909 0.888 0.930

Association des probabilités prédites et des réponses observées

Percent Concordant 47.8 Somers' D 0.087

Percent Discordant 39.1 Gamma 0.100
Percent Tied 13.0 Tau-a 0.061
Pairs 161 c 0.543

13