TP3 Régression linéaire simple : intervalles de confiance, tests et diagnostics

Objectifs du TP

• Savoir déterminer un intervalle de confiance pour \(\beta_{1}\), \(\beta_{2}\), \(\sigma^2\).
• Savoir déterminer un intervalle de confiance pour une valeur moyenne, un intervalle de prévision, et comprendre la différence entre les deux.
• Savoir effectuer des tests sur \(\beta_{1}\), \(\beta_{2}\).
• Interpréter les sorties de la fonction lm, en particulier les tests et le coefficient de détermination.
• Faire un diagnostic sur les résidus, et utiliser le résultat pour construire un meilleur modèle.

Intervalles de confiance, intervalle de prévision

library(tidyverse)
data <- readr::read_delim(
  here::here("data", "ozone_simple.txt"), delim = ";", 
  col_types = c(col_double(), col_double())
)

1. Charger le jeu de données ozone_simple.txt et ajuster un modèle linéaire expliquant \(O3\) à l’aide de \(T12\).

Avant tout, on trace les données pour vérifier qu’il n’y ait pas des valeurs aberrantes.

df <- data.frame(data)
n <- nrow(df)

tibble(df) %>%
  ggplot(mapping = aes(x = O3)) +
    geom_histogram(
      aes(y = stat(density)),
      bins = 10,
      colour = "black", fill = "white"
    ) +
      geom_density(alpha = .2, fill = "blue")

tibble(df) %>%
  ggplot(mapping = aes(x = T12)) +
    geom_histogram(
      aes(y = stat(density)),
      bins = 10,
      colour = "black", fill = "white"
    ) +
      geom_density(alpha = .2, fill = "blue")

tibble(df) %>% 
  ggplot(mapping = aes(x = T12, y = O3)) +
    geom_point()

summary(df)

##        O3             T12       
##  Min.   : 41.8   Min.   : 7.90  
##  1st Qu.: 66.6   1st Qu.:17.40  
##  Median : 83.9   Median :19.75  
##  Mean   : 86.3   Mean   :20.32  
##  3rd Qu.:102.2   3rd Qu.:24.32  
##  Max.   :139.0   Max.   :29.50

var(df$O3) 

## [1] 571.2494

var(df$T12) 

## [1] 21.85224

Nous modélisons \(O3_i = \beta_1 + \beta_2T12_i + \varepsilon_i\)

reg <- lm(O3 ~ T12, data = df)
summary(reg)

## 
## Call:
## lm(formula = O3 ~ T12, data = df)
## 
## 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 *  
## 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

Les coefficients \(\beta_1\) et \(\beta_2\) sont significatifs au seuil de significativité \(\alpha = 5\%\), avec une \(p\)-value respective de \(0.02\) et \(8.04\times10^{-5}\). Le coefficient de détermination \(\mathfrak{R}^2 = 0.2791\). Cela veut dire que le modèle estimé n’est peut-être pas complet.

2. En utilisant les formules du cours, déterminer des intervalles de confiance de probabilité de couverture \(0.95\) pour \(\beta_1\), puis pour \(\beta_2\), et enfin pour \(\sigma^2\). Faire de même pour les coefficients \(\beta_1\) et \(\beta_2\) avec la fonction confint.

hat_beta1 <- reg$coefficients[1]
hat_beta2 <- reg$coefficients[2]

hat_sigma2 <- summary(reg)[["sigma"]]^2

var_beta1 <- (hat_sigma2/n) * ((sum(df$T12^2) / n) / (sum((df$T12 - mean(df$T12))^2) /n ))
var_beta2 <- (hat_sigma2/n) * (1 / (sum((df$T12 - mean(df$T12))^2) / n))

alpha <- 0.05
t.val <- qt(1 - alpha/2, n - 2)

IC1 <- c(hat_beta1 - t.val * sqrt(var_beta1), reg$coefficients[1] + t.val * sqrt(var_beta1))
IC2 <- c(hat_beta2 - t.val * sqrt(var_beta2), reg$coefficients[2] + t.val * sqrt(var_beta2))

n.ic.s2 <- (n - 2) * hat_sigma2
d.ic.low <-  qchisq(1 - alpha/2, df = n-2)
d.ic.up <-  qchisq(alpha/2, df = n-2)
IC3 <- n.ic.s2 * c(1/d.ic.low, 1/d.ic.up)

IC.confint <- confint(reg, level = 0.95)

Paramètre	Borne Inférieure	Borne Supérieure
\(\beta_{1}\)	5.1592319	57.6707146
\(\beta_{2}\)	1.4411798	3.9608898
\(\sigma^2\)	292.3593056	656.1443531

Les intervalles de confiance calculés à partir des formules discutées en cours sont égaux à ceux obtenus à l’aide de la fonction confint. À partir des intervalles de confiance à \(95\%\), nous remarquons que ceux associés aux coefficients \(\beta\) ne comprennent pas zéro. Cela n’est pas étonnant puisque nous avons vu précédemment que ces coefficients étaient significatifs. Cependant, \(\beta_1\) montre une grande variabilité ce qui est expliqué par l’ordre de magnitude des résidus.

