Une ébauche de correction pour les 50 questions à choix multiples de la semaine passée est en ligne (avec les statistiques de réponses). Comme je le rappelle en introduction, il s’agit d’un extrait d’un vrai examen donné en début d’année. Merci de me dire rapidement si vous n’êtes pas d’accord avec ma correction.
Category Archives: 2014-E
Variance of the Average of a Sequence
In the case where are i.i.d. random variables, then
Now, what if are identically distributed, but no longer independent. What if we have an autoregressive process? Assume that
Then
can be written
Here, we will express the variance as a function of and
, but it is possible to use also
, since, in the context of an
,
Now, since we get
which can be simplified, since
i.e.
So, the variance of the mean can be writen as
Observe that if is large enough,
This asymptotic relationship is well known actually. A simple way to get it is the following. One can can write
or equivalently
But actually, the first relationship is probably more interesting to get an asymptotic approximation,
In the context of an process, this can be writen
Thus, we get the following well-known relationship
In the case where is an i.i.d. sequence, i.e.
, then we get the relationship mentioned initially. And in the case of a random walk… unfortunately, we cannot use that relationship. But observe that
i.e.
which can be written
If we compare the true value and the approximation, we get the following graph,
> V=function(phi,s2=1,n=100){
+ g0=s2/(1-phi^2)
+ if(phi<1){
+ if(phi==0){v1=g0/n}
+ if(phi>0){v1=g0/n^2*(n+2*((n-1)*
+ phi^(-1)-n+phi^(n-1))/(phi^(-1)-1)^2)}
+ v2=g0/n*(1+phi)/(1-phi)
+ }
+ if(phi==1){
+ v1=(2*n+1)*(n+1)*s2/(6*n)
+ v2=NA
+ }
+ return(c(v1,v2))}
>
> Vphi=function(phi) V(phi,1,100)
> x=seq(.01,1,by=.02)
> M=matrix(unlist(lapply(x,V)),nrow=2)
> plot(x,M[1,],type="l",col="red",log="y",
+ ylab="Variance of the average (log scale)",
+ xlab="Autoregressive coefficient")
> lines(x,M[2,],col="blue")
Tukey and Mosteller’s Bulging Rule (and Ladder of Powers)
When discussing transformations in regression models, I usually briefly introduce the Box-Cox transform (see e.g. an old post on that topic) and I also mention local regressions and nonparametric estimators (see e.g. another post). But while I was working on my ACT6420 course (on predictive modeling, which is a VEE for the SOA), I read something about a “Ladder of Powers Rule” also called “Tukey and Mosteller’s Bulging Rule“. To be honest, I never heard about this rule before. But that won’t be the first time I learn something while working on my notes for a course !
The point here is that, in a standard linear regression model, we have
But sometimes, a linear relationship is not appropriate. One idea can be to transform the variable we would like to model, , and to consider
This is what we usually do with the Box-Cox transform. Another idea can be to transform the explanatory variable, , and now, consider,
For instance, this year in the course, we considered – at some point – a continuous piecewise linear functions,
It is also possible to consider some polynomial regression. The “Tukey and Mosteller’s Bulging Rule” is based on the following figure.
and the idea is that it might be interesting to transform and
at the same time, using some power functions. To be more specific, we will consider some linear model
for some (positive) parameters and
. Depending on the shape of the regression function (the four curves mentioned on the graph above, in the four quadrant) different powers will be considered.
To be more specific, let us generate different models, and let us look at the associate scatterplot,
> fakedataMT=function(p=1,q=1,n=99,s=.1){
+ set.seed(1)
+ X=seq(1/(n+1),1-1/(n+1),length=n)
+ Y=(5+2*X^p+rnorm(n,sd=s))^(1/q)
+ return(data.frame(x=X,y=Y))}
> par(mfrow=c(2,2))
> plot(fakedataMT(p=.5,q=2),main="(p=1/2,q=2)")
> plot(fakedataMT(p=3,q=-5),main="(p=3,q=-5)")
> plot(fakedataMT(p=.5,q=-1),main="(p=1/2,q=-1)")
> plot(fakedataMT(p=3,q=5),main="(p=3,q=5)")
If we consider the South-West part of the graph, to get such a pattern, we can consider
or more generally
where and
are both larger than 1. And the larger
and/or
, the more convex the regression curve.
Let us visualize that double transformation on a dataset, say the cars dataset.
> base=cars
> names(base)=c("x","y")
> MostellerTukey=function(p=1,q=1){
+ regpq=lm(I(y^q)~I(x^p),data=base)
+ u=seq(min(min(base$x)-2,.1),max(base$x)+2,length=501)
+ par(mfrow=c(1,2))
+ plot(base$x,base$y,xlab="X",ylab="Y",col="white")
+ vic=predict(regpq,newdata=data.frame(x=u),interval="prediction")
+ vic[vic<=0]=.1
+ polygon(c(u,rev(u)),c(vic[,2],rev(vic[,3]))^(1/q),col="light blue",density=40,border=NA)
+ lines(u,vic[,2]^(1/q),col="blue")
+ lines(u,vic[,3]^(1/q),col="blue")
+ v=predict(regpq,newdata=data.frame(x=u))^(1/q)
+ lines(u,v,col="blue")
+ points(base$x,base$y)
+
+ plot(base$x^p,base$y^q,xlab=paste("X^",p,sep=""),ylab=paste("Y^",q,sep=""),col="white")
+ polygon(c(u,rev(u))^p,c(vic[,2],rev(vic[,3])),col="light blue",density=40,border=NA)
+ lines(u^p,vic[,2],col="blue")
+ lines(u^p,vic[,3],col="blue")
+ abline(regpq,col="blue")
+ points(base$x^p,base$y^q)
+ }
For instance, if we call
> MostellerTukey(2,1)
we get the following graph,
On the left, we have the original dataset, and on the right, the transformed one,
, with two possible transformations. Here, we did only consider the square of the speed of the car (and only one component was transformed, here). On that transformed dataset, we run a standard linear regression. We add, here, a confidence tube. And then, we consider the inverse transformation of the prediction. This line is plotted on the left. The problem is that it should not be considered as our optimal prediction, since it is clearly biased because
. But quantiles associated with a monotone transformation are the transformed quantiles. So confidence tubes can still be considered as confidence tubes.
Note that here, it could have be possible to consider another transformation, with the same shape, but quite different
> MostellerTukey(1,.5)
Of course, there is no reason to consider a simple power function, and the Box-Cox transform can also be used. The interesting point is that the logarithm can be obtained as a particular case. Furthermore, it is also possible to seek optimal transformations, seen here as a pair of parameters. Consider
> p=.1 > bc=boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$y > for(p in seq(.2,3,by=.1)) bc=cbind(bc,boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$y) > vp=boxcox(y~I(x^p),data=base,lambda=seq(.1,3,by=.1))$x > vq=seq(.1,3,by=.1) > library(RColorBrewer) > blues=colorRampPalette(brewer.pal(9,"Blues"))(100) > image(vp,vq,bc,col=blues) > contour(vp,vq,bc,levels=seq(-60,-40,by=1),col="white",add=TRUE)
The darker, the better (here the log-likelihood is considered). The optimal pair is here
> bc=function(a){p=a[1];q=a[2]; as.numeric(-boxcox(y~I(x^p),data=base,lambda=q)$y[50])}
> optim(c(1,1), bc,method="L-BFGS-B",lower=c(0,0),upper=c(3,3))
$par
[1] 0.5758362 0.3541601
$value
[1] 47.27395
and indeed, the model we get is not bad,
Fun, ins’t it?
Google, et prévision
Pour le second devoir du cours ACT6420 (modèles de prévisions), le but est de prévoir des recherches sur Google, via https://www.google.com/trends/. Soit vous avez un mot clé qui vous intéresse, soit vous allez chercher une série extraite dans /ACT6420-TS/. Dans les bases qui ont été mises en ligne, il y a le mot clé gym, par exemple.
> report=read.table( + "http://freakonometrics.free.fr/ACT6420-TS/report-GYM.csv", + skip=4,header=TRUE,sep=",",nrows=548)
On va nettoyer un peu la base (en particulier les valeurs manquantes de la fin)
> tail(report)
Semaine gym
543 2014-05-25 - 2014-05-31 80
544 2014-06-01 - 2014-06-07 80
545 2014-06-08 - 2014-06-14 78
546 2014-06-15 - 2014-06-21 NA
547 2014-06-22 - 2014-06-28 NA
548 2014-06-29 - 2014-07-05 NA
> report=report[!is.na(report[,2]),]
> tail(report)
Semaine gym
540 2014-05-04 - 2014-05-10 79
541 2014-05-11 - 2014-05-17 80
542 2014-05-18 - 2014-05-24 79
543 2014-05-25 - 2014-05-31 80
544 2014-06-01 - 2014-06-07 80
545 2014-06-08 - 2014-06-14 78
Les données sont ici hebdomadaires, comme on peut le voir graphiquement
> hebdo=ts(report$gym,start=2004,frequency=52) > hebdo Time Series: Start = c(2004, 1) End = c(2014, 25) Frequency = 52 [1] 68 60 60 53 55 50 49 49 51 48 48 45 45 [14] 47 42 48 46 47 46 47 47 46 47 46 45 46 [27] 46 50 48 48 52 51 57 55 53 56 55 50 48 [40] 50 46 49 46 48 49 48 46 50 47 46 43 54 [53] 69 64 63 60 57 57 53 54 55 54 50 53 54 [66] 46 50 50 49 49 49 47 49 48 49 50 49 49 [79] 47 47 50 49 52 51 55 55 55 54 54 52 53 [92] 54 53 52 51 51 50 51 48 52 50 47 45 56 [105] 76 72 66 64 63 59 53 56 57 58 54 55 54 [118] 53 52 52 50 53 50 51 49 51 51 50 50 48 [131] 48 53 52 56 58 57 60 62 62 63 59 58 58 [144] 56 54 54 53 53 54 53 50 55 54 53 51 56 [157] 77 73 68 68 67 66 64 67 64 63 63 63 62 [170] 62 61 62 61 63 62 63 63 63 63 59 58 59 [183] 61 60 60 58 61 61 62 62 64 68 66 61 58 [196] 58 55 54 51 55 54 55 53 55 53 52 50 55 [209] 76 77 68 67 64 64 58 61 59 59 57 55 57 [222] 59 57 54 56 55 54 52 52 53 54 53 55 55 [235] 55 57 54 56 55 58 65 63 64 67 66 63 62 [248] 60 60 57 55 56 56 58 58 53 56 55 54 52 [261] 69 77 71 68 66 63 60 60 62 59 59 57 57 [274] 60 58 58 56 54 58 57 56 57 57 57 57 54 [287] 54 55 57 57 56 64 60 59 62 62 64 59 58 [300] 57 57 54 53 52 53 53 55 53 55 53 50 49 [313] 63 76 73 68 65 66 60 61 61 60 58 58 61 [326] 61 62 57 57 58 55 58 57 58 59 57 55 55 [339] 57 57 58 59 59 60 60 63 63 63 66 65 62 [352] 60 59 58 57 56 58 59 56 54 57 54 54 53 [365] 66 87 77 74 72 69 68 64 65 65 68 63 65 [378] 65 65 62 61 62 62 63 63 61 65 63 64 63 [391] 61 62 62 61 62 65 63 67 67 71 74 71 70 [404] 68 68 65 66 64 65 68 64 64 65 62 61 60 [417] 69 91 88 83 81 78 75 71 73 74 73 70 69 [430] 66 68 69 66 68 68 65 69 66 69 70 69 70 [443] 69 72 72 71 69 76 74 72 77 77 82 78 72 [456] 72 69 68 67 67 64 66 66 63 65 64 62 61 [469] 67 88 90 83 81 83 77 76 76 74 75 74 74 [482] 77 77 77 73 72 76 72 71 72 74 72 74 73 [495] 73 76 73 73 71 76 76 79 79 83 83 81 78 [508] 78 76 78 80 74 73 75 74 75 72 71 70 69 [521] 73 92 100 94 93 91 86 84 84 85 85 83 82 [534] 83 83 78 79 80 80 79 80 79 80 80 78 > plot(hebdo)
Pour avoir des modèles plus simples, on peut agréger les données, pour les rendre mensuelles (par interpolation linéaire). La fonction à utiliser est ici
H2M=function(BASE){
X=BASE[,2]
Y=BASE[,1]
date1=substr(as.character(Y),1,10)
date2=substr(as.character(Y),14,23)
D1=as.Date(date1,"%Y-%m-%d")
D2=as.Date(date2,"%Y-%m-%d")
vm=vy=N=NA
for(t in 1:length(D1)){
mois=seq(D1[t],D2[t],length=7)
vm=c(vm,as.POSIXlt(mois)$mon+1)
vy=c(vy,as.POSIXlt(mois)$year+1900)
N=c(N,rep(X[t],7))}
N=N[-1]; vm=vm[-1]; vy=vy[-1]
YM=vy*100+vm
Z=tapply(N,as.factor(YM),mean)
Zts=ts(as.numeric(Z),start=c(2004,1),frequency=12)
return(Zts)}
Si on utilise cette fonction sur nos données, on obtient
> mensuel=H2M(report)
> mensuel
Jan Feb Mar Apr May Jun
2004 60.25000 50.75862 47.51613 45.66667 46.67742 46.00000
2005 63.22581 55.10714 53.03226 49.10000 48.45161 49.10000
2006 68.87097 57.10714 55.51613 51.86667 50.70968 49.93333
2007 70.74194 65.57143 62.87097 61.60000 62.77419 59.96667
2008 70.45161 60.79310 57.19355 56.13333 52.96774 54.30000
2009 70.35484 61.25000 58.19355 57.13333 56.80645 56.00000
2010 69.87097 61.78571 59.45161 58.76667 57.16129 56.40000
2011 76.58065 66.21429 65.22581 62.66667 62.51613 63.16667
2012 85.00000 73.82759 69.93548 67.76667 67.25806 69.46667
2013 84.93548 76.39286 75.00000 74.80000 72.67742 73.13333
2014 94.29032 84.96429 83.29032 79.80000 79.54839 79.00000
Jul Aug Sep Oct Nov Dec
2004 47.80645 53.67742 52.63333 47.77419 47.96667 47.61290
2005 48.41935 53.51613 53.43333 52.41935 50.20000 49.74194
2006 52.48387 59.77419 59.66667 54.12903 52.86667 54.35484
2007 59.87097 62.03226 63.10000 54.45161 54.30000 54.09677
2008 55.45161 62.16129 63.96667 57.41935 56.23333 56.09677
2009 55.96774 61.12903 60.16667 54.29032 53.60000 53.35484
2010 58.12903 61.64516 63.43333 57.67742 56.73333 56.48387
2011 61.80645 66.22581 70.86667 65.77419 65.23333 63.19355
2012 71.48387 75.06452 75.80000 67.22581 64.90000 65.12903
2013 73.51613 78.93548 79.73333 76.41935 73.93333 72.80645
2014
> ts.plot(mensuel)
Cette base est maintenant utilisable pour le devoir. Le but est ici de faire de la prévision pour les 2 prochaines années, avec un intervalle de confiance. Mais j’en reparlerais par courriel.
Prévision de séries chronologiques
Dans la seconde partie du cours de modèles de prévision, on quittera (un peu) les données individuelles pour parler de données chronologiques. Les slides de cette semaine (et probablement la semaine prochaine) sont en ligne.
J’ai mis en ligne, en parallèle, quelques notes de cours sur les séries temporelles, qui pourront peut être servir de complément. Comme le disait Doug Martin “Time series is the worst subject to teach. First, you have to teach the standard theory. Then, if you are beging honest, you have to tell the students ‘none of this stuff works, and this is what people really do'” (cité la semaine passée par David J. Thomson, au congrès de la SSC). On essayera de garder ça en mémoire tout au long du cours !
Modèles de prévision, fourre-tout
Dans le dernier cours de modèle de prévision, la semaine passée, nous avions passé un peu de temps sur l’étude des points aberrants et des points influents. Tout est expliqué dans les slides (avec les codes) donc je ne reviendrais pas dessus. Je pourrais juste évoquer quelques lignes de codes utilisées pour voir l’impact d’une observation sur la régression, en enlevant l’observation de la base, et en regardant ce que ça donne sur la prévision
> plot(cars) > text(cars[c(23,49),1],cars[c(23,49),2]+4,c(23,49)) > abline(lm(dist~speed,data=cars)) > abline(lm(dist~speed,data=cars[-23,]),col="red") > abline(lm(dist~speed,data=cars[-49,]),col="blue")
On a ensuite vu les problèmes de choix de modèles, et les méthodes de pénalisation (pénalisation du R2, et de la log-vraisemblance). On a ensuite vu les méthodes automatiques de choix de modèle,
> US=read.table("http://freakonometrics.free.fr/US.txt",sep=";") > reg=lm(Murder~.,data=US) > regs=step(reg) Start: AIC=63.01 Murder ~ Population + Income + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - Income 1 0.236 128.27 61.105 - HS.Grad 1 0.973 129.01 61.392 <none> 128.03 63.013 - Area 1 7.514 135.55 63.865 - Illiteracy 1 8.299 136.33 64.154 - Frost 1 9.260 137.29 64.505 - Population 1 25.719 153.75 70.166 - Life.Exp 1 127.175 255.21 95.503 Step: AIC=61.11 Murder ~ Population + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - HS.Grad 1 0.763 129.03 59.402 <none> 128.27 61.105 - Area 1 7.310 135.58 61.877 - Illiteracy 1 8.715 136.98 62.392 - Frost 1 9.345 137.61 62.621 - Population 1 27.142 155.41 68.702 - Life.Exp 1 127.500 255.77 93.613 Step: AIC=59.4 Murder ~ Population + Illiteracy + Life.Exp + Frost + Area Df Sum of Sq RSS AIC <none> 129.03 59.402 - Illiteracy 1 8.723 137.75 60.672 - Frost 1 11.030 140.06 61.503 - Area 1 15.937 144.97 63.225 - Population 1 26.415 155.45 66.714 - Life.Exp 1 140.391 269.42 94.213
Enfin, on a aussi vu le passage au logarithme sur la variable d’intérêt, avec un modèle log-normal,
> plot(cars$speed,cars$dist) > regln=lm(log(dist)~speed,data=cars) > summary(regln) Call: lm(formula = log(dist) ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -1.46604 -0.20800 -0.01683 0.24080 1.01519 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.67612 0.19614 8.546 3.34e-11 *** speed 0.12077 0.01206 10.015 2.41e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4463 on 48 degrees of freedom Multiple R-squared: 0.6763, Adjusted R-squared: 0.6696 F-statistic: 100.3 on 1 and 48 DF, p-value: 2.413e-13 > beta=coefficients(regln) > for(x in seq(2,30,by=.1)){ + zhat=beta[1]+beta[2]*x + yhat=exp(zhat+1/2*0.4463^2) + points(x,exp(zhat),col="red") + points(x,yhat,col="blue") + } > abline(lm(dist~speed,data=cars))
Pour comparer les modèles (linéaire versus log-linéaire), on avait suggéré de comparer les résidus des deux modèles,
> U=cars$dist - exp(predict(regln)+1/2*0.4463^2) > sd(U) [1] 17.12949 > sd(residuals(lm(dist~speed,data=cars))) [1] 15.22184
Ce soir, on fini cette première partie du cours, sur les données individuelles.
Régression linéaire, quelques codes
Un rapide billet pour mettre en ligne les codes utilisés la semaine passée, complétant les codes des transparents. On travaille toujours sur la même base, ou on cherche à prévoir une distance de freinage d’un véhicule, tenant compte de la vitesse du véhicule.
> plot(cars) > reg=lm(dist~speed,data=cars) > summary(reg) Call: lm(formula = dist ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -29.069 -9.525 -2.272 9.215 43.201 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 1.49e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.38 on 48 degrees of freedom Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438 F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12
Pour faire plusieurs prévisions, à la main, on peut utiliser le code suivant (la boucle permet de faire des prévisions pour plusieurs valeurs)
> for(x in seq(3,30,by=.25)){ + b0=coef(reg)[1] + b1=coef(reg)[2] + Yx=b0+b1*x + V=vcov(reg) + Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2 + IC1=Yx+c(-1,+1)*1.96*sqrt(Vx) + s=summary(reg)$sigma + IC2=Yx+c(-1,+1)*1.96*s + points(x,Yx,pch=19,col="red") + points(c(x,x),IC1,pch=3,col="blue") + points(c(x,x),IC2,pch=3,col="purple")}
On avait ensuite fait une régression linéaire sur une sous-base, avec 20 observations tirées au hasard
> I=sample(1:50,size=20) > reg=lm(dist~speed,data=cars[I,])
Le but était de visualiser l’impact du nombre d’observation sur la qualité de la régression
> summary(reg) Call: lm(formula = dist ~ speed, data = cars[I, ]) Residuals: Min 1Q Median 3Q Max -23.529 -7.998 -5.394 11.634 39.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -20.7408 9.4639 -2.192 0.0418 * speed 4.2247 0.6129 6.893 1.91e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 16.62 on 18 degrees of freedom Multiple R-squared: 0.7252, Adjusted R-squared: 0.71 F-statistic: 47.51 on 1 and 18 DF, p-value: 1.91e-06 > for(x in seq(3,30,by=.25)){ + b0=coef(reg)[1] + b1=coef(reg)[2] + Yx=b0+b1*x + V=vcov(reg) + Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2 + IC=Yx+c(-1,+1)*1.96*sqrt(Vx) + points(x,Yx,pch=19,col="purple") + points(c(x,x),IC,pch=3,col="green")}
Notons qu’il est possible d’utiliser des fonctions de R pour faire des prévisions, avec des intervalles de confiance
> predict(reg, + newdata=data.frame(speed=c(15,25)),interval= "confidence") fit lwr upr 1 42.62976 34.75450 50.50502 2 84.87677 68.92746 100.82607 > predict(reg, + newdata=data.frame(speed=15),interval= "prediction") fit lwr upr 1 42.62976 6.836077 78.42344
Quand on a plus d’une variable explicative, c’est plus compliqué de “visualiser” la régression
> chicago=read.table("http://freakonometrics.free.fr/chicago.txt", + header=TRUE,sep=";") > Y=chicago$Fire > X1=chicago$X_1 > X2=chicago$X_2 > X3=chicago$X_3 > base=data.frame(Y,X1,X2,X3) > plot(X2,X3) > reg=lm(Y~X2+X3,data=base) > y=function(x2,x3) predict(reg,newdata=data.frame(X2=x2,X3=x3)) > VX2=seq(0,80) > VX3=seq(5,25) > VY=outer(VX2,VX3,y) > image(VX2,VX3,VY) > contour(VX2,VX3,VY,add=TRUE)
qui correspond à un plan de régression
> persp(VX2,VX3,VY,theta=30,ticktype=detailed)
On reviendra plus en détails sur ce point, mais il est possible de faire des régressions non linéaires assez facilement, à partir de ce modèle linéaire. On avait commencé par un modèle linéaire sur le logarithme de la distance
> plot(cars$speed,log(cars$dist)) > reg1=lm(log(dist)~speed,data=cars) > abline(reg1,col="red")
(on le verra, ce n’est pas fini, car on n’a pas ici de prévision sur la distance, juste sur son logarithme… mais promis, on en reparlera) ou sur la racine carrée
> plot(cars$speed,sqrt(cars$dist)) > reg1=lm(sqrt(dist)~speed,data=cars) > abline(reg1,col="red")
Au lieu de transformer la variable d’intérêt, on peut aussi transformer la variable explicative. On peut pendre des puissances, ou des fonctions simples, mais aussi mettre des ruptures. On avait commencé par une variable indicatrice,
> plot(cars$speed,cars$dist) > s=10 > abline(v=s,col="green") > regs=lm(dist~speed+I(speed>s),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + I(speed > s), data = cars) Residuals: Min 1Q Median 3Q Max -29.472 -9.559 -2.088 7.456 44.412 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.2964 6.7709 -2.555 0.0139 * speed 4.3140 0.5762 7.487 1.5e-09 *** I(speed > s)TRUE -7.5116 7.8511 -0.957 0.3436 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.39 on 47 degrees of freedom Multiple R-squared: 0.6577, Adjusted R-squared: 0.6432 F-statistic: 45.16 on 2 and 47 DF, p-value: 1.141e-11
Mais on peut aussi mettre des fonctions afin d’avoir un modèle linéaire par morceaux, tout en étant continu
> plot(cars) > s=15 > abline(v=s,col="green") > positive=function(x) ifelse(x>0,x,0) > regs=lm(dist~speed+positive(speed-s),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + positive(speed - s), data = cars) Residuals: Min 1Q Median 3Q Max -29.502 -9.513 -2.413 5.195 45.391 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.6519 10.6254 -0.720 0.47500 speed 3.0186 0.8627 3.499 0.00103 ** positive(speed - s) 1.7562 1.4551 1.207 0.23350 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.31 on 47 degrees of freedom Multiple R-squared: 0.6616, Adjusted R-squared: 0.6472 F-statistic: 45.94 on 2 and 47 DF, p-value: 8.761e-12
On a ici une rupture, mais on pourrait imaginer en avoir plusieurs
> nouvellebase=data.frame(speed=5:25) > y=predict(regs,newdata=nouvellebase) > lines(5:25,y,col="red") > > plot(cars$speed,cars$dist) > s1=10 > s2=20 > abline(v=c(s1,s2),col="green") > positive=function(x) ifelse(x>0,x,0) > regs=lm(dist~speed+positive(speed-s1)+positive(speed-s2),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + positive(speed - s1) + positive(speed - s2), data = cars) Residuals: Min 1Q Median 3Q Max -24.374 -9.475 -2.625 6.639 43.914 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.6305 16.2941 -0.468 0.6418 speed 3.0630 1.8238 1.679 0.0998 . positive(speed - s1) 0.2087 2.2453 0.093 0.9263 positive(speed - s2) 4.2812 2.2843 1.874 0.0673 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15 on 46 degrees of freedom Multiple R-squared: 0.6821, Adjusted R-squared: 0.6613 F-statistic: 32.89 on 3 and 46 DF, p-value: 1.643e-11
Comme vu en cours, le test de significativité des deux derniers coefficients ne veut pas dire que la pente est nulle, mais qu’elle est significativement différente de cette obtenue sur la zone de gauche (avant les deux seuils).
Modèles de Prévision, suite
La seconde partie de transparents pour le cours ACT6420 est en ligne. Il s’agit d’un document qui permettra à tout le monde garder des traces, mais je n’entends pas forcément les suivre… Par exemple, ce soir, je propose de jouer un peu avec les données d’incendies à Chicago
> chicago=read.table("http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
ou une base sur les états américains (cf un ancien billet)
> US=read.table("http://freakonometrics.free.fr/US.txt", + header=TRUE,sep=";")
Sinon, la semaine passée, on avait évoqué la matrice chapeau (dont on reparlera davantage lorsque l’on évoquera la phase de diagnostique). L’idée était de noter que
avec
On peut visualiser cette matrice avec le code suivant
> X=cbind(1,cars$speed)
> M=X %*% solve( t(X) %*% X) %*% t(X)
> Q=matrix(as.numeric(cut(M,c(-10,quantile(M,(1:5)/6),10),labels=1:6)),50,50)
> plot(1:50,1:50,col="white",axes=FALSE,xlab="",ylab="",ylim=c(50,0))
> library(RColorBrewer)
> couleur=rev(brewer.pal(6,"RdBu"))
> for(i in 1:50){for(j in 1:50){ points(i,j,col=couleur[Q[i,j]],pch=15,cex=1)}}
> axis(1)
> axis(2)
avec des cases rouges pour les valeurs importantes (et positive) et des case bleues pour des valeurs importantes (mais négatives)
Modèles de Prévision
Cet session d’été, je redonne le cours ACT6420, modèles de prévision. Le plan de cours est maintenant en ligne (comme souvent, les salles de cours ne sont pas encore attribuées définitivement mais elles seront disponibles sur le site de l’UQAM). Les transparents pour la première partie de cours sont en ligne (il s’agit de rappels de statistique, et d’une introduction aux modèles de régression).
Des exercices et des codes informatiques seront mis en ligne sur /courses/modeles-de-prevision.
Durant le cours, j’essayerais d’illustrer avec le maximum d’exemples sur des données réelles, à partir de codes R. Quelques séances de démonstrations se feront en salle informatique, afin de découvrir les rudiments de R. Pour aller plus loin, je recommande
- “R pour les débutants”d’Emmanuel Paradis (PDF)
- “Brise Glace-R” d’Andrew Robinson er Arnaud Schloesing (PDF)
- “Introduction à la programmation en R” de Vincent Goulet (PDF)
ainsi que quelques tutoriels en ligne, par exemple, dans un ancien billet, mais aussi sur https://youtube.com/.