Tag Archives: ACT6420

Test, valeur critique et p-value

Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la 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), \beta
• la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) \alpha
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. \alpha de l’ordre de 5%), et on chercher a un test qui, à \alpha 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é \alpha, qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse H_0 est vraie, l’alternative étant une hypothèse H_1.

H_0 vraie H_1 vraie
accepter H_0 OK erreur
type 2
rejeter 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 H_0 avec un seuil de significativité donné.

  • La p-value

La 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 p-value est la probabilité, sous 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é \alpha donné, on compare p et \alpha, afin d’accepter, ou de rejeter H_0,
• si p\leq\alpha, on va rejeter l’hypothèse H_0 (en faveur de H_1)
• si p>\alpha, on va rejeter H_1 (en faveur de H_0).
On peut alors interpréter la 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 p-value.

  • Valeur critique versus p-value

Si on formalise un peu, on peut vouloir tester H_0:\theta=\theta_0 contre H_1:\theta>theta_0 (par exemple). De manière très générale, on dispose d’une statistique de test T qui a pour loi, sous H_0, F_{\theta_0}(\cdot) (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme H_1:\theta\neq\theta_0, c’est juste plus pénible parce qu’il faut travailler sur \vert T\vert, 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 \alpha petit (disons 5%), et on cherche une valeur critique T_{1-alpha} telle que

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

\pi(\theta\vert \alpha)=\mathbb{P}(T\geq T_{1-\alpha}\vert \theta)=1-F_{\theta}(T_{1-\alpha})

Formellement, la p-value associée au test T est la variable aléatoire P définie par
P=1-F_{\theta_0}(T).
Donc effectivement, la p-value et la puissance sont liées, puisque

\mathbb{P}(P\leq \alpha\vert \theta)=\pi(\theta\vert \alpha)

autrement dit, la puissance peut-être vue comme la fonction de répartition de la p-value.

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

D’un point de vue computationnel, la 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 H_0:\mu_X=\mu_Y contre H_1:\mu_X>\mu_Y 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)

image manquante

La statistique de test est ici

T = \frac{\overline{X} - \overline{Y}}{\displaystyle{ \sqrt{ {s_X^2 \over n_X} + {s_Y^2 \over n_Y} }}}

et sous H_0, T va suivre une loi de Student à \nu degrés de liberté, où \nu est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)),

\nu = {{\left( {s_X^2 \over n_X} + {s_Y^2 \over n_Y}\right)^2 } \over {{s_X^4 \over n_X^2 \cdot \left({n_X-1}\right)}+{s_Y^4 \over n_Y^2 \cdot \left({n_Y-1}\right)}}}

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

image manquante

