2.2 Modèles linéaires généralisés (GLM)
Les GLM, comme pour la régression linéaire, sont utilisés pour déterminer et quantifier la relation entre une variable réponse et des variables explicatives incluses dans le paramètre de moyenne de la distribution.
Cette modélisation diffère de la régression linéaire classique en deux points majeurs:
- La variable réponse est un membre de la famille exponentielle linéaire;
- On peut utiliser une transformation pour inclure les variables explicatives dans le paramètre de moyenne.
Une conséquence de l’utilisation des distributions de la famille exponentielle linéaire est la présence d’hétéroscédasticité, signifiant que la variance n’est pas nécessairement constante en fonction de la moyenne.
L’utilisation des GLM est importante en assurance non-vie où, comme expliqué précedemment, les données sont rarement normales. De plus, la relation entre la sinistralité et les caractéristiques du risque est rarement additive, mais plutôt multiplicative.
2.2.1 Caractéristiques des GLM
Supposons une variable réponse Y, ainsi la fonction de densité (ou de probabilité) s’exprime comme:
f(y)=c(y,ϕ)exp[yθ−a(θ)ϕ],
avec:
g(μ)=XTβ,
appelé la fonction de lien, qui est une transformation de la moyenne μ. La fonction de lien g(μ) est, quant à elle, reliée linéairement avec les variables explicatives exprimées dans le vecteur X.
Pour complètement définir un modèle GLM, deux éléments sont ainsi nécessaires:
Le choix de a(θ) déterminera la distribution de la variable réponse;
Le choix de g(μ) indiquera comment la moyenne est reliée avec les variables explicatives.
Finalement, pour travailler avec la théorie des GLM lorsque nous travaillons avec un échantillon de m observations, ces observations seront supposées indépendantes.
2.2.2 Fonction de lien et lien canonique
Si g(μ)=θ, alors g est appelé le lien canonique correspondant à a(θ). Dans une telle situation, θ=XTβ. Comme nous le verrons plus tard, choisir le lien canonique dans un modèle GLM comporte plusieurs avantages.
Exemple 2.10 Trouvez le lien canonique g(μ) pour une loi de Poisson.
Nous savons que:
Pr[Y=y]=λye−λy!
et que:
c(y,ϕ)=−ln(y!)θ=ln(λ)a(θ)=λ=exp(θ)ϕ=1
Ainsi, puisque la moyenne μ=λ pour une loi de Poisson:
θ=ln(λ)=ln(μ)
et donc g(μ)=ln(μ). On appelle ce lien, le lien logarithmique.
Exemple 2.11 Montrez les liens canoniques g(μ) suivants:
Distribution | g(μ) | Fonction de lien canonique |
---|---|---|
Binomiale | lnμ1−μ | logit |
Poisson | lnμ | logarithmique |
Normale | μ | Identité |
Gamma | μ−1 | puissance (-1) |
Inv.-gaus. | μ−2 | puissance (-2) |
En exercice à la maison.
2.2.3 Estimation par maximum de vraisemblance
Les estimateurs du maximum de vraisemblance des paramètres $ {}$ se calculent de manière similaire à ce qu’on avait vu plus tôt. La fonction de logvraisemblance s’exprime comme:
ℓ(β,ϕ)=n∑i=1lnf(yi;β,ϕ)=n∑i=1(ln(c(yi,ϕ))+[yiθi−a(θi)ϕ])=1ϕn∑i=1[yiθi−a(θi)]+n∑i=1ln(c(yi,ϕ))
avec ηi=XTiβ=g(μi) pour des variables yi,i=1,...,n indépendantes. En utilisant la règle de dérivation en chaîne, on peut estimer le MLE des βj,j=1,...,p comme:
δℓδβ=n∑i=1δlnf(yi;β,ϕ)δθiδθiδβ=n∑i=11ϕ[yi−a′(θi)]δθiδβ.
Puisque nous savons que E[Yi]=μi=a′(θi), et que Var[Yi]=ϕa″(θi), nous obtenons:
δμiδβ=a″(θi)δθiδβ=1ϕVar[Yi]δθiδβ⇒δθiδβ=δμiδβϕVar[Yi]
Ainsi:
δℓδβ=n∑i=1(yi−a′(θi))Var[Yi]δμiδβ.
Le dernière dérivée peut se développer encore une fois selon le principe des dérivations en chaîne:
δμiδβ=δμiδηiδηiδβ
où:
δηiδβ=Xiδμiδηi=(δηiδμi)−1=(g′(μi))−1.
Ainsi, la condition de premier ordre s’exprime comme:
δℓδβj=n∑i=1(yi−μi)g′(μi)Var[Yi]Xij=0
et puisque Var[Yi]=ϕa″(θi), on peut réécrire la formule précédente par:
δℓδβj=n∑i=1(yi−μi)g′(μi)a″(θi)Xij=0
et ainsi montrer que les équations de vraisemblance relatives à β peuvent être résolues sans se préoccuper du paramètre ϕ.
Dans le cas où le modèle utilise le lien canonique pour g(μ), il est possible de simplifier l’équation de premier ordre. En effet, dans ce cas, on sait que θi=g(μi), ainsi:
a″(θi)=δa′(θi)δθi=δa′(θi)δg(μi)=δa′(θi)δμiδμiδg(μi)=δa′(θi)δa′(θi)δμiδg(μi)=δ(μi)δg(μi)=1/g′(μi).
En conséquence, nous avons:
δℓδβj=n∑i=1(yi−μi)Xij=0.
On peut remarquer que cette relation traduit l’orthogonalité entre les variables explicatives et les résidus (ce qui est similaire à ce qu’on obtient lors d’une régression linéaire classique).
Exemple 2.12 On suppose que Yi∼Poisson(λi), avec λi=exp(XTiβ). Trouvez la condition de premier ordre pour l’estimation par maximum de vraisemblance de $ {}$.
Nous savons que le lien canonique pour une loi de Poisson est le lien logarithmique. Ainsi, XTiβ=log(μi)=log(λi), et donc λi=exp(XTiβ).
Ainsi, la condition de premier ordre est simplement:
δℓδβj=n∑i=1(yi−μi)Xij=n∑i=1(yi−λi)Xij=0
ou encore, en notation vectorielle:
δℓδβ=n∑i=1(yi−λi)Xi=0
On peut aussi utiliser l’autre forme, plus générale, de la condition de premier ordre:
δℓδβj=n∑i=1(yi−μi)g′(μi)Var[Yi]Xij=n∑i=1(yi−λi)1λiλiXij=n∑i=1(yi−λi)Xi=0
Exemple 2.13 On suppose encore que Yi∼Poisson(λi), avec λi=exp(XTiβ). Trouvez la condition de premier ordre en utilisant directement la fonction de vraisemblance de la Poisson.
L(β)=n∏i=1λyiie−λiyi!ℓ(β)=n∑i=1[yiln(λi)−λi−ln(yi!)]ℓ(β)=n∑i=1[yiln(exp(XTiβ))−exp(XTiβ)−ln(yi!)]
Pour trouver les estimateurs, on dérive par rapport à β et on pose la dérivée égale à zéro pour obtenir le maximum:
δℓ(λi)δβ0=n∑i=1(yi−λi)Xi=0
Nous reverrons plus en détails la loi de Poisson et plusieurs généralisations dans un prochain chapitre.
Exemple 2.14 On suppose encore que Yi∼Gamma(μi,ν). Trouvez la condition de premier ordre en supposant:
- μi=(XTiβ)−1;
- μi=exp(XTiβ).
Regardons chacune des situations séparemment:
- μi=(XTiβ)−1 est la fonction de lien canonique. Ainsi, on peut directement utiliser la fonction de premier ordre de résolution des paramètres:
δℓδβ=n∑i=1(yi−μi)Xi=0
On peut aussi utiliser l’autre forme, plus générale, de la condition de premier ordre:
δℓδβj=n∑i=1(yi−μi)g′(μi)Var[Yi]Xij=n∑i=1(yi−μi)−1μ2iμ2i/νXij=n∑i=1(yi−λi)Xi=0
Finalement, en utilisant directement la fonction de vraisemblance de la gamma, nous avons:
L(β)=n∏i=11Γ(ν)(νyiμi)νexp(−yiνμi)1yiℓ(β)=n∑i=1−log(Γ(ν))+νlog(νyi)−νlog(μi)−yiνμi−log(yi)δℓ(β)δβ=n∑i=1νμi(μi)2Xi−yiνμ2iμ2iXi=n∑i=1(yi−μi)Xi=0
- μi=exp(XTiβ) n’est pas la fonction de lien canonique:
On peut toutefois utiliser la forme plus générale de la condition de premier ordre:
δℓδβj=n∑i=1(yi−μi)g′(μi)Var[Yi]Xij=n∑i=1(yi−μi)1μiμ2i/νXij=n∑i=1(1−yiμi)Xi=0
En utilisant directement la fonction de vraisemblance de la gamma, nous avons:
L(β)=n∏i=11Γ(ν)(νyiμi)νexp(−yiνμi)1yiℓ(β)=n∑i=1−log(Γ(ν))+νlog(νyi)−νlog(μi)−yiνμi−log(yi)δℓ(β)δβ=n∑i=1νμi(μi)Xi−yiνμ2iμiXi=n∑i=1(1−yiμi)Xi=0
2.2.4 Algorithme d’estimation
Nous avons vu que les estimateurs du maximum de vraisemblance ˆβj des paramètres βj sont les solutions à l’équation:
δℓδβj=n∑i=1(yi−μi)g′(μi)Var[Yi]Xij=0.
De façon générale, cette équation ne possède pas de solutions explicites et doit donc être résolue numériquement. Plusieurs techniques peuvent être utilisées pour résoudre le système d’équations à p inconnues.
2.2.4.1 Newton-Raphson
À l’origine, la méthode de Newton-Raphson a été développée pour estimer la racine d’une fonction (continue et dérivable) à une seule dimension.
En choisissant une valeur de départ initiale arbitraire (x0), l’algorithme trouve la fonction f(x0).
L’idée de l’algorithme est de calculer la tangente de la fonction estimée à ce point (pouvant être estimée par la première dérivée de f(.)), afin de longer la tangente et de trouver la racine de la droite engendrée par cette tangente, générant ainsi une solution à l’algorithme (x1).
Si la fonction générée à ce nouveau point (f(x1)) n’est pas assez proche de 0, la procédure recommence avec ce nouveau point de départ.
L’algorithme est illustré dans la figure suivante:

