{"word_count": 789, "character_count": 7297}
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.&#62; plot(cars)&#62; reg=lm(dist~speed,data=cars)&#62; 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(&#62;|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 ‘ ’ 1Residual standard error: 15.38 on 48 degrees of freedomMultiple R-squared:  0.6511,Adjusted R-squared:  0.6438 F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12Pour faire plusieurs prévisions, à la main, on peut utiliser le code suivant (la boucle permet de faire des prévisions pour plusieurs valeurs)&#62; 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&#62; I=sample(1:50,size=20)&#62; 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&#62; 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(&#62;|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 ‘ ’ 1Residual standard error: 16.62 on 18 degrees of freedomMultiple R-squared:  0.7252,Adjusted R-squared:   0.71 F-statistic: 47.51 on 1 and 18 DF,  p-value: 1.91e-06&#62; 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&#62; predict(reg,+ newdata=data.frame(speed=c(15,25)),interval= "confidence")       fit      lwr       upr1 42.62976 34.75450  50.505022 84.87677 68.92746 100.82607&#62; predict(reg,+ newdata=data.frame(speed=15),interval= "prediction")       fit      lwr      upr1 42.62976 6.836077 78.42344Quand on a plus d'une variable explicative, c'est plus compliqué de "visualiser" la régression&#62;  chicago=read.table("http://freakonometrics.free.fr/chicago.txt",+  header=TRUE,sep=";")&#62;  Y=chicago$Fire&#62;  X1=chicago$X_1&#62;  X2=chicago$X_2&#62;  X3=chicago$X_3&#62;  base=data.frame(Y,X1,X2,X3)&#62; plot(X2,X3)&#62; reg=lm(Y~X2+X3,data=base)&#62; y=function(x2,x3) predict(reg,newdata=data.frame(X2=x2,X3=x3))&#62; VX2=seq(0,80)&#62; VX3=seq(5,25)&#62; VY=outer(VX2,VX3,y)&#62; image(VX2,VX3,VY)&#62; contour(VX2,VX3,VY,add=TRUE)qui correspond à un plan de régression&#62; 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&#62; plot(cars$speed,log(cars$dist))&#62; reg1=lm(log(dist)~speed,data=cars)&#62; 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&#62; plot(cars$speed,sqrt(cars$dist))&#62; reg1=lm(sqrt(dist)~speed,data=cars)&#62; 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,&#62; plot(cars$speed,cars$dist)&#62; s=10&#62; abline(v=s,col="green")&#62; regs=lm(dist~speed+I(speed&#62;s),data=cars)&#62; summary(regs)Call:lm(formula = dist ~ speed + I(speed &#62; 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(&#62;|t|)    (Intercept)      -17.2964     6.7709  -2.555   0.0139 *  speed              4.3140     0.5762   7.487  1.5e-09 ***I(speed &#62; s)TRUE  -7.5116     7.8511  -0.957   0.3436    ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 15.39 on 47 degrees of freedomMultiple R-squared:  0.6577,Adjusted R-squared:  0.6432 F-statistic: 45.16 on 2 and 47 DF,  p-value: 1.141e-11Mais on peut aussi mettre des fonctions afin d'avoir un modèle linéaire par morceaux, tout en étant continu&#62; plot(cars)&#62; s=15&#62; abline(v=s,col="green")&#62; positive=function(x) ifelse(x&#62;0,x,0)&#62; regs=lm(dist~speed+positive(speed-s),data=cars)&#62; 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(&#62;|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 ‘ ’ 1Residual standard error: 15.31 on 47 degrees of freedomMultiple R-squared:  0.6616,Adjusted R-squared:  0.6472 F-statistic: 45.94 on 2 and 47 DF,  p-value: 8.761e-12On a ici une rupture, mais on pourrait imaginer en avoir plusieurs&#62; nouvellebase=data.frame(speed=5:25)&#62; y=predict(regs,newdata=nouvellebase)&#62; lines(5:25,y,col="red")&#62; &#62; plot(cars$speed,cars$dist)&#62; s1=10&#62; s2=20&#62; abline(v=c(s1,s2),col="green")&#62; positive=function(x) ifelse(x&#62;0,x,0)&#62; regs=lm(dist~speed+positive(speed-s1)+positive(speed-s2),data=cars)&#62; 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(&#62;|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 ‘ ’ 1Residual standard error: 15 on 46 degrees of freedomMultiple R-squared:  0.6821,Adjusted R-squared:  0.6613 F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11Comme 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).