On peut aussi calculer la 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 p-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une 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 \mathcal{N}(60,15^2) serait possible (par exemple, ça sera notre hypothèse 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

T=\sup\{\vert \widehat{F}_n(x)-F_0(x)\vert ,x\in\mathbb{R}\}

F_0(\cdot) est la fonction de répartition de la loi \mathcal{N}(60,15^2), et \widehat{F}_n(\cdot) est la fonction de répartition empirique

\widehat{F}_n(x)=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(x_i\leq x)

La loi de T n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les 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")

image manquante


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.

Régression: données individuelles ou temporelles ?

Avant de parler davantage du premier devoir, un mot pour expliquer pourquoi dans la première partie il est important d’éviter d’utiliser des données temporelles (le cours présente des outils pour modéliser des données individuelles). Considérons un exemple simple (évoqué voilà plusieurs mois ici même), où on cherche à prédire un tour de poitrine conditionnellement à un tour de taille (pour une personne de sexe féminin, cf. le mannequin ci-contre). L’idée intuitive est qu’il doit exister une relation positive, possiblement linéaire entre les deux mensurations. A partir de données, on devrait pouvoir quantifier cette relation. On dispose pour cela de données, en ligne sur

> mensurations<-read.table("
+ http://freakonometrics.free.fr/mensurations.csv",
+ header=TRUE,sep=";")
> head(mensurations[,-1])
POITRINE TAILLE HANCHE HAUTEUR POIDS      AGE
1       NA     NA     NA      NA    NA       NA
2       94     58     89     165    51 22.24230
3       NA     NA     NA      NA    NA       NA
4       94     58     89     165    51 22.40383
5       91     58     97     165    55       NA
6       NA     NA     NA      NA    NA       NA

Si on regarde rapidement la corrélation entre la taille et la poitrine, cette dernière est faible,

> Z=mensurations[,c("TAILLE","POITRINE")]
> I=apply(Z,1,function(x) sum(is.na(x)))
> cor(Z[I==0,])
TAILLE   POITRINE
TAILLE   1.00000000 0.09700002
POITRINE 0.09700002 1.00000000

… mais significativement non nulle,

> cor.test(Z[,1],Z[,2],,method = "pearson")

Pearsons product-moment correlation

data:  Z[, 1] and Z[, 2]
t = 2.4559, df = 635, p-value = 0.01432
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.0194634 0.1733769
sample estimates:
cor
0.09700002

Avec une https://latex.codecogs.com/gif.latex?p-value de l’ordre de 1.5%, on va rejeter l’hypothèse https://latex.codecogs.com/gif.latex?H_0 de nullité de la corrélation (c’est ici un test basé sur la transformation de Fisher). On peut aussi faire une régression simple,

> reg=lm(POITRINE~TAILLE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7253 -3.4424 -0.2091  1.5576 13.8579

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.67798    2.82953  29.220   <2e-16 ***
TAILLE       0.11663    0.04749   2.456   0.0143 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.133 on 635 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.009409,	Adjusted R-squared: 0.007849
F-statistic: 6.031 on 1 and 635 DF,  p-value: 0.01432

(on notera que le test sur la corrélation, le test de significativité de la – seule – variable explicative, et le test de significativité globale donnent ici exactement la même https://latex.codecogs.com/gif.latex?p-value de l’ordre de 1.5%). Là encore, on voit que le tour de taille a un effet significativement non nul… mais marginal. Graphiquement, la prédiction avec le tour de taille est le trait rouge,

La valeur moyenne est en trait pointillé. On peut voir que la différence est infime,

> base=data.frame(TAILLE=seq(50,70,by=5))
> base$Pred.POITRINE=predict(reg,newdata=base)
> base$Moy.POITRINE=mean(POITRINE,na.rm=TRUE)
> base
TAILLE Pred.POITRINE Moy.POITRINE
1     50      88.50935     89.61538
2     55      89.09248     89.61538
3     60      89.67562     89.61538
4     65      90.25876     89.61538
5     70      90.84189     89.61538

Le soucis est que cette base est en fait une base sur données individuelles. En fait, il s’agit des mensurations des demoiselles qui ont posé comme playmate dans la version U.S. du magasine Playboy. Si on se restreint sur des sous-période, par exemple les années 70, on a des prédictions relativement différentes,

> submensurations=mensurations[abs(DATE-1975)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> summary(subreg)

Call:
lm(formula = POITRINE ~ TAILLE, data = submensurations)

Residuals:
Min     1Q Median     3Q    Max
-8.655 -1.871 -0.655  1.750 13.155

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  66.1485     7.3839   8.958 6.12e-15 ***
TAILLE        0.4053     0.1230   3.296   0.0013 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.574 on 117 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.08496,	Adjusted R-squared: 0.07714
F-statistic: 10.86 on 1 and 117 DF,  p-value: 0.001299

> base$Pred70s.POITRINE=predict(subreg,newdata=base)
> base
TAILLE Pred.POITRINE Moy.POITRINE Pred70s.POITRINE
1     50      88.50935     89.61538         86.41274
2     55      89.09248     89.61538         88.43917
3     60      89.67562     89.61538         90.46559
4     65      90.25876     89.61538         92.49201
5     70      90.84189     89.61538         94.51844

que l’on peut visualiser ci-dessous,

L’explication est que la morphologie des femmes (ou en tous les cas de celles qui ont posé nues dans Playboy) a beaucoup changé en 60 ans, comme évoquéprécédemment. On peut le voir en regardant l’évolution de la pente, par date (https://latex.codecogs.com/gif.latex?\pm5ans)

(le trait horizontal est la pente obtenue sur l’ensemble de la période). On peut d’ailleurs intégrer l’année dans notre modèle, en rajoutant le temps comme variable explicative dans la régression,

> reg=lm(POITRINE~TAILLE+DATE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE + DATE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7304 -2.2184 -0.5405  2.0742 12.1041

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 368.296266  18.008455  20.451  < 2e-16 ***
TAILLE        0.323319   0.042140   7.673 6.39e-14 ***
DATE         -0.150301   0.009393 -16.002  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.491 on 634 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.2944,	Adjusted R-squared: 0.2922
F-statistic: 132.3 on 2 and 634 DF,  p-value: < 2.2e-16

Visuellement, on aurait la prédiction suivante,

Si on se place autour de 1970,

> submensurations=mensurations[abs(DATE-1970)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
89.14361 95.75476

alors qu’autour de 1990,

> submensurations=mensurations[abs(DATE-1990)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
87.69853 92.70228

Notons qu’on fait ici des régressions locales (dans ces deux derniers calculs), ce qui n’est pas ce qui est représenté visuellement… Si on veut une meilleure modélisation locale, il faudrait quitter le monde des modèles linéaires pour aller un peu plus loin. On verra en fin de cours si on a le temps d’en parler

Cela dit, globalement, on avait un effet positif, mais il aurait été tout à fait possible d’avoir un effet (global négatif). C’est ce qui arrive dans le paradoxe de Simpson.

La recherche, l’enseignement et les élèves

Depuis quelques temps, lorsqu’on me demande ce que je fais dans la vie, je réponds “enseignant-chercheur” et je trouve que le terme me convient bien…

I.− Part. prés. de enseigner*: a) [L’obj. désigne le destinataire] Transmettre un savoir à. Synon.

Prononc. et Orth. : [ʃ ε ʀ ʃ œ:ʀ], fém. [-ø:z]. Ds Ac. 1694-1932. Étymol. et Hist. 1. le plus souvent avec compl. introduit par de a) 1538 cercheur « personne qui cherche, qui s’informe » (Est., s.v. conquisitores); 1636 chercheur (Monet); b) 1840 adj. « qui cherche (qqc.) » (Proudhon, Qu’est-ce que la propriété, p. 133); 2. emploi abs. 1835 adj. chercheuse « avide de connaissances, qui essaie de découvrir » (Balzac, loc. cit.);

Si on regarde les définitions, on voit qu’un chercheur s’informe (et reçoit du savoir) alors que l’enseignant le transmet (et communique son savoir). D’aucuns prétendront que le chercheur aussi transmet, mais quand on voit l’impact des articles académiques, on peut se demander à qui on transmet. Sans chercher à polémiquer, je trouve d’ailleurs surprenant qu’en France on ait pu séparer les deux, en créant des postes de chercheurs qui n’étaient pas tenu de partager leur savoir (mais j’ouvre une parenthèse que je ne tiens pas à ouvrir davantage)…. Bref, je suis très content d’enseigner car c’est cette activité qui me donne une utilité sociale.
Sauf que comme tous les enseignants, je déteste faire des examens. Pire encore, les corriger…
Lors de mon dernier examen, j’avais fait un questionnaire, avec des choix multiples (4 réponses, une seule bonne). Ça demande un temps fou à préparer, mais la correction est beaucoup plus rapide que lire des dizaines de pages (activité que je me laisse pour les heures à venir, car en plus des questionnaires, il y avait deux gros devoirs sur des bases de données). Et en plus, avec les choix multiples la correction est impartiale !
Dans un monde idéal, j’aimerais pouvoir dire “si vous pensez avoir réussi, passez au cours suivant, sinon reprenez le…” Alors j’ai tenté l’expérience. J’ai demandé aux étudiants d’écrire sur la copie le nombre de bonnes réponses qu’ils pensaient avoir. L’exercice était libre, et sans influence sur la note finale. Enfin, presque. Suite à une suggestion, j’ai accepté de donner des points à toute personne qui prédirait exactement son score !

> ACT6420=read.table("ACT6420-H2012.csv",
> base=ACT6420[ACT6420$obs>0,]
> attach(base)

Tout d’abord, comme l’exercice n’était pas obligatoire, certains (plusieurs) n’ont pas fait de prédiction… et comme souvent, la “non-réponse” est (très) informative,

> reponse=(is.na(pred)==FALSE)
> boxplot(obs~reponse,horizontal=TRUE)

Ceux qui n’ont pas répondu ont moins bien réussi que ceux qui ont répondu. Et sans vouloir revenir sur un point vu en cours (corrélation n’est pas causalité), je suis toujours surpris de fait que les personnes qui ne répondent pas aient un comportement aussi clairement différent de ceux qui répondent.

Maintenant regardons ce qui a été répondu par ceux qui ont joué le jeu, et qui ont fait une prédiction,

> plot(pred,obs,xlim=c(0,40),ylim=c(0,40),
+ xlab="Nombre de bonne réponses prédit",
+ ylab="Nombre de bonne réponses observé")
> abline(a=0,b=1,lty=1,col="red")
> abline(lm(obs~pred),lty=2)

En abscisse la prédiction, et en ordonnée, la valeur observée (un peu comme dans le cours de prévision justement). La courbe en rouge est une prédiction parfaite. Et la courbe en pointillé, la droite de régression. Les meilleurs se sous-estiment (en moyenne), alors que les moins bons se surestiment. Comme le note @AlbertineP on voit que les moins bons sont en fait un peu meilleurs qu’ils ne le pensent, alors que les meilleurs se surestiment ! Un effet Dunning-Kruger inversé, pour reprendre la terminologie proposée par @Tournyol. Damned, je ne pourrais donc pas demander aux étudiants de s’auto-évaluer, il va falloir que je continue à corriger mes copies ! (et je vais renouveler l’expérience, pour voir sur plus de copies ce que ça donne !)

ACT6420 compléments de lecture

Histoire de revoir, d’approfondir certaines notions abordées dans la première partie du cours, quelques liens vers de documents, ici ou là… Parmi les ouvrages cités dans le plan de cours, il y a le site de Jed Frees qui propose des codes R et des bases de données, utilisés dans le livre. De même pour le livre de William Greene, un site permet d’avoir accès aux données par exemple. Enfin, Russell Davidson et James MacKinnon ont mis en ligne l’intégralité de leur livre d’estimation et d’inférence en économétrie. Bien que le livre soit plus poussé que ce que l’on devrait faire, il apportera de nombreux éclaircissements, et une autre manière de voir le cours. Je peux aussi renvoyer aux notes de cours de Bruno Crepon et de Nicolas Jacquemet, d’économétrie linéaire appliquée.

  • statistique mathématique
La section 3. des notes de cours de Didier Concordet rappelle les principaux résultats à connaître sur l’inférence statistique, ou sinon, je peux renvoyer vers les notes de Jérôme Dupuis.
  • tests statistiques
Une référence un peu générale, le livre de Gaël Millot, en ligne sur Google Books . Sinon, les sections 4. et 5. des notes de cours de Didier Concordet reprennent l’essentiel des résultats à connaître sur les tests pour suivre le cours. Je renverrais aussi vers les transparents de Wissem Maazoun sur les tests usuels. Enfin, sur la théorie des tests, des compléments peuvent être trouvés dans les notes de cours écrites par Paul Doukhan, ou les notes de Nicolas Chopin, Emmanuelle Cretois et Jean Diebolt. (ces deux documents vont beaucoup plus loin que ce qui est demandé comme prérequis).
  • optimisation et aspects numériques
On reviendra plus longuement sur ces aspects là sur la section de séries temporelles, mais les transparents de Grégory Vial peuvent être un bon rappel (en particulier sur la descente de gradient). Je peux aussi mentionner une revue d’algorithme dans un document de Joël Pothier (qui va beaucoup plus loin que ce qui pourrait servir dans ce cours, avec page 10 la descente de gradient).
  • calcul différentiel et calcul matriciel
Sur les rappels d’algèbre matricielle, je peux renvoyer à un court document synthétique de Xavier Joutard. Et sur le calcul différentiel appliqué à des matrices (que l’on utilisera lorsque l’on abordera la régression multiple) je renvoie à la première section du document de Thomas Minka.
  • espérance conditionnelle et géométrie
Pour des rappels sur l’espérance conditionnelle, et l’interprétation géométrique en terme de projection, je revois vers des transparents de Samy Tindel, qui reprennent tout ce qu’il faut savoir. Sinon Massimiliano Gubinelli a mis en ligne un document de 4 pages, un cours document de Vincent Bansaye, ou la Section 1.4. des notes de cours d’Arnaud Guyader permettent d’apporter des éclaircissements.
  • un peu de programmation (en R)
Quelques liens vers des livres de R, d’introduction à la programmation, dont l’introduction d’Emmanuel Paradis, le brise glace d’Andrew Robinson et Arnaud Schloesing, ou le livre de Vincent Goulet. Bien qu’aucun prérequis en informatique ne soit demandé pour le cours, je m’attends à ce que l’on connaissent les principes généraux de programmation. Par exemple connaitre la différence entre une variable locale et une variable globale dans un code.
> y=cars$dist 
> f=function(z) {return(mean(log(z)))} 
> meanlog=f(y)
D’autres compléments seront mis en ligne, en particulier sur les processus, dès que nous aborderons les séries temporelles.

J’ai probablement raté des documents intéressants, les commentaires sont ouverts si quelqu’un souhaite partager des documents…

Fin de session !

La session d’hiver est enfin terminée (pour les cours en tous cas). L’examen du cours ACT6420 méthodes de prévisions est maintenant en ligne, avec l’énoncé, les annexes (qui correspondent à celles mises en ligne en début de semaine, avec quelques nombres entourés, qu’il fallait commenter expliquer) et quelques éléments (succincts) de correction (qui seront en ligne à partir de 15 heures).

Il reste encore le deuxième devoir à rendre, pour vendredi prochain dernier délai (comme indiqué mardi, après le cours). Merci aux étudiants qui ont joué le jeu, et qui ont mis une prévision de leur propre note.

Prévision des vente de voitures au Québec

Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu.

Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne,

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/car-sales-quebec.csv",
header=TRUE,sep=";",nrows=108)
Xt=ts(X[,2],start=c(1960,1),frequency=12)

On peut regarder l’évolution de cette série temporelles (car les dessins, et la visualisation, c’est important),

plot(Xt)

On note une tendance linéaire (ou qu’on pourrait supposer linéaire), que l’on va estimer,

X=as.numeric(Xt)
temps=1:length(X)
base=data.frame(temps,X)
reg=lm(X~temps)

et que l’on pourra représenter

plot(temps,X,type="l",xlim=c(0,120))
abline(reg,col="red")

Si on veut une série stationnaire (et c’est effectivement ce que l’on cherche), on va retrancher cette tendance à notre série, et travailler sur la série résiduelle.

La série résiduelle est ici

Y=X-predict(reg)
plot(temps,Y,type="l")

que l’on devrait pouvoir supposer stationnaire. Classiquement, on regarde les autocorrélations,

acf(Y,36,lwd=4)

où on repère un beau cycle annuel, mais on va surtout utiliser les autocorrélations partielles pour identifier l’ordre de la composante autorégressive.

pacf(Y,36,lwd=4)

La 12ème est non-nulle, et tous les autres ensuite sont significativement nulles (ou presque). On va donc tenter un https://latex.codecogs.com/gif.latex?AR(12).

> fit.ar12=arima(Y,order=c(12,0,0))
> fit.ar12
Series: Y
ARIMA(12,0,0) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
      0.1975  0.0832  -0.1062  -0.1212  0.1437  -0.1051  0.0319
s.e.  0.0838  0.0809   0.0826   0.0843  0.0850   0.0833  0.0854
ar9     ar10    ar11    ar12  intercept
      -0.0332  -0.0616  0.2635  0.4913  -148.3180
s.e.   0.0853   0.0840  0.0840  0.0841   384.5095

sigma^2 estimated as 2177974:  log likelihood=-946.75
AIC=1921.51   AICc=1926.03   BIC=1959.06

La douzième composante https://latex.codecogs.com/gif.latex?\phi_{12} semble significative, si on repense au test de Student,

u=seq(-6,6,by=.1)
plot(u,dnorm(u),type="l")
points(0.4913/0.0841,0,pch=19,col="red")

En revanche, la huitième https://latex.codecogs.com/gif.latex?\phi_{8} ne l’est pas

points(-0.1018/0.0847,0,pch=19,col="blue")

(voilà pour le petit retour sur la lecture des tests). Vérifions maintenant que le bruit résiduel est bien blanc… Si on regarde les autocorrélations

acf(residuals(fit.ar12),36,lwd=4)

Mais ce n’est pas un test de bruit blanc. En particulier, si une autocorrélation semblerait significative, on pourrait accepter que – globalement – les autocorrélations soient très proches de 0. On va alors faire des tests de Box-Pierce

> Box.test(residuals(fit.ar12),lag=12,
+  type='Box-Pierce')

Box-Pierce test

data:  residuals(fit.ar12)
X-squared = 7.7883, df = 12, p-value = 0.8014

Pour l’interpétation du test, c’est toujours pareil: on a une statistique de l’ordre de 7, et la loi sous-jacente (la somme des carrés des 12 premières autocorrélations du bruit) est une loi du chi-deux, à 12 degrés de liberté,

u=seq(0,30,by=.1)
plot(u,dchisq(u,df=12),type="l")
points(7.7883,0,pch=19,
col="red")

On accepte l’hypothèse que les 12 premières autocorrélations soient nulles. Si on va un peu plus loin,

> Box.test(residuals(fit.ar12),lag=18,
+  type='Box-Pierce')

Box-Pierce test

data:  residuals(fit.ar12)
X-squared = 20.3861, df = 18, p-value = 0.3115

on a la même conclusion,

Pour rappels, on peut retrouver à la main la p-value,

> 1-pchisq(20.3861,df=18)
[1] 0.3115071

Faisons notre petite dessin de toutes les p-value, pour s’assurer que le bruit est bien un bruit blanc,

 BP=function(h) Box.test(residuals(fit.ar12),lag=h,
type='Box-Pierce')$p.value
plot(1:24,Vectorize(BP)(1:24),type='b',
ylim=c(0,1),col="red")
abline(h=.05,lty=2,col="blue")

Cette fois c’est bon, on tient notre premier modèle, un https://latex.codecogs.com/gif.latex?AR(12). Si on veut aller plus loin, on peut regarder ce que donnerait de la sélection automatique,

> library(caschrono)
> armaselect(Y,nbmod=5)
p q      sbc
[1,] 14 1 1635.214
[2,] 12 1 1635.645
[3,] 15 1 1638.178
[4,] 12 3 1638.297
[5,] 12 4 1639.232

On peut être tenté d’aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant,

> fit.arma12.1=arima(Y,order=c(12,0,1))
> fit.arma12.1
Series: Y
ARIMA(12,0,1) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
      0.0301  0.1558  -0.0941  -0.1461  0.1063  -0.0688  -0.002
s.e.  0.1235  0.0854   0.0757   0.0784  0.0807   0.0774   0.080
ar9     ar10    ar11    ar12     ma1  intercept
      -0.0646  -0.0798  0.2538  0.5786  0.2231  -131.3495
s.e.   0.0802   0.0766  0.0751  0.0861  0.1393   368.8156

sigma^2 estimated as 2127759:  log likelihood=-945.65
AIC=1921.31   AICc=1926.52   BIC=1961.54

Le nouveau coefficient, https://latex.codecogs.com/gif.latex?\theta_1 est à peine significatif. Éventuellement avec un seuil à 10%… Pourquoi pas? Si on regarde les résidus, sans grande surprise, on a toujours un bruit blanc, encore plus blanc qu’auparavant (en violet sur le dessin ci-dessous)

BP=function(h) Box.test(
residuals(fit.arma12.1),lag=h,
type='Box-Pierce')$p.value
plot(1:24,Vectorize(BP)(1:24),type='b',
ylim=c(0,1),col="red")
abline(h=.05,lty=2,col="blue")

BP=function(h) Box.test(residuals(fit.ar12),lag=h,
type='Ljung-Box')$p.value
lines(1:24,Vectorize(BP)(1:24),col="purple",
type="b")

On peut aussi aller voir parmi les modèles proposés, en particulier le modèle https://latex.codecogs.com/gif.latex?AR(14).

> fit.ar14=
+ arima(Y,order=c(14,0,0),method="CSS")
> fit.ar14
Series: Y
ARIMA(14,0,0) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
      0.2495  0.2105  -0.0584  -0.1569  0.1282  -0.1152  0.0268
s.e.  0.0956  0.0972   0.0854   0.0830  0.0838   0.0840  0.0847
ar9     ar10    ar11    ar12     ar13     ar14  intercept
      -0.0327  -0.1116  0.2649  0.5887  -0.1575  -0.1572    80.5
s.e.   0.0855   0.0851  0.0853  0.0886   0.1031   0.0999   338.9

sigma^2 estimated as 2218612:  part log likelihood=-942.31

C’est un peu tiré par les cheveux, mais on pourrait accepter l’hypothèse que https://latex.codecogs.com/gif.latex?\phi_{14}soit significativement non-nuls. Mais on est encore limite…  Allez, on l’accepte aussi dans notre gang de modèle.

On a finalement trois modèles. Si on fait un peu de backtesting, sur les 12 derniers mois,

T=length(Y)
backtest=12
subY=Y[-((T-backtest+1):T)]
subtemps=1:(T-backtest)
plot(temps,Y,type="l")
lines(subtemps,subY,lwd=4)
fit.ar12.s=arima(subY,
order=c(p=12,d=0,q=0),method="CSS")
fit.arma12.1.s=arima(subY,
order=c(p=12,d=0,q=1),method="CSS")
fit.ar14.s=arima(subY,
order=c(p=14,d=0,q=0),method="CSS")
p.ar12=predict(fit.ar12.s,12)
pred.ar12=as.numeric(p.ar12$pred)
p.arma12.1=predict(fit.arma12.1.s,12)
pred.arma12.1=as.numeric(p.arma12.1$pred)
p.ar14=predict(fit.ar14.s,12)
pred.ar14=as.numeric(p.ar14$pred)

on obtient les prévisions suivantes (avec les valeurs observées dans la première colonne)

> (M=cbind(observé=Y[((T-backtest+1):T)],
+  modèle1=pred.ar12,
+  modèle2=pred.arma12.1,
+  modèle3=pred.ar14))
observé    modèle1    modèle2    modèle3
97  -4836.2174 -5689.3331 -5885.4486 -6364.2471
98  -3876.4199 -4274.0391 -4287.2193 -4773.8116
99   1930.3776  1817.8411  2127.9915  2290.1460
100  3435.1751  4089.3598  3736.1110  4039.4150
101  7727.9726  6998.9829  7391.6694  7281.4797
102  2631.7701  3456.8819  3397.5478  4230.5324
103  -509.4324 -2128.6315 -2268.9672 -2258.7216
104 -1892.6349 -3877.7609 -3694.9409 -3620.4798
105 -4310.8374 -3384.0905 -3430.4090 -2881.4942
106  2564.9600  -504.6883  -242.5018   183.2891
107 -1678.2425 -1540.9904 -1607.5996  -855.7677
108 -4362.4450 -3927.4772 -3928.0626 -3718.3922
>  sum((M[,1]-M[,2])^2)
[1] 19590931
>  sum((M[,1]-M[,3])^2)
[1] 17293716
>  sum((M[,1]-M[,4])^2)
[1] 21242230

I.e. on aurait envie de retenir le second modèle, https://latex.codecogs.com/gif.latex?ARMA(12,1). On va maintenant l’utiliser pour faire un peu de prévision,

library(forecast)
fit.arma12.1=
arima(Y,order=c(12,0,1))
fit.arma12.1
PREDARMA=forecast(fit.arma12.1,12)
plot(Y,xlim=c(1,120),type="l")
temps=T+1:12
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2],
rev(PREDARMA$upper[,2])),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1],
rev(PREDARMA$upper[,1])),col="orange",border=NA)
lines(temps,PREDARMA$mean,col="red")

