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

Introduction

• L’objet de ce chapitre est faire une introduction aux modèles linéaires généralisés d’un point de vue theorique et pratique avec une mise en oeuvre sous R.
• Cette notion qui regroupe tous les modèles (linéaire gaussien, logit, log-linéaire).
• Ils se basent sur une façon pour exprimer l’espérance d’une variable réponse \(y\) en fonction d’une combinaison linéaire des variables explicatives.

Composantes du modèle

Composante 1: la distribution de la variable réponse.

• La composante aléatoire identifie la distribution de probabilités de la variable à expliquer.
• On suppose que l’échantillon statistique est constitué de \(n\) variables aléatoires \(\{y_i;i=1,\ldots,n\}\) indépendantes admettant des distributions issues d’une famille exponentielle.
• Cela signifie que les lois de ces variables ont une densité de la forme :

\[f(y_i;\theta_i,\phi)=\exp \left\{\displaystyle\frac{y_i\theta_i-v(\theta_i)}{u(\phi)}+w(y_i,\phi)\right\}\]

• Cette formulation inclut la plupart des lois usuelles comportant un ou deux paramètres : gaussienne, gaussienne inverse, gamma, Poisson, binomiale….(Exercice)

• Le paramètre \(\theta_i\) est appelé paramètre naturel de la famille exponentielle.

• Calcul de \(E(y_i)\) : Comme \[\displaystyle\int f(y_i;\theta_i,\phi)dy_i=1.\] Donc \[\displaystyle\frac{\partial}{\partial \theta_i} \displaystyle\int f(y_i;\theta_i,\phi)dy_i=0\]

