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.