Cela dit, on peut aussi aller beaucoup plus loin dans la prévision,

PREDARMA=forecast(fit.arma12.1,120)
plot(Y,xlim=c(1,210),type="l")
temps=T+1:120
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2],
rev(PREDARMA$upper[,2])),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1],
rev(PREDARMA$upper[,1])),col="orange",border=NA)
lines(temps,PREDARMA$mean,col="red")

Bon, on y est presque, car on a modélisé https://latex.codecogs.com/gif.latex?(Y_t), la série obtenue en enlevant la tendance linéaire. Pour remonter sur https://latex.codecogs.com/gif.latex?(X_t), on va rajouter la tendance à la prévision faite auparavant,

X=as.numeric(Xt)
temps=1:length(X)
plot(temps,X,type="l",xlim=c(0,210),
ylim=c(5000,30000))
base=data.frame(temps,X)
reg=lm(X~temps)
abline(reg,col="red")
PREDTENDANCE=predict(reg,newdata=
data.frame(temps=T+1:120))
temps=T+1:120
polygon(c(temps,rev(temps)),c(PREDTENDANCE+
PREDARMA$lower[,2],rev(PREDTENDANCE+PREDARMA$upper[,2])),
col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDTENDANCE+
PREDARMA$lower[,1],rev(PREDTENDANCE+PREDARMA$upper[,1])),
col="orange",border=NA)
lines(temps,PREDTENDANCE+PREDARMA$mean,col="red")

