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).
L'examen de mercredi est en ligne ici (les données et les annexes étaient dans un précédant billet), et je poste des éléments de correction. Une correction avec les statistiques de réponses sera en ligne dans les jours qui viennent (le temps de corriger). Si quelqu'un voit des coquilles dans la correction, voire des réponses erronées, merci de me le dire avant que je ne valide les notes.
L'examen final du cours de modèles de prévision aura lieu dans 15 jours. Comme annoncé ce matin, la forme sera proche de celle de l'examen intra,
Pour la fin de cours , quelques transparents sur l'identification et l'estimation des modèles SARIMA, quelques compléments sur les tests (racines unités et non-stationnarité, ainsi que saisonnalité), et enfin, quelques pistes pour construire des prédictions (avec une quantification de l'incertitude). Les transparents sont en ligne 
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 
Nous allons commencer demain la modélisation des séries temporelles. Les transparents de la séance sont en ligne 
On célèbre cette année les 35 ans d'un livre qui a révolutionné la
statistique, tant sur les aspects computationnels que graphiques: Exploratory Data Analysis par John Tukey. La légende prétend (on pourra relire 
Pour les étudiants en groupes de deux, merci de mettre le coéquipier en copie du courriel. Je rappelle que des instructions sont toujours en ligne sur un vieux 
..
