# 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.

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

on the variable of interest. Then assume that

As mentioned in Chapter 14 of Davidson & MacKinnon (1993) – in French – the log-likelihood of this model (assuming that observations are independent, with distribution ) can be written

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 matrix (since we assume that  is a full-rank matrix). More precisely,  where  is a  matrix,  is a  orthonornal matrix, and is a  upper triangle matrix. It might be convenient to use this matrix since, for instance, .  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, > 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
 0.430591

$objective  -197.6966 > abline(v=maxf$maximum,lty=2) The good point is that the optimal value of  is the same as the one we got before. The only problem is that the -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)  0.1383507 > (upper=uniroot(f, c(maxf$maximum,1))$root)  0.780573 > segments(lower,ic,upper,ic,lwd=2,col="red") Actually, it possible to rewrite the log-likelihood as (let us just get rid of the constant), where 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  0.430591$objective
 -47.73436

with again the same optimal value for , 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 . # Qui sort en avance lors des examens?

Tous les profs qui surveillent leur examen se posent la même question: que signifie un départ au bout d’une heure et demi, si l’examen est supposé durer trois heures? Que l’examen était trop facile? trop dur?

Loin de moi l’idée de blâmer les étudiants qui partent tôt ! Tout d’abord parce que j’ai toujours souvent été dans les premiers à quitter la salle lors d’examens, ou de concours. Et parce que je ne pense pas qu’il y ait un lien évident entre le fait de partir tôt, ou tard, et la note obtenue.

Mais on peut regarder (et répondre au commentaire de @kl suite à l’examen de la session passée). Car l’examen du cours ACT6420 avait lieu la semaine passée, et il était court. Autrement dit, 95% des étudiants étaient sortis avant que j’annonce la fin de l’épreuve. Et j’ai noté sur toutes les copies l’heure de sortie (pour ceux qui veulent plus d’information sur le protocole de l’expérience, les étudiants ne savaient pas que je notais l’heure).

Il s’agissait d’un examen à choix multiples, ce qui garantit une forme d’impartialité dans la correction (sinon on pourrait objecter que si je note en prenant le paquet de copie tel qu’il est arrivé, les notes pour les élèves partis tôt – corrigés en dernier – et ceux partis à la fin – corrigés en premier – pourraient ne pas être vraiment comparables). Et comme à la dernière session, je demandais aux étudiants de prédire leur nombre de réponses. La question était facultative, mais les élèves pouvaient gagner un point en répondant (à condition de prédire exactement le bon nombre).

Les données sont (partiellement) en ligne ici

> act6420=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/baseact6420a1.txt",
> names(act6420)=c("n","heure","obs","pred")

On peut alors comparer les notes en fonction de l’heure de sortie, ainsi que la note prédites (les notes sont sur 30, ou plutôt les nombres de bonnes réponses), en faisant des régressions linéaires simples,

> summary(lm(obs~heure,data=act6420))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1519     6.0546   0.521   0.6039
heure         1.3209     0.5391   2.450   0.0162 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.819 on 93 degrees of freedom
Multiple R-squared: 0.06063,	Adjusted R-squared: 0.05053
F-statistic: 6.002 on 1 and 93 DF,  p-value: 0.01616

> summary(lm(pred~heure,data=act6420))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.1178     8.8373   0.126   0.8996
heure         1.7234     0.7889   2.185   0.0317 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.999 on 85 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.05316,	Adjusted R-squared: 0.04202
F-statistic: 4.773 on 1 and 85 DF,  p-value: 0.03167 

On observe qu’à la fois les bons élèves (ceux avec une bonne note) et les élèves qui se pensent bons (ceux qui ont prédit une note élevée) sortent plutôt tard. Avec dans les deux cas,  une significativité de la variable heure de sortie sur la prédiction de la note. Mais peut-être est-ce plus complexe que ça, et qu’il se cache un effet non-linéaire. On essayer de lisser un peu, avec des splines

> D=data.frame(Y=note,X=heure)
> library(splines)
> rg=lm(Y~bs(X,df=5),data=D)
> newheure=seq(10.2,12.1,by=.01)
> notep=predict(rg,newdata=data.frame(X=
+ newheure),interval="confidence")

On obtient les régressions (non-linéaires) suivantes

ou, si on souhaite rajouter des intervalles de confiance autour de la prédiction faite par nos modèles,

(avec en abscisse l’heure de sortie, l’examen finissant un peu après midi). L’analyse est un peu plus complexe qu’il n’y paraissait, avec quelques bons élèves qui finissent très tôt. Si on souhaite d’ailleurs mieux comprendre le lien entre la note obtenue (ou plutôt, encore une fois le nombre de bonnes réponses données, car aucun bonus n’est pris en compte ici), on peut faire un graphique proche de celui fait à la session passée

Le trait bleu est ici la régression obtenue à la session passée. On retrouve exactement la même chose sur les deux sessions, étonnant, non?