Examen final ACT6420

Chose promise, chose due: les deux bases de données qui seront utilisées pour l’examen sont désormais en ligne. 10 questions porteront sur l’analyse de la série temporelle suivante

serie=read.table("http://freakonometrics.free.fr/exam-ts.txt") 
X=serie$x XTS=serie=ts(X,start=c(1999,9),frequency=12)

(qui peut être modélisée soit comme un objet série temporel, soit comme un vecteur numérique).

Enfin, 10 questions porteront sur la base de données suivante

examen=read.table( "http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv", header=TRUE,sep=";")

(si l’importation ne marche pas, veillez à ce que l’adresse du fichier tienne sur une ligne). Cette base devait servir initialement pour l’examen intra. Elle contient des matchs de baskets universitaires (aux États-Unis). Les variables sont les suivantes

  • EQUIPE le nom de l’équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
  • WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps

La variable d’intérêt sera le nombre de points d’écart à la fin du match.

Enfin, 20 questions supplémentaires porteront sur la compréhension générale du cours. Les calculatrices seront autorisées. Maintenant, pour tous ceux qui ont encore du mal avec R, 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)

(les questions d’analyse de sortie ne porteront pas sur la programmation en R, mais davantage sur l’analyse critique des sorties, et leur interprétation).

Ordres d’un processus ARMA

