Tag Archives: ACT6420

Google, et prévision

Pour le second devoir du cours ACT6420 (modèles de prévisions), le but est de prévoir des recherches sur Google, via https://www.google.com/trends/. Soit vous avez un mot clé qui vous intéresse, soit vous allez chercher une série extraite dans /ACT6420-TS/. Dans les bases qui ont été mises en ligne, il y a le mot clé gym, par exemple.

> report=read.table(
+ "http://freakonometrics.free.fr/ACT6420-TS/report-GYM.csv",
+ skip=4,header=TRUE,sep=",",nrows=548)

On va nettoyer un peu la base (en particulier les valeurs manquantes de la fin)

> tail(report)
                    Semaine gym
543 2014-05-25 - 2014-05-31  80
544 2014-06-01 - 2014-06-07  80
545 2014-06-08 - 2014-06-14  78
546 2014-06-15 - 2014-06-21  NA
547 2014-06-22 - 2014-06-28  NA
548 2014-06-29 - 2014-07-05  NA
> report=report[!is.na(report[,2]),]
> tail(report)
                    Semaine gym
540 2014-05-04 - 2014-05-10  79
541 2014-05-11 - 2014-05-17  80
542 2014-05-18 - 2014-05-24  79
543 2014-05-25 - 2014-05-31  80
544 2014-06-01 - 2014-06-07  80
545 2014-06-08 - 2014-06-14  78

Les données sont ici hebdomadaires, comme on peut le voir graphiquement

> hebdo=ts(report$gym,start=2004,frequency=52)
> hebdo
Time Series:
Start = c(2004, 1) 
End = c(2014, 25) 
Frequency = 52 
  [1]  68  60  60  53  55  50  49  49  51  48  48  45  45
 [14]  47  42  48  46  47  46  47  47  46  47  46  45  46
 [27]  46  50  48  48  52  51  57  55  53  56  55  50  48
 [40]  50  46  49  46  48  49  48  46  50  47  46  43  54
 [53]  69  64  63  60  57  57  53  54  55  54  50  53  54
 [66]  46  50  50  49  49  49  47  49  48  49  50  49  49
 [79]  47  47  50  49  52  51  55  55  55  54  54  52  53
 [92]  54  53  52  51  51  50  51  48  52  50  47  45  56
[105]  76  72  66  64  63  59  53  56  57  58  54  55  54
[118]  53  52  52  50  53  50  51  49  51  51  50  50  48
[131]  48  53  52  56  58  57  60  62  62  63  59  58  58
[144]  56  54  54  53  53  54  53  50  55  54  53  51  56
[157]  77  73  68  68  67  66  64  67  64  63  63  63  62
[170]  62  61  62  61  63  62  63  63  63  63  59  58  59
[183]  61  60  60  58  61  61  62  62  64  68  66  61  58
[196]  58  55  54  51  55  54  55  53  55  53  52  50  55
[209]  76  77  68  67  64  64  58  61  59  59  57  55  57
[222]  59  57  54  56  55  54  52  52  53  54  53  55  55
[235]  55  57  54  56  55  58  65  63  64  67  66  63  62
[248]  60  60  57  55  56  56  58  58  53  56  55  54  52
[261]  69  77  71  68  66  63  60  60  62  59  59  57  57
[274]  60  58  58  56  54  58  57  56  57  57  57  57  54
[287]  54  55  57  57  56  64  60  59  62  62  64  59  58
[300]  57  57  54  53  52  53  53  55  53  55  53  50  49
[313]  63  76  73  68  65  66  60  61  61  60  58  58  61
[326]  61  62  57  57  58  55  58  57  58  59  57  55  55
[339]  57  57  58  59  59  60  60  63  63  63  66  65  62
[352]  60  59  58  57  56  58  59  56  54  57  54  54  53
[365]  66  87  77  74  72  69  68  64  65  65  68  63  65
[378]  65  65  62  61  62  62  63  63  61  65  63  64  63
[391]  61  62  62  61  62  65  63  67  67  71  74  71  70
[404]  68  68  65  66  64  65  68  64  64  65  62  61  60
[417]  69  91  88  83  81  78  75  71  73  74  73  70  69
[430]  66  68  69  66  68  68  65  69  66  69  70  69  70
[443]  69  72  72  71  69  76  74  72  77  77  82  78  72
[456]  72  69  68  67  67  64  66  66  63  65  64  62  61
[469]  67  88  90  83  81  83  77  76  76  74  75  74  74
[482]  77  77  77  73  72  76  72  71  72  74  72  74  73
[495]  73  76  73  73  71  76  76  79  79  83  83  81  78
[508]  78  76  78  80  74  73  75  74  75  72  71  70  69
[521]  73  92 100  94  93  91  86  84  84  85  85  83  82
[534]  83  83  78  79  80  80  79  80  79  80  80  78

> plot(hebdo)

Pour avoir des modèles plus simples, on peut agréger les données, pour les rendre mensuelles (par interpolation linéaire). La fonction à utiliser est ici

        H2M=function(BASE){
 	X=BASE[,2]
 	Y=BASE[,1]
 	date1=substr(as.character(Y),1,10)
 	date2=substr(as.character(Y),14,23)
 	D1=as.Date(date1,"%Y-%m-%d")
 	D2=as.Date(date2,"%Y-%m-%d")
 	vm=vy=N=NA
 	for(t in 1:length(D1)){
 		mois=seq(D1[t],D2[t],length=7)
 		vm=c(vm,as.POSIXlt(mois)$mon+1)
 		vy=c(vy,as.POSIXlt(mois)$year+1900)
 		N=c(N,rep(X[t],7))}
 	N=N[-1]; vm=vm[-1]; vy=vy[-1]
 	YM=vy*100+vm
 	Z=tapply(N,as.factor(YM),mean)
 	Zts=ts(as.numeric(Z),start=c(2004,1),frequency=12)
 	return(Zts)}

Si on utilise cette fonction sur nos données, on obtient

> mensuel=H2M(report)
> mensuel
          Jan      Feb      Mar      Apr      May      Jun