Sous des conditions de régularités, on a \[ \displaystyle\int \displaystyle\frac{y_i-v'(\theta_i)}{u(\phi)} \exp \left\{\displaystyle\frac{y_i\theta_i-v(\theta_i)}{u(\phi)}+w(y_i,\phi)\right\} dy_i=0\] Donc \[ E(y_i)=v'(\theta_i)\]

• Exercice Calculer de la même façon \(\mbox{Var}(y_i)\)

• L’expression de la structure exponentielle se met alors sous la forme canonique en posant :

\[Q(\theta)=\displaystyle\frac{\theta}{u(\phi)}\] \[ a(\theta)=\exp\left\{\displaystyle\frac{-v(\theta)}{u(\phi)}\right\}\] \[b(y)=\exp\{w(y,\phi)\}\] on alors \[f(y_i,\theta_i)=a(\theta_i)b(y_i)\exp\left\{y_i Q(\theta_i)\right\}\]

Composante 2 : les prédicteurs

• Les observations planifiées des variables explicatives sont organisées dans la matrice \(X\in \mathbb{R}^{n\times p}\) de planification d’expérience (design matrix)

• Soit \(\beta\in \mathbb{R}^{p\times 1}\) un vecteur de \(p\) paramètres, le prédicteur linéaire, composante déterministe du modèle, est le vecteur à \(n\) composantes

\[\eta= X\beta\in \mathbb{R}^{n \times 1}\]

Composante 3 : la fonction lien.

• Elle exprime une relation fonctionnelle entre la composante aléatoire et le prédicteur linéaire.

• Soit \(\{\mu_i = E(y_i); i = 1,\dots, n\}\), on pose

\[\eta_i=g(\mu_i)=g(v'(\theta_i)),\;i=1,\ldots,n\]

où \(g,\) appelée fonction lien, est supposée monotone et différentiable.

• Ceci revient donc à écrire un modèle dans lequel une fonction de la moyenne appartient au sous-espace engendré par les variables explicatives :

\[g(\mu_i)=x_i'\beta,\;i=1,\ldots,n\]

• La fonction lien qui associe la moyenne \(\mu_i\) au paramètre naturel est appelée fonction lien canonique. Dans ce cas,

\[g(\mu_i) = \theta_i = x_i'\beta.\]

Donc \(g=(v')^{-1}\).

Exemples

La loi Gaussienne

• Dans le cas d’un échantillon gaussien, les densités d’une famille de lois \(\mathcal{N}(\mu_i,\sigma^2)\) : \[ f(y_i,\mu_i) = \displaystyle\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\displaystyle\frac{(y_i-\mu_i)^2}{2\sigma^2}\right\}\]
• Exercice Calculer les éléments \(Q\), \(a\) et \(b\) de l’expression canonique de la densité.

La loi Bernoulli

• Considérons \(n\) variables \(y_i\) indépendantes de probabilité de succès \(\pi_i\), donc \(E(y_i)=\pi_i\).

• Les fonctions densités sont donc sous la forme :

\[f(y_i,\pi_i)=\pi_i^{y_i}(1-\pi_i)^{(1-y_i)}=(1-\pi_i)\exp\left\{y_i\log \displaystyle\frac{\pi_i}{1-\pi_i}\right\}\]

• Donc le paramètre naturel \[\theta_i=\log \displaystyle\frac{\pi_i}{1-\pi_i}\]

• La fonction lien est donc la fonction logit

La loi de Poisson

• Les \(n\) variables \(y_i\) sont considérées suivre la loi de Poisson de paramètre \(\mu_i=E(y_i)\)

• Les densités : \[f(y_i,\mu_i)=\displaystyle\frac{\mu_i^{y_i} e^{-\mu_i}}{y_i!}=\exp \{-\mu_i\}\displaystyle\frac{1}{y_i!}\exp\left\{y_i\log \mu_i\right\}\]

• Donc le paramètre naturel \[\theta_i=\log \mu_i\]

• La fonction lien est donc la fonction logarithmique

Estimation

L’estimation des paramètres \(\beta_j\) est calculée en maximisant la \(\log\)-vraisemblance du modèle linéaire généralisé.

Expression des moments

• Soit \(l(\theta_i,\phi;y_i)=\log f(y_i;\theta_i,\phi)\) le \(\log\)-vraisemblance de la \(i-\)ième observation.

• L’estimation du maximum du \(\log\)-vraisemblance nécessite la connaissance des dérivées :

\[\begin{array}{rcl} \displaystyle\frac{\partial l}{\partial \theta_i}&=&\displaystyle\frac{y_i-v'(\theta_i)}{u(\phi)}\\ \displaystyle\frac{\partial^2 l}{\partial \theta_i^2}&=&\displaystyle\frac{-v''(\theta_i)}{u(\phi)} \end{array} \]

• Avec les conditions de régularité on a : \[E\left(\displaystyle\frac{\partial l}{\partial \theta_i}\right)=0\;\mbox{ et }\; E\left(\displaystyle\frac{\partial^2 l}{\partial \theta_i^2}\right)=-E\left[\left(\displaystyle\frac{\partial l}{\partial \theta_i}\right)^2\right]\]

• Proposition On a \[E(y_i)=\mu_i=v'(\theta_i)\] et \[ \mbox{Var}(y_i)=v''(\theta_i)u(\phi)\]

Equations de vraisemblance

• Considérons
  • les \(p\) variables explicatives dont les observations sont rangées dans la matrice de plan d’expérience \(X\),
  • \(\beta\in \mathbb{R}^p\)
    
  • le prédicteur linéaire à \(n\) composantes : \[\eta= X\beta\]

• La fonction lien \(g\) est supposée monotone différentiable telle que \(\eta_i=g(\mu_i)\)

• c’est la fonction lien canonique si \(\theta_i=g(\mu_i)\)

• Pour \(n\) observations supposées indépendantes et en tenant compte que \(\theta\) dépend de \(\beta\) , la \(\log-\)vraisemblance s’éccri : \[\mathcal{L}(\beta)=\displaystyle\sum_{i=1}^n \log f(y_i;\theta_i,\phi)= \displaystyle\sum_{i=1}^n l(\theta_i,\phi;y_i)\]

• On a

\[ \displaystyle\frac{\partial l}{\partial \beta_j}= \displaystyle\frac{\partial l}{\partial \theta_i} \displaystyle\frac{\partial \theta_i}{\partial \mu_i} \displaystyle\frac{\partial \mu_i}{\partial \eta_i} \displaystyle\frac{\partial \eta_i}{\partial \beta_j}\]

• Or \[\begin{array}{rcl} \displaystyle\frac{\partial l}{\partial \theta_i}& =& \displaystyle\frac{y_i-v'(\theta_i)}{u(\phi)}= \displaystyle\frac{y_i-\mu_i}{u(\phi)}\\ \displaystyle\frac{\partial \mu_i}{\partial \theta_i}&=& v''(\theta_i)=\mbox{Var}(y_i)/u(\phi)\\ \displaystyle\frac{\partial \eta_i}{\partial \beta_j}&=& x_{i}^j \mbox{ car }\eta_i=x_i'\beta\\ \displaystyle\frac{\partial \mu_i}{\partial \eta_i} &=&(g^{-1})'(\mu_i) \end{array} \]

• Les équation de vraisemblances, pour tout \(j=1,\ldots,p\)

\[ \displaystyle\sum_{i=1}^n \displaystyle\frac{(y_i-\mu_i)x_i^j}{\mbox{Var}(y_i)} \displaystyle\frac{\partial \mu_i}{\partial \eta_i}=0\]

• Exercice Ecrire les équation de vraisemblances dans les cas des modèles Gaussiens et Bernoulli

• Ces équations non-linéaires en \(\beta\) dont la résolution requiert des méthodes itératives dans lesquelles interviennent le Hessien (pour Newton- Raphson) ou la matrice d’information (pour les Scores de Fisher).

• La matrice d’information est la matrice \[\mathcal{I}=X'WX\] où \[ \left[\mathcal{I}\right]_{jk}=E \displaystyle\frac{\partial ^2\mathcal{L}(\beta)}{\partial \beta_j \partial \beta_k}= \displaystyle\sum_{i=1}^n \displaystyle\frac{x_i^jx_j^k}{\mbox{Var}(y_i)} \left(\displaystyle\frac{\partial \mu_i}{\partial \eta_i}\right)^2 \]

• \(W\) est appele la matrice diagonale des pondérations : \[\left[W\right]_{ii}=\displaystyle\frac{1}{\mbox{Var}(y_i)} \left(\displaystyle\frac{\partial \mu_i}{\partial \eta_i}\right)^2\]

• Exercice Calculer \(W\) dans le cas Gaussien et Bernoulli.

Cas de la fonction de lien canonique.

• Dans le cas particulier où la fonction lien du modèle linéaire généralisé utilisée est la fonction lien canonique associée à la structure exponentielle (\(g=(v')^{-1}\)). Alors \[\eta_i =\theta_i=x'_i\beta\] \[\displaystyle\frac{\partial \mu_i}{\partial \eta_i}= \displaystyle\frac{\partial \mu_i}{\partial \theta_i}=v''(\theta_i)\]

• Ainsi

\[ \displaystyle\frac{\partial l_i}{\partial \beta_j}= \displaystyle\frac{y_i-\mu_i}{\mbox{Var}(y_i)}v''(\theta_i)x_i^j= \displaystyle\frac{y_i-\mu_i}{u(\phi)}x_i^j \]

• Donc le terme \(\displaystyle\frac{\partial^2 \mathcal{L}(\beta)}{\partial \beta_j\partial \beta_k}\) ne dépendent plus de \(y_i\), on montre que le Hessien est égal à la matrice d’information et donc les méthodes de résolution du score de Fisher et de Newton-Raphson coïncident.

Qualité d’ajustement

Déviance

• Le modèle estimé est comparé avec le modèle saturé, c’est-à-dire le modèle possédant autant de paramètres que d’observations et estimant donc exactement les données.

• Cette comparaison est basée sur l’expression de la déviance \(D\) des \(\log\)-vraisemblances \(\mathcal{L}\) et \(\mathcal{L}_{sat}\) :

\[D=-2(\mathcal{L}-\mathcal{L}_{sat})\]

• On montre qu’asymptotiquement, \(D\) suit une loi du \(\chi^2\) à \(n-p\) degrés de liberté ce qui permet de construire un test de rejet ou d’acceptation du modèle selon que la déviance est jugée significativement ou non importante.

Test de Pearson

• Un test du \(\chi^2\) est également utilisé pour comparer les valeurs observées \(y_i\) à leur prévision par le modèle. La statistique du test est définie par

\[ X^2=\displaystyle\sum_{i=1}^n \displaystyle\frac{(y_i-\widehat{\mu_i})}{\widehat{\mbox{Var}}(\widehat{\mu_i})}\]

• On montre que asymptotiquement \(X^2\) a la même loi que la déviance.

• Remarque: Quand les deux approches conduisent à des résultats peu différents et, dans le cas contraire, c’est une indication de mauvaise approximation de la loi asymptotique.

• Sachant que l’espérance d’une loi du \(\chi^2\) est son nombre de degrés de liberté et, connaissant les aspects approximatifs des tests construits, l’usage est souvent de comparer les statistiques avec le nombre de degrés de liberté. le modèle peut être jugé satisfaisant pour un rapport \(D/ddl\) plus petit que 1.

Tests

On propose deux critères pour le choix du modèle

Rapport de vraisemblance

• Ce test permet de comparer un modèle avec un modèle réduit (deux modèles emboîtés)

• Le rapport de vraisemblance ou la différence de déviance est une évaluation de l’apport des variables explicatives supplémentaires dans l’ajustement du modèle.

• La différence des déviances entre deux modèles emboîtés respectivement à \(q_1\) et \(q_2(q_2> q_1)\) variables explicatives

\[\begin{array}{rcl} D_2-D_1 &=& 2(\mathcal{L}_1-\mathcal{L}_{sat})-2(\mathcal{L}_2-\mathcal{L}_{sat}) \\ &=& 2(\mathcal{L}_1-\mathcal{L}_{2}) \end{array} \]

• \(D_2-D_1\approx \chi^2(q_2-q_1)\) pour les lois à 1 paramètre (Binomial, Poisson)

• \(D_2-D_1\approx \mbox{Fisher}(q_2,q_1)\) pour les lois à 2 paramètre (Gaussienne)

• Ce test permet de tester la significativité de la diminution de la déviance par l’ajout de variables explicatives ou la prise en compte d’interaction.

Test de Wald

• Ce test est basé sur la forme quadratique faisant intervenir la matrice de covariance des paramètres, l’inverse de la matrice d’information observée \((X'WX)^{-1}\).

• Cette matrice est calculée à partir du Hessien approché par l’algorithme de maximisation. Elle généralise la matrice \((X'X)^{-1}\) 1utilisée dans le cas du modèle linéaire gaussien en faisant intervenir une matrice \(W\) de pondération.

• Le test de Wald et test de Fisher sont équivalents dans le cas particulier du modèle gaussien.

Mise en oeuvre sous R

Estimation du modèle

• On considére le jeu de données suivant

> mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")

• Les données sont des informations sur un groupe d’étudiants concernant leur parcours scolaire (Scores obtenus au cours de leur éducation) en fonction d’une variable binaire (admis ou pas) dans une institution universitaire.

• Résumé des données

> head(mydata)
  admit gre  gpa rank
1     0 380 3.61    3
2     1 660 3.67    3
3     1 800 4.00    1
4     1 640 3.19    4
5     0 520 2.93    4
6     1 760 3.00    2

> summary(mydata)
     admit             gre             gpa             rank      
 Min.   :0.0000   Min.   :220.0   Min.   :2.260   Min.   :1.000  
 1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   1st Qu.:2.000  
 Median :0.0000   Median :580.0   Median :3.395   Median :2.000  
 Mean   :0.3175   Mean   :587.7   Mean   :3.390   Mean   :2.485  
 3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670   3rd Qu.:3.000  
 Max.   :1.0000   Max.   :800.0   Max.   :4.000   Max.   :4.000  

• Tableau de contingence entre les deux variables qualitative admit et rank

> xtabs(~admit + rank, data = mydata)
     rank
admit  1  2  3  4
    0 28 97 93 55
    1 33 54 28 12

• Estimation du modèle logistique

> mydata$rank <- factor(mydata$rank)
> mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

• Résumé du modèle

> summary(mylogit)

Call:
glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

• Dans la sortie R on nous rappelle l’équation du modèle
• Ensuite on donne un résumé des résidus de la déviance.

• Le tableau donne une estimation des paramètres, leur écart-type, la Z-statistique du test de Wald et la p-valeur du test.

• Les estimations des coefficients donnent en effet une estimation de la variation du log odds ratio pour chaque unité de la variable correspondante. Par exemple,
  • Pour chaque changement d’unité dans la variable gre le log odds de l’admission (versus non-admission) croit de 0.0023
  • Pour chaque changement d’unité dans la variable gpa le log odds de l’admission (versus non-admission) croit de 0.804
  • Les indicateurs de la variable rank s’interprètent différement. En effet, ayant ete accepte dans une institution avec un rang 2, versus être accepté avec un rang 1 (la modalité de référence) change le log du odds ratio de -0.6754.
• Rappelons que le log du odds ratio s’exprime \[ \mbox{log du odds ratio}=\log \left(\displaystyle\frac{\pi_1}{1-\pi_1}\right)\] Alors un log du odds ratio
  • égal à 0 les deux groupes ont les mêmes chances
  • positif donc le groupe 1 est plus en faveur que l’autre groupe 2.
  • negatif donc le groupe 2 est plus en faveur que l’autre groupe 1.

• Dans R il y a deux facçons pour calculer les Intervalles de confiances des coefficients.

calcul basé sur le profil \(\log\)-vraisemblance.

> confint(mylogit)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -6.2716202334 -1.792547080
gre          0.0001375921  0.004435874
gpa          0.1602959439  1.464142727
rank2       -1.3008888002 -0.056745722
rank3       -2.0276713127 -0.670372346
rank4       -2.4000265384 -0.753542605

calcul classique

> confint.default(mylogit)
                    2.5 %       97.5 %
(Intercept) -6.2242418514 -1.755716295
gre          0.0001202298  0.004408622
gpa          0.1536836760  1.454391423
rank2       -1.2957512650 -0.055134591
rank3       -2.0169920597 -0.663415773
rank4       -2.3703986294 -0.732528724

• On peut aussi tester l’effect global de la variable rank

> library(aod)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
Wald test:
----------

Chi-squared test:
X2 = 20.9, df = 3, P(> X2) = 0.00011

• On doit préciser dans R, b les coefficients du modèle logit, Sigma la matrice covariance des erreurs et Terms les rangs de ces modélités dans l’objet b.

• On peut de la même facçon tester une fonction des coefficients \(\beta\) : Si \(K\in \mathbb{R}^{q\times p}\) est une matrice donnée. On teste \[H_0\,:\, K\beta=0\] Par exemple, si on pose \(K=(0,0,1,-1,0)\), le test décrit ici est équivalent \[H_0\,:\, \beta_3-\beta_4=0\]

> l <- cbind(0, 0, 0, 1, -1, 0)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 5.5, df = 1, P(> X2) = 0.019

• Estimation des odds ratios

> exp(coef(mylogit))
(Intercept)         gre         gpa       rank2       rank3       rank4 
  0.0185001   1.0022670   2.2345448   0.5089310   0.2617923   0.2119375 

Et les intervalles de confiance

> exp(cbind(OR = coef(mylogit), confint(mylogit)))
Waiting for profiling to be done...
                   OR       2.5 %    97.5 %
(Intercept) 0.0185001 0.001889165 0.1665354
gre         1.0022670 1.000137602 1.0044457
gpa         2.2345448 1.173858216 4.3238349
rank2       0.5089310 0.272289674 0.9448343
rank3       0.2617923 0.131641717 0.5115181
rank4       0.2119375 0.090715546 0.4706961

Prédiction.

• Montrons maintenant comment on peut faire de la prédiction avec le modèle déjà estimé.

• On crée une nouvelle base ou les nouvelles observations concernées par la prédiction : Prédire les probabilités pour les valeurs moyennes des variables quantitatives et les rangs.

> newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
> newdata1
    gre    gpa rank
1 587.7 3.3899    1
2 587.7 3.3899    2
3 587.7 3.3899    3
4 587.7 3.3899    4

• Les probabilités prédictes.

> newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")
> newdata1
    gre    gpa rank     rankP
1 587.7 3.3899    1 0.5166016
2 587.7 3.3899    2 0.3522846
3 587.7 3.3899    3 0.2186120
4 587.7 3.3899    4 0.1846684

• Montrons comment les probablités dans newdata1$rankP sont calculées.

> coef(mylogit)
 (Intercept)          gre          gpa        rank2        rank3 
-3.989979073  0.002264426  0.804037549 -0.675442928 -1.340203916 
       rank4 
-1.551463677 
>  x=sum(coef(mylogit)[1:3]*c(1,newdata1$gre[1],newdata1$gpa[1]))
>  exp(x)/(1+exp(x)) ## Pour le rang 1
[1] 0.5166016
>  x=sum(coef(mylogit)[1:4]*c(1,newdata1$gre[1],newdata1$gpa[1],1))
>  exp(x)/(1+exp(x)) ## Pour le rang 1
[1] 0.3522846
>   x=sum(coef(mylogit)[c(1:3,5)]*c(1,newdata1$gre[1],newdata1$gpa[1],1))
>  exp(x)/(1+exp(x)) ## Pour le rang 1
[1] 0.218612

• On peut aussi prédire les \(\nu_i=x'_i\beta\).

> newdata1$rankNU <- predict(mylogit, newdata = newdata1, type = "link")
> newdata1
    gre    gpa rank     rankP      rankNU
1 587.7 3.3899    1 0.5166016  0.06643085
2 587.7 3.3899    2 0.3522846 -0.60901208
3 587.7 3.3899    3 0.2186120 -1.27377307
4 587.7 3.3899    4 0.1846684 -1.48503283

Ceci se calcule de la facon suivante

> sum(coef(mylogit)[1:3]*c(1,newdata1$gre[1],newdata1$gpa[1]))
[1] 0.06643085
> sum(coef(mylogit)[1:4]*c(1,newdata1$gre[1],newdata1$gpa[1],1))
[1] -0.6090121
> sum(coef(mylogit)[c(1:3,5)]*c(1,newdata1$gre[1],newdata1$gpa[1],1))
[1] -1.273773

• Faisons un calcul de l’erreur de prédiction avec la méthode leave-one-out.

• On écrit une fonction une fonction qui fait de la prédiction.

>       PredictLogit=function(x,y,rk,data){
+         mylogit <- glm(admit ~ gre + gpa + rank, data = data, family = "binomial")
+         newdata1 <- with(data, data.frame(gre = x, gpa = y, rank = factor(rk)))
+         Z=predict(mylogit, newdata = newdata1, type = "response")
+         return(Z)
+       }

• Calcul des prédictions

> n=nrow(mydata)
> df=list()
> for(i in 1:n){
+   dt=mydata[-i,]
+   x=mydata$gre[i]
+   y=mydata$gpa[i]
+   rk=mydata$rank[i]
+   df[[i]]=PredictLogit(x=x,y=y,rk = rk,data = dt)
+ }
> df=do.call(rbind,df)

• Calculons les prédictions sur l’échantillon

> pr=predict(mylogit, newdata = mydata, type = "response")

• Courbe de ROC, PresenceAbsence package

> df1=data.frame(mydata$admit,df[,1],pr)
> colnames(df1)=c("OBSERVED","LOO","SUB")

• Dans tout probléme de prediction on a le tableau suivant : si \(y\) est la variable binaire à predire \(1/0\) ou malade/sain (notre objectif est de prédire les malades.)

• On a alors le tableau suivant

	Malade	Non Malade
Test positif	TP \((y=1,\widehat{y}=1)\)	FP \((y=0,\widehat{y}=1)\)
Test négatif	FN \((y=1,\widehat{y}=0)\)	TN \((y=0,\widehat{y}=0)\)

• Les positifs \(P=|\{\widehat{y}=1|\) et les négatifs \(N=|\{\widehat{y}=0\}|\)

• Donc quand on prédit correctement un malade, on parle de vrais positifs, \(TP=|\{y=1,\widehat{y}=1\}|\) et les faux positifs (erreur de type 1) \(FP=|\{y=0,\widehat{y}=1\}|\)

• Et quand on détecte correctement un sain, on parle de vrais négatifs, \(TN=|\{y=0,\widehat{y}=0\}|\) et les faux négatifs (erreur de type 2) \(FN=|\{y=1,\widehat{y}=0\}|\)

• On mesure les statistiques suivantes

  • Sensitivité ou le taux de vrais positifs \[\mbox{Sensitivité}=\mbox{TPR}=\displaystyle\frac{TP}{TP+FN}\]
  • Spécificité ou le taux de vrais négatifs \[\mbox{Spécificité}=\mbox{TNR}=\displaystyle\frac{TN}{TN+FP}\]

• Calcul sous R

> library(PresenceAbsence)
> ## Il faut mettre les donnees sous le format lisible par le 
> ## package PresenceAbsence
> df1$ID=rep(1,nrow(df1))
> df1=df1[,c(4,1,2,3)]
> presence.absence.accuracy(DATA = df1)
  model threshold    PCC sensitivity specificity     Kappa       AUC
1   LOO       0.5 0.7025   0.2283465   0.9230769 0.1807229 0.6654841
2   SUB       0.5 0.7100   0.2362205   0.9304029 0.1993650 0.6928413
      PCC.sd sensitivity.sd specificity.sd   Kappa.sd     AUC.sd
1 0.02288654     0.03739582     0.01615708 0.04739614 0.02936373
2 0.02271652     0.03784056     0.01542931 0.04745029 0.02829281

• Voyons certains calculs ?

> y=mydata$admit
> yhat.loo=1*(df1$LOO>.5)
> tt=xtabs(~y+yhat.loo)
> tt
   yhat.loo
y     0   1
  0 252  21
  1  98  29

• Donc

> P=sum(tt[2,]);P
[1] 127
> N=sum(tt[1,]);N
[1] 273
> sensit=tt[2,2]/P;sensit
[1] 0.2283465
> speci=tt[1,1]/N;speci
[1] 0.9230769

• On calcule aussi le
  • PCC pourcentage des correctement classés ou prédits
  • Kappa
  • AUC l’aire sous la courbe de ROC. C’est une valeur comprise entre 0.5 et 1.
  • FPR (false postive rate) = erreur de type 1= 1- spécificité
  • FNR (false negtive rate) = erreur de type 2= 1- sensitivité
  • Puissance (Power) = sensitivité
• Calcul du PCC

> tt
   yhat.loo
y     0   1
  0 252  21
  1  98  29
> pcc(tt)
     PCC     PCC.sd
1 0.7025 0.02288654
> (tt[1,1]+tt[2,2])/sum(tt)
[1] 0.7025

• On change le seuil threshold

> ## Pour SUB
> presence.absence.accuracy(DATA = df1,which.model = 2,threshold = 10)
   model threshold    PCC sensitivity specificity      Kappa       AUC
1    SUB 0.0000000 0.3175  1.00000000  0.00000000 0.00000000 0.6928413
2    SUB 0.1111111 0.3500  1.00000000  0.04761905 0.03077296 0.6928413
3    SUB 0.2222222 0.5175  0.81102362  0.38095238 0.14707442 0.6928413
4    SUB 0.3333333 0.6425  0.59842520  0.66300366 0.24065421 0.6928413
5    SUB 0.4444444 0.6950  0.32283465  0.86813187 0.21343606 0.6928413
6    SUB 0.5555556 0.7050  0.17322835  0.95238095 0.15578608 0.6928413
7    SUB 0.6666667 0.6900  0.05511811  0.98534799 0.05354349 0.6928413
8    SUB 0.7777778 0.6825  0.00000000  1.00000000 0.00000000 0.6928413
9    SUB 0.8888889 0.6825  0.00000000  1.00000000 0.00000000 0.6928413
10   SUB 1.0000000 0.6825  0.00000000  1.00000000 0.00000000 0.6928413
       PCC.sd sensitivity.sd specificity.sd    Kappa.sd     AUC.sd
1  0.02330434     0.00000000    0.000000000 0.000000000 0.02829281
2  0.02387835     0.00000000    0.012912527 0.008753977 0.02829281
3  0.02501597     0.03487669    0.029445092 0.036146840 0.02829281
4  0.02399320     0.04367198    0.028660654 0.048407426 0.02829281
5  0.02304920     0.04165356    0.020515321 0.050717125 0.02829281
6  0.02283069     0.03371452    0.012912527 0.043647108 0.02829281
7  0.02315362     0.02033062    0.007285495 0.028254764 0.02829281
8  0.02330434     0.00000000    0.000000000 0.000000000 0.02829281
9  0.02330434     0.00000000    0.000000000 0.000000000 0.02829281
10 0.02330434     0.00000000    0.000000000 0.000000000 0.02829281
> ## Pour le leave-one-out
> presence.absence.accuracy(DATA = df1,which.model = 1,threshold = 10)
   model threshold    PCC sensitivity specificity      Kappa       AUC
1    LOO 0.0000000 0.3175  1.00000000  0.00000000 0.00000000 0.6654841
2    LOO 0.1111111 0.3450  0.98425197  0.04761905 0.02065228 0.6654841
3    LOO 0.2222222 0.5050  0.80314961  0.36630037 0.12919186 0.6654841
4    LOO 0.3333333 0.6325  0.59055118  0.65201465 0.22242793 0.6654841
5    LOO 0.4444444 0.6875  0.30708661  0.86446886 0.19219336 0.6654841
6    LOO 0.5555556 0.7000  0.15748031  0.95238095 0.13697005 0.6654841
7    LOO 0.6666667 0.6875  0.05511811  0.98168498 0.04856142 0.6654841
8    LOO 0.7777778 0.6825  0.00000000  1.00000000 0.00000000 0.6654841
9    LOO 0.8888889 0.6825  0.00000000  1.00000000 0.00000000 0.6654841
10   LOO 1.0000000 0.6825  0.00000000  1.00000000 0.00000000 0.6654841
       PCC.sd sensitivity.sd specificity.sd   Kappa.sd     AUC.sd
1  0.02330434     0.00000000    0.000000000 0.00000000 0.02936373
2  0.02379818     0.01109126    0.012912527 0.01114689 0.02936373
3  0.02503006     0.03542265    0.029212986 0.03591733 0.02936373
4  0.02413640     0.04380698    0.028881830 0.04839847 0.02936373
5  0.02320464     0.04109459    0.020754380 0.05056271 0.02936373
6  0.02294157     0.03245024    0.012912527 0.04259904 0.02936373
7  0.02320464     0.02033062    0.008130276 0.02861126 0.02936373
8  0.02330434     0.00000000    0.000000000 0.00000000 0.02936373
9  0.02330434     0.00000000    0.000000000 0.00000000 0.02936373
10 0.02330434     0.00000000    0.000000000 0.00000000 0.02936373

• Tracer les sensitivité, spécificité et Kappa

> par(oma = c(0, 5, 0, 0), mar = c(3, 3, 3, 1), mfrow = c(1, 2),
+    cex = 0.7, cex.lab = 1.4, mgp = c(2, 0.5, 0))
> error.threshold.plot(df1, which.model = 2, color = TRUE,
+     add.legend = T, legend.cex = .8)
> error.threshold.plot(df1, which.model = 1, color = TRUE,
+     add.legend = F)

• La courbe de ROC

>  auc.roc.plot(df1, color = TRUE, legend.cex = 1.4, main = "")

• Calcul du seuil optimal: dans le package PresenceAbsence on a plusieurs critères pour calculer le seuil optimal

  • ID=1, par défaut, seuil=0.5
  • ID=2, Sens=Spec, sensitivité=spécificité
  • ID=3, MasSens+Spec, le seuil qui maximise la somme sensitivité+spécificité.
  • ID=4, MaxKappaa, le seuil qui maximise le Kappa
  • ID=5, MaxPACC, le seuil qui maximise le PCC
  • ID=6, PredPrev=Obs, La prévalence à partir de \(\widehat{y}\) égale à l’observée, \(|\{y=1\}|=|\{\widehat{y}=1\}|\).
  • ID=7, ObsPrev, le seuil est égal la prévalence à partir de \(\widehat{y}\).
  • ID=8, MeanProb, le seuil est égal à la moyenne des probabilités prédites.
  • ID=9, MinROCdist, le seuil qui rend la courbe de ROC très proche du point (0,1)
  • ID=10, ReqSens, Le plus grand seuil où la sensitivité est proche de la sensitivité spécifié par l’utilisateur.
  • ID=11, ReqSpec, même chose qu’un ID=9 mais pour la spécificité
  • ID=12, Cost, le seuil base sur un coût relatif spécifié par l’utilisateur.

> optimal.thresholds(df1,opt.methods = 1:12,req.sens = .9,req.spec = .9)
         Method       LOO    SUB
1       Default 0.5000000 0.5000
2     Sens=Spec 0.3200000 0.3200
3  MaxSens+Spec 0.3600000 0.3600
4      MaxKappa 0.3600000 0.3600
5        MaxPCC 0.4800000 0.4900
6  PredPrev=Obs 0.3700000 0.3700
7       ObsPrev 0.3175000 0.3175
8      MeanProb 0.3174214 0.3175
9    MinROCdist 0.3200000 0.3100
10      ReqSens 0.1700000 0.1800
11      ReqSpec 0.4800000 0.4700
12         Cost 0.4800000 0.4900
> mean(df1$LOO)
[1] 0.3174214
> mean(df1$SUB)
[1] 0.3175

• Représentation des seuils optimaux sur la courbe de ROC.

>  auc.roc.plot(df1, color = TRUE, legend.cex = 1.4, main = "",opt.thresholds = T,opt.methods = 1:4)