Dans la méthodologie de Box & Jenkins, une étape qui arrive très rapidement est le choix des ordres du processus http://freakonometrics.blog.free.fr/public/perso6/ARMA__latex.gif, une fois que l’on validé l’hypothèse de stationnarité de la série, comme on l’a vu en cours la semaine passée. Considérons la série du trafic autoroutier,

> autoroute=read.table(
+"http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)

Un outils pratique de sélections des ordres http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif dans un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMAlatex.gif est la fonction d’autocorrélation étendue. La définition est donnée dans les notes de cours (Def. 223) à partir des statistiques proposées par Tsay & Tiao (1984),

>  EACF=eacf(A7)
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o x x x x x o o x  x  x  o
1 x x o o x x x o o o x  x  x  x
2 o x o o o o x o o o o  x  x  o
3 o x o x o o o o o o o  x  x  x
4 x x x x o o o o o o o  x  o  o
5 x x o o o x o o o o o  x  o  x
6 x x o o o x o o o o o  x  o  o
7 o x o o o x o o o o o  x  o  o
>  EACF
$eacf
           [,1]       [,2]        [,3]        [,4]
[1,]  0.6476234  0.2124105 -0.02413173 -0.24234535
[2,]  0.4889076  0.2797152 -0.02494135 -0.06094037
[3,]  0.1006541 -0.2285771  0.03514148  0.08000588
[4,]  0.1390240 -0.2788742  0.04386746  0.28307260
[5,] -0.5091680  0.3144899 -0.34572269  0.31450865
[6,] -0.4224571 -0.4877505  0.16054232 -0.09130728
[7,] -0.4731353 -0.4324857 -0.04847184 -0.10500350
[8,] -0.2129591 -0.4072901  0.09487899 -0.06493243
             [,5]        [,6]        [,7]        [,8]
[1,] -0.514330187 -0.61634046 -0.52314403 -0.28008661
[2,] -0.254912957 -0.28966664 -0.33963243 -0.21863077
[3,] -0.156624357 -0.01199786 -0.25116738  0.13079231
[4,] -0.183283544  0.03651508 -0.08711829  0.11626377
[5,] -0.190885091 -0.09786463 -0.09182557  0.08818875
[6,] -0.072904071  0.29271777 -0.09334712  0.01972648
[7,] -0.009873289  0.36909726  0.01698660 -0.03317456
[8,]  0.020485930  0.38342158  0.16981715 -0.02592442
             [,9]       [,10]       [,11]     [,12]
[1,] -0.082607058  0.15178887  0.56583179 0.8368975
[2,] -0.085026400  0.02731460  0.26158357 0.7844748
[3,]  0.001866994 -0.11658312  0.03038621 0.7026361
[4,]  0.025183793 -0.21608692  0.05660781 0.6674301
[5,] -0.006831894 -0.02514440 -0.07390257 0.4762563
[6,]  0.010058718 -0.03888613 -0.04382043 0.5338091
[7,]  0.032124905 -0.07022090 -0.04427400 0.4674165
[8,]  0.024189179 -0.20818201  0.01459933 0.4830369
          [,13]       [,14]
[1,]  0.5637439  0.17862571
[2,]  0.4530716  0.24413569
[3,]  0.2534178 -0.20160890
[4,]  0.2409861 -0.29462510
[5,] -0.1517324 -0.14763294
[6,] -0.1701182 -0.28771495
[7,] -0.1651515  0.05466457
[8,] -0.1403198 -0.04095030

$ar.max
[1] 8

$ma.ma
[1] 14

Les ronds dans la matrice désignent des valeurs non-significatives. Par défaut, le nombre de retards, pris en compte pour la composante autorégressive est faible, mais on peut l’augmenter.

> EACF=eacf(A7,13,13)
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12 13
0  x o o x x x x x o o x  x  x  o
1  x x o o x x x o o o x  x  x  x
2  o x o o o o x o o o o  x  x  o
3  o x o x o o o o o o o  x  x  x
4  x x x x o o o o o o o  x  o  o
5  x x o o o x o o o o o  x  o  x
6  x x o o o x o o o o o  x  o  o
7  o x o o o x o o o o o  x  o  o
8  o x o o o x o o o o o  x  o  o
9  x x o o x o x o o o o  x  o  o
10 x x o o x o x o o o o  x  o  o
11 x x o o x x x o o o x  x  o  o
12 x x o o o o o o o o o  o  o  o
13 x x o o o o o o o o o  o  o  o
> EACF$eacf
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,]  0.64  0.21 -0.02 -0.24 -0.51 -0.61 -0.52 -0.28 -0.08  0.15
[2,]  0.48  0.27 -0.02 -0.06 -0.25 -0.28 -0.33 -0.21 -0.08  0.02
[3,]  0.10 -0.22  0.03  0.08 -0.15 -0.01 -0.25  0.13  0.00 -0.11
[4,]  0.13 -0.27  0.04  0.28 -0.18  0.03 -0.08  0.11  0.02 -0.21
[5,] -0.50  0.31 -0.34  0.31 -0.19 -0.09 -0.09  0.08  0.00 -0.02
[6,] -0.42 -0.48  0.16 -0.09 -0.07  0.29 -0.09  0.01  0.01 -0.03
[7,] -0.47 -0.43 -0.04 -0.10  0.00  0.36  0.01 -0.03  0.03 -0.07
[8,] -0.21 -0.40  0.09 -0.06  0.02  0.38  0.16 -0.02  0.02 -0.20
[9,] -0.14 -0.50 -0.10 -0.06  0.11  0.38 -0.01 -0.02 -0.08 -0.20
[10,] -0.29  0.48  0.00  0.18  0.27 -0.12  0.41  0.00  0.16  0.15
[11,]  0.24  0.48 -0.04  0.02  0.46  0.10  0.37 -0.10  0.05  0.16
[12,] -0.59  0.49 -0.18  0.05  0.31 -0.32  0.34 -0.16  0.07  0.16
[13,] -0.31  0.31 -0.12  0.16 -0.11  0.04 -0.04  0.12 -0.09  0.00
[14,]  0.47  0.26  0.11  0.13 -0.01 -0.01  0.00  0.07 -0.03  0.02