2004 60.25000 50.75862 47.51613 45.66667 46.67742 46.00000
2005 63.22581 55.10714 53.03226 49.10000 48.45161 49.10000
2006 68.87097 57.10714 55.51613 51.86667 50.70968 49.93333
2007 70.74194 65.57143 62.87097 61.60000 62.77419 59.96667
2008 70.45161 60.79310 57.19355 56.13333 52.96774 54.30000
2009 70.35484 61.25000 58.19355 57.13333 56.80645 56.00000
2010 69.87097 61.78571 59.45161 58.76667 57.16129 56.40000
2011 76.58065 66.21429 65.22581 62.66667 62.51613 63.16667
2012 85.00000 73.82759 69.93548 67.76667 67.25806 69.46667
2013 84.93548 76.39286 75.00000 74.80000 72.67742 73.13333
2014 94.29032 84.96429 83.29032 79.80000 79.54839 79.00000
          Jul      Aug      Sep      Oct      Nov      Dec
2004 47.80645 53.67742 52.63333 47.77419 47.96667 47.61290
2005 48.41935 53.51613 53.43333 52.41935 50.20000 49.74194
2006 52.48387 59.77419 59.66667 54.12903 52.86667 54.35484
2007 59.87097 62.03226 63.10000 54.45161 54.30000 54.09677
2008 55.45161 62.16129 63.96667 57.41935 56.23333 56.09677
2009 55.96774 61.12903 60.16667 54.29032 53.60000 53.35484
2010 58.12903 61.64516 63.43333 57.67742 56.73333 56.48387
2011 61.80645 66.22581 70.86667 65.77419 65.23333 63.19355
2012 71.48387 75.06452 75.80000 67.22581 64.90000 65.12903
2013 73.51613 78.93548 79.73333 76.41935 73.93333 72.80645
2014                                                      
> ts.plot(mensuel)

Cette base est maintenant utilisable pour le devoir. Le but est ici de faire de la prévision pour les 2 prochaines années, avec un intervalle de confiance. Mais j’en reparlerais par courriel.

Prévision de séries chronologiques

Dans la seconde partie du cours de modèles de prévision, on quittera (un peu) les données individuelles pour parler de données chronologiques. Les slides de cette semaine (et probablement la semaine prochaine) sont en ligne.

J’ai mis en ligne, en parallèle, quelques notes de cours sur les séries temporelles, qui pourront peut être servir de complément. Comme le disait Doug MartinTime series is the worst subject to teach. First, you have to teach the standard theory. Then, if you are beging honest, you have to tell the students ‘none of this stuff works, and this is what people really do'” (cité la semaine passée par David J. Thomson, au congrès de la SSC). On essayera de garder ça en mémoire tout au long du cours !

Modèles de prévision, fourre-tout

Dans le dernier cours de modèle de prévision, la semaine passée, nous avions passé un peu de temps sur l’étude des points aberrants et des points influents. Tout est expliqué dans les slides (avec les codes) donc je ne reviendrais pas dessus. Je pourrais juste évoquer quelques lignes de codes utilisées pour voir l’impact d’une observation sur la régression, en enlevant l’observation de la base, et en regardant ce que ça donne sur la prévision

> plot(cars)
> text(cars[c(23,49),1],cars[c(23,49),2]+4,c(23,49))
> abline(lm(dist~speed,data=cars))
> abline(lm(dist~speed,data=cars[-23,]),col="red")
> abline(lm(dist~speed,data=cars[-49,]),col="blue")

On a ensuite vu les problèmes de choix de modèles, et les méthodes de pénalisation (pénalisation du R2, et de la log-vraisemblance). On a ensuite vu les méthodes automatiques de choix de modèle,

> US=read.table("http://freakonometrics.free.fr/US.txt",sep=";")
> reg=lm(Murder~.,data=US)
> regs=step(reg)

Start:  AIC=63.01
Murder ~ Population + Income + Illiteracy + Life.Exp + HS.Grad + 
    Frost + Area

             Df Sum of Sq    RSS    AIC
- Income      1     0.236 128.27 61.105
- HS.Grad     1     0.973 129.01 61.392
<none>                    128.03 63.013
- Area        1     7.514 135.55 63.865
- Illiteracy  1     8.299 136.33 64.154
- Frost       1     9.260 137.29 64.505
- Population  1    25.719 153.75 70.166
- Life.Exp    1   127.175 255.21 95.503

Step:  AIC=61.11
Murder ~ Population + Illiteracy + Life.Exp + HS.Grad + Frost + 
    Area

             Df Sum of Sq    RSS    AIC
- HS.Grad     1     0.763 129.03 59.402
<none>                    128.27 61.105
- Area        1     7.310 135.58 61.877
- Illiteracy  1     8.715 136.98 62.392
- Frost       1     9.345 137.61 62.621
- Population  1    27.142 155.41 68.702
- Life.Exp    1   127.500 255.77 93.613

Step:  AIC=59.4
Murder ~ Population + Illiteracy + Life.Exp + Frost + Area

             Df Sum of Sq    RSS    AIC
<none>                    129.03 59.402
- Illiteracy  1     8.723 137.75 60.672
- Frost       1    11.030 140.06 61.503
- Area        1    15.937 144.97 63.225
- Population  1    26.415 155.45 66.714
- Life.Exp    1   140.391 269.42 94.213

Enfin, on a aussi vu le passage au logarithme sur la variable d’intérêt, avec un modèle log-normal,

> plot(cars$speed,cars$dist)
> regln=lm(log(dist)~speed,data=cars)
> summary(regln)

Call:
lm(formula = log(dist) ~ speed, data = cars)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.46604 -0.20800 -0.01683  0.24080  1.01519 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.67612    0.19614   8.546 3.34e-11 ***
speed        0.12077    0.01206  10.015 2.41e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4463 on 48 degrees of freedom
Multiple R-squared:  0.6763,	Adjusted R-squared:  0.6696 
F-statistic: 100.3 on 1 and 48 DF,  p-value: 2.413e-13

> beta=coefficients(regln)
> for(x in seq(2,30,by=.1)){
+ zhat=beta[1]+beta[2]*x
+ yhat=exp(zhat+1/2*0.4463^2)
+ points(x,exp(zhat),col="red")
+ points(x,yhat,col="blue")
+ }
> abline(lm(dist~speed,data=cars))

Pour comparer les modèles (linéaire versus log-linéaire), on avait suggéré de comparer les résidus des deux modèles,

> U=cars$dist - exp(predict(regln)+1/2*0.4463^2)

> sd(U)
[1] 17.12949
> sd(residuals(lm(dist~speed,data=cars)))
[1] 15.22184

Ce soir, on fini cette première partie du cours, sur les données individuelles.

Régression linéaire, quelques codes

Un rapide billet pour mettre en ligne les codes utilisés la semaine passée, complétant les codes des transparents. On travaille toujours sur la même base, ou on cherche à prévoir une distance de freinage d’un véhicule, tenant compte de la vitesse du véhicule.

