Bachelor regression Linéaire
S. Jaubert
10 avril 2020
Rappel sur le modèle linéaire simple
la régression linéaire simple est une méthode statistique classique, qui est employée pour évaluer la significativité du lien linéaire entre deux variables numériques continues.
Exemple : Une machine produit automatiquement des pièces cylindriques. Réglée initialement pour un diamètre de 8 mm, elle se dérègle en cours d’utilisation. Afin de contrôler la fabrication et de procéder aux réglages éventuellement nécessaires, on mesure le diamètre de la dernière pièce dans chaque série de dix pièces produites. Les résultats sont les suivants :
``
X1<-seq(10,100,10)
Y<-c(8.03,8,8.01,8.01,8.02,8.03,8.03,8.04,8.05,8.06)
plot(x = X1,y = Y)
Le modèle linéaire est visiblement le plus approprié :
\(Y_{i}=aX1_{i}+b+\epsilon_{i}\) où chaque \(\epsilon_{i}\) est un terme résiduel aléatoire suivant la loi normale \(N(0,\sigma^{2})\) (on dit que les \(\epsilon_{i}\) sont des variables iid) a et b sont des paramètres inconnus.
L’idée est de minimiser (méthode des moindres carrés) la somme des \(\epsilon_{i}^{2}\) :
\[ \sum_{i=1}^{n} \epsilon_{i}^{2}=\sum_{i=1}^{n}(Y_{i}-(aX1_{i}+b))^{2}\] avec \(Y=(Y_{1},Y_{2},...,Y_{n})\) et \(X1=(X1_{1},X1_{2},...,X1_{n})\)
Où Y et X1 sont les vecteurs de nos données observées
Le plus simple est de construire un système matriciel :
(X<-matrix(cbind(rep(1,length(X1)),X1),ncol=2) )## [,1] [,2]
## [1,] 1 10
## [2,] 1 20
## [3,] 1 30
## [4,] 1 40
## [5,] 1 50
## [6,] 1 60
## [7,] 1 70
## [8,] 1 80
## [9,] 1 90
## [10,] 1 100(Y<-matrix(Y,ncol=1))## [,1]
## [1,] 8.03
## [2,] 8.00
## [3,] 8.01
## [4,] 8.01
## [5,] 8.02
## [6,] 8.03
## [7,] 8.03
## [8,] 8.04
## [9,] 8.05
## [10,] 8.06On peut alors écrire :
\(Y=X\beta+\epsilon\) où \(\beta\) est la matrice (2,1) : \(\begin{pmatrix}b\\a\end{pmatrix}\) On a :
\(\epsilon=Y-X\beta\) et on démontre que l’estimateur de \(\beta\) est :
\(\hat{\beta} = (X'X)^{-1}X'Y\) (\(X'\) représente la transposée de \(X\))
Effectuons les calculs avec R :
(beta<-solve(t(X)%*%X)%*%t(X)%*%Y)## [,1]
## [1,] 8.0000000000
## [2,] 0.0005090909Nous pouvons calculer à présent le vecteur \(\hat{Y}\) prédit par notre modèle :
Y_hat<-X%*%betapuis la somme des carrés des résidus :
(Scres<-sum((Y-Y_hat)^{2}))## [1] 0.001021818Nous verrons plus loin que l’estimateur sans biais de \(\sigma^{2}\) est donné par \(CMres=\frac{Scres}{n-p}\) (ici n=10 et p=2, n-p est le nombre de ddl associé à Scres)
n<-10 #nombre de lignes de la matrice X
p<-2 # " " colonnes " " " XSoit :
(CMres<-Scres/(n-p))## [1] 0.0001277273La somme des carrés expliqués par la régression :
(Screg<-sum((Y_hat-mean(Y))^{2}))## [1] 0.002138182CMreg<-Screg/(p-1)Retrouvons avec R les valeurs critiques permettant la prise de décision (Test de Fisher) :
(Fobs<-CMreg/CMres)## [1] 16.74021(Fth<-qf(p = 1-0.05,df1 = p-1,df2 = n-p))## [1] 5.317655pf(Fobs,df1 = p-1,df2 = n-p,lower.tail = F)## [1] 0.003478466\(F_{obs}\) est bien supérieur au \(F_{th}\), on rejette donc l’hypothèse nulle \(H_{0} : \beta_{1}=0\)
Pour rappel (voir cours sur l’ANOVA) :
\(\begin{matrix}\sum_{i}(Y_{i}-\bar{Y})^{2} &=&\sum_{i}(Y_{i}-\hat{Y_{i}})^{2} &+& \sum_{i} (\hat{Y_{i}}-\bar{Y})^{2} \\ SCT &=& SCres &+& SCreg \\ \end{matrix}\)
La capacité du modèle à expliquer les variations de Y est donné par le fameux \(R^{2}=\frac{Screg}{Screg+Scres}\) (c’est le pourcentage de Y qui est expliqué par la régression)
Soit ici :
(R2<-Screg/(Screg+Scres))## [1] 0.6766398Sous R il est très facile d’obtenir tous ces résultats très simplement avec l’instruction lm :
(model<-lm(Y~X1))##
## Call:
## lm(formula = Y ~ X1)
##
## Coefficients:
## (Intercept) X1
## 8.0000000 0.0005091Les coefficients de notre modèle sont donnés avec :
coefficients(model)## (Intercept) X1
## 8.0000000000 0.0005090909le vecteur \(\hat{Y}\) prédit par notre modèle sera donné par :
fitted.values(model)## 1 2 3 4 5 6 7 8
## 8.005091 8.010182 8.015273 8.020364 8.025455 8.030545 8.035636 8.040727
## 9 10
## 8.045818 8.050909et la SCres est obtenue très simplement par :
sum(residuals(model)^2)## [1] 0.001021818Que l’on retrouve aussi par une ANOVA classique :
anova(model)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 0.0021382 0.00213818 16.74 0.003478 **
## Residuals 8 0.0010218 0.00012773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Donnons un peu plus d’explications sur les résultats fournis par le summary(model) :
summary(model)##
## Call:
## lm(formula = Y ~ X1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.010364 -0.005591 -0.003000 0.003000 0.024909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.0000000 0.0077205 1036.203 < 2e-16 ***
## X1 0.0005091 0.0001244 4.091 0.00348 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0113 on 8 degrees of freedom
## Multiple R-squared: 0.6766, Adjusted R-squared: 0.6362
## F-statistic: 16.74 on 1 and 8 DF, p-value: 0.003478Ligne (Intercept) : nous avons le premier coefficient de notre modèle :8.00000 puis son écart-type 0.0077205 et enfin la t value pour le test de Student de nullité des coefficients, donnée par 8/0.0077205 = 1036.203. Ici la probabilité que le hasard puisse expliquer cette valeur est très faible 2e-16, on rejette donc l’hypothèse de nullité de ce coefficient
Ligne X1 : Nous avons le coefficient directeur de notre modèle : 0.0005091, son écart-type et la t-value, comme précédemment on rejettera donc l’hypothèse de nullité de ce coefficient
Les codes ‘***’ ‘**’ ‘*’ nous donne le risque \(\alpha\) de 1ère espèce de rejeter l’hypothèse de nullité des coefficients.
Residual standard error : est l’écart-type résiduel \(\sigma_{X}=\sqrt{CMres}\)
CMres<-Scres/(n-p)
sqrt(CMres)## [1] 0.01130165Multiple R-squared : est le fameux \(R^{2}=\frac{Screg}{Screg+Scres}\) que nous avions déjà calculé plus haut par :
(R2<-Screg/(Screg+Scres))## [1] 0.6766398F-statistic : Nous donne la F-value du test de Fisher de significativité de notre modèle. Que nous avions obtenu aussi par une ANOVA ou par :
(Fobs<-CMreg/CMres)## [1] 16.74021Traçons à présent la droite prédite par le modèle :
plot(x = X1,y = Y,type="p")
abline(model)
Cas du modèle multiple linéaire (on se restreint au cas où le nombre d’observations est supérieur au nombre de caractéristiques)
Un exemple : En fonction des paramètres de coupe : V (vitesse de coupe m/min), f (avance mm/rev), D (profondeur de coupe mm) nous obtenons une rugosité Ra (µm) et des vibrations radiales Ve (mm/s)
Pouvons nous modéliser par une régression linéaire la réponse Ra ou Ve en fonction de V, f et D ?
Commençons par entrer les données (téléchargeables ici : https://sjaubert.github.io/SPCR/rugo_vib.csv)
setwd("D:/OneDrive - CFAI Centre/Bachelor/Regression-Analysis-with-R-master") #à adapter
donnees<-read.csv2("rugo_vib.csv")
donnees## D V f Ra Ve
## 1 0.15 140 0.08 0.87 0.20
## 2 0.15 140 0.16 1.37 0.21
## 3 0.15 140 0.22 1.98 0.24
## 4 0.15 210 0.08 1.24 0.37
## 5 0.15 210 0.16 1.29 0.41
## 6 0.15 210 0.22 1.82 0.50
## 7 0.15 280 0.08 1.03 0.50
## 8 0.15 280 0.16 1.18 0.52
## 9 0.15 280 0.22 1.87 0.55
## 10 0.33 140 0.08 0.93 0.20
## 11 0.33 140 0.16 1.58 0.24
## 12 0.33 140 0.22 2.04 0.25
## 13 0.33 210 0.08 1.43 0.47
## 14 0.33 210 0.16 1.75 0.48
## 15 0.33 210 0.22 2.21 0.54
## 16 0.33 280 0.08 1.24 0.55
## 17 0.33 280 0.16 1.43 0.55
## 18 0.33 280 0.22 2.15 0.57
## 19 0.50 140 0.08 1.15 0.33
## 20 0.50 140 0.16 1.66 0.55
## 21 0.50 140 0.22 2.47 0.69
## 22 0.50 210 0.08 1.09 0.48
## 23 0.50 210 0.16 1.57 0.58
## 24 0.50 210 0.22 2.44 0.70
## 25 0.50 280 0.08 0.98 0.47
## 26 0.50 280 0.16 1.45 0.53
## 27 0.50 280 0.22 2.27 0.57Construisons les matrices de notre modèle linéaire
attach(donnees)
(X<-matrix(cbind(rep(1,27),D,V,f),ncol = 4))## [,1] [,2] [,3] [,4]
## [1,] 1 0.15 140 0.08
## [2,] 1 0.15 140 0.16
## [3,] 1 0.15 140 0.22
## [4,] 1 0.15 210 0.08
## [5,] 1 0.15 210 0.16
## [6,] 1 0.15 210 0.22
## [7,] 1 0.15 280 0.08
## [8,] 1 0.15 280 0.16
## [9,] 1 0.15 280 0.22
## [10,] 1 0.33 140 0.08
## [11,] 1 0.33 140 0.16
## [12,] 1 0.33 140 0.22
## [13,] 1 0.33 210 0.08
## [14,] 1 0.33 210 0.16
## [15,] 1 0.33 210 0.22
## [16,] 1 0.33 280 0.08
## [17,] 1 0.33 280 0.16
## [18,] 1 0.33 280 0.22
## [19,] 1 0.50 140 0.08
## [20,] 1 0.50 140 0.16
## [21,] 1 0.50 140 0.22
## [22,] 1 0.50 210 0.08
## [23,] 1 0.50 210 0.16
## [24,] 1 0.50 210 0.22
## [25,] 1 0.50 280 0.08
## [26,] 1 0.50 280 0.16
## [27,] 1 0.50 280 0.22(Y<-matrix(Ra,ncol=1))## [,1]
## [1,] 0.87
## [2,] 1.37
## [3,] 1.98
## [4,] 1.24
## [5,] 1.29
## [6,] 1.82
## [7,] 1.03
## [8,] 1.18
## [9,] 1.87
## [10,] 0.93
## [11,] 1.58
## [12,] 2.04
## [13,] 1.43
## [14,] 1.75
## [15,] 2.21
## [16,] 1.24
## [17,] 1.43
## [18,] 2.15
## [19,] 1.15
## [20,] 1.66
## [21,] 2.47
## [22,] 1.09
## [23,] 1.57
## [24,] 2.44
## [25,] 0.98
## [26,] 1.45
## [27,] 2.27Nous cherchons donc à estimer le vecteur \(a\) tel que : \[Y=aX+\epsilon\] (\(X\) matrice à n=27 lignes et p+1=4 colonnes) \[y_{i}=a_{0}+a_{1}x_{i,1}+a_{2}x_{i,2}+...+a_{p}x_{i,p}+\epsilon_{i}\]
les hypothèses sont :
Les X sont observés sans erreur
\(E(\epsilon)=0\) en moyenne le modèle est bien spécifié
\(\sigma_{\epsilon}^{2}=E(\epsilon^{2})-E(\epsilon)^{2}=E(\epsilon^{2})\) la variance de l’erreur est constante
\(\epsilon \hookrightarrow N(0,\sigma_{\epsilon})\)
On rappelle, comme pour le cas simple, que l’idée est de minimiser \(\left \| \epsilon \right \|^{2}=\epsilon'\epsilon=(Y-aX)'(Y-aX)\) (méthode des moindres carrés)
et on trouve de façon classique un estimateur de \(a\) par : \[\hat{a}=(X'X)^{-1}X'Y\] (\(X'X\) est définie positive donc inversible car n>p+1)
(a<-solve(t(X)%*%X)%*%t(X)%*%Y)## [,1]
## [1,] 0.2873671985
## [2,] 0.7766291863
## [3,] -0.0003571429
## [4,] 7.2237237237Biais de \(\hat{a}\)
\[\hat{a}=(X'X)^{-1}X'Y\] \[\hat{a}=(X'X)^{-1}X'(Xa+\epsilon)\] \[\hat{a}=a+(X'X)^{-1}X'\epsilon\]
Or comme X est non aléatoire
\[E(\hat{a})=a \] donc \(\hat{a}\) est sans biais
Matrice de variance covariance de \(\hat{a}\)
\(\hat{a}-a=(X'X)^{-1}X'\epsilon\) alors
\(\Omega_{\hat{a}}= E((\hat{a}-a)(\hat{a}-a)')=(X'X)^{-1}X'E(\epsilon\epsilon')X(X'X)^{-1}\)
et comme \(E(\epsilon_{i}\epsilon_{j})=0 \,\ i\neq j\)
On a : \(\Omega_{\hat{a}}=\sigma_{\epsilon}^{2}(X'X)^{-1}\)
Faisons les calculs :
\(Y-aX=\epsilon\)
Y_hat<-X%*%a
(SCres<-sum((Y-Y_hat)^2))## [1] 0.8950231\(\sigma_{\epsilon}^{2}=\frac{Scres}{(27-3)}=Cmres\)
(Cmres<-SCres/24)## [1] 0.03729263et on obtient la matrice de covariance :
(Omega<-round(Cmres*(solve(t(X)%*%X)),3))## [,1] [,2] [,3] [,4]
## [1,] 0.037 -0.022 0 -0.064
## [2,] -0.022 0.068 0 0.000
## [3,] 0.000 0.000 0 0.000
## [4,] -0.064 0.000 0 0.420Les écart-types des estimateurs de \(\hat{a}\) sont donnés par les racines carrées des éléments diagonaux :
(sqrt(diag(Omega)))## [1] 0.1923538 0.2607681 0.0000000 0.6480741Calculons pour terminer, le pourcentage de variation de la rugosité expliquée par D, V et f :
Pour cela nous avons besoin du \(SCreg=\sum_{i} (\hat{Y_{i}}-\bar{Y})^{2}\)
(SCreg<-sum((Y_hat-mean(Y))^2))## [1] 4.977606\(R^{2}=\frac{SCreg}{SCreg+SCres}\)
SCreg/(SCreg+SCres)## [1] 0.8475941Calcul direct avec R
modele<-lm(Ra~D+V+f)summary(modele)##
## Call:
## lm(formula = Ra ~ D + V + f)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28148 -0.14061 -0.05358 0.13689 0.38345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2873672 0.1968052 1.460 0.15777
## D 0.7766292 0.2656561 2.923 0.00764 **
## V -0.0003571 0.0006642 -0.538 0.59596
## f 7.2237237 0.6619828 10.912 1.44e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1973 on 23 degrees of freedom
## Multiple R-squared: 0.8476, Adjusted R-squared: 0.8277
## F-statistic: 42.64 on 3 and 23 DF, p-value: 1.473e-09pairs(donnees, lower.panel = NULL)
Terminons par une ANOVA
anova(modele)## Analysis of Variance Table
##
## Response: Ra
## Df Sum Sq Mean Sq F value Pr(>F)
## D 1 0.3326 0.3326 8.5465 0.007643 **
## V 1 0.0112 0.0112 0.2891 0.595963
## f 1 4.6338 4.6338 119.0772 1.436e-10 ***
## Residuals 23 0.8950 0.0389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nous pouvons remarquer que la variable V a très peu d’influence pour notre modèle (qu’en pense les spécialistes ?!)
Refaire la même étude pour cette fois rechercher un modèle de régression de Ve en fonction de V, f et D
TP
L’entreprise CITRON fabrique un matériau en matière plastique qui est utilisé dans la fabrication de jouets. Le département de contrôle qualité de l’entreprise a effectué une étude qui a pour but d’établir dans quelle mesure la résistance à la rupture (en kg/cm2) de cette matière plastique pouvait être affectée par l’épaisseur du matériau ainsi que la densité de ce matériau. Douze essais ont été effectués et les résultats sont présentés dans le tableau ci-dessous.
Y<-c(37.8,22.5,17.1,10.8,7.2,42.3,30.2,19.4,14.8,9.5,32.4,21.6)
X1<-c(4,4,3,2,1,6,4,4,1,1,3,4)
X2<-c(4.0,3.6,3.1,3.2,3.0,3.8,3.8,2.9,3.8,2.8,3.4,2.8)effectuer les analyses suivantes avec R et retrouver les résultats dans les tableaux.
Quel pourcentage de variation dans la résistance à la rupture est expliqué Pour chacune de ces régressions ?
Entre l’épaisseur du matériau et la densité qu’est-ce qui a le plus d’influence sur la résistance à la rupture ?
Corrigé (sans explication ! ce qui n’a aucun intérêt… A vous d’interpréter les résultats !!)
model<-lm(formula = Y~X1)
summary(model)##
## Call:
## lm(formula = Y ~ X1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.266 -4.887 -1.208 3.232 10.770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.523 4.383 0.804 0.440237
## X1 6.036 1.279 4.721 0.000816 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.633 on 10 degrees of freedom
## Multiple R-squared: 0.6903, Adjusted R-squared: 0.6593
## F-statistic: 22.29 on 1 and 10 DF, p-value: 0.0008155anova(model)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 980.63 980.63 22.285 0.0008155 ***
## Residuals 10 440.03 44.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model<-lm(formula = Y~X2)
summary(model)##
## Call:
## lm(formula = Y ~ X2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.1923 -5.1780 -0.2298 6.1123 12.3077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.373 20.489 -1.775 0.1062
## X2 17.464 6.069 2.878 0.0164 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.815 on 10 degrees of freedom
## Multiple R-squared: 0.453, Adjusted R-squared: 0.3983
## F-statistic: 8.282 on 1 and 10 DF, p-value: 0.01645anova(model)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X2 1 643.57 643.57 8.2816 0.01645 *
## Residuals 10 777.10 77.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model<-lm(formula = Y~X1+X2)
summary(model)##
## Call:
## lm(formula = Y ~ X1 + X2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.897 -2.135 -1.126 1.714 10.122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -30.081 11.455 -2.626 0.027542 *
## X1 4.905 1.014 4.838 0.000923 ***
## X2 11.072 3.621 3.058 0.013617 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.897 on 9 degrees of freedom
## Multiple R-squared: 0.8481, Adjusted R-squared: 0.8143
## F-statistic: 25.12 on 2 and 9 DF, p-value: 0.0002075anova(model)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 980.63 980.63 40.8959 0.000126 ***
## X2 1 224.22 224.22 9.3509 0.013617 *
## Residuals 9 215.81 23.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1