On peut aussi visualiser graphiquement les différentes valeurs, avec en ordonnée les ordres autoréregressifs (http://freakonometrics.blog.free.fr/public/perso6/platex.gif) et en abscisse les ordres moyenne mobile (http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif).

> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> ceacf=matrix(as.numeric(cut(EACF$eacf,
+1,nrow(EACF$eacf),
+ ncol(EACF$eacf))
> for(i in 1:ncol(EACF$eacf)){
+ for(j in 1:nrow(EACF$eacf)){
+ polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
+ col=CL[ceacf[j,i]])
+ }}


Un bruit blanc est en bas à gauche (http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif nuls). Dans cette méthode, dite méthode des coins, on cherche un coin telle que dans le quadrant supérieur, les valeurs soient non-significatives (ternes sur le dessin ci-dessus). La figure ci-desssous correspond à l’analyse d’un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMA126latex.gif

Les fortes valeurs positives sont en bleu foncé, les fortes valeurs négatives sont en rouge foncé. On peut aussi regarder la fonction suivante, qui utilise une identification de modèle par MINIC (Minimum Information Criterion) à l’aide du critère de Shcwarz (BIC, ou SBC Schwarz’s Bayesian Criterion)

> armaselect(A7,nbmod=5)
p q      sbc
[1,] 12 1 1441.798
[2,] 13 0 1442.628
[3,] 12 0 1443.188
[4,] 12 2 1443.362
[5,] 14 0 1445.069

Enfin, une dernière fonction possible est évoquée dans la section 6.5. du livre deCryer & Chan (2008),

> ARMA.SELECTION=
+ armasubsets(A7,nar=14,nma=14,ar.method='ols')
> plot(ARMA.SELECTION)

basé sur le critère de Schwarz.
Je vais me répéter, mais ces méthodes ne sont que des outils, histoire d’avoir des pistes si on ne sait trop dans quelle direction partir. Et compte tenu de la saisonnalité de la série, je serais pour ma part parti sur un modèle avec http://freakonometrics.blog.free.fr/public/perso6/platex.gif=12, histoire de voir si on ne pourrait pas avoir un modèle simple, et facilement interprétable en plus…

  1. -3):3)/3,labels=1:6 []

Unit root, or not ? is it a big deal ?

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X[1]=0
for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

The idea is to assume that an autoregressive model can be considered, but we don’t know the value of the parameter. More precisely, we can’t choose if the parameter is either one (and the series is integrated), or some value strictly smaller than 1 (and the series is stationary). Based on past observations, the higher the autocorrelation, the lower the variance of the noise.

rhoest=0.9; H=260
u=241:(240+H)
P=X[240]*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)

Now that we have a model, consider the following forecast, including a confidence interval,

 
plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")
Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series,

As we can see above, assuming either that the series is stationary (parameter lower – strictly – than 1) or integrated (parameter equal to 1), the shape of the prediction can be quite different. So yes, assuming an integrated model is a big deal, since it has a strong impact on predictions.

Les tests de non-stationnarité (racine unité)

Pour commencer la modélisation d’une série temporelle, la première étape est de savoir si on peut la considérer comme stationnaire. L’alternative étant qu’elle soit intégrée. Le test le plus classique est le test de Dickey-Fuller.De manière générale, considérons un modèle autorégressive à l’ordre 1,

http://freakonometrics.blog.free.fr/public/perso6/adf-1.gif

On peut éventuellement réécrire ce processus sous la forme

http://freakonometrics.blog.free.fr/public/perso6/acf-18.gif

http://freakonometrics.blog.free.fr/public/perso6/adf-21.gif. On aura un processus non-stationnaire dès lors que http://freakonometrics.blog.free.fr/public/perso6/adf-20.gif ou encore http://freakonometrics.blog.free.fr/public/perso6/adf-16.gif. Si on n’a pas de tendance linéaire, tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif peut se faire avec un test de Student. L’idée de Dickey-Fuller est de généraliser le test de Student, en rajoutant cette tendance linéaire (pour commencer).

Afin d’illustrer l’utilisation de ce test, commençons par un cas “simple”: on va intégrer un bruit blanc simulé.

>  set.seed(1)
>  E=rnorm(240)
>  X=cumsum(E)
>  plot(X,type="l")

On peut regarder s’il y a une constante non nulle (on parlera de modèle avec drift) voire une tendance linéaire (on parlera de trend),

La première étape pourrait être de tenter un modèle sans rien.