Une version multivariée de cet algorithme pour estimer les paramètres de la fonction. Une autre manière de voir cet algorithme d’estimation est l’approximation linéaire/quadratique de la condition de premier ordre par une série de Taylor.
Avec les GLM, nous devons trouver les racines de la première dérivée de ℓ par rapport à β, ce qui correspond à la fonction f(.) de l’illustration de l’algorithme de Newton-Raphson. Ainsi, nous devons trouver la première dérivée de f(.) ou, finalement, la seconde dérivée de ℓ.
Notons
- U(β) le vecteur gradient de la logvraisemblance, dont la composante j (pour le je paramètre) est:
U(β)=δℓ(β)δβj,
- H(β) la matrice hessienne de ℓ(β), celle dont l’élément (j,k) (la seconde dérivée par rapport aux je et ke paramètres) est:
H(β)=δ2ℓ(β)δβjδβk.
Partant d’une valeur initiale de ˆβ(0) que l’on espère proche de l’esimtateur du maximum de vraisemblance ˆβ, on définit la (r+1)-ème valeur aprochée de ˆβ(r+1) de ˆβ à partir de la r-ème ˆβ(r) par:
ˆβ(r+1)=ˆβ(r)−H−1(ˆβ(r))U(ˆβ(r)).
On arrête l’algorithme lorsque :
ℓ(ˆβ(r+1))−ℓ(ˆβ(r))ℓ(ˆβ(r))≤ϵ
pour un petit ϵ.
Exemple 2.15 On suppose que Yi∼Poisson(λi), avec λi=exp(XTiβ). Trouvez U(β) et H(β).
Nous savons que:
ℓ(β)=n∑i=1−yilog(λi)−λi−log(yi!)
et que:
U(β)=n∑i=1Xi(yi−λi), une vecteur p×1
La matrice hessienne se calcule comme:
H(β)=δ2δβ2ℓ(β)=δδβn∑i=1Xi(yi−λi)=−n∑i=1XiXTiλi, une matrice p×p
Il faut faire attention aux calculs de la matrice hessienne, qui est la seconde dérivée du logvraisemblance et pas seulement la dérivée de l’équation de premier ordre.
2.2.4.2 Exemple numérique: Estimation
Exemple 2.16 On suppose que Yi∼Poisson(λi), avec λi=exp(XTiβ). Trouvez les estimateurs ^β0,^β1,...,^βp si nous avons les données suivantes:
Observations (i) | Territoire | Etat civil | Nb. de sin. (y) |
---|---|---|---|
1 | Urbain | Celibataire | 1 |
2 | Urbain | Divorcé | 0 |
3 | Urbain | Marié | 2 |
4 | Rural | Marié | 3 |
5 | Rural | Divorcé | 1 |
6 | Urbain | Marié | 0 |
7 | Rural | Celibataire | 1 |
8 | Urbain | Celibataire | 2 |
La première étape est de codifier les caractéristiques.
Bien qu’il existe d’autres possibilités, nous codifierons en variables binaires.
X1={1si le territoire est urbain0si le territoire est rural,X2={1si l'assuré est marié0sinon,X3={1si l'assuré est célibataire0sinon
Remarques:
Le territoire a 2 modalités: urbain et rural. Ainsi, on introduit 2-1=1 variable binaire X pour codifier l’information.
L’état civil a 3 modalités: marié, célibataire et divorcé. Dans ce cas, on introduit 3-1=2 variables binaires X pour codifier toutes les possibilités.
- si l’assuré est marié: {X2=1,X3=0};
- si l’assuré est célibataire: {X2=0,X3=1};
- si l’assuré est divorcé: {X2=0,X3=0}.
Ainsi, en incluant une valeur de X0=1, qui sera utilisée pour l’intercept β0, notre tableau de données correspond à:
i | X0 | X1 | X2 | X3 | y |
---|---|---|---|---|---|
1 | 1 | 1 | 0 | 1 | 1 |
2 | 1 | 1 | 0 | 0 | 0 |
3 | 1 | 1 | 1 | 1 | 2 |
4 | 1 | 0 | 1 | 1 | 3 |
5 | 1 | 1 | 0 | 0 | 1 |
6 | 1 | 1 | 1 | 1 | 0 |
7 | 1 | 0 | 0 | 1 | 1 |
8 | 1 | 1 | 0 | 1 | 2 |
Notre modèle est ainsi, pour l’assuré i:
Yi∼Poisson(λi=exp(β0X0,i+β1X1,i+β2X2,i+β3X3,i)).
Nous devons ensuite trouver les paramètres {β0,β1,β2,β3} afin que:
U(ˆβ)=n∑i=1Xi(yi−λi)=0
avec Xi={X0,i,X1,i,X2,i,X3,i}.
Utilisons l’algorithme de Newton-Raphson pour trouver les MLE des β:
ˆβ(r+1)=ˆβ(r)−H−1(ˆβ(r))U(ˆβ(r)).
Pour débuter l’algorithme, nous devons commencer avec des valeurs initiales de β, i.e. {β(0)0,β(0)1,β(0)2,β(0)3}.
Ce choix initial est arbitraire. Mais il est plus simple de débuter avec un effet nul pour {X1,i,X2,i,X3,i}, et choisir une valeur de β(0)0 qui fait en sorte que λ(0)i=¯Y,∀i∈{1,…,8}.
Puisque nous avons une fonction de lien logarithmique, nous avons ainsi:
(ˆβ(0))T={log(¯Y),0,0,0}.
Ainsi, tous les λ(0)i=exp(Xi′β)=¯Y=1.25, pour i=1,...,8.
On peut donc estimer la valeur de U pour la première itération:
U(ˆβ(0))=8∑i=1Xi(ni−ˆλ(0)i), un vecteur 4×1=[1101]×(1−1.25)+[1100]×(0−1.25)+[1111]×(2−1.25)+[1011]×(3−1.25)+[1100]×(1−1.25)+[1111]×(0−1.25)+[1001]×(1−1.25)+[1101]×(2−1.25)=[0−1.501.251.50]
alors que:
H(ˆβ(0))=−8∑i=1XiXTiˆλ(0)i, une matrice 4×4=−[1101]×[1101]×1.25−[1100]×[1100]×1.25−[1111]×[1111]×1.25−[1011]×[1011]×1.25−[1100]×[1100]×1.25−[1111]×[1111]×1.25−[1001]×[1001]×1.25−[1101]×[1101]×1.25=[−10.00−7.5−3.75−7.50−7.50−7.5−2.50−5.00−3.75−2.5−3.75−3.75−7.50−5.0−3.75−7.50]
Ainsi:
\begin{align*} \hat{\beta}^{(1)} =& \hat{\beta}^{(0)} - H^{-1}(\hat{\beta}^{(0)}) U(\hat{\beta}^{(0)}) \\ =& \begin{bmatrix} log(1.25) \\ 0 \\ 0 \\ 0 \end{bmatrix} - \begin{bmatrix} -10.00 & -7.5 & -3.75 & -7.50 \\ -7.50 & -7.5 & -2.50 & -5.00 \\ -3.75 & -2.5 & -3.75 & -3.75 \\ -7.50 & -5.0 & -3.75 &-7.50 \end{bmatrix}^{-1} \begin{bmatrix} 0 \\ -1.50 \\ 1.25 \\ 1.50 \end{bmatrix} \\ =& \begin{bmatrix} 0.2231 \\ 0 \\ 0 \\ 0 \end{bmatrix} - \begin{bmatrix} 0 \\ 0.6000 \\ -0.2667 \\ -0.4667 \end{bmatrix} = \begin{bmatrix} 0.2231 \\ -0.6000 \\ 0.2667 \\ 0.4667 \end{bmatrix} \end{align*}
On peut ainsi recommencer l’algorithme avec ces nouvelles valeurs de \hat{\beta}^{(1)}. Attention, car pour cette nouvelle itération, la valeur des \lambda^{(1)} ne sera pas la même pour tous les assurés car \lambda_i^{(1)} = \exp(\mathbf{X}_i^T \beta^{(1)}).
i | X_0 | X_1 | X_2 | X_3 | y | \lambda^{(1)} |
---|---|---|---|---|---|---|
1 | 1 | 1 | 0 | 1 | 1 | 1.093 |
2 | 1 | 1 | 0 | 0 | 0 | 0.686 |
3 | 1 | 1 | 1 | 1 | 2 | 1.428 |
4 | 1 | 0 | 1 | 1 | 3 | 2.602 |
5 | 1 | 1 | 0 | 0 | 1 | 0.686 |
6 | 1 | 1 | 1 | 1 | 0 | 1.428 |
7 | 1 | 0 | 0 | 1 | 1 | 1.993 |
8 | 1 | 1 | 0 | 1 | 2 | 1.094 |
On aurait ainsi:
On peut donc estimer la valeur de U pour la première itération:
\begin{align*} U(\hat{\beta}^{(0)}) =& \sum_{i=1}^8 \mathbf{X}_i (n_i - \hat{\lambda}_i^{(0)}) \text{, un vecteur } 4 \times 1 \\ =& \begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 1 - 1.093 \end{array} \right) + \begin{bmatrix} 1 \\ 1 \\ 0 \\ 0 \end{bmatrix} \times \left( \begin{array}{c} 0 - 0.686 \end{array} \right) + \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 2 - 1.428 \end{array} \right) + \begin{bmatrix} 1 \\ 0 \\ 1 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 3 - 2.602 \end{array} \right) + \\ & \begin{bmatrix} 1 \\ 1 \\ 0 \\ 0 \end{bmatrix} \times \left( \begin{array}{c} 1 - 0.686 \end{array} \right) + \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 0 - 1.428 \end{array} \right) + \begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 1 - 1.993 \end{array} \right) + \begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} \times \left( \begin{array}{c} 2 - 1.094 \end{array} \right) \\ =& ... \end{align*}
alors que:
\begin{align*} H(\hat{\beta}^{(0)}) =& - \sum_{i=1}^8 \mathbf{X}_i \mathbf{X}_i^T \hat{\lambda}_i^{(0)} \text{, une matrice } 4 \times 4 \\ =& -\begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 0 & 1 \end{bmatrix} \times 1.093 -\begin{bmatrix} 1 \\ 1 \\ 0 \\ 0 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 0 & 0 \end{bmatrix} \times 0.686\\ & -\begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 1 & 1 \end{bmatrix} \times 1.428 -\begin{bmatrix} 1 \\ 0 \\ 1 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 1 & 1 \end{bmatrix} \times 2.602\\ & -\begin{bmatrix} 1 \\ 1 \\ 0 \\ 0 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 0 & 0 \end{bmatrix} \times 0.686 -\begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 1 & 1 \end{bmatrix} \times 1.428\\ & -\begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 0 & 0 & 1 \end{bmatrix} \times 1.993 -\begin{bmatrix} 1 \\ 1 \\ 0 \\ 1 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 0 & 1 \end{bmatrix} \times 1.094 \\ =& ... \end{align*}
Ainsi, après convergence, on aurait
\begin{align*} \hat{\beta} = \begin{bmatrix} -0.2231 \\ -0.4700 \\ 0.2231 \\ 0.7985 \end{bmatrix} \end{align*}
Et donc, on obtiendrait les primes suivantes pour les 6 profils possibles:
profil | X_0 | X_1 | X_2 | X_3 | \hat{\lambda} |
---|---|---|---|---|---|
1 | 1 | 1 | 1 | 0 | 0.625 |
2 | 1 | 1 | 0 | 1 | 1.111 |
3 | 1 | 1 | 0 | 0 | 0.500 |
4 | 0 | 1 | 1 | 0 | 1.000 |
5 | 0 | 1 | 0 | 1 | 1.778 |
6 | 0 | 1 | 0 | 0 | 0.800 |
Il est assez simple de refaire ce dernier exemple en utilisant R.
### entrées les données
<- matrix(c(1,1,1,1,1,1,1,1,
x 1,1,1,0,1,1,0,1,
0,0,1,1,0,1,0,0,
1,0,1,1,0,1,1,1),
8, 4)
<- matrix(c(1,0,2,3,1,0,1,2), 8, 1)
y
####
<- mean(y)
ybar <- matrix(c(log(ybar), 0, 0, 0), 4, 1)
beta
### 10 itérations
for(iter in 0:9){
print(c(iter, beta))
<-0
U for(i in 1:8){
<- as.numeric(exp(x[i,] %*% beta))
lambda <- U + t(x[i,])*(y[i] - lambda)
U
}
U
<-0
H for(i in 1:8){
<- as.numeric(exp(x[i,] %*% beta))
lambda <- H - lambda*(x[i,] %*% t(x[i,]))
H
}
H
<- beta - solve(H) %*% t(U)
beta }
## [1] 0.0000000 0.2231436 0.0000000 0.0000000 0.0000000
## [1] 1.0000000 0.2231436 -0.6000000 0.2666667 0.4666667
## [1] 2.0000000 -0.1688347 -0.4791742 0.2259235 0.7520675
## [1] 3.0000000 -0.2220958 -0.4700478 0.2231551 0.7974987
## [1] 4.0000000 -0.2231430 -0.4700036 0.2231436 0.7985072
## [1] 5.0000000 -0.2231436 -0.4700036 0.2231436 0.7985077
## [1] 6.0000000 -0.2231436 -0.4700036 0.2231436 0.7985077
## [1] 7.0000000 -0.2231436 -0.4700036 0.2231436 0.7985077
## [1] 8.0000000 -0.2231436 -0.4700036 0.2231436 0.7985077
## [1] 9.0000000 -0.2231436 -0.4700036 0.2231436 0.7985077
2.2.5 Fisher scoring
Au lieu d’utiliser \mathbf{H}(\mathbf{\beta}) dans l’algorithme de Newton-Raphson, il a été suggéré par Ronald Fisher d’utiliser la matrice d’information de Fisher \mathbf{I(\beta)}. En effet, on peut montrer que:
\mathbf{H}(\mathbf{\beta}) \approx E[\mathbf{H}(\mathbf{\beta})] = - \mathbf{I(\beta)},
Ainsi, nous avons:
\mathbf{\hat{\beta}}^{(r+1)} = \mathbf{\hat{\beta}}^{(r)} + \mathbf{I}^{-1}(\mathbf{\hat{\beta}}^{(r)}) \mathbf{U}(\mathbf{\hat{\beta}}^{(r)}) .
La substitution de la matrice hessienne par l’information de Fisher dans l’algorithme d’estimation s’appelle le Fisher scoring, ou l’algorithme de score.
Dans le but d’éviter le calcul d’une deuxième dérivée du logvraisemblance, et sachant qu’il peut être montré que:
E\left[\frac{\delta^2 \ell(\mathbf{\beta})}{\delta \mathbf{\beta}^2}\right] = - E\left[\left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right) \left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right)^T \right],
il est aussi possible d’utiliser l’espérance du produit croisé de la première dérivée, au lieu de l’espérance de la seconde dérivée dans l’algorithme de score, dans l’algorithme d’estimation.
Exemple 2.17 On suppose que Y_i \sim Poisson(\lambda_i), avec \lambda_i = \exp(\mathbf{X_i}' \mathbf{\beta}).
Montrez que:
E\left[\frac{\delta^2 \ell(\mathbf{\beta})}{\delta \mathbf{\beta}^2}\right] = - E\left[\left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right) \left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right)^T \right],
\begin{eqnarray*} E\left[\left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right) \left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right)^T \right] &=& E\left[\sum_{i=1}^n \left(\mathbf{X}_{i} \left(y_i - \lambda_i \right) \right) \left(\mathbf{X}_{i} \left(y_i - \lambda_i \right)\right)^T \right]\\ &=& \sum_{i=1}^n \mathbf{X}_{i} \mathbf{X}_{i}^T E\left[\left(y_i - \lambda_i \right)^2 \right]\\ &=& \sum_{i=1}^n \mathbf{X}_{i} \mathbf{X}_{i}^T Var[Y_i] \\ &=& \sum_{i=1}^n \mathbf{X}_{i} \mathbf{X}_{i}^T \lambda_i \end{eqnarray*}
2.2.6 Variance des estimateurs
Un autre résultat majeur des statistiques mathématiques est qu’il peut être prouvé que la deuxième dérivée du logvraisemblance peut être utilisée pour estimer la variance des MLE. Nous y reviendrons.
Ainsi, nous pouvons donc utiliser l’information de Fisher, ou l’espérance du produit croisé de la première dérivée pour estimer la matrice variance-covariance des MLE. Plus précisément, nous avons:
\begin{eqnarray*} \sigma^2(\hat{\beta}) &=& - E\left[\frac{\delta^2 \ell(\mathbf{\beta})}{\delta \mathbf{\beta}^2}\right]^{-1} \end{eqnarray*}
ou encore:
\begin{eqnarray*} \sigma^2(\hat{\beta}) &=& E\left[\left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right) \left(\frac{\delta \ell(\mathbf{\beta})}{\delta \mathbf{\beta}}\right)^T \right]^{-1} \end{eqnarray*}
2.2.6.1 Exemple numérique: Variance
Exemple 2.18 En reprenant le modèle et les données indiquées dans l’exemple numérique précédent, obtenez la matrice variance-covariance des estimateurs \hat{\beta}.
Nous avions montré que les MLE de \beta étaient:
\begin{align*} \hat{\beta} = \begin{bmatrix} -0.2231 \\ -0.4700 \\ 0.2231 \\ 0.7985 \end{bmatrix} \end{align*}
Comme lors de l’exemple précédent, on peut donc estimer la matrice hessienne:
\begin{align*} H(\hat{\beta}) =& \begin{bmatrix} -10.00 & -6.00 & -5.00 & -9.00\\ -6.00 & -6.00 & -2.78 & -5.00\\ -5.00 & -2.78 & -5.00 & -5.00\\ -9.00 & -5.00 & -5.00 & -9.00\\ \end{bmatrix} \end{align*}
Ainsi la matrice variance-covariance peut s’exprimer comme:
\begin{eqnarray*} \sigma^2(\hat{\beta}) &=& - H(\hat{\beta})^{-1} = \begin{bmatrix} 1.450& -0.450& 0.000& -1.200\\ -0.450& 0.450& 0.000& 0.200 \\ 0.000& 0.000& 0.450& -0.250 \\ -1.200& 0.200& -0.250& 1.339 \\ \end{bmatrix} \end{eqnarray*}
Exemple 2.19 En reprenant les résultats de l’exemple précédent, trouvez:
- La variance de \hat{\beta}_0;
- La variance de \hat{\beta}_2;
- La covariance entre \hat{\beta}_1 et \hat{\beta}_2.
Il suffit de voir que la matrice variance-covariance représente:
et donc:
- Var[\hat{\beta}_0] = 1.450;
- Var[\hat{\beta}_2] = 0.4500;
- Cov[\hat{\beta}_1, \hat{\beta}_2] = 0.000.
2.2.7 Déviance
Définition 2.2 Le modèle saturé correspond au modèle contenant autant de paramètres que d’observations et fournissant ainsi une description parfaite des données. Le modèle saturé est caractérisé par \hat{\mu}_i = y_i, pour i=1,\ldots,n.
La fonction de vraisemblance du modèle saturé est notée L(y|y) et la fonction de logvraisemblance \ell(y|y). En pratique, ce modèle n’est pas intéressant car il se limite à reproduire exactement les observations passées, sans les résumer.
Définition 2.3 La déviance d’un modèle ajusté par \hat{\beta}, notée D, est définie comme une mesure de distance entre l’ajustement de ce modèle et le modèle saturé. Ainsi:
D = 2\left(\ell(y|y) - \ell(\hat{\mu}|y)\right),
où \ell(y|y) est la fonction logvraisemblance du modèle saturé et \ell(\hat{\mu}|y) la fonction logvraisemblance du modèle ajusté par \hat{\beta}.
Une grande valeur de D laisse penser que le modèle ajusté est de piètre qualité. D est aussi appelé la déviance réduite dans le cadre des modèles linéaires généralisés alors que la déviance non réduite est elle donnée par D^* = \phi D.
En pratique, on jugera le modèle de mauvaise qualité si:
D_{\text{observé}} > \chi^2_{n-p-1; 1-\alpha}
où \chi^2_{n-p-1; 1-\alpha} est le quantile d’ordre 1-\alpha de la loi chi-deux à n-p-1 degrés de liberté, n le nombre d’observations et $p +1 $ le nombre de paramètres dans le modèle (\beta_0, \ldots, \beta_p).
2.2.7.1 Exemple numérique: Déviance
Exemple 2.20 Trouvez la déviance dans un modèle de régression Poisson.
Nous savons que la fonction de probabilité d’une loi de Poisson s’exprime comme:
\Pr(Y =y) = \frac{\lambda^y e^{-\lambda}}{y!}
Ainsi, la fonction logvraisemblance d’un modèle staturé peut s’exprimer comme:
\begin{eqnarray*} \ell(y|y) &=& \ln\left( \prod_{i=1}^n \frac{\lambda_i^{y_i} e^{-\lambda_i}}{y_i!} \right), \text{ avec } \lambda_i = y_i \\ &=& \ln\left( \prod_{i=1}^n \frac{y_i^{y_i} e^{-y_i}}{y_i!} \right) \end{eqnarray*}
Ainsi:
\begin{eqnarray*} D(\boldsymbol{y}, \hat{\beta}) &=& 2 \ln\left( \prod_{i=1}^n \frac{y_i^{y_i} e^{-y_i}}{y_i!} \right) - 2 \ln\left( \prod_{i=1}^n \frac{\widehat{\lambda}_i^{y_i} e^{-\widehat{\lambda_i}}}{y_i!} \right) \\ &=& 2 \left( \left(\sum_{i=1}^n y_i \ln y_i -y_i - \ln(y_i!)\right) - \left(\sum_{i=1}^n y_i \ln \widehat{\lambda}_i - \widehat{\lambda}_i - \ln(y_i!)\right) \right) \\ &=& 2 \sum_{i=1}^n \left( y_i \ln y_i - y_i \ln \widehat{\lambda}_i + y_i - \widehat{\lambda}_i \right) \\ &=& 2 \sum_{i=1}^n \left( y_i \ln\left(\frac{y_i}{\widehat{\lambda}_i}\right) + (y_i - \widehat{\lambda}_i) \right)\\ &=& 2 \sum_{i=1}^n \left( y_i \ln\left(\frac{y_i}{\widehat{\lambda}_i}\right) \right) \end{eqnarray*}
car nous savons que le modèle de Poisson suppose un équilibre financier:
\begin{eqnarray*} \frac{\delta \ell}{\delta \beta_)} &=& \sum_{i=1}^n (y_i - \widehat{\lambda}_i) X_{i0} = 0\\ &=& \sum_{i=1}^n (y_i -\widehat{\lambda}_i) \\ &\Rightarrow& \sum_{i=1}^n y_i = \sum_{i=1}^n \widehat{\lambda}_i \end{eqnarray*} ::::
Exemple 2.21 Trouvez la déviance dans les modèles de régression normal, binomial, gamma et inverse-gaussien.
À faire en exercice.
Exemple 2.22 En reprenant les résultats du modèle de Poisson de l’exemple précédentt, et en supposant que les paramètres ont été estimés par maximum de vraisemblance, trouvez la déviance du modèle.
Nous savons que:
\begin{eqnarray*} D(y, \hat{\beta}) &=& 2 \sum_{i=1}^n \left( y_i \ln\left(\frac{y_i}{\widehat{\lambda}_i}\right) \right) = \sum_{i=1}^n D_i(y, \hat{\beta}) \end{eqnarray*}
i | X_0 | X_1 | X_2 | X_3 | y | \widehat{\lambda} | D_i |
---|---|---|---|---|---|---|---|
1 | 1 | 1 | 0 | 1 | 1 | 1.093 | -0.211 |
2 | 1 | 1 | 0 | 0 | 0 | 0.686 | 0.000 |
3 | 1 | 1 | 1 | 1 | 2 | 1.428 | 1.459 |
4 | 1 | 0 | 1 | 1 | 3 | 2.602 | 1.800 |
5 | 1 | 1 | 0 | 0 | 1 | 0.686 | 1.386 |
6 | 1 | 1 | 1 | 1 | 0 | 1.428 | 0.000 |
7 | 1 | 0 | 0 | 1 | 1 | 1.993 | -1.151 |
8 | 1 | 1 | 0 | 1 | 2 | 1.094 | 2.351 |
La déviance totale est égale à 5.635192.
Pour tester l’ajustement du modèle, il suffit ainsi de comparer cette valeur de la déviance (5.635192) à une chi-carré avec (8-4 = 4) degrés de liberté.
A un niveau 95%, un table de chi-carré à 4 degrés de liberté nous donne 9.488. On ne peut donc pas rejeter le modèle Poisson à ce niveau.
Autre remarque: on peut aussi voir que \sum_{i=1}^8 \widehat{\lambda}_i = \sum_{i=1}^8 y_i = 8, ce qu’on nommait un équilibre financier.
2.2.8 Autres éléments
Classiquement, dans les ouvrages sur les GLM, beaucoup d’autres éléments sont couverts. Puisqu’ils sont un peu moins importants pour la suite du cours, je laisse la liberté à chacun de lire sur le sujet.
Par contre, au moins deux autres thématiques doivent être couvertes:
1- L’estimation du paramètre de dispersion \phi;
2- Les tests statistiques sur les paramètres.
Nous reviendrons un peu plus en détails sur ces deux thèmes dans un prochain chapitre.