> plot(cars)
> reg=lm(dist~speed,data=cars)
> summary(reg)

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.069  -9.525  -2.272   9.215  43.201 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791     6.7584  -2.601   0.0123 *  
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared:  0.6511,	Adjusted R-squared:  0.6438 
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Pour faire plusieurs prévisions, à la main, on peut utiliser le code suivant (la boucle permet de faire des prévisions pour plusieurs valeurs)

> for(x in seq(3,30,by=.25)){
+ b0=coef(reg)[1]
+ b1=coef(reg)[2]
+ Yx=b0+b1*x
+ V=vcov(reg)
+ Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2
+ IC1=Yx+c(-1,+1)*1.96*sqrt(Vx)
+ s=summary(reg)$sigma
+ IC2=Yx+c(-1,+1)*1.96*s
+ points(x,Yx,pch=19,col="red")
+ points(c(x,x),IC1,pch=3,col="blue")
+ points(c(x,x),IC2,pch=3,col="purple")}

On avait ensuite fait une régression linéaire sur une sous-base, avec 20 observations tirées au hasard

> I=sample(1:50,size=20)
> reg=lm(dist~speed,data=cars[I,])

Le but était de visualiser l’impact du nombre d’observation sur la qualité de la régression

> summary(reg)

Call:
lm(formula = dist ~ speed, data = cars[I, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-23.529  -7.998  -5.394  11.634  39.348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -20.7408     9.4639  -2.192   0.0418 *  
speed         4.2247     0.6129   6.893 1.91e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.62 on 18 degrees of freedom
Multiple R-squared:  0.7252,	Adjusted R-squared:   0.71 
F-statistic: 47.51 on 1 and 18 DF,  p-value: 1.91e-06

> for(x in seq(3,30,by=.25)){
+   b0=coef(reg)[1]
+   b1=coef(reg)[2]
+   Yx=b0+b1*x
+   V=vcov(reg)
+   Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2
+   IC=Yx+c(-1,+1)*1.96*sqrt(Vx)
+   points(x,Yx,pch=19,col="purple")
+   points(c(x,x),IC,pch=3,col="green")}

Notons qu’il est possible d’utiliser des fonctions de R pour faire des prévisions, avec des intervalles de confiance

> predict(reg,
+ newdata=data.frame(speed=c(15,25)),interval= "confidence")
       fit      lwr       upr
1 42.62976 34.75450  50.50502
2 84.87677 68.92746 100.82607
> predict(reg,
+ newdata=data.frame(speed=15),interval= "prediction")
       fit      lwr      upr
1 42.62976 6.836077 78.42344

Quand on a plus d’une variable explicative, c’est plus compliqué de “visualiser” la régression

>  chicago=read.table("http://freakonometrics.free.fr/chicago.txt",
+  header=TRUE,sep=";")
>  Y=chicago$Fire
>  X1=chicago$X_1
>  X2=chicago$X_2
>  X3=chicago$X_3
>  base=data.frame(Y,X1,X2,X3)
> plot(X2,X3)
> reg=lm(Y~X2+X3,data=base)
> y=function(x2,x3) predict(reg,newdata=data.frame(X2=x2,X3=x3))
> VX2=seq(0,80)
> VX3=seq(5,25)
> VY=outer(VX2,VX3,y)
> image(VX2,VX3,VY)
> contour(VX2,VX3,VY,add=TRUE)

qui correspond à un plan de régression

> persp(VX2,VX3,VY,theta=30,ticktype=detailed)

On reviendra plus en détails sur ce point, mais il est possible de faire des régressions non linéaires assez facilement, à partir de ce modèle linéaire. On avait commencé par un modèle linéaire sur le logarithme de la distance

> plot(cars$speed,log(cars$dist))
> reg1=lm(log(dist)~speed,data=cars)
> abline(reg1,col="red")

(on le verra, ce n’est pas fini, car on n’a pas ici de prévision sur la distance, juste sur son logarithme… mais promis, on en reparlera) ou sur la racine carrée

> plot(cars$speed,sqrt(cars$dist))
> reg1=lm(sqrt(dist)~speed,data=cars)
> abline(reg1,col="red")

Au lieu de transformer la variable d’intérêt, on peut aussi transformer la variable explicative. On peut pendre des puissances, ou des fonctions simples, mais aussi mettre des ruptures. On avait commencé par une variable indicatrice,

> plot(cars$speed,cars$dist)
> s=10
> abline(v=s,col="green")
> regs=lm(dist~speed+I(speed>s),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + I(speed > s), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.472  -9.559  -2.088   7.456  44.412 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -17.2964     6.7709  -2.555   0.0139 *  
speed              4.3140     0.5762   7.487  1.5e-09 ***
I(speed > s)TRUE  -7.5116     7.8511  -0.957   0.3436    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.39 on 47 degrees of freedom
Multiple R-squared:  0.6577,	Adjusted R-squared:  0.6432 
F-statistic: 45.16 on 2 and 47 DF,  p-value: 1.141e-11

Mais on peut aussi mettre des fonctions afin d’avoir un modèle linéaire par morceaux, tout en étant continu

> plot(cars)
> s=15
> abline(v=s,col="green")
> positive=function(x) ifelse(x>0,x,0)
> regs=lm(dist~speed+positive(speed-s),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + positive(speed - s), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.502  -9.513  -2.413   5.195  45.391 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)          -7.6519    10.6254  -0.720  0.47500   
speed                 3.0186     0.8627   3.499  0.00103 **
positive(speed - s)   1.7562     1.4551   1.207  0.23350   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.31 on 47 degrees of freedom
Multiple R-squared:  0.6616,	Adjusted R-squared:  0.6472 
F-statistic: 45.94 on 2 and 47 DF,  p-value: 8.761e-12

On a ici une rupture, mais on pourrait imaginer en avoir plusieurs

> nouvellebase=data.frame(speed=5:25)
> y=predict(regs,newdata=nouvellebase)
> lines(5:25,y,col="red")
> 
> plot(cars$speed,cars$dist)
> s1=10
> s2=20
> abline(v=c(s1,s2),col="green")
> positive=function(x) ifelse(x>0,x,0)
> regs=lm(dist~speed+positive(speed-s1)+positive(speed-s2),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + positive(speed - s1) + positive(speed - s2), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.374  -9.475  -2.625   6.639  43.914 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)  
(Intercept)           -7.6305    16.2941  -0.468   0.6418  
speed                  3.0630     1.8238   1.679   0.0998 .
positive(speed - s1)   0.2087     2.2453   0.093   0.9263  
positive(speed - s2)   4.2812     2.2843   1.874   0.0673 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11

Comme vu en cours, le test de significativité des deux derniers coefficients ne veut pas dire que la pente est nulle, mais qu’elle est significativement différente de cette obtenue sur la zone de gauche (avant les deux seuils).

Modèles de Prévision

Cet session d’été, je redonne le cours ACT6420, modèles de prévision. Le plan de cours est maintenant en ligne (comme souvent, les salles de cours ne sont pas encore attribuées définitivement mais elles seront disponibles sur le site de l’UQAM). Les transparents pour la première partie de cours sont en ligne (il s’agit de rappels de statistique, et d’une introduction aux modèles de régression).

Des exercices et des codes informatiques seront mis en ligne sur /courses/modeles-de-prevision.

Durant le cours, j’essayerais d’illustrer avec le maximum d’exemples sur des données réelles, à partir de codes R. Quelques séances de démonstrations se feront en salle informatique, afin de découvrir les rudiments de R. Pour aller plus loin, je recommande

  • “R pour les débutants”d’Emmanuel Paradis (PDF)
  • “Brise Glace-R” d’Andrew Robinson er Arnaud Schloesing (PDF)
  • “Introduction à la programmation en R” de Vincent Goulet (PDF)

ainsi que quelques tutoriels en ligne, par exemple, dans un ancien billet, mais aussi sur https://youtube.com/.

ACT6420 examen final

Mercredi prochain, c’est l’examen final (qui compte pour 30%). Au programme, comme annoncé ce matin, la forme sera proche de celle de l’examen intra, avec 33 questions à choix multiple

  • quelques questions de compréhension générales sur la modélisation des séries temporelles,
  • quelques questions portant sur de l’analyse de sorties obtenues suite à une modélisation d’une série.

Cette session, la série à étudier sera celle obtenue sur la fréquentation d’un aéroport, sur une quinzaine d’années. Les données sont mensuelles, et sont en ligne via le code suivant

> base=read.table(
"http://freakonometrics.blog.free.fr/public/data/TS-examen.txt",
+ sep=";",header=TRUE)
> X=ts(base$X,start=c(base$A[1],base$M[1]),frequency=12)
> plot(X)

Les annexes qu’il faudra discuter à l’examen sont en ligne. Est-il utile de préciser que je ne répondrais aps aux questions sur ce document d’ici mercredi ?

Bon courage.

Save R objects, and other stuff

Yesterday, Christopher asked me how to store an R object, in order to save some time, when working on the project.

First, download the csv file for searches related to some keyword, via http://www.google.com/trends/, for instance “sunglasses“. Recall that csv files store tabular data (numbers and text) in plain-text form, with comma-separated values (where csv term comes from). Even if it is not a single and well-defined format, it is a text file, not an excel file!

Recherche sur le Web : intérêt pour sunglasses
Dans tous les pays; De 2004 à ce jour

Intérêt dans le temps
Semaine,sunglasses
2004-01-04 - 2004-01-10,48
2004-01-11 - 2004-01-17,47
2004-01-18 - 2004-01-24,51
2004-01-25 - 2004-01-31,50
2004-02-01 - 2004-02-07,52

(etc) The file can be downloaded from the blog,

> report=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/glasses.csv",
+ skip=4,header=TRUE,sep=",",nrows=464)

Then, we have run a function on this data frame, to transform it. It can be found from a source file

> source(“http://freakonometrics.blog.free.fr/public/code/H2M.R”)

In this source file, there is function that transforms a weekly series into a monthly one. The output is either a time series, or a numeric vector,

> sunglasses=H2M(report,lang="FR",type="ts")

Here, we asked for a time series,

> sunglasses
          Jan      Feb      Mar      Apr
2004 49.00000 54.27586 66.38710 80.10000
2005 48.45161 58.25000 69.93548 80.06667
2006 49.70968 57.21429 67.41935 82.10000
2007 47.32258 55.92857 70.87097 84.36667
2008 47.19355 54.20690 64.03226 79.36667
2009 45.16129 50.75000 63.58065 76.90000
2010 32.67742 44.35714 58.19355 70.00000
2011 44.38710 49.75000 59.16129 71.60000
2012 43.64516 48.75862 64.06452 70.13333
          May      Jun      Jul      Aug
2004 83.77419 89.10000 84.67742 73.51613
2005 83.06452 91.36667 89.16129 76.32258
2006 86.00000 92.90000 93.00000 72.29032
2007 86.83871 88.63333 84.61290 72.93548
2008 80.70968 80.30000 78.29032 64.58065
2009 77.93548 70.40000 62.22581 51.58065
2010 71.06452 73.66667 76.90323 61.77419
2011 74.00000 79.66667 79.12903 66.29032
2012 79.74194 82.90000 79.96774 69.80645
          Sep      Oct      Nov      Dec
2004 56.20000 46.25806 44.63333 53.96774
2005 56.53333 47.54839 47.60000 54.38710
2006 51.23333 46.70968 45.23333 54.22581
2007 56.33333 46.38710 44.40000 51.12903
2008 51.50000 44.61290 40.93333 47.74194
2009 37.90000 30.38710 28.43333 31.67742
2010 50.16667 46.54839 42.36667 45.90323
2011 52.23333 45.32258 42.60000 47.35484
2012 54.03333 46.09677 43.45833

that we can plot using

> plot(sunglasses)

Now we would like to store this time series. This is done easily using

> save(sunglasses,file="sunglasses.RData")

Next time we open R, we just have to use

> load("/Users/UQAM/sunglasses.RData")

to load the time series in R memory. So saving objects is not difficult.

Last be not least, for the part on seasonal models, we will be using some functions from an old package. Unfortunately, on the CRAN website, we see that

but nicely, files can still be found on some archive page. On Linux, one can easily install the package using (in R)

> install.packages(
+ "/Users/UQAM/uroot_1.4-1.tar.gz",
+ type="source")

With a Mac, it is slightly more complicated (see e.g. Jon’s blog): one has to open a Terminal and to type

R CMD INSTALL /Users/UQAM/uroot_1.4-1.tar.gz

On Windows, it is possible to install a package from a zipped file: one has to download the file from archive page, and then to spot it from R.

The package is now installed, we just have to load it to play with it, and use functions it contains to tests for cycles and seasonal behavior,

> library(uroot)
> CH.test
function (wts, frec = NULL, f0 = 1, DetTr = FALSE, ltrunc = NULL)
{
s <- frequency(wts)
t0 <- start(wts)
N <- length(wts)
if (class(frec) == "NULL")
frec <- rep(1, s/2)
if (class(ltrunc) == "NULL")
ltrunc <- round(s * (N/100)^0.25)
R1 <- SeasDummy(wts, "trg")
VFEalg <- SeasDummy(wts, "alg")

(etc)

Introduction aux processus SARIMA

Quelques transparents en plus, qui devraient correspondre aux deux prochains cours de séries temporelles, sur les processus autorégressifs (AR) et moyennes mobiles (MA), les ARMA, les ARIMA (intégrés) et les SARIMA (saisonniers). J’ai mis des notes sur les tests de racine unités, je rajouterais quelques transparents la semaine prochaine sur les tests de saisonnalité, et quelques exemples pratiques de prévision. Les transparents sont en ligne ici,

Introduction aux séries temporelles

Nous allons commencer demain la modélisation des séries temporelles. Les transparents de la séance sont en ligne ici. Je rappelle que les notes de cours (complètes) sont également en ligne, . Je continuerai à poster régulièrement des billets contenant des commandes R.

Au programme cette semaine, l’utilisation des méthodes de régression pour extraire tendance et cycle, le lissage exponentielle (simple, double et saisonnier), et une présentation des notions importantes dans le cours (stationnarité, autocorrélations, bruit blanc, etc).

sur les transformations dans un modèle linéaire

Je voulais prendre 5 minutes pour reprendre une question posée par courriel, qui me permettra de poursuivre sur des choses évoquées mercredi dernier en cours. Je vais reformuler la question, mais en gros, cela disait: la méthode de Box-Cox (évoquée en début de semaine, ici) avait pour objet de choisir entre un modèle linéaire (que l’on étudié dans tous les sens depuis le début du cours) et un modèle log-linéaire (évoqué ici, par exemple). L’idée était d’associer le cas linéaire à  (dans la transformée de Box-Cox) et à  pour le cas multiplicatif (ou log-linéaire). Mais que se passe-t-il si la valeur optimale rejette ces deux cas, et est proche (disons) de  ?
Considérons le cas suivant (exemple que je ne me lasse d’utiliser)

> reglm=lm(dist~speed,data=cars)
> library(MASS)
> boxcox(reglm)

La valeur optimale est effectivement proche de . En posant , on aurait envie de considérer un modèle de la forme

(comme le suggère la transformation de Box-Cox). On peut faire la régression,

> regsqrt=lm(sqrt(dist)~speed,data=cars)
+ summary(regsqrt)

Call:
lm(formula = sqrt(dist) ~ speed, data = cars)

Residuals:
Min      1Q  Median      3Q     Max
-2.0684 -0.6983 -0.1799  0.5909  3.1534

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.27705    0.48444   2.636   0.0113 *
speed        0.32241    0.02978  10.825 1.77e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.102 on 48 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.7034
F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14

et effectivement, le résultat est concluant. Mais on ne peut pas en rester là, comme pour pour le passage au logarithme, notre variable d’intéret reste .
La prédiction sur le modèle transformé était

et

Comme avec le logarithme, on pourrait prendre comme prédiction pour  mais

Bref, on a un estimateur biaisé. Et là encore, comme on est sur des variables positives, l’inégalité de Jensen doit meme nous garantir que l’une des quantités domine toujours l’autre (je laisse les amateurs de convexité l’écrire dans le bon sens). Comment faire ?
La première solution est de noter que

Donc de manière très simple, on a notre prédiction (en prenant le carré, comme intuité, mais en rajoutant une terme – positif – lié à la variance).
Une seconde consiste à noter que Z (conditionnellement à la variable explicative) était gaussienne. Donc en prenant le carré (moyennant quelques changement d’échelle), on devrait tomber sur une loi du chi-deux qui est une loi que l’on connait bien (sinon je peux renvoyer ici).
Essayons de creuser un peu ces idées (en particulier, pour dériver des intervalles de confiance). Si on visualise la prédiction de ce modèle sur la racine carrée de la distance de freinage, on obtient

> plot(speed,sqrt(dist))
> x=seq(0,30,by=.2)
> distsqrtp=predict(regsqrt,newdata=
+ data.frame(speed=x),interval="prediction")
> polygon(c(x,rev(x)),c(distsqrtp[,2],
+ rev(distsqrtp[,3])),col="yellow",border=NA)
> lines(x,distsqrtp[,1],lwd=2,col="red")

Pour passer d’une loi normale à une loi du chi-deux, il faut retrancher la moyenne (pour centrer la variable) et diviser par l’écart-type (pour avoir une variance unitaire),

> s=summary(regsqrt)$sigma
> mu=predict(regsqrt)
> distsqrtp01=(sqrt(dist)-mu)/s

On a ainsi une variable qui, conditionnellement à la variable explicative est supposé suivre une loi ,

> plot(speed,distsqrtp01,ylim=c(-3,3))
> abline(h=qnorm(c(.025,.975),0,1),lty=2,col="red")

On prend ensuite le carré pour obtenir notre loi , on on trace la bande de confiance (il faut que la valeur à la deux soit faible),

> plot(speed,distsqrtp01^2)
> abline(h=qchisq(.95,df=1),lty=2,col="red")

On a maitenant nos intervalles de confiance, en utilisant cette loi du chi-deux, en inversant notre transformation,

> mu=predict(regsqrt,newdata=data.frame(speed=x))
> distsup=(mu+s*sqrt(qchisq(.95,df=1)))^2
> distinf=(mu-s*sqrt(qchisq(.95,df=1)))^2

que l’on peut visualiser simplement (sur les données de base, avec en ordonnée la distance)

> plot(cars)
> lines(x,distsup,lty=2,col="red")
> lines(x,distinf,lty=2,col="red")

Et la valeur prédite ? On va utiliser notre relation liant la variance, l’espérance du carré, et le carré de l’espérance,

> distesp=mu^2+s^2
> lines(x,distesp,lwd=2,col="red")

Nice, n’est-ce pas ?
Sinon, comme je le disais en cours, si la transformation optimale sur  est de prendre , peut-etre pourrait-on envisager de prendre  (de manière un peu duale) comme variable explicative. Encore une fois, c’est une idée, car rien ne le garantit. En effet

En particulier, le second modèle est un modèle linéaire classique, avec de l’homoscédasticité: on aura une dispersion uniforme autour de notre parabole. En revanche, avec le premier modèle, on a une espèce d’hétéroscédasticité (comme tenu du double produit), avec une variance du terme d’erreur qui croit avec  (car les coefficients sont positifs). Bref, en terme d’intervalles de confiance, on devrait avoir des choses assez différentes.
Regardons la régression de la distance (cette fois) sur la vitesse, et le carré de la vitesse,

> reglm=lm(dist~speed+I(speed^2),data=cars)
> distp=predict(reglm,newdata=
+ data.frame(speed=x),interval="prediction")

Si on visualise la prédiction, avec un intervalle de confiance, on obtient

> plot(cars)
> polygon(c(x,rev(x)),c(distp[,2],
+ rev(distp[,3])),col="yellow",border=NA)
> lines(x,distp[,1],lwd=2,col="red")

Autrement dit, on suppose que notre modèle est homoscédastique. Les régions de confiance entre les deux approches sont clairement différentes… meme si les prédictions sont presque superposées (oui oui, il y a deux courbes, une rouge et une bleue).

On Box-Cox transform in regression models

A few days ago, a former student of mine, David, came back to me about Box-Cox tests in linear models. It made me look more carefully at the test, and I do not understand what is computed, to be honest. Let us start with something simple, like a linear simple regression, i.e.

https://latex.codecogs.com/gif.latex?Y_i=\beta_0+\beta_1%20X_i+\varepsilon_i

Let us introduced – as suggested in Box & Cox (1964) – the following family of (power) transformations

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}%20=%20\begin{cases}%20\dfrac{Y_i^\lambda-1}{\lambda}%20&\text{%20%20}%20(\lambda%20\neq%200)\\[8pt]%20\log{(Y_i)}%20%20&\text{%20}%20(\lambda%20=%200)%20\end{cases}

on the variable of interest. Then assume that

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}=\beta_0+\beta_1%20X_i+\varepsilon_i

As mentioned in Chapter 14 of Davidson & MacKinnon (1993) – in French – the log-likelihood of this model (assuming that observations are independent, with distribution https://latex.codecogs.com/gif.latex?\varepsilon_i\sim\mathcal{N}(0,\sigma^2)) can be written

https://latex.codecogs.com/gif.latex?\log%20\mathcal{L}=-\frac{n}{2}\log(2\pi)-n\log\sigma%20\\-\frac{1}{2\sigma^2}\sum_{i=1}^n\left[Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)\right]^2+(\lambda-1)\sum_{i=1}^n\log%20Y_i

We can then use profile-likelihood techniques (see here) to derive the optimal transformation.

This can be done in R extremely simply,

> library(MASS)
> boxcox(lm(dist~speed,data=cars),lambda=seq(0,1,by=.1))

we then get the following graph,

If we look at the code of the function, it is based on the QR decomposition of the https://latex.codecogs.com/gif.latex?\boldsymbol{X}matrix (since we assume that https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a full-rank matrix). More precisely, https://latex.codecogs.com/gif.latex?\boldsymbol{X}=QR where https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a https://latex.codecogs.com/gif.latex?n\times%202 matrix, https://latex.codecogs.com/gif.latex?Q is a https://latex.codecogs.com/gif.latex?n\times%202 orthonornal matrix, andhttps://latex.codecogs.com/gif.latex?R is a https://latex.codecogs.com/gif.latex?2\times2 upper triangle matrix. It might be convenient to use this matrix since, for instance, https://latex.codecogs.com/gif.latex?R\widehat{\boldsymbol{\beta}}=Q%27Y.  Thus, we do have an upper triangle system of equations.

> X=lm(dist~speed,data=cars)$qr

The code used to get the previous graph is (more or less) the following,

> g=function(x,lambda){
+ y=NA
+ if(lambda!=0){y=(x^lambda-1)/lambda}
+ if(lambda==0){y=log(x)}
+ return(y)} 
> n=nrow(cars)
> X=lm(dist~speed,data=cars)$qr
> Y=cars$dist
> logv=function(lambda){
+ -n/2*log(sum(qr.resid(X, g(Y,lambda)/
+ exp(mean(log(Y)))^(lambda-1))^2))}
> L=seq(0,1,by=.05)
> LV=Vectorize(logv)(L)
> points(L,LV,pch=19,cex=.85,col="red")

As we can see (with those red dots) we can reproduce the R graph. But it might not be consistent with other techniques (and functions described above). For instance, we can plot the profile likelihood function, https://latex.codecogs.com/gif.latex?\lambda\mapsto\log\mathcal{L}

> logv=function(lambda){
+ s=summary(lm(g(dist,lambda)~speed,
+ data=cars))$sigma
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ -n/2*log(2 * pi)-n*log(s)-.5/s^2*(sum(e^2))+
+ (lambda-1)*sum(log(Y))
+ }
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> (maxf=optimize(logv,0:1,maximum=TRUE))
$maximum
[1] 0.430591

$objective
[1] -197.6966

> abline(v=maxf$maximum,lty=2)

The good point is that the optimal value of https://latex.codecogs.com/gif.latex?\lambda is the same as the one we got before. The only problem is that the https://latex.codecogs.com/gif.latex?y-axis has a different scale. And using profile likelihood techniques to derive a confidence interval will give us different results (with a larger confidence interval than the one given by the standard function),

> ic=maxf$objective-qchisq(.95,1)
> #install.packages("rootSolve")
> library(rootSolve)
> f=function(x)(logv(x)-ic)
> (lower=uniroot(f, c(0,maxf$maximum))$root)
[1] 0.1383507
> (upper=uniroot(f, c(maxf$maximum,1))$root)
[1] 0.780573
> segments(lower,ic,upper,ic,lwd=2,col="red")

Actually, it possible to rewrite the log-likelihood as

https://latex.codecogs.com/gif.latex?\mathcal{L}=\star-\frac{n}{2}\log\left[\sum_{i=1}^n\left(\frac{Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)}{\dot{Y}^\lambda}\right)^2\right]

(let us just get rid of the constant), where

https://latex.codecogs.com/gif.latex?\dot{Y}=\exp\left[\frac{1}{n}\sum_{i=1}^n%20\log%20Y_i\right]

Here, it becomes

> logv=function(lambda){
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ elY=(exp(mean(log(Y))))
+ -n/2*log(sum((e/elY^lambda)^2))
+ }
>
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> optimize(logv,0:1,maximum=TRUE)
$maximum
[1] 0.430591

$objective
[1] -47.73436

with again the same optimal value for https://latex.codecogs.com/gif.latex?\lambda, and the same confidence interval, since the function is the same, up to some additive constant.

So we have been able to derive the optimal transformation according to Box-Cox transformation, but so far, the confidence interval is not the same (it might come from the fact that here we substituted an estimator to the unknown parameter https://latex.codecogs.com/gif.latex?\sigma.

Le passage au log dans les modèles linéaires

Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l’on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à  suit une loi normale. Pour rappel, si on a une loi normale, , alors  et . Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant),

Dans un modèle Gaussien avec homoscédasiticité,  i.e.  alors que . On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire,

Bon, maintenant, que se passe-t-il si on prend l’exponentiel ? Pour la loi normale, rappelons que l’on obtient une loi lognormale, i.e. , les deux paramètres étant liés à la loi normale sous jacente, car désormais

alors que

Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est   alors que le point bleu est l’espérance de la loi lognormale.

On notera que le quantile de la loi log-normale est l’exponentiel du quantile de la loi normale. En effet, si https://latex.codecogs.com/gif.latex?\mathbb{P}(Y\leq%20q)=\alpha alors https://latex.codecogs.com/gif.latex?\mathbb{P}(\exp(Y)\leq%20\exp(q))=\alpha. En particulier,  n’est pas la moyenne de , mais la médiane (puisque  était la médiane de ).

Mais il n’est pas rare de voir utilisé un intervalle de confiance de la forme

qui est la forme classique de l’intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants

Notons qu’il n’y a aucune raison ici d’avoir une probabilité https://latex.codecogs.com/gif.latex?1-\alpha d’être dans l’intervalle de confiance obtenu avec les quantiles https://latex.codecogs.com/gif.latex?q_{1-\alpha/2} de la loi normale.

Maintenant, si on prend l’exponentiel d’un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a

avec une variance (conditionnelle) qui dépend de la variable explicative

Là encore, le plus naturel est d’utiliser comme bornes de l’intervalle de confiance des quantiles associés à la loi lognormale,

mais il n’est pas rare de voir utilisé des intervalles de type Gaussiens,

On perd là encore en interprétation car les bornes n’ont plus rien à voir avec les quantiles.

Test, valeur critique et p-value

Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la https://latex.codecogs.com/gif.latex?p-value dans la lecture de la sortie d’un test.

  • Les erreurs dans un test statistique

Mais avant, rappelons qu’un test est une prise de décision: accepter ou rejeter une hypothèse. Et qu’on peut commettre une erreur. Ou pour être plus précis, on peut commettre deux types d’erreur,
• accepter l’hypothèse alors que cette dernière est fausse
• rejeter l’hypothèse alors que cette dernière était vraie
Pour reprendre une terminologie plus médicale, un test de grossesse peut dire à une femme qu’elle n’est pas enceinte, alors qu’elle l’est; ou dire qu’elle l’est, alors qu’elle ne l’est pas (voir tous les exemples dans les exercices de probabilités de l’examen P de la SOA, ou le cours ACT2121).
Formellement, on a deux probabilités,
• la probabilité d’accepter à tort notre hypothèse (on parlera d’erreur de second espèce), 
• la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) 
Dans un monde idéal on voudrait que les deux probabilités soient aussi petites que possibles…  Mais c’est impossible, et le plus souvent, baisser une des probabilités se fait en augmentant l’autre. Les cas extrêmes étant
• avoir un test de grossesse qui déclare tout le monde enceinte: on ne rejette alors jamais à tort (on ne rejette jamais tout court en fait), mais on a un fort taux d’acceptation à tort,
• avoir un test de grossesse qui ne déclare personne enceinte: on n’accepte jamais à tort (car on n’accepte jamais) mais on a un fort taux de rejet à tort.
Bref, on a un arbitrage à faire entre deux types d’erreurs. Souvent, en pratique on va demander à contrôler l’erreur de première espèce (i.e.  de l’ordre de 5%), et on chercher a un test qui, à  donné, possède la plus faible erreur de première espèce. Voilà en gros pour la théorie: on se donne un seuil de significativité , qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse https://latex.codecogs.com/gif.latex?H_0 est vraie, l’alternative étant une hypothèse https://latex.codecogs.com/gif.latex?H_1.

https://latex.codecogs.com/gif.latex?H_0 vraie https://latex.codecogs.com/gif.latex?H_1 vraie
accepter https://latex.codecogs.com/gif.latex?H_0 OK erreur
type 2
rejeter https://latex.codecogs.com/gif.latex?H_0 erreur
type 1
OK
  • La “valeur critique”

La notion de valeur critique a été introduite dans Neyman & Pearson (1928). Cette valeur dépend de la forme de l’hypothèse alternative, en particulier savoir si le test est bilatéral, unilatéral à gauche, ou unilatéral à droite. Pour un test donné, la valeur critique peut-être vue comme la valeur limite a partir de laquelle on pourra rejeter https://latex.codecogs.com/gif.latex?H_0 avec un seuil de significativité donné.

  • La p-value

La https://latex.codecogs.com/gif.latex?p-value a été introduite dans Gibbons & Pratt (1975), meme si on peut retrouve l’idée beaucoup plus tôt, comme Pearson (1900), qui propose de calculer “the probability that the observed value of the chi-square statistic would be exceeded under the null hypothesis“. La https://latex.codecogs.com/gif.latex?p-value est la probabilité, sous https://latex.codecogs.com/gif.latex?H_0, d’obtenir une statistique aussi extrême (pour ne pas dire aussi grande) que la valeur observée sur l’échantillon. Aussi, pour un seuil de significativité  donné, on compare https://latex.codecogs.com/gif.latex?p et , afin d’accepter, ou de rejeter https://latex.codecogs.com/gif.latex?H_0,
• si , on va rejeter l’hypothèse https://latex.codecogs.com/gif.latex?H_0 (en faveur de )
• si , on va rejeter  (en faveur de https://latex.codecogs.com/gif.latex?H_0).
On peut alors interpréter la https://latex.codecogs.com/gif.latex?p-value comme le plus petit seuil de significativité pour lequel l’hypothèse nulle est acceptée. Gibbons & Pratt (1975) reviennent longuement sur les interprétations, et surtout les mauvaises interprétations, de cette https://latex.codecogs.com/gif.latex?p-value.

  • Valeur critique versus p-value

Si on formalise un peu, on peut vouloir tester  contre  (par exemple). De manière très générale, on dispose d’une statistique de test  qui a pour loi, sous  (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme , c’est juste plus pénible parce qu’il faut travailler sur , et calculer des probabilités à gauche, ou à droite. Donc pour notre exemple, on va prendre un test unilatéral.
Dans l’approche classique (telle que présentée dans tous les cours de statistiques), on se donne un seul d’acceptation  petit (disons 5%), et on cherche une valeur critique  telle que

Pour ceux qui se souviennent de leur cours de stats, cela peut faire penser à la puissance du test, définie par

Formellement, la https://latex.codecogs.com/gif.latex?p-value associée au test  est la variable aléatoire  définie par

Donc effectivement, la https://latex.codecogs.com/gif.latex?p-value et la puissance sont liées, puisque

autrement dit, la puissance peut-être vue comme la fonction de répartition de la https://latex.codecogs.com/gif.latex?p-value.

  • Intérêt computationnel de la p-value

D’un point de vue computationnel, la https://latex.codecogs.com/gif.latex?p-value est l’outil le plus important pour interpréter la sortie d’un test. Commençons par un test simple, comme une comparaison de moyennes. On cherche ici à tester  contre  pour des moyennes calculées sur deux groupes. Pour reprendre l’exemple abordé dans un précédant billet, on a les notes obtenues en ACT6420 par deux groupes différents. Et on veut savoir s’ils sont vraiment différents (ci-dessous le nombre de bonnes réponses, sur 40 questions, on travaillera ensuite sur la note sur 100)

La statistique de test est ici

https://latex.codecogs.com/gif.latex?%20%20%20%20%20T%20=%20\frac{\overline{X}%20-%20\overline{Y}}{\displaystyle{%20\sqrt{%20{s_X^2%20\over%20n_X}%20+%20{s_Y^2%20\over%20n_Y}%20}}}

et sous https://latex.codecogs.com/gif.latex?H_0https://latex.codecogs.com/gif.latex?%20%20%20%20%20T va suivre une loi de Student à https://latex.codecogs.com/gif.latex?%20%20%20%20\nu degrés de liberté, où https://latex.codecogs.com/gif.latex?%20%20%20%20\nu est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)),

https://latex.codecogs.com/gif.latex?%20%20\nu%20=%20{{\left(%20{s_X^2%20\over%20n_X}%20+%20{s_Y^2%20\over%20n_Y}\right)^2%20}%20%20\over%20{{s_X^4%20\over%20n_X^2%20\cdot%20\left({n_X-1}\right)}+{s_Y^4%20\over%20n_Y^2%20\cdot%20\left({n_Y-1}\right)}}}%20\,

Numériquement, on a ici

> Xbar=mean(X)
> Ybar=mean(Y)
> Sx2=var(X)
> Sy2=var(Y)
> nX=length(X)
> nY=length(Y)
> (T=(Xbar-Ybar)/sqrt(Sx2/nX+Sy2/nY))
[1] -2.155754

et pour les degrés de liberté

> (nu=(Sx2/nX+Sy2/nY)^2/(Sx2^2/nX^2/(nX-1)+
+ Sy2^2/nY^2/(nY-1)))
[1] 36.35279

La valeur critique est obtenue en lisant dans les tables,

(car ici on a des probabilité pour un test bilatéral dans la table) comme on apprenait dans les cours de statistique au siècle passé. D’un point de vue informatique, on cherche à savoir si on est à gauche, ou à droite de la valeur critique

> qt(.05,df=nu)
[1] -1.687865

On peut aussi calculer la https://latex.codecogs.com/gif.latex?p-value,

> pt(T,df=nu)
[1] 0.01889768

Si on regarde, sous R, il existe des fonctions de tests, pour comparer des moyennes. Et dans ce cas, la sortie est

> t.test(X,Y,alternative = "less")

Welch Two Sample t-test

data:  X and Y
t = -2.1558, df = 36.353, p-value = 0.0189
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf -1.772507
sample estimates:
mean of x mean of y
48.75000  56.91667

Autrement dit, on a automatiquement la https://latex.codecogs.com/gif.latex?p-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une https://latex.codecogs.com/gif.latex?p-value (et que l’on vérifié au préalable les conditions d’application d’un test), on peut faire tous les tests que l’on veut !
Si on veut faire un peu plus compliqué, on peut regarder la distribution des notes, et se demander si une loi  serait possible (par exemple, ça sera notre hypothèse https://latex.codecogs.com/gif.latex?H_0, l’hypothèse alternative étant que ce n’est pas cette loi). Pour faire ce test, il existe le test de Kolmogorov-Smirnov. La statistique de test est ici

où  est la fonction de répartition de la loi , et  est la fonction de répartition empirique

La loi de https://latex.codecogs.com/gif.latex?%20%20%20%20%20T n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les https://latex.codecogs.com/gif.latex?p-values automatiquement,

> ks.test(Y, "pnorm", 60, 15)

One-sample Kolmogorov-Smirnov test

data:  Y
D = 0.1421, p-value = 0.5796
alternative hypothesis: two-sided

Aussi, on peut accepter ici l’hypothèse nulle. On peut d’ailleurs faire un petit dessin pour s’en convaincre,

> Femp=function(x) mean(Y<=x)
> plot(0:100,Vectorize(Femp)(0:100),type="s")
> lines(0:100,pnorm(0:100,60,15),col="red")

 


Et ça va nous servir dans ce cours ? A priori oui… parce qu’on parlera du test de Student (pour tester si une variable dans une régression est significative), du test de Fisher (pour tester si plusieurs variables dans une régression sont significatives, ou plus généralement si une contrainte – linéaire – sur les coefficients peut être acceptée), du test de Chow (pour tester des ruptures dans un modèle linéaire, mais c’est un test de Fisher un peu déguisé), du test d’Anderson-Darling (pour tester si des résidus sont Gaussiens), du test de Breuch-Pagan voire le test de White (pour tester si les résidus peuvent être considérés de variance constante), du test de Durbin-Watson (pour tester s’il n’y a pas d’auto-corrélation dans la série des résidus), du test de Dickey-Fuller (pour tester si une série temporelle est – ou n’est pas – stationnaire), des tests de Franses (pour tester si une série peut être considérée comme saisonnière, ou pas), du test de Ljung-Box (pour tester si un bruit est un bruit blanc)… Et j’en oublie un paquet. Donc quand il est dit (dans le plan de cours) que le cours de statistique est un prérequis, il ne s’agit pas de l’avoir suivi, mais bel et bien de l’avoir compris, car on passera notre temps à utiliser des notions entrevues dans ce cours.