Une courte histoire des expériences randomisées

On a beaucoup parlé ces derniers mois des expériences randomisées. Cet automne pour commencer, lorsque le prix Nobel d’économie a été attribué à Esther Duflo, Michael Kremer et Abhijit Banerjee. Plus récemment, lors de la polémique, en France, où, dans une conférence donnée le 13 février 2020, intitulée contre la méthode, un médecin marseillais affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ».

Les expériences randomisées

Quand on se demande l’effet d’un traitement (par rapport à un autre), la bonne méthode consiste à se poser la question « que ce serait-il passé si la personne avait choisi l’autre option ? » – ce que John Stuart Mill parlait de la « méthode des différences ». On parle ici d’inférence causale car on ne peut pas savoir ce qu’aurait donné l’autre option, qui n’a pas été retenue. L’effet causal d’une option, ou d’un traitement, n’est jamais observé car il s’agit de la différence potentielle entre deux alternatives. La méthode la plus simple consiste à prendre deux patients proches (à défaut d’être identiques), et de voir comment les deux répondent aux deux options, la proximité se mesurant en fonction de covariables (comme l’âge, le genre, etc), et les patients n’ayant pas choisi les options proposées. Il y a trois points importants dans cette approche : (1) la notion de groupe de contrôle, autrement dit il faut qu’il y ait des personnes « non traitées » (2) la notion d’aveuglement, ou d’absence d’information, autrement dit, les patients ne doivent pas savoir quel traitement ils ont eu (3) la notion d’affectation aléatoire, autrement dit, le traitement doit être choisi de manière totalement indépendante des patients. Ce sont des points qui opposent fondamentalement les expériences randomisées aux données observationnelles, comme les oppose par exemple Paul Rosenbaum dans Observation and Experiment. Si la méthode d’expérience randomisée a été formalisée par Ronald Fisher en 1934, nous allons revenir ici sur deux moments importants dans l’histoire de cette méthodologie. La première est l’utilisation de la randomisation pour faire des tests autour de 1890, en particulier avec les travaux de Peirce puis Richet, sur le spiritualisme et la transmission de pensée (qui ont très probablement influencé Fisher). Et la seconde est l’expérience qui a été menée en 1954 aux États-Unis, pour tester l’efficacité d’un vaccin contre la polio (et qui montre au passage que faire des expériences randomisées dans un contexte de maladies infectieuses a tout son sens).

Continue reading Une courte histoire des expériences randomisées

STT5100, quiz (régression de Poisson #1)

Pour la fin de la semaine, j’avais donné un petit quiz, basé sur la base suivante, qui donne le nombre de cyclistes à une intersection, à New York City,

download.file("http://freakonometrics.free.fr/NYCVelo.RData","velo.RData")
load("velo.RData")
str(base)
'data.frame':	214 obs. of  7 variables:
 $ Date    : chr  "1-Apr-17" "2-Apr-17" "3-Apr-17" "4-Apr-17" ...
 $ HIGH_T  : num  46 62.1 63 51.1 63 48.9 48 55.9 66 73.9 ...
 $ LOW_T   : num  37 41 50 46 46 41 43 39.9 45 55 ...
 $ PRECIP  : num  0 0 0.03 1.18 0 0.73 0.01 0 0 0 ...
 $ BB_COUNT: int  606 2021 2470 723 2807 461 1222 1674 2375 3324 ...
 $ DAY     : chr  "Sam" "Dim" "Lun" "Mar" ...
 $ DIFF_T  : num  9 21.1 13 5.1 17 7.9 5 16 21 18.9 ......

En utiisant une régression de Poisson, il fallait prédire combien de cyclistes passeront un dimanche, s’il fait une température maximale de 85F, minimale de 70F, et s’il ne pleut pas. Et il fallait voir ce que donnerait la prévision pour un lundi.

newbase = data.frame(DAY=as.factor(c("Lun","Dim")),
 HIGH_T=c(85,85),LOW_T=c(70,70),
 PRECIP=c(0,0))

Faisons un modèle avec toutes les variables explicatives.

reg0 = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+DAY,data=base,family=poisson)

Rajoutons aussi une indicatrice pour indiquer les jours où il ne pleut pas,

reg = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+I(PRECIP==0)+DAY,data=base,family=poisson)
summary(reg) 
 
Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         6.8844970  0.0110463 623.241   <2e-16 ***
HIGH_T              0.0210950  0.0003133  67.328   <2e-16 ***
LOW_T              -0.0114006  0.0003351 -34.024   <2e-16 ***
PRECIP             -0.6570450  0.0071899 -91.384   <2e-16 ***
I(PRECIP == 0)TRUE  0.1303908  0.0033283  39.176   <2e-16 ***
DAYJeu              0.1683475  0.0049690  33.880   <2e-16 ***
DAYLun              0.1144129  0.0050480  22.665   <2e-16 ***
DAYMar              0.1655886  0.0049936  33.160   <2e-16 ***
DAYMer              0.1688035  0.0049190  34.317   <2e-16 ***
DAYSam              0.0466447  0.0051838   8.998   <2e-16 ***
DAYVen              0.1003536  0.0050919  19.709   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for poisson family taken to be 1)
 
    Null deviance: 70021  on 213  degrees of freedom
Residual deviance: 26493  on 203  degrees of freedom
AIC: 28580
 
Number of Fisher Scoring iterations: 4

Tout semble vraiment significatif. Et même avoir du sens (si on regarde les signes, par exemple). Et si on a peur de rater un effet non-linéaire, on peut mettre des splines sur toutes les variables continues

library(gam)
reggam = gam(BB_COUNT~bs(HIGH_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson)
plot(reggam, se=TRUE)

pour la température maximale, ou pour la température minimale

et la courbe suivante pour la précipitation, avec un lissage linéaire entre l’observation maximale (3) et celle juste avant (vers 1.8)

On peut aussi régresser sur la température minimale, et l’écart entre le maximum et le minimum (dans un modèle linéaire les modèles sont équivalents, mais avec de transformation nonlinéaire, la différence pourrait donner un modèle plus simple)

library(gam)
reggam2 = gam(BB_COUNT~bs(DIFF_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson)
plot(reggam2, se=TRUE)

On peut maintenant comparer les quatre modèles, et leurs prévision. Par exemple, pour le modèle linéaire (avec la variable indicatrice indiquant qu’il ne pleut pas),

P = predict(reg,newdata=newbase,type="response",se.fit=TRUE)

pour le lundi, on obtient l’intervalle de confiance à 95% pour \widehat{\lambda}

P$fit[1]+c(-2,2)*P$se.fit[1]
[1] 3349.842 3401.395

et pour le dimanche l’intervalle de confiance à 95% est

P$fit[2]+c(-2,2)*P$se.fit[2]
[1] 2987.497 3033.861

On peut visualiser les quatre intervalles de confiance (avec, de haut en bas, le second modèle gam, le premier, le modèle linéaire avec l’indicatrice d’absence de pluie, puis le premier modèle linéaire)

alors que pour le dimanche, on a

Autrement dit, en changeant de modèle, on change beaucoup les intervalles de confiance sur la prévision (et les intervalles sont parfois complètement disjoints). Ce qui n’est pas forcément une bonne nouvelle.