Category Archives: STT5100

GLM, STT5100

Dernière ligne droite dans le cours STT5100 de modèles linéaires appliqués. Les supports de cours sont en ligne sur https://github.com/freakonometrics/STT5100. Cette session étant en distanciel, le cours est asynchrone, et je poste régulièrement des capsules vidéos. Les capsules en lien avec les modèles linéaires généralisés (GLM) sont maintenant en ligne,

  1. introduction générale video pdf (10:15)
  2. lois de Bernoulli, binomiale, multnomiale video pdf (29:23)
  3. régression logistique (Bernoulli) video pdf (23:04)
  4. régression multinomiale video pdf (21:45)
  5. régression logistique sur variables catégorielles video pdf (30:24)
  6. régression logistique sur variables continues video pdf (21:35)
  7. analyse discriminante et courbe ROC video pdf (56:53)
  8. modèles de comptage et loi de Poisson video pdf (19:28)
  9. régression de Poisson video pdf (25:38)
  10. régression de Poisson et interprétations video pdf (40:15)
  11. régression de Poisson et méthode des marges video pdf (25:29)
  12. régression de Poisson et application en assurance video pdf (36:23)
  13. famille exponentielle video pdf (30:01)
  14. famille exponentielle et GLM video pdf (41:02)
  15. loi et lien video video pdf (36:05)
  16. déviance et résidus video pdf (15:10)
  17. modèle Tweedie et poids video pdf (29:33)
  18. surdispersion video pdf (21:01)
  19. tests et GLM video pdf (19:41)
  20. GLM en petite dimension video pdf (23:34)
  21. méthode stepwise video pdf (17:57)
  22. Poisson vs. Binomiale, application en démographie video pdf (22:25)
  23. Exemple (1) video + pdf

Moindres carrés, STT5100

On avance dans le cours STT5100 de modèles linéaires appliqués. Les supports de cours sont en ligne sur https://github.com/freakonometrics/STT5100. Cette session (mais aussi celle d’hiver) étant en distanciel, le cours est asynchrone, et je poste régulièrement des capsules vidéos. Les capsules en lien avec les moindres carrés sont en ligne

  1. Régressions video + pdf (40:01)
  2. Régression simple sur variable continue (1) video + pdf (41:33)
  3. Régression simple sur variable continue (2) video + pdf (30:26)
  4. Régression simple sur variable continue (3) video + pdf (30:39)
  5. Régression simple sur variable factorielle video + pdf (39:59)
  6. Régression multiple (1) video + pdf (45:56)
  7. Régression multiple (2) video + pdf (37:40)
  8. Modèle Gaussien et tests video + pdf (55:43)
  9. Bootstrap video + pdf (45:37)
  10. Diagnostique video + pdf (50:08)
  11. ANOVA (1) video + pdf (37:52)
  12. ANOVA (2) video + pdf (28:28)
  13. ANOVA (3) video + pdf (16:42)
  14. Régression pondérée video + pdf (16:33)
  15. Multicolinéarité (1) video + pdf (43:25)
  16. Sélection de variables video + pdf (44:13)
  17. MCG & hétéroscédasticité video + pdf (21:30)
  18. Multicolinéarité (2) video + pdf (25:05)
  19. ANOVA (4) pdf (hors programme)
  20. Transformations video + pdf (13:20)
  21. Lissage, non-linéarités video + pdf (34:51)
  22. Discontinuités video + pdf (21:36)
  23. Exemple video + pdf (40:24)
  24. Exemple video + pdf (28:17)

STT5100, introduction et rappels

Reprise des cours ce vendredi pour STT5100, modèles linéaires appliqués. Les supports de cours sont en ligne sur https://github.com/freakonometrics/STT5100. Avec la session en distanciel, le cours sera asynchrone, et je posterais régulièrement des capsules vidéos. Les premières, présentant le cours, et reprenant les principaux concepts déjà vus dans les prérequis, sont en ligne

1. Syllabus video pdf (36:44)
2. Epistémologie video pdf (29:50)

1. Rmarkdown video pdf (30:48)
2. Probabilités video pdf (47:24)
3. Loi normale et vecteur Gaussien video pdf (36:05)
4. Maximum de vraisemblance video pdf (39:56)
5. Estimation paramétrique video pdf (27:32)
6. Estimation nonparamétrique video pdf (37:12)
7. Simulations, Bootstrap video pdf (40:50)
8. Tests (1) video pdf (32:29)
9. Tests (2) video pdf (32:30)
10. Tests (3) video pdf (25:44)
11. Optimisation video pdf (49:57)
12. Algèbre linéaire et matrices video pdf (53:54)

 