> df=ur.df(X,type="none",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.87492 -0.53977 -0.00688  0.64481  2.47556

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared: 0.003292,	Adjusted R-squared: -0.005155
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

Value of test-statistic is: -0.7328

Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

La région critique du test est ici (pour un seuil à 5%) l’ensemble des valeurs inférieures à -1.95. Or ici la statistique de test est ici -0.73, on est alors dans la région d’acceptation du test, et on va retenir l’hypothèse de racine unité, i.e. la série n’est pas stationnaire. Mais peut-être faudrait-il juste prendre en compte une constante ?

> df=ur.df(X,type="drift",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.91930 -0.56731 -0.00548  0.62932  2.45178

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared: 0.02313,	Adjusted R-squared: 0.01482
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

Value of test-statistic is: -2.3039 2.7329

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Deux statistiques de test sont calculées, ici: la première relative à la racine unité, la seconde relative à la constante.  On observe ici que la statistique de test pour le test de racine unité (-2.3) est ici supérieure à toutes les valeurs critiques associées (données sur la première ligne). On va encore accepter l’hypothèse de racine unité. Mais le modèle était peut-être faux, et peut-être avait-on en fait en tendance linéaire ?

> df=ur.df(X,type="trend",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.87727 -0.58802 -0.00175  0.60359  2.47789

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
tt          -0.0004194  0.0009767  -0.429   0.6680
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared: 0.0239,	Adjusted R-squared: 0.01139
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

Value of test-statistic is: -1.9828 1.8771 2.7371

Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

On obtient cette fois trois statistique, la première relative au test de racine unité, et les deux suivantes aux tests sur la constante, et sur la tendance (la pente de l’ajustement linéaire). Là encore, la valeur de test (-1.98) excède les valeurs critiques: la p-value serait ici supérieure à 10%. On va là encore accepter l’hypothèse de racine unité.
Mais la modélisation autorégressive à l’ordre 1 était peut-être elle aussi fausse. Aussi, il existe un test de Dickey-Fuller augmenté. L’idée est de considérer, de manière beaucoup plus générale, un processus autorégressif à un ordre supérieur,
http://freakonometrics.blog.free.fr/public/perso6/adf-2.gif
Comme auparavant, on peut réécrire ce modèle en fonction des variations (on parlera de modèle à correction d’erreur)
http://freakonometrics.blog.free.fr/public/perso6/acf-29.gif
où on construit les nouveaux coefficients par une relation de récurrence, du genre http://freakonometrics.blog.free.fr/public/perso6/adf-6.gif. Bref, à nouveau, on est ramené à tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif . Et là encore, on a le choix sur la tendance.
Si on commence par supposer qu’il n’y a pas de tendance,

> df=ur.df(X,type="none",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.86738 -0.53887 -0.02009  0.67058  2.45970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1     -0.005168   0.007389  -0.699    0.485
z.diff.lag1 -0.029619   0.065289  -0.454    0.650
z.diff.lag2 -0.037856   0.065436  -0.579    0.563

Residual standard error: 0.9685 on 234 degrees of freedom
Multiple R-squared: 0.004704,	Adjusted R-squared: -0.008056
F-statistic: 0.3687 on 3 and 234 DF,  p-value: 0.7757

Value of test-statistic is: -0.6994

Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

la valeur de la statistique de test excède là encore les valeurs critiques, i.e. la p-value dépasse 10%. Si on rajoute une constante,

> df=ur.df(X,type="drift",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.91609 -0.56702  0.01025  0.62109  2.43970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.31105    0.13332   2.333   0.0205 *
z.lag.1     -0.03744    0.01565  -2.392   0.0175 *
z.diff.lag1 -0.01917    0.06483  -0.296   0.7677
z.diff.lag2 -0.02794    0.06496  -0.430   0.6676
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9594 on 233 degrees of freedom
Multiple R-squared: 0.02664,	Adjusted R-squared: 0.0141
F-statistic: 2.125 on 3 and 233 DF,  p-value: 0.09774

Value of test-statistic is: -2.3924 2.9709

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

on valide encore l’hypothèse de racine unité et de même avec une tendance linéaire,

> df=ur.df(X,type="trend",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.85835 -0.58826 -0.00867  0.61407  2.47280

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3530591  0.1522373   2.319   0.0213 *
z.lag.1     -0.0338831  0.0168523  -2.011   0.0455 *
tt          -0.0005674  0.0009879  -0.574   0.5663
z.diff.lag1 -0.0237595  0.0654158  -0.363   0.7168
z.diff.lag2 -0.0328215  0.0656102  -0.500   0.6174
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9608 on 232 degrees of freedom
Multiple R-squared: 0.02802,	Adjusted R-squared: 0.01126
F-statistic: 1.672 on 4 and 232 DF,  p-value: 0.1573

Value of test-statistic is: -2.0106 2.0849 3.0185

Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

On notera toutefois que ces trois modèles nous suggèrent de ne pas retenir autant de retard, qui ne semblent pas significatifs.
On peut bien entendu faire ces tests sur de vraies données par exemple le niveau du Nil

> library(datasets)
> NILE=Nile

ou encore les taux d’intérêt américains,

> base=read.table(
+ "http://freakonometrics.free.fr/basedata.txt",
+ header=TRUE)
> Y=base[,"R"]
> Y=Y[(base$yr>=1960)&(base$yr<=1996.25)]
> TAUX=ts(Y,frequency = 4, start = c(1960, 1))

Par exemple, sur cette dernière, on rejette l’hypothèse de stationnarité

>  df=ur.df(y=TAUX,lags=3,type="drift")
>  summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min      1Q  Median      3Q     Max
-3.1982 -0.2947 -0.0629  0.3524  3.1899

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.40609    0.16706   2.431 0.016362 *
z.lag.1     -0.06339    0.02500  -2.535 0.012354 *
z.diff.lag1  0.32613    0.08145   4.004 0.000102 ***
z.diff.lag2 -0.31102    0.08027  -3.875 0.000165 ***
z.diff.lag3  0.28712    0.08103   3.543 0.000541 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7805 on 137 degrees of freedom
Multiple R-squared: 0.2051,	Adjusted R-squared: 0.1819
F-statistic: 8.839 on 4 and 137 DF,  p-value: 2.227e-06

Value of test-statistic is: -2.5354 3.2457

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Mais toutes sortes d’autres tests (plus robustes) peuvent être faits. Les plus connus sont le test de Philips-Perron et le test dit KPSS. Pour ce dernier, il faut spécifier s’il l’on suppose que la série est de moyenne constante, ou si une tendance doit être prise en compte. Si l’on suppose une constante non-nulle, mais pas de tendance, on utilise

> summary(ur.kpss(X,type="mu"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 4 lags.

Value of test-statistic is: 0.972

Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

alors que pour prendre en compte une tendance

> summary(ur.kpss(X,type="tau"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 4 lags.

Value of test-statistic is: 0.5057

Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

Derrière se cache un teste du multiplicateur de Lagrange. L’hypothèse nulle correspond à l’absence de racine unité: plus la statistique de test est grande, plus on s’éloigne de la stationnarité (hypothèse nulle). Avec ces deux tests, on rejette là encore l’hypothèse de stationnarité de notre marche aléatoire simulée.

Pour le test de Philips-Perron, on a un test de type Dickey-Fuller,

> PP.test(X)

Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller = -2.0116, Trunc lag parameter = 4, p-value = 0.571

Là encore, la p-value nous recommande de valider l’hypothèse de non-stationnarité. A titre de comparaison, si on avait travaillé sur la série différenciée, on aurait accepté l’hypothèse de stationnarité

> PP.test(diff(X))

Phillips-Perron Unit Root Test

data:  diff(X)
Dickey-Fuller = -15.9522, Trunc lag parameter = 4, p-value = 0.01

Notes de cours sur les séries temporelles

La session d’hiver n’étant pas terminée, je vais poster mes notes de cours sur la dernière section (sur la modélisation de séries temporelles) pour le cours ACT2040. Il s’agit – comme je l’avais dit en cours – d’une remise au goût du jour de notes tapées il y a une dizaine d’années. J’ai également rajouté du code R, mais il doit resté un certain nombre de coquilles et de fautes de frappe. Je profiterais des jours qui viennent pour réviser cette version.

La normalité ? avant tout une histoire de taille

Un billet rapide pour préciser le fond de ma pensée lorsque j’explique aux étudiants (qui travaillent sur leur devoir) que sur les bases de données avec beaucoup d’observations, il est très dur d’accepter l’hypothèse de normalité.

Pour ça, simulons des jeux de données: on a une seule variable explicative, un bruit gaussien, homoscédastique. Mais on va mal spécifier le modèle. En particulier, le vrai modèle ne sera pas linéaire. Voilà ce qu’on simule, et que l’on va chercher à ajuster.

n=50
X1=runif(n)
epsilon=rnorm(n)
Y=-2+8*X1-4*(X1<.2)+epsilon
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)

Autrement dit, les résidus (du modèle linéaire) ne devraient pas être gaussiens, car on n’ajuste pas le bon modèle. Maintenant, mon point est le suivant: avec peu de données, la plupart des tests de normalité ne vont pas rejeter l’hypothèse de normalité. Alors qu’avec beaucoup de données, l’hypothèse de normalité sera presque toujours rejetée.
Prenons par exemple le test Lillie, évoqué dans un vieux billet, qui est une adaptation du test de Kolmogorov-Smirnov quand on compare la distribution empirique avec la loi normale dont les paramètres ont été obtenues par maximum de vraisemblance. Sur le graphique ci-dessous, à gauche, on a en abscisse la taille de l’échantillon (qui va de 10 à 1000), et en ordonnée la p-value du test Lillie. Ou plus précisément, la courbe bleue est la région où s’est trouvée la p-value dans 90% des cas (la courbe supérieure étant le quantile à 95% de la p-value sur 10,000 échantillons simulés). La courbe rouge est la valeur moyenne de la p-value. A droite, on a – toujours en fonction de la taille de l’échantillon – la proportion des échantillons simulés pour lesquels la p-value était supérieure à 5%.

Autrement dit, plus la taille de l’échantillon est grande, plus on a de chance de rejeter (fort justement d’ailleurs) l’hypothèse de normalité. Et ce, uniquement sur un effet taille. Ou pour être plus rigoureux, sur les petits échantillons, on accepte beaucoup trop souvent (à tort cette fois) l’hypothèse de normalité.
On peut aussi regarder le test de Shapiro-Wilk

ou d’Anderson-Darling,

ou encore de Cramér-von Mises

Le code R pour faire ces tests est relativement simple,

X1=runif(n)
Y=-2+8*X1-4*(X1<.2)+rnorm(n)
base=data.frame(Y,X1)
reg=lm(Y~X1,data=base)
AD[n,i]=ad.test(residuals(reg))$p.value
SW[n,i]=shapiro.test(residuals(reg))$p.value
CV[n,i]=cvm.test(residuals(reg))$p.value
LL[n,i]=lillie.test(residuals(reg))$p.value

Pour ceux qui veulent aller plus loin, se cache derrière la notion de puissance du test…

Examen intra, modèles de régression

L’examen intra comportera deux exercices proches de ceux faits en démonstration, et une analyse de sorties de régression. Cette dernière partie portera sur la base suivante (ou peut être une extraction de cette base)

> examen=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv",
+ header=TRUE,sep=";")
> tail(examen)
EQUIPE halfdiff winner
9917                 MARYLAND       -1      0
9918             RHODE ISLAND      -15      1
9919            ROBERT MORRIS        4      1
9920 NORTH CAROLINA-ASHEVILLE       11      1
9921                    BROWN       -2      1
9922          MOUNT ST. MARYS        1      0
totalpointsdiff secondhalfpoints
9917             -10               26
9918               1               44
9919               6               43
9920              19               39
9921               6               48
9922              -9               26
nodiff_pts_1000_left hometeam halfwinner
9917                    9        1          0
9918                   18        1          0
9919                   17        1          1
9920                   20        1          1
9921                   16        1          0
9922                   13        1          1
halflooser secondhalfpointsother
9917          1                    35
9918          1                    28
9919          0                    41
9920          0                    31
9921          1                    40
9922          0                    36

qui contient des matchs de baskets universitaires.

  • EQUIPE le nom de l’équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
  • WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps

Visualisation des résidus pour des données spatiales

Mardi, nous allons travailler un peu sur la modélisation du nombre d’homicides aux États-Unis, à partir de la base

> US=read.table("http://freakonometrics.free.fr/US.txt", 
+ header=TRUE,sep=";")

(je renvoie sur le précédant billet pour un descriptif précis). Idéalement, ça serait parfait si on pouvait visualiser sur une carte les variables. Pour cela, il faut rajouter une colonne à notre base, avec le nom complet des états,

> abreviation=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/etatus.csv", 
> header=TRUE,sep=",") 
> US$USPS=rownames(US) 
> US=merge(US,abreviation) 
> US$nom=tolower(US$NOM)

Cette fois, on va pouvoir faire de la cartographie, les noms de nos états étant (presque) les mêmes que ceux des cartes de R,

> library(maps) 
> VL0=strsplit(map("state")$names,":") 
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,US$nom)

Cette fois-ci, on a toutes les informations pour faire une carte, avec une couleur fonction de la variable d’intérêt (espérance de vie à la naissance, taux d’homicides, taux illettrisme, etc).

> library(RColorBrewer) 
> carte=function(V=US$Murder,titre= 
+ "Taux d'homicides aux Etats-Unis"){ 
+ variable=as.numeric(as.character(cut(V, 
+ quantile(V,seq(0,1,by=1/6)),labels=1:6))) 
+ niveau=variable[ETAT] 
+ couleur=rev(brewer.pal(6,"RdBu")) 
+ noml=levels(cut(V,quantile(V,seq(0,1,by=1/6)))) 
+ map("state", fill = TRUE, col=couleur[niveau]); 
+ legend(-78,34,legend=noml,fill=couleur, + cex=1,bty="n"); 
+ title(titre)}

Commençons par l’analyse du nombre de jours durant lesquels la température passe en dessous de 0°C, par an, en moyenne, dans la capitale (ou la plus grande ville) de l’état, afin de tester notre fonction,

> carte(US$Frost, + titre="Nombre de jours de gel par an")

Pour le taux d’homicide (qui est la variable par défaut) on a

> carte()

Ça sera notre variable d’intérêt lors de la modélisation de mardi. Enfin, on peut lancer une régression, et représenter spatialement les résidus,

> reg=lm(Murder~.-NOM-USPS-nom,data=US) 
> regs=step(reg) 
> carte(residuals(regs), 
+ titre="Résidus de la régression")

Nous voila équipés pour commencer l’économétrie spatiale…

Pour aller maintenant un peu plus loin dans la modélisation, je vais rajouter une variable qualitative, par exemple l’appartenance politique du gouverneur en 1977 (les données datent de cette époque). Les données ont été extraites de wikipedia, suite aux élections de 1974, 1975, 1976 et 1977,

> GV=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/governor.csv", 
+ header=TRUE,sep=";") 
> etat=strsplit(as.character(GV$State),"-") 
> listeetat=rep(NA,nrow(GV)) 
> for(i in 1:nrow(GV)){ 
+ listeetat[i]=etat[[i]][1] 
+ } 
> indice=which(is.na(listeetat)==FALSE) 
> basegv=data.frame(state=tolower(listeetat[indice]), 
+ party=GV$Party[indice])

On a la visualisation suivante

> library(maps) 
> library(RColorBrewer) 
> couleur=rev(brewer.pal(6, "RdBu")) 
> Z=rep(6,length(basegv$party)) 
> Z[basegv$party=="Democratic"]=1 
> VL0=strsplit(map("state")$names,":")
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,basegv$state)
> niveau=Z[ETAT] 
> map("state", fill = TRUE, col=couleur[niveau])

For those who think more variates should be added to the dataset, some can be found e.g. on http://www.statemaster.com/, like the total executions since 1930, or the date the state joint the U.S.A.