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)
(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).
soit stationnaire, de moyenne 
. Alors en prenant l’espérance de part et d’autre dans l’expression précédente, 
, i.e. 
, ou 
. Plus généralement, avec un ARMA, 
, en prenant la encore l’espérance des deux cotés, on en déduit 
.
(afin de se débarrasser de cette tendance). Alors
ou encore 
suit alors un processus ARIMA(1,1,0) sans constante cette fois.
avec
dans le modèle de
intégrable, alors pour tout 
et 
, 
. La conséquence de cette propriété (qui est juste) est que 
dès lors que la transformation n’est pas linéaire (ou au moins, il n’y a aucune raison pour qu’il y ait égalité).
(définies sur le même espace), 
, pour tout 
, et que l’on simule les couples 
, alors
sera grande. Donc effectivement, l’approximation pourrait être relativement mauvaise si on est sur des couples présentant une forte dépendance.