Reprise des cours pour l’automne

Les cours de la session d’automne reprennent cette semaine, avec ACT6100 (analyse de données en assurance) jeudi et STT5100 (modèles linéaires appliqués). Les slides sont en ligne sur github, ainsi que les liens vers les vidéos (ces derniers seront mis à jour au fur et à mesure…)

Suite aux soucis sur viméo l’hiver dernier, toutes les vidéos sont maintenant sur youtube, avec une chaîne pour STT5100, et une pour ACT6100.

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.

On the “correlation” between a continuous and a categorical variable

Let us get back on the Titanic dataset,

loc_fichier = "http://freakonometrics.free.fr/titanic.RData"
download.file(loc_fichier, "titanic.RData")
load("titanic.RData")
base = base[!is.na(base$Age),]

On consider two variables, the age x (the continuous one) and the survivor indicator y (the qualitative one)

X = base$Age
Y = base$Survived

It looks like the age might be a valid explanatory variable in the logistic regression,

summary(glm(Survived~Age,data=base,family=binomial))
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.05672    0.17358  -0.327   0.7438  
Age         -0.01096    0.00533  -2.057   0.0397 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 960.23  on 712  degrees of freedom
AIC: 964.23

The significance test here has a p-value just below 4%. Actually, one can relate it with the value of the deviance (the null deviance and the residual deviance). Recall thatD=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\widehat{\boldsymbol{\mu}})\big)whileD_0=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\overline{y})\big)Under the assumption that x is worthless, D_0-D tends to a \chi^2 distribution with 1 degree of freedom. And we can compute the p-value dof that likelihood ratio test,

1-pchisq(964.52-960.23,1)
[1] 0.03833717

(which is consistent with a Gaussian test). But if we consider a nonlinear transformation

summary(glm(Survived~bs(Age),data=base,family=binomial))
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8648     0.3460   2.500 0.012433 *  
bs(Age)1     -3.6772     1.0458  -3.516 0.000438 ***
bs(Age)2      1.7430     1.1068   1.575 0.115299    
bs(Age)3     -3.9251     1.4544  -2.699 0.006961 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 948.69  on 710  degrees of freedom

which seems to be “more significant”

1-pchisq(964.52-948.69,3)
[1] 0.001228712

So it looks like the variable x is interesting here.

To visualize the non-null correlation, one can consider the condition distribution of x given y=1, and compare it with the condition distribution of x given y=0,

ks.test(X[Y==0],X[Y==1])
 
	Two-sample Kolmogorov-Smirnov test
 
data:  X[Y == 0] and X[Y == 1]
D = 0.088777, p-value = 0.1324
alternative hypothesis: two-sided

i.e. with a p-value above 10%, the two distributions are not significatly different.

F0 = function(x) mean(X[Y==0]<=x)
F1 = function(x) mean(X[Y==1]<=x)
vx = seq(0,80,by=.1)
vy0 = Vectorize(F0)(vx)
vy1 = Vectorize(F1)(vx)
plot(vx,vy0,col="red",type="s")
lines(vx,vy1,col="blue",type="s")

(we can also look at the density, but it looks like that there is not much to see)

An alternative is discretize variable x and to use Pearson’s independence test,

k=5
LV = quantile(X,(0:k)/k)
LV[1] = 0
Xc = cut(X,LV)
table(Xc,Y)
           Y
Xc           0  1
  (0,19]    85 79
  (19,25]   92 45
  (25,31.8] 77 50
  (31.8,41] 81 63
  (41,80]   89 53
chisq.test(table(Xc,Y))
 
	Pearson's Chi-squared test
 
data:  table(Xc, Y)
X-squared = 8.6155, df = 4, p-value = 0.07146

The p-value is here 7%, with five categories for the age. And actually, we can compare the p-value

pvalue = function(k=5){
LV = quantile(X,(0:k)/k)
LV[1] = 0
Xc = cut(X,LV)
chisq.test(table(Xc,Y))$p.value}
vk = 2:20
vp = Vectorize(pvalue)(vk)
plot(vk,vp,type="l")
abline(h=.05,col="red",lty=2)

which gives a p-value close to 5%, as soon as we have enough categories. In the slides of the course (STT5100), I claim that actually, the age is an important variable when trying to predict if a passenger survived. Test mentioned here are not as conclusive, nevertheless…