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


Leave a Reply

Your email address will not be published. Required fields are marked *