Suite de l’article Le Juge est-il raciste ? Si vous ne l’avez pas lu, je vous conseille de le faire avant de regarder ce qui suit.
Dans l’épisode précédent, nous avons utilisé un modèle logit pour estimer le modèle suivant :
\begin{aligned}
Pr(prison=1) = f(\beta_0 &+ \beta_1 \times afroamericain\\
&+\beta_2 \times unSeulParent\\
&+\beta_3 \times moinsDe14Ans)
\end{aligned}Le coefficient associé à la variable “afroamericain” était statistiquement significatif ; autrement dit la vraie valeur de ce coefficient avait peu de chance d’être égale à zéro. Nous en avons conclu que le juge était probablement raciste.
Mais si Alicia — l’héroïne-avocate de notre histoire — présentait cet argument devant un tribunal, remporterait-elle le procès ? Ce n’est pas sûr…
Statistiquement significatif ?
Pour déterminer la significativité de ce coefficient, nous avons utilisé un test de Student. On commence par calculer la statistique suivante :
t =\dfrac{\text{coefficient estime}}{\text{ecart-type estimé du coefficient}}Sous certaines conditions, la loi de cette statistique est connue et l’on peut alors déterminer la probabilité d’observer une valeur aussi éloignée de 0 si la vraie valeur du coefficient était nulle. Si cette probabilité est suffisamment faible, on conclut que le coefficient est statistiquement significatif.
Le problème est que pour que ce test soit valide, il faut que l’estimateur utilisé suive une loi normale et que l’estimation de sa variance ne soit pas biaisée. Dans le cas du modèle logit, ces deux conditions tendent à être vérifiées lorsqu’on a un très grand nombre de données. Sauf qu’ici on n’a que 16 observations ; par conséquent on ne peut avoir aucune certitude sur la validité de ces deux hypothèses et l’on n’a pas trop d’idée de comment se comporte l’estimateur.
Bootstrap et simulation
Lorsqu’on ne connait pas le comportement d’un estimateur, une bonne pratique consiste à utiliser la méthode du bootstrap ou à réaliser une simulation.
Le bootstrap consiste à simuler de nouveaux échantillons en tirant aléatoirement avec remise des observations dans l’échantillon d’origine. La simulation consiste quant à elle à créer de nouveaux échantillons à partir du modèle qu’on a estimé. Dans les deux cas, pour chaque nouvel échantillon on recalcule le coefficient d’intérêt et l’on obtient ainsi une estimation non paramétrique de sa distribution.
Dans notre cas, le bootstrap donne cette distribution :
Distribution des estimations du coefficient « afroamericain » réalisées sur 2000 échantillons obtenus par bootstrap. On est loin d’une distribution gaussienne !
et la simulation cette distribution-là :
Où est le problème ?
Le modèle logit est estimé par maximum de vraisemblance. C’est-à-dire que les coefficients sont choisis de manière à maximiser la fonction de vraisemblance :
\mathcal{L}(\beta)=\prod_{i}Pr(prison_{i}|X_{i},\beta)où \beta représente le vecteur des coefficients à estimer et X_i les variables explicatives pour l’individu i.
Un premier problème est qu’il n’existe pas de méthode permettant de trouver à coup sûr le maximum de cette fonction. Les algorithmes utilisés sont généralement itératifs et s’arrêtent au bout d’un certain nombre d’itérations, qu’ils aient convergé ou non sur un maximum local. Néanmoins, dans les simulations précédentes, ce problème de convergence concerne moins de 2 % des échantillons simulés.
Le principal problème vient du fait que pour trouver le maximum de la fonction de vraisemblance, il faut que cette fonction ait un maximum ! ce qui n’est pas toujours le cas. Supposons que dans un échantillon, tous les Afro-Américains aient été condamnés. Dans ce cas, plus le coefficient associé à la variable “afroaméricain” est élevé, plus la fonction de vraisemblance augmente : on ne peut donc jamais en atteindre le maximum. De manière générale, si la variable dépendante est toujours égale à 0 ou à 1 lorsqu’une variable explicative dichotomique est égale à 1 (ou à 0), alors la fonction de vraisemblance n’a pas de maximum.
Dans l’échantillon d’origine, il n’y a qu’un Afro-Américain qui n’a pas été condamné et un seul non Afro-Américain qui a été condamné. Donc, lorsqu’on fait du bootstrap ou de la simulation, la probabilité que ce problème apparaisse est très élevée.
Que faire ?
On ne peut pas utiliser le test de Student et en plus le bootstrap et les simulations donnent des résultats aberrants inutilisables pour construire un intervalle de confiance. Que peut-on faire ? La seule solution que je vois est d’essayer de changer de modèle. Un modèle probit a exactement les mêmes problèmes que le modèle logit, donc ce n’est pas la peine d’y songer. En revanche, ce n’est pas le cas d’un modèle linéaire.
Ce modèle revient à estimer l’équation suivante :
\begin{aligned}
Pr(prison)=\beta_0&+\beta_1 \times afroamericain\\
&+\beta_2 \times unSeulParent \\
&+\beta_3\times moinsDe14Ans+\varepsilon
\end{aligned}Sur l’échantillon initial, les probabilités estimées par ce modèle sont quasi identiques à celles obtenues avec le modèle logit. En revanche, lorsqu’on rééchantillonne ou qu’on simule des échantillons, les coefficients estimés sont beaucoup plus stables, comme on peut le voir sur les graphiques suivants.
À partir de ces résultats, on peut construire des intervalles de confiance et déterminer si le coefficient de la variable “afroaméricain” est significatif ou non. D’après les résultats du bootstrap ce coefficient a plus de 90 % de chances d’être non nul et d’après les simulations, cette probabilité est supérieure à 95 %. Il y a donc de fortes chances pour que le juge soit raciste.
Néanmoins, ces résultats sont eux aussi contestables. Lorsqu’on cherche à estimer une probabilité, le modèle linéaire est peu adapté car il conduit souvent à des résultats aberrants : il arrive régulièrement que les probabilités estimées soient supérieures à 1 ou inférieures à 0. Dans le cas présent, 40 % des échantillons créés par bootstrap ont conduit à des probabilités estimées supérieures à 1. On peut donc sérieusement s’interroger sur la validité de cette méthode.
Morale de l’histoire
Lorsqu’on cherche à identifier les causes de phénomènes, il est rare que les statistiques permettent de conclure de manière certaine. Elles sont rarement des preuves absolues, mais plutôt des indices qui renforcent ou affaiblissent nos croyances.
Dans le cas présent, les données ne démontrent pas clairement que le juge est raciste, mais elles ne permettent pas non plus de conclure qu’il ne l’est pas. Elles doivent inciter à rester vigilant et à rechercher d’autres indices qui confirment que le juge est ou n’est pas raciste.
Le Code R
Le code suivant permet de générer les données utilisées :
data <- data.frame(id = 1:23,
jail = c(rep(TRUE, 12), rep(FALSE, 11)),
afroAmerican = c(1,1,1,1,1,1,1,0,1,1,1,1,0,0,0,
0,0,0,0,1,0,0,1),
singleParent = c(0,NA,1,1,NA,1,1,1,NA,0,1,NA,0,
1,0,0,0,NA,0,1,0,NA,NA),
underFourteen = c(1,NA,0,0,NA,1,0,0,NA,0,0,NA,
0,1,1,1,1,NA,1,0,1,NA,NA)
)
Pour utiliser le bootstrap dans R il faut installer la librairie bootstrap. Si vous ne l’avez pas, vous pouvez l’installer avec la ligne de code suivante :
install.packages("bootstrap")
Le code suivant génère 2000 échantillons par bootstrap ou par simulation et estime le modèle sur chacun d’eux. Pour utiliser un modèle linéaire, il suffit de supprimer l’argument « family=binomial(« logit ») ».
library(bootstrap)
model <- jail~ afroAmerican+singleParent+underFourteen
reg <- lm(model, data=data, family=binomial("logit"))
# Bootstrap
#----------
bt <- bootstrap(1:16, 2000, function(x) {
dt <- reg$model[x,]
reg <- lm(model, data=dt, family=binomial("logit"))
c(coef(reg), reg$converged)
})
# Simulation
#-----------
fit <- fitted(reg)
bt <- bootstrap(1:16, 2000, function(x) {
dt <- reg$model
dt$jail <- rbinom(nrow(dt), 1, fit)
reg <- lm(model, data=dt, family=binomial("logit"))
c(coef(reg), reg$converged)
})
Les résultats sont contenus dans « bt$thetastar ». Il s’agit d’une matrice dont les lignes représentent les coefficients estimés et les colonnes les échantillons. La dernière ligne de ce tableau indique si l’algorithme d’optimisation a convergé ou non. Pour calculer un intervalle de confiance, on peut utiliser la fonction « quantile » :
quantile(bt$thetastar[2,] , c(0.025, 0.975))
Pour visualiser les résultats sous forme d’histogramme, on utilise la fonction « hist » :
hist(bt$thetastar[2,], 30)