ozo = read.csv("~/data/ozone.txt", header = TRUE, sep =";" )
summary(ozo)
## O3 T12 T15 Ne12
## Min. : 41.8 Min. : 7.90 Min. :10.10 Min. :0.00
## 1st Qu.: 66.6 1st Qu.:17.40 1st Qu.:18.15 1st Qu.:3.00
## Median : 83.9 Median :19.75 Median :21.05 Median :6.00
## Mean : 86.3 Mean :20.32 Mean :21.39 Mean :5.02
## 3rd Qu.:102.2 3rd Qu.:24.32 3rd Qu.:25.70 3rd Qu.:7.00
## Max. :139.0 Max. :29.50 Max. :30.60 Max. :8.00
## N12 S12 E12 W12 VENT12
## Min. :0.00 Min. :0.00 Min. :0.00 Min. :0.00 E:16
## 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.00 1st Qu.:0.00 N: 9
## Median :0.00 Median :0.00 Median :0.00 Median :0.00 S: 7
## Mean :0.94 Mean :0.58 Mean :1.36 Mean :1.54 W:18
## 3rd Qu.:0.00 3rd Qu.:0.00 3rd Qu.:3.00 3rd Qu.:2.75
## Max. :7.00 Max. :6.00 Max. :8.00 Max. :8.00
## Vx maxO3v
## Min. :-27.0600 Min. : 38.00
## 1st Qu.:-10.8000 1st Qu.: 63.40
## Median : -3.2550 Median : 83.40
## Mean : -0.8312 Mean : 84.28
## 3rd Qu.: 9.3500 3rd Qu.:102.75
## Max. : 28.3600 Max. :142.80
ozo
## O3 T12 T15 Ne12 N12 S12 E12 W12 VENT12 Vx maxO3v
## 19960422 63.6 13.4 15.0 7 0 0 3 0 E 9.35 95.6
## 19960429 89.6 15.0 15.7 4 3 0 0 0 N 5.40 100.2
## 19960506 79.0 7.9 10.1 8 0 0 7 0 E 19.30 105.6
## 19960514 81.2 13.1 11.7 7 7 0 0 0 N 12.60 95.2
## 19960521 88.0 14.1 16.0 6 0 0 0 6 W -20.30 82.8
## 19960528 68.4 16.7 18.1 7 0 3 0 0 S -3.69 71.4
## 19960605 139.0 26.8 28.2 1 0 0 3 0 E 8.27 90.0
## 19960612 78.2 18.4 20.7 7 4 0 0 0 N 4.93 60.0
## 19960619 113.8 27.2 27.7 6 0 4 0 0 S -4.93 125.8
## 19960627 41.8 20.6 19.7 8 0 0 0 1 W -3.38 62.6
## 19960704 65.0 21.0 21.1 6 0 0 0 7 W -23.68 38.0
## 19960711 73.0 17.4 22.8 8 0 0 0 2 W -6.24 70.8
## 19960719 126.2 26.9 29.5 2 0 0 4 0 E 14.18 119.8
## 19960726 127.8 25.5 27.8 3 0 0 5 0 E 13.79 103.6
## 19960802 61.6 19.4 21.5 7 6 0 0 0 N -7.39 69.2
## 19960810 63.6 20.8 21.4 7 0 0 0 5 W -13.79 48.0
## 19960817 134.2 29.5 30.6 2 0 3 0 0 S 1.88 118.6
## 19960824 67.2 21.7 20.3 7 0 0 0 7 W -24.82 60.0
## 19960901 87.8 19.7 21.7 5 0 0 3 0 E 9.35 74.4
## 19960908 96.8 19.0 21.0 6 0 0 8 0 E 28.36 103.8
## 19960915 89.6 20.7 22.9 1 0 0 4 0 E 12.47 78.8
## 19960923 66.4 18.0 18.5 7 0 0 0 2 W -5.52 72.2
## 19960930 60.0 17.4 16.4 8 0 6 0 0 S -10.80 53.4
## 19970414 90.8 16.3 18.1 0 0 0 5 0 E 18.00 89.0
## 19970422 104.2 13.6 14.4 1 0 0 1 0 E 3.55 97.8
## 19970429 70.0 15.8 16.7 7 7 0 0 0 N -12.60 61.4
## 19970708 96.2 26.0 27.3 2 0 0 5 0 E 16.91 87.4
## 19970715 65.6 23.5 23.7 7 0 0 0 3 W -9.35 67.8
## 19970722 109.2 26.3 27.3 4 0 0 5 0 E 16.91 98.6
## 19970730 86.2 21.8 23.6 6 4 0 0 0 N 2.50 112.0
## 19970806 87.4 24.8 26.6 3 0 0 0 2 W -7.09 49.8
## 19970813 84.0 25.2 27.5 3 0 0 0 3 W -10.15 131.8
## 19970821 83.0 24.6 27.9 3 0 0 0 2 W -5.52 113.8
## 19970828 59.6 16.8 19.0 7 0 0 0 8 W -27.06 55.8
## 19970904 52.0 17.1 18.3 8 5 0 0 0 N -3.13 65.8
## 19970912 73.8 18.0 18.3 7 0 5 0 0 S -11.57 90.4
## 19970919 129.0 28.9 30.0 1 0 0 3 0 E 8.27 111.4
## 19970926 122.4 23.4 25.4 0 0 0 2 0 E 5.52 118.6
## 19980504 106.6 13.0 14.3 3 7 0 0 0 N 12.60 84.0
## 19980511 121.8 26.0 28.0 2 0 4 0 0 S 2.50 109.8
## 19980518 116.2 24.9 25.8 2 0 0 5 0 E 18.00 142.8
## 19980526 81.4 18.4 16.8 7 0 0 0 4 W -14.40 80.8
## 19980602 88.6 18.7 19.6 5 0 0 0 5 W -15.59 60.4
## 19980609 63.0 20.4 16.6 7 0 0 0 8 W -22.06 79.8
## 19980617 104.0 19.6 21.2 6 0 0 0 3 W -10.80 84.6
## 19980624 88.4 23.2 23.9 4 0 4 0 0 S -7.20 92.6
## 19980701 83.8 19.8 20.3 8 0 0 5 0 E 17.73 40.2
## 19980709 56.4 18.9 19.3 8 0 0 0 4 W -14.40 73.6
## 19980716 50.4 19.7 19.3 7 0 0 0 5 W -17.73 59.0
## 19980724 79.2 21.1 21.9 3 4 0 0 0 N 9.26 55.2
X = ozo$T12
X2 = ozo$T15
X3 = ozo$Vx
Y = ozo$O3
reg = lm(Y~ozo$T12)
reg2 = lm(Y~X2)
reg3 = lm(Y~X3)
plot(ozo$T12,Y, lwd = 3, ylab = 'O3', xlab = 'T12', main ="Pollution/ temperature 12")
abline(c(reg[1],reg[2]), lwd = 2, col = 'red')
legend("topleft", c("O2~T12","Reg. lineaire"), col= c("black", "red"), lty = c(3,1), lwd = 3)
plot(X2,Y, lwd = 3, ylab = 'O3', xlab = 'T15', main ="Pollution/ temperature 15")
abline(c(reg2[1],reg2[2]), lwd = 2, col = 'red')
legend("topleft", c("O2~T15","Reg. lineaire"), col= c("black", "red"), lty = c(3,1), lwd = 3)
plot(X3,Y, lwd = 3, ylab = 'O3', xlab = 'Vx')
abline(c(reg3[1],reg3[2]), lwd = 2, col = 'red', main ="Pollution/ Vx")
legend("topleft", c("O2~Vx","Reg. lineaire"), col= c("black", "red"), lty = c(3,1), lwd = 3)
Modelisation de la variable cible par une modelisation lineaire simple
Ainsi le modèle de régression linéaire simple semble pertinent pour modéliser la relation entre O3 et T12
\[ Y_i = a.1 + b*X_i + u \] où \(u\) est l’erreur tel que :
\[ \begin{aligned} &1-\mathbb E[U_i] =0 \\ & 2-Var(U_i) = \sigma^2 \\ & 3- \mathbb E[U_iU_j] = 0; \; \forall i \neq j \\ \end{aligned} \]
lm
lm
summary(reg)
##
## Call:
## lm(formula = Y ~ ozo$T12)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.256 -15.326 -3.461 17.634 40.072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.4150 13.0584 2.406 0.02 *
## ozo$T12 2.7010 0.6266 4.311 8.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.5 on 48 degrees of freedom
## Multiple R-squared: 0.2791, Adjusted R-squared: 0.2641
## F-statistic: 18.58 on 1 and 48 DF, p-value: 8.041e-05
shapiro.test(reg$residuals)
##
## Shapiro-Wilk normality test
##
## data: reg$residuals
## W = 0.97852, p-value = 0.4912
mean(reg$residuals)
## [1] 8.300131e-16
La moyenne est quasi nulle \(\mathbb E[U] = 8.300131e-16\). La p-value est supérieure à 5%, que nous pouvons poser comme notre seuil \(\alpha\) d’erreur. La statistique de test W étant très proche de 1, on en conclut que les résidus suivent une loi normale.
Comme \(\hat b\) est une combinaison linéaire de vecteurs gaussiens sous l’hypothèse que notre vecteur d’erreur \(U\) est gaussien. Comme la variance est inconnue, nous estimerons la variance e qui implique que nous changerons la loi de modélisation \(\hat {Var(\hat b)}\). Autrement dit :\(\hat b\) ~ \(\mathcal N(b,\sigma^2\frac{n}{\sum(X_i -\bar X)^2})\)
Par méconnaissance de la variance de \(\hat b\), nous devrons alors estimer la variance \(\sigma^2\) par \(\hat{V(\hat b)}\). Alors \(\hat b\)~\(t _{(n-2=48)}\).
Par l’estimation de la variance, l’intervalle de confiance (de \(95\)%) se construira à partir de quantile de la loi de Student. Donc les quantiles choisis sont \(q_{0.025}^t\), \(q_{0.975}^t\)
reg$coefficients[2] + qt(0.975,48) * 0.6
## ozo$T12
## 3.907416
reg$coefficients[2] + qt(0.025,48) * 0.6
## ozo$T12
## 1.494654
Ci-dessus les bornes de l’intervalle de confiance. Que l’on retrouve avec la commande confint
confint(reg, level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 5.159232 57.67071
## ozo$T12 1.441180 3.96089
Pb. de test : la variable est-elle significative ?
\(H_0=\) : n’est pas signficatif
\(H_1=\) : est signficatif
Ici la statistique de test est \(\mathcal T \not = 0\), tq :
\[ \mathcal T = \frac{\hat b}{\sqrt{\hat {V(\hat b)}}} \]
ainsi nous allons rejeter \(H_0\) car la p-value; donnée dans la régression linéaire lm
, est de \(9* 10^{-5}\), et la t-value est différente de 0.
plot(X,Y, main = "Pollution/Temperature", lwd =3)
abline(reg[1],reg[2], lwd = 2, col = "red")
legend("topleft",c("Pollution ~ Temperature12", "Droite ajustée"), lty = c(3,1), lwd = c(3,2), col = c("black","red"))
Rappel :
\(SCR + SCE = SCT\\\)
\(\mathcal R^2 = \frac{SCE}{SCT} \Leftrightarrow \mathcal R^2 = 1 - \frac{SCR}{SCT}\)
Yhat = reg$coefficients[1] + reg$coefficients[2]*X
var(Yhat)/var(Y) ## SCE/SCT
## [1] 0.2790813
1 - var(reg$residuals)/var(Y) ##1 - SCR/SCT
## [1] 0.2790813
On retrouve bien la valeur obtenue lors de lm
pour la variable statistique \(\mathcal R^2\) que l’on retrouve ici par \(\mathcal R^2 = \frac{SCE}{SCT}\). ; On retrouve bien \(\mathcal R^2\) par la relation \(\mathcal R^2 = 1 - \frac{SCR}{SCT}\). Malgré la modélisation qui semblait satisfaisante, la part de variance expliquée de la pollution par la température est de 27% seulement. Nous allons donc vérifier que notre régression est bien utile.
Pb de test : la regression est-elle utile ?
\(H_0\) : la régression est inutile.
\(H_1\) : la régression utile.
ou en termes mathématiques : si \(H_0\) est vraie alors \(\mathcal F\) est toute petite et si \(H_1\) est vraie alors \(\mathcal F\) est grande. Ici la statistique de test est \(\mathcal F\), tq :
\[ \mathcal F = \frac{R^2/1}{(1-R^2)/48} \]
Il est donné auparavant par la commande lm
et vaut environ \(18,58\). Ainsi nous rejetons l’hypothèse \(H_0\) d’inutilité de la régression et acceptons \(H_1\) qui montre que la régression est utile. La faible explication du regresseur sur notre variable d’intérêt n’est pas tellement étonnante. Effectivement, on peut se douter que la pollution ne s’explique pas que par une variable et qu’il en existe une certaine quantité expliquant la variation de la pollution plus ou moins signifivative.