Mais j’arrête là pour l’analyse de ces notes. Pour ceux qui veulent en savoir plus, je mettrais bientôt en ligne une correction, avec des statistiques par réponse. On en profitera pour faire une discussion sur l’utilisation de règles démocratiques, à savoir, la majorité a-t-elle toujours raison? Et pour ceux qui n’auraient pas encore compris, on a fait ici une analyse de corrélations entre les heures de sortie, et les notes, et non pas une analyse de causalité. Donc inutile de se forcer à attendre la fin de l’examen pour sortir la prochaine fois, ça ne fera a pas (a priori) monter la note.

# 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",
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 -value de l’ordre de 1.5%, on va rejeter l’hypothèse  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 -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 (5ans) (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. # 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 Cretois 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…

# Transformation logarithmique de séries temporelles

Pour poursuivre une discussion amorcée en fin de cours, dans certains cas, on peut avoir l’impression que modéliser une série pourrait être compliqué,

plot(X,xlim=c(1,length(X)+20)) Mais on peut avoir l’intuition que modéliser le logarithme de la série pourrait être plus simple,

> X=log(Y)
> plot(X,xlim=c(1,length(X)+20)) On va alors tenter une modélisation par un processus ARMA de cette dernière série,

> md=arima(X,c(12,0,1))
> P=predict(md,24)
> E=P$pred > V=P$se^2

On peut alors faire une prévision sur cette série plus simple à modéliser, et visualiser cette prévision.

> temps=length(X)+1:24
> ciu=(E+2*sqrt(V))
> cil=(E-2*sqrt(V))
> polygon(c(temps,rev(temps)),c(ciu,
+ rev(cil)),col="yellow",border=NA)
> lines(temps,E,col="red",lwd=2)
> lines(temps,ciu,col="red",lty=2)
> lines(temps,cil,col="red",lty=2) Maintenant, on va devoir remonter. On va utiliser un résultat que l’on a vu sur la transformation logarithmique dans une régression: si après avoir pris le logarithme, on a un modèle simple, Gaussien, c’est que le modèle initial était log-normal. On peut alors utiliser les propriétés de la loi lognormale, dont on connait les moments à partir de ceux de la loi Gaussienne sous-jacente. Pour la prévision, on n’a pas trop le choix,

> mu=exp(E+.5*V)

Par contre, pour construire un intervalle de confiance, soit on utilise la variance de notre loi lognormale pour avoir la variance de notre processus, et on oublie cette histoire de loi lognormale pour construire un intervalle Gaussien,

> sig2=(exp(V)-1)*exp(2*E+V)
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)

ou alors on utilise le fait que comme la transformation est monotone, l’intervalle de confiance peut etre vu comme une transformation du précédant intervalle de confiance,

> ci2u=exp(E+2*sqrt(V))
> ci2l=exp(E-2*sqrt(V))

Si on compare visuellement les deux, on a dans le premier cas,

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> polygon(c(temps,rev(temps)),c(ci1u,
> rev(ci1l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci1u,col="red",lty=2)
> lines(temps,ci1l,col="red",lty=2) (qui est symétrique et centré sur notre prévision) et dans le second > plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)
> polygon(c(temps,rev(temps)),c(ci2u,
> rev(ci2l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci2u,col="red",lty=2)
> lines(temps,ci2l,col="red",lty=2) Pour le devoir de série temporelles, les données fournies par Google Insight sont hebdomadaires, ce qui peut rendre la modélisation compliquée. Comme on l’a évoqué en fin de cours, il peut être plus simple de travailler sur des données mensuelles. La petite fonction suivante permet de transformer les données pour avoir des données mensuelles (avec des moyennes par mois, ce qui fait qu’un mois de 28 jours et un mois de 31 jours sont comparables). Histoire d’illustrer, prenons un mot au hasard, disons…. épilation. On commence par sauver le fichier csv, et on le lit sous R,

EPILATION=read.table("/home/charpentier/report-epilation.csv", skip=4,header=TRUE,sep=",",nrows=426)

La petite fonction suivant va nous aider à convertir la base en une série temporelle mensuelle,

H2M=function(BASE){
X=BASE[,2]
date1=substr(as.character(BASE$Semaine),1,10) date2=substr(as.character(BASE$Semaine),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(Z,start=c(2004,1),frequency=12)
return(Zts)}

Si on travaille sur la série hebdomadaire, on a la série suivante

hebdo=ts(EPILATION\$épilation,start=2004, frequency=52) La fonction d’autocorrélation est alors

acf(hebdo,160) Maintenant, on peut travailler sur les données mensuelles. On utilise alors

mensuel=H2M(EPILATION)

Le graphique est alors le suivant La fonction d’autocorrélation est cette fois la suivante

acf(mensuel,40) On retrouve le comportement cyclique, avec une saisonnalité annuelle. Mais avec 12 retards, on devrait avoir des modèles plus simples qu’avec 52 retards. Bref, il peut être plus simple de travailler sur des données mensuelles qu’hebdomadaires. Mais chacun est libre de la série qu’il ou elle analysera…