3. Sur un même graphe, représenter les données, tracer la droite de régression, les intervalles de confiance à \(95\%\) pour \(E[y]\) et les intervalles de prévision à \(95\%\). On pourra utiliser la fonction predict. Commenter.

plot(O3 ~ T12, data = df, ylim = c(min(df$O3), max(df$O3)))

seq.x <- seq(min(df$T12), max(df$T12), length=3*n)
grid.x <- data.frame(seq.x)
dimnames(grid.x)[[2]] <- "T12"

ICconf <- predict(reg, new=grid.x, interval="confidence", level=0.95)
ICprev <- predict(reg, new=grid.x, interval="prediction", level=0.95)

matlines(
  grid.x, cbind(ICconf,ICprev[,-1]),
  lty=c(1,2,2,3,3), col= c(1, 2, 2, 3, 3)
)

legend(
  "topleft", lty = 1:3,
  c("Régression","Int. confiance", "Int. prevision")
)

L’écart des points par rapport à la ligne de régression est en accord avec la discussion précédente sur les résidus et la variance. Les intervalles des prévision sont toujours plus amples que ceux de confiance. Ceci est interprétable comme suit : Les intervalles de confiance concernent l’espérance de la concentration \(O3\), pour une température \(T12\) donnée, alors que les intervalles de prévision essaient de prendre en compte toute concentration possible dans le seuil fixé.

Tests, analyse de la variance

4. Tester les paramètres du modèle (en précisant bien \(H_{0}\) et \(H_{1}\) pour chaque test). Commenter.

Paramètre	\(H_{0}\)	\(H_{1}\)
\(\beta_{1}\)	\(\beta_{1} = 0\)	\(\beta_{1} \ne 0\)
\(\beta_{2}\)	\(\beta_{2} = 0\)	\(\beta_{2} \ne 0\)

t.val <- qt(0.975, n - 2)
Test1 <- hat_beta1[[1]] / sqrt(var_beta1)
Test1 > t.val

## [1] TRUE

Le coefficient \(\beta_{2}\) est significatif au seuil de \(\alpha = 0.05\).

Test2 <- hat_beta2[[1]] / sqrt(var_beta2)
Test2 > t.val

## [1] TRUE

Le coefficient \(\beta_{1}\) est significatif au seuil de \(\alpha = 0.05\).

Pour les deux paramètres, on rejette l’hypothèse nulle.

Les valeurs sont \(T1 =\frac{\hat{\beta}_1}{\hat{\sigma}_{\beta_1}} =\) 2.4057229 et \(T2 =\frac{\hat{\beta}_2}{\hat{\sigma}_{\beta_2}} =\) 4.3106503.

5. Calculer à la main la valeur du coefficient de détermination et la valeur de la statistique de Fisher. Commenter.

\(H_{0}\)	\(Y_{i} = \beta_1 + \varepsilon_i\) est vrai
\(H_{1}\)	\(Y_{i} = \beta_1 + \beta_2 x_{i} + \varepsilon_i\) est vrai

SCR <- sum(residuals(reg)^2)
SCM <- sum((fitted(reg) - mean(df$O3))^2)
SCT <- sum((df$O3 - mean(df$O3))^2)

r2.hat <- SCM /SCT

F.hat <- (SCM / 1) / (SCR / (n - 2))
f.test <- qf(0.95, 1, n - 2)
F.hat > f.test

## [1] TRUE

La valeur estimé du test Fisher est \(\mathcal{F}_{0.95}(1, 48) =\) 18.5817059. On rejette \(H_{0}\) : Le modèle \(Y_{i} = \beta_1 + \beta_2 x_{i} + \varepsilon_i\) explique mieux la relation entre la concentration d’ozone et la température que \(Y_{i} = \beta_1 + \varepsilon_i\).

Diagnostics

6. Diagnostiquer les résidus. On pourra notamment utiliser la commande plot appliquée à la sortie de la fonction lm.

Les résidus ne semblent pas suivre une loi normale d’après le qqplot puisque les queues de la densité des résidus standardardisés sont plus faibles que celles de la loi normale \(\mathcal{N}(0, 1)\).

L’hypothèse d’homoscedasticité ne semble pas être vérifiée car les \(\hat{\epsilon}_{i}\) n’ont pas tous la même variance. En effet, \(\epsilon_{1}, \cdots, \epsilon_{n}\) dépend de \(x_{1}, \cdots, x_{n}\) \(\rightarrow\) hétéroscedasticité.

hat_sigma <- sqrt(hat_sigma2)
residus <- tibble(
  n = 1:50,
  residus = reg$residuals,
  T12 = df$T12,
  O3 = df$O3,
  e.standardized = reg$residuals / hat_sigma,
  adjusted = predict(reg)
)

residus %>%
  ggplot(mapping = aes(sample = e.standardized)) +
    stat_qq() +
      stat_qq_line()

residus %>%
  ggplot(mapping = aes(x = T12, y = e.standardized)) +
    geom_point() + 
      geom_smooth() +
        ggtitle("epsilon ~ temperature")

residus %>%
  ggplot(mapping = aes(x = O3, y = e.standardized)) +
    geom_point() + 
      geom_smooth() +
        ggtitle("epsilon ~ O3")

residus %>%
  ggplot(mapping = aes(x = e.standardized)) +
    geom_histogram(
      aes(y = stat(density)),
      bins = 10,
      colour = "black", fill = "white"
    ) +
      geom_density(alpha = .2, fill = "blue") +
      geom_vline(xintercept = 2, colour="red") +
      geom_vline(xintercept = -2, colour="red")

residus %>% filter( abs(e.standardized) > 2)

On ne peut pas conclure qu’il y a des valeurs aberrantes à partir de ce graphique. En effet, si l’on filtre les erreurs standardisées dont la valeur absolue est supérieure à 2, on constate qu’il s’agit d’une seule observation \(\hat{\epsilon}/\sigma = -2.20\) qui est plausible.

7. Proposer une transformation des données susceptible de mieux expliquer la relation entre concentration maximale \(O3\) et température à midi \(T12\). Interpréter le modèle final (i.e. refaire les questions précédentes avec ce nouveau modèle).

D’après le graphique \(\left(x_i,\hat{\varepsilon}_i\right)\), nous remarquons une rupture entre les températures inférieures et supérieures à 20-22 dégrés. Nous pouvons donc proposer deux modèles univariés: modèle (1) \(O3_{i} = \beta_1^\star + \beta_2^\star \times T12_i + \varepsilon_i \quad \forall\; T_i < 21\) et modèle (2) \(O3_{i} = \beta_1^\ast + \beta_2^\ast \times T12_i + \varepsilon_i \quad \forall\; T_i > 21\).

reg1 <- lm(O3 ~ T12, data = df[df$T12<21,])
summary(reg1)

## 
## Call:
## lm(formula = O3 ~ T12, data = df[df$T12 < 21, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.661 -11.422  -1.862  13.147  31.950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  103.192     18.112   5.698 4.14e-06 ***
## T12           -1.589      1.035  -1.535    0.136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.54 on 28 degrees of freedom
## Multiple R-squared:  0.07765,    Adjusted R-squared:  0.04471 
## F-statistic: 2.357 on 1 and 28 DF,  p-value: 0.1359

reg2 <- lm(O3 ~ T12, data = df[df$T12>21,])
summary(reg2)

## 
## Call:
## lm(formula = O3 ~ T12, data = df[df$T12 > 21, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.103 -12.500  -3.174   9.842  31.458 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -87.176     41.529  -2.099 0.051042 .  
## T12            7.612      1.647   4.623 0.000243 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.09 on 17 degrees of freedom
## Multiple R-squared:  0.5569, Adjusted R-squared:  0.5309 
## F-statistic: 21.37 on 1 and 17 DF,  p-value: 0.000243

Dans le modèle (1), nous remarquons que le coefficient \(\beta_2^\star\) n’est pas significatif. En effet, la \(p-\text{value}\) est de \(0.136\) donc il n’y a pas de significativité même au seuil de \(\alpha = 10\%\). Concernant le modèle (2), le coefficient \(\beta_2^\ast\) est très significatif avec une \(p-\text{value} \approx 0\). Nous remarquons également un coefficient de détermination \(\mathfrak{R}^2 = 0.56\). Ce coefficient est très supérieur au coefficient \(\mathfrak{R}^2 = 0.28\) obtenu à la question (1). En revanche, le coefficient de détermination est très faible dans le modèle (1) car la concentration en \(O3\) n’est pas significativement expliquée par la température en deçà de 21 degrés celsius.

reg3 <- lm(O3 ~ 1, data = df[df$T12<21,])
summary(reg3)

## 
## Call:
## lm(formula = O3 ~ 1, data = df[df$T12 < 21, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.973 -12.623   0.227  12.677  30.827 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   75.773      3.089   24.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.92 on 29 degrees of freedom

df.mod1 <-  df[df$T12 < 21, ]
df.mod2 <-  df[df$T12 > 21, ]

tibble(df.mod1) %>%
  ggplot(mapping = aes(x = T12, y = O3)) +
    geom_point() +
      geom_abline(
        slope = reg1$coefficients[["T12"]],
        intercept = reg1$coefficients[["(Intercept)"]],
        colour = "blue", 
      ) + 
      geom_abline(
        slope = 0,
        intercept = reg3$coefficients[["(Intercept)"]],
        colour = "red"
      ) +
        ggtitle("Modele T12 < 21, bleu : avec Beta2, rouge sans Beta2")

tibble(df.mod2) %>%
  ggplot(mapping = aes(x = T12, y = O3)) +
    geom_point() +
      geom_abline(
        slope = reg2$coefficients[["T12"]],
        intercept = reg2$coefficients[["(Intercept)"]],
        colour = "green"
      ) +
        ggtitle("Modele T12 > 21")