Florent univarie

1 R?gression lin?aire simple

On cherche une relation entre la variable observ?e \(Y\) et la variable explicative \(x\).

\[ Y \sim x \]

Les objectifs peuvent ?tre multiples :

• mod?le explicatif

• recherche d'une corr?lation significative

• mod?les de pr?diction

Domaines d'application : chimie, astronomie, sismologie, ?conomie, d?mographie, ?pid?miologie, g?ographie, ?cologie, g?ologie, physique, m?decine...

Historiquement, le mod?le lin?aire a ?t? d?velopp? par Fisher, avec applications en g?n?tique et en agronomie.

---

1.1 Pr?sentation du mod?le

Consid?rons \(x\) et \(Y\) deux variables. Y est expliqu?e (mod?lis?e) par la variable explicative \(x\).

Important, dans le cadre des mod?les lin?aire Y est quantitatives (Ex : taille, poids, age...).

Dans le cas de la regression lin?aire simple, \(x\) est ?galement quantitatives et la relation entre \(Y\) et \(x\) est une fonction affine (ou lin?aire) de \(x\) (Cf. programme du brevet des coll?ges).

\[\begin{eqnarray} \text{f : } \mathbb{R} & \rightarrow & \mathbb{R} \\ x & \rightarrow & f(x) = ax + b \\ \end{eqnarray}\]

Attention, la fonction abline de R utilise la formule \(bx+a\) ! Dans ce cours nous utliserons la notation \(f(x) = \beta_0 + \beta_1 x\)

## Warning: Missing column names filled in: 'X1' [1]
## 
## -- Column specification --------------------------------------------------------
## cols(
##   X1 = col_double(),
##   sexe = col_character(),
##   situation = col_character(),
##   the = col_double(),
##   cafe = col_double(),
##   taille = col_double(),
##   poids = col_double(),
##   age = col_double(),
##   viande = col_character(),
##   poisson = col_character(),
##   fruit_crus = col_character(),
##   fruit_legume_cuits = col_character(),
##   chocol = col_character(),
##   matgras = col_character()
## )

Les \(n\) observations vont permettre :

• de construire le mod?le (d?finir \(\beta_0\) et \(\beta_1\)),
• de quantifier l'ad?quation du mod?le aux donn?es.

Remarques :

• On utilise un mod?le lin?aire Gaussien pour des observations pouvant ?tre mod?lis?es par une loi normales. Pour d'autres distributions (Poisson, Bernoulli...), on utilisera un mod?le lin?aire g?n?ralis? (GLM).
• La r?gression lin?aire se caract?rise par des variables explicatives quantitatives (e.g., \(poids~taille\)). L'ANOVA se caract?rise par des variables explicatives qualitatives (e.g., \(poids~sexe\)).

---

Exemple d'une relation d?terministe

La temp?rature \(x\) en degr?s Celsius explique enti?rement la temp?rature \(y\) en degr?s Farenheit par la relation : \(y=32 + 9/5 x\).

\[\begin{eqnarray} \text{f : } \mathbb{R} & \rightarrow & \mathbb{R} \\ x & \rightarrow & f(x) = \beta_0 + \beta_1 x = 32 + 9/5 x \\ 0 & \rightarrow & f(0) = 32 \\ 10 & \rightarrow & f(10) = 50 \\ \end{eqnarray}\]

---

Exemple d'une relation stochastique

Dans le cas d'une variable al?atoire \(Y\), il faut chercher la droite qui ajuste le mieux les observations et on ajoute un terme qui portera la partira non d?terministe des observations.

\[Y = \beta_0 + \beta_1x + \epsilon \]

\[Y_i = \beta_0 + \beta_1x_i + \epsilon_i \ \ \ \ \ \forall i \in {1, ..., n}\]

o? \(\beta_0\) et \(\beta_1\) sont des r?els fixes, mais inconnus, et \(\epsilon_i\) une valeur caract?risant le comportement individuel de l'observation \(i\) (r?sidus).

Mod?le : \[\mathbb{E} (Y|x) = \widehat Y = f(x) = \beta_0 + \beta_1 x\]

Par exemple \(x\) est la taille, \(Y\) le poids, ? une m?me taille, plusieurs poids peuvent correspondrent. Les donn?es ne sont plus align?es, mais s'ajustent autour de la droite de r?gression laissant apparaitre les r?sidus.

---

1.2 Estimation des param?tres

Question : comment estimer les param?tres de la droite de regression ?

R?ponse : Par la m?thode des moindres carr?s.

---

Droite de regression, estimation et r?sidus

En analyse de r?gression lin?aire, pour un individu \(i\) :

• \(x_i\) est la valeur de la varible explicative
• \(Y_i\) est la valeur observ?e de la variable al?atoire
• la composante al?atoire de \(Y_i\) est le \(\epsilon_i\) correspondant
• la valeur de \(Y_i\) estim?e par le mod?le est \(\widehat Y_i\).

Ainsi :

\[Y_i = \widehat Y_i + \epsilon_i\]

Lorsque \(x = x_i\), alors \(\widehat Y(x_i) = \widehat Y_i\) , c'est-?-dire :

\[ \widehat Y_i=\widehat{\beta}_0+ \widehat{\beta}_1x_i \] \(\widehat Y_i\) est appel?e la valeur estim?e par le mod?le.

Les quantit?s inobservables : \[\epsilon_i = Y_i - \beta_0 - \beta_1x_i\] sont estim?es par les quantit?s observables : \[ e_i= Y_i-\widehat Y_i\]

• Les quantit?s \(e_i\) sont appel?es les r?sidus du mod?le.
• La plupart des m?thodes d'estimation cherchent ? d?finir une droite qui minimise une fonction de r?sidus.
• La plus connue : la m?thode des moindres carr?s ordinaires (MCO).

---

Intutition des moindres carr?s

Pour l'instant, la droite de r?gression est inconnue. Tout le probl?me est d'estimer \(\beta_0\) et \(\beta_1\) ? partir d'un ?chantillon de donn?es.

On essaie de d?terminer la droite de regression qui approche le mieux les donn?es :

\[\mathbb{E} (Y|x)) = \widehat Y(x) = \widehat{\beta}_0 + \widehat{\beta}_1x\]

Avec \(\widehat Y(x)\) un estimateur de la moyenne th?orique \(\mu_Y(x)\) (\(\widehat Y(x)\) correspond ? la moyenne de \(Y\) m?sur?e sur tous les individus pour lesquels \(x\) prend une valeur donn?e).

\(\mu_Y(x)\) n'est ni observable, ni calculable (il faudrait alors recenser tous les individus de la population).

Remarque : la moyenne empirique de \(Y\), d?finie par \(\overline{Y}_n(x) = \frac{1}{n}\sum^n_{i=1}Y_i(x)\) est un autre estimateur de \(\mu_Y(x)\) . Si le mod?le est bon, \(\widehat Y(x)\) est plus pr?cis que \(\overline{Y}_n(x)\) qui lui correspondrait au mod?le nul, i.e, \(Y\) ne depend pas de \(x\).

layout(1, respect=TRUE)
beta_0 = mean(d$poids)
beta_1 = 0
res =  d$poids - (beta_0 + beta_1*d$taille)
eqm = sum(res^2)
plot(d$taille, d$poids, main=paste0("poids~taille, EQM=", round(eqm)), ylab="poids", xlab="taille")
abline(a=beta_0, b=beta_1, col=2) # /!\ y = b.x + a
arrows(d$taille, d$poids, d$taille, d$poids-res, length=0.1, col=adjustcolor(4, alpha.f=0.5))

layout(matrix(1:2, 1), respect=TRUE)
beta_0 = -100
beta_1 = 1
res =  d$poids - (beta_0 + beta_1*d$taille)
eqm = sum(res^2)
plot(d$taille, d$poids, main=paste0("poids=", signif(beta_0, 3), "+", signif(beta_1, 3), " taille, EQM=", round(eqm)), ylab="poids", xlab="taille")
abline(a=beta_0, b=beta_1, col=2) # /!\ y = b.x + a
arrows(d$taille, d$poids, d$taille, d$poids-res, length=0.1, col=adjustcolor(4, alpha.f=0.5))
beta_0 = -16
beta_1 = 0.5
res =  d$poids - (beta_0 + beta_1*d$taille)
eqm = sum(res^2)
plot(d$taille, d$poids, main=paste0("poids=", signif(beta_0, 3), "+", signif(beta_1, 3), " taille, EQM=", round(eqm)), ylab="poids", xlab="taille")
abline(a=beta_0, b=beta_1, col=2) # /!\ y = b.x + a
arrows(d$taille, d$poids, d$taille, d$poids-res, length=0.1, col=adjustcolor(4, alpha.f=0.5))

layout(1, respect=TRUE)
results = list()
for(beta_0 in 0:-260) {
  for(beta_1 in seq(0,2, length.out=21)) {
    res =  d$poids - (beta_0 + beta_1*d$taille)
    eqm = sum(res^2)
    results[[length(results)+1]] = c(beta_0=beta_0, beta_1=beta_1, eqm=eqm)
  }
}
results = as.data.frame(do.call(rbind, results))
# plot(results$beta_0, results$beta_1, pch=16, col=colorRampPalette(c("cyan", "black", "red"))(100)[as.numeric(cut(results$eqm, breaks=quantile(results$eqm, probs=(0:100)/100), include.lowest=TRUE))])
r = results[results$eqm==min(results$eqm),]
beta_0 = r[["beta_0"]]
beta_1 = r[["beta_1"]]
res =  d$poids - (beta_0 + beta_1*d$taille)
eqm = sum(res^2)
plot(d$taille, d$poids, main=paste0("poids=", signif(beta_0, 3), "+", signif(beta_1, 3), " taille, EQM=", round(eqm)), ylab="poids", xlab="taille")
abline(a=beta_0, b=beta_1, col=2) # /!\ y = b.x + a
suppressWarnings(arrows(d$taille, d$poids, d$taille, d$poids-res, length=0.1, col=adjustcolor(4, alpha.f=0.5)))

---

M?thode et d?monstration

L'objectifest de d?finir des estimateurs qui minimisent la somme des carr?s des r?sidus (ou erreur quadratique moyenne EQM).

Les estimateurs sont donc les coordonn?es du minimum de la fonction ? deux variables :

\[EQM(\beta_0, \beta_1) = \sum^{n}_{i=1} e^2_i = \sum^{n}_{i=1} (Y_i - \widehat Y_i)^2 = \sum^{n}_{i=1} (Y_i - \widehat\beta_0 - \widehat\beta_1x_i)^2\]

Cette fonction est appel?e la fonction objectif.

Les estimateurs correspondent aux valeurs annulant les d?riv?es partielles de cette fonction :

\[\frac{\delta EQM(\beta_0, \beta_1)}{\delta \beta_0} = -2 \sum (Y_i - \beta_0 - \beta_1x_i)\] \[\frac{\delta EQM(\beta_0, \beta_1)}{\delta \beta_1} = -2 \sum x_i (Y_i - \beta_0 - \beta_1x_i)\]

Les estimateurs \((\widehat\beta_0^{MCO},\widehat\beta_1^{MCO})\) sont les solutions du syst?me d'?quation : \[ \Bigg\{ \begin{align} -2 \sum (Y_i - \beta_0 - \beta_1x_i) = 0 \\ -2 \sum x_i (Y_i - \beta_0 - \beta_1x_i) = 0 \end{align} \]

Soient les equations suivantes :

\[(1) \sum{Y}_i = n\widehat\beta_0 + \widehat\beta_1 \sum{x_i} \]

\[(2) \sum{x_iY_i} = \widehat\beta_0 \sum{x_i} + \widehat\beta_1 \sum{x^2_i}\] On utilise les notations usuelles des moyennes empiriques \[\overline{x}_n = \frac{1}{n}\sum{x_i} ,\ \ \ \ \ \ \overline{Y}_n = \frac{1}{n}\sum{Y_i}\]

D'apr?s (1): \[\widehat\beta_0 = \overline{Y}_n - \widehat\beta_1 \overline{x}_n\]

D'apr?s (2) :

\[\begin{eqnarray} \widehat\beta_1 \sum{x^2_i} & = & \sum{x_iY_i} - \widehat\beta_0 \overline{x}_n \\ & = & \sum{x_iY_i} - n\overline{x}_n\overline{Y}_n + \widehat\beta_1 (\overline{x}_n)^2 \end{eqnarray}\]

Ainsi \[\widehat\beta_1 = \frac{\sum{x_iY_i} - n\overline{x}_n\overline{Y}_n}{\sum{x_i^2} - n(\overline{x}_n)^2}\]

Et comme nous avons :

\[\begin{eqnarray} \sum{(x_i-\overline{x}_n)(Y_i-\overline{Y}_n)} & = & \sum{x_iY_i} - n\overline{x}_n\overline{Y}_n \\ \sum{(x_i-\overline{x}_n)} & = & \sum{x_i^2} - n(\overline{x}_n)^2 \end{eqnarray}\]

Nous obtenons : \[\widehat\beta_1 = \frac{\sum{(x_i-\overline{x}_n)(Y_i-\overline{Y}_n)}}{\sum{(x_i-\overline{x}_n)}}\]

Dans la pratique, nous calculons \(\widehat\beta_1\) puis \(\widehat\beta_0\).

Nous obtenons ainsi une estimation de la droite de r?gression, appel?e la droite des moindres carr?s ordinaires :

\[\widehat Y(x) = \widehat\beta_0 + \widehat\beta_1x\]

---

Exemple

\(\widehat\beta_0\) correspond ? l'ordonn?e ? l'origine. Ce param?tre n'est pas toujours interpr?table (il d?pend de la signifaction de \(X\) et du fait que \(X\) soit centr? ou non).

\(\widehat\beta_1\) correspond ? la pente.

Attention, la technique des MCO cr?e des estimateurs sensibles aux valeurs atypiques.

layout(1, respect=TRUE)
plot(d$taille, d$poids, main="poids~taille")
## Model
# Y~X
# E(Y) = b.X
# E(Y) = b_0 + b_1.X
# Y_i = b_0 + b_1.X_i + e_i
m = lm(d$poids~d$taille)
m$coefficients
## (Intercept)    d$taille 
## -64.6134804   0.8006312
abline(a=m$coef[[1]], b=m$coef[[2]], col=2) # /!\ y = b.x + a
# residuals
# head(m$residuals)
arrows(d$taille, d$poids, d$taille, d$poids-m$residuals, col=adjustcolor(4, alpha.f=0.5), length=0.1)
legend("topleft",c("regression line", "residuals"), col=c(2,4), lty=1, cex=.8)

---

1.3 Qualit? du mod?le

1.3.1 D?composition de la variance

La variation de \(Y\) vient du fait de sa d?pendance ? la variable explicative \(X\). C'est la variation expliqu?e par le mod?le.

Dans l'exemple \(poids \sim taille\), nous avons remarqu? que lorsque nous mesurons \(Y\) avec une m?me valeur de \(X\), nous observons une certaine variation sur \(Y\). C'est la variation inexpliqu?e par le mod?le.

On dit que la variation totale de Y se d?compose en Variation expliqu?e par le mod?le et variation inexpliqu?e par le mod?le

Soit :

\[\begin{eqnarray} Y_i - \overline{Y}_n &=& Y_i - \widehat Y_i + \widehat Y_i - \overline{Y}_n \\ &=& (Y_i - \widehat Y_i) + (\widehat Y_i - \overline{Y}_n) \end{eqnarray}\]

avec \(\widehat Y_i - \overline{Y}_n\) la diff?rence expliqu?e par le mod?le et \(Y_i - \widehat Y_i\) la diff?rence inexpliqu?e (ou r?siduelle).

La m?thode des moindres carr?s, en consid?rant la somme des carr?s de ces diff?rences, nous assure la d?composition de la variance totale en variance expliqu?e et variance r?siduelle :

\[\begin{eqnarray} \sum(Y_i - \overline{Y}_n)^2 &=& \sum((Y_i - \widehat Y_i) + (\widehat Y_i - \overline{Y}_n))^2 \\ &=& \sum(Y_i - \widehat Y_i)^2 + \sum(\widehat Y_i - \overline{Y}_n)^2 + 2 \sum(Y_i - \widehat Y_i) (\widehat Y_i - \overline{Y}_n) \\ &=& \sum(\widehat Y_i - \overline{Y}_n)^2 + \sum(Y_i - \widehat Y_i)^2 \end{eqnarray}\]

\[ \]

Soit :

\[ SC_{tot} = SC_{reg} + SC_{res}\]

Avec \(SC_{tot}\) la somme des carr?s totale, \(SC_{reg}\) la somme des carr?s due ? la r?gression et \(SC_{res}\) la somme des carr?s des r?sidus.

Notons que \(2 \sum(Y_i - \widehat Y_i) (\widehat Y_i - \overline{Y}_n)=0\) en raison des moindres carr?s [demo].

1.3.2 Coefficient de determination \(R^2\)

La mesure du pourcentage de la variation totale expliqu?e par le mod?le se fait par le coefficient de d?termination.

\[ R^2 = \frac{SC_{reg}}{SC_{tot}}\]

• \(R^2\) est compris entre 0 et 1.
• \(R^2 = 1\) : cas o? les donn?es sont parfaitement align?es (comme c'est le cas pour un mod?le d?terministe).
• \(R^2 = 0\) : cas o? la variation de \(Y\) n'est pas due ? la variation de \(X\). Les donn?es ne sont pas du tout align?es.
• Plus \(R^2\) est proche de 1, plus les donn?es sont align?es sur la droite de r?gression.

1.3.3 Travaux pratiques

Jeu de donn?es CO2

data(CO2)
head(CO2)
##   Plant   Type  Treatment conc uptake
## 1   Qn1 Quebec nonchilled   95   16.0
## 2   Qn1 Quebec nonchilled  175   30.4
## 3   Qn1 Quebec nonchilled  250   34.8
## 4   Qn1 Quebec nonchilled  350   37.2
## 5   Qn1 Quebec nonchilled  500   35.3
## 6   Qn1 Quebec nonchilled  675   39.2
layout(matrix(1:2, 1), respect=TRUE)
plot(
  CO2[CO2$Type=="Quebec", ]$conc, 
  CO2[CO2$Type=="Quebec", ]$uptake, 
  xlab="food", ylab="growth", main="growth~food"
)
plot(
  log10(CO2[CO2$Type=="Quebec", ]$conc), 
  CO2[CO2$Type=="Quebec", ]$uptake, 
  xlab="log(food)", ylab="growth", main="growth~log(food)"
)

On s'int?resse uniquement aux plantes dont le Type est Quebec.

• Calculer les mod?les de regression lin?aire m1: growth~food et m2: growth~log(food)
• Tracer les nuages de points et les droites de regression correspondantes.
• Calculer, comparer et commenter le pourcentage de variance expliqu?e de chaque mod?le.
• Comparer les co?ficiants de chacun des mod?les.

Jeu de donn?es anscombe

data(anscombe)
print(anscombe)
##    x1 x2 x3 x4    y1   y2    y3    y4
## 1  10 10 10  8  8.04 9.14  7.46  6.58
## 2   8  8  8  8  6.95 8.14  6.77  5.76
## 3  13 13 13  8  7.58 8.74 12.74  7.71
## 4   9  9  9  8  8.81 8.77  7.11  8.84
## 5  11 11 11  8  8.33 9.26  7.81  8.47
## 6  14 14 14  8  9.96 8.10  8.84  7.04
## 7   6  6  6  8  7.24 6.13  6.08  5.25
## 8   4  4  4 19  4.26 3.10  5.39 12.50
## 9  12 12 12  8 10.84 9.13  8.15  5.56
## 10  7  7  7  8  4.82 7.26  6.42  7.91
## 11  5  5  5  8  5.68 4.74  5.73  6.89
layout(matrix(1:4,2), respect=TRUE)
plot(anscombe$x1, anscombe$y1)
plot(anscombe$x2, anscombe$y2)
plot(anscombe$x3, anscombe$y3)
plot(anscombe$x4, anscombe$y4)

• Calculer les mod?les de regression lin?aire m1 : y1~x1, m2 : y2~x2, m3 : y3~x3, m4 : y4~x4.
• Tracer les nuages de points et les droites de regression correspondantes.
• Calculer, comparer et commenter les co?ficiants de chacun des mod?les.
• Calculer, comparer et commenter le pourcentage de variance expliqu?e de chaque mod?le.

---

1.4 Test statistique

Il est important de noter que la construction du mod?le de r?gression et l'estimation des param?tres par MCO ne fait pas appel aux hypoth?ses de distribution.

Les hypoth?ses de distribution sont essentielles lorsqu'il s'agit de construire des tests et des intervalles de confiance et de pr?diction.

---

1.4.1 Test de Student

Pour rappel, l'objectif d'un test statitsique est de :

1. valider une hypoth?ses,
2. quantifier le risque associ? ? cette d?cision.

Souvent, on cherche ? invalider l'hypoth?se d?te nulle (\(\mathcal{H}_0\)) et donc valider l'hypoth?se d?te alternative (\(\mathcal{H}_1\)) au risque de premi?re esp?ce \(\alpha\) de 5%.

Conditions d'applications du test

• Les observations sont ind?pendantes
• La variance des erreurs est constante
• La loi des erreurs est une loi normale

Hypoth?se nulle et hypoth?se alternative

On se pose la question de l'effet de la variable \(x\) :

\(\mathcal{H}_0 : {\beta_1 = 0}\) contre \(\mathcal{H}_1 : {\beta_1 \neq 0}\)

Risque de premi?re esp?ce \(\alpha\)

Le risque de premi?re esp?ce ou risque \(\alpha\) est le risque de rejeter l'hypoth?se nulle \(\mathcal{H}_0\) alors que celle-ci est vraie.

Notons que le risque de seconde esp?ce ou risque \(\beta\) est le risque de conserver \(\mathcal{H}_0\) alors que celle-ci est fausse. Attention, le calcul du risque \(\beta\) est plus compliqu?. Souvent en minimisant \(\alpha\) on augmente \(\beta\). Pour s'assurer que \(\beta\) reste faible il faut g?n?ralement augmenter \(n\).

Statistique du test

Sous \(\mathcal{H}_0\),

\[t_{\widehat\beta_1,n-2} = \frac{\widehat\beta_1}{s_{\widehat\beta_1}} \sim \mathcal{T}_{n-2}\]

avec la variance r?siduelle estim?e :

\[ \widehat\sigma^2_n = \frac{1}{n-2}\sum^n_{i=1} (Y_i - \widehat\beta_0 - \widehat\beta_1x_i)^2 = \frac{1}{n-2}\sum^n_{i=1} \widehat\epsilon_i^2 \]

et l'erreur standard sur \(\beta_1\) estim?e :

\[ s^2_{\widehat\beta_1} = \frac{\widehat\sigma^2_n}{ns^2_x} = \frac{\widehat\sigma^2_n}{\sum^n_{i=1} (x_i - \overline x)^2} \]

D?cision et conclusion du test

La valeur critique du test, not?e \(c_\alpha\) est lue dans une table de la loi de Student.

• Si la valeur absolue de la valeur de la statistique calcul?e sur l'echantillon, not?e \(t_{\widehat\beta_1,n-2}\) est sup?rieure ou ?gale ? \(c_\alpha\), alors le test est significatif. Vous rejetez \(\mathcal{H}_0\) et vous d?cidez que \(\mathcal{H}_1\) est vraie avec un risque d'erreur de premi?re esp?ce \(\alpha\).

• Si la valeur absolue de la valeur de la statistique calcul?e sur l'echantillon, not?e \(t_{\widehat\beta_1,n-2}\) est strictement inf?rieure ? \(c_\alpha\), alors le test n'est pas significatif. Vous conservez \(\mathcal{H}_0\) avec un risque d'erreur de deuxi?me esp?ce \(\beta\) ? ?valuer.

Mais alors pourquoi ?

Si l'hypoth?se nulle \(\mathcal{H}_0\) est v?rifi?e, i.e. \(\beta_1=0\), alors la variable al?atoire \(t_{\widehat\beta_1,n-2} = \frac{\widehat\beta_1}{s_{\widehat\beta_1}}\), i.e, l'effet de la variable \(x\) extim? ? partir des observations "r?duit" de la "variabilit?" observ?e sur cet estimateur, suit une la loi de Student ? \(n-2\) d?gr? de libert? not?e \(\mathcal{T}_{n-2}\) en raison des conditions d'applications du test.

Exemple 0

m = lm(poids~taille, d)
# Student test
summary(m)
## 
## Call:
## lm(formula = poids ~ taille, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.488  -5.836   1.911   4.613  22.509 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -64.6135    39.2021  -1.648  0.11665   
## taille        0.8006     0.2387   3.354  0.00354 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.921 on 18 degrees of freedom
## Multiple R-squared:  0.3846, Adjusted R-squared:  0.3504 
## F-statistic: 11.25 on 1 and 18 DF,  p-value: 0.003535
n = length(m$residuals)
xi = m$model$taille
s_beta1 = sqrt( sum(m$residuals^2)/(n-2) / sum((xi-mean(xi))^2))
s_beta1
## [1] 0.238727
stat_t =  m$coefficient[2]/s_beta1
stat_t
##   taille 
## 3.353753
pval = pt(-stat_t, n-2) + 1-pt(stat_t, n-2)
pval
##      taille 
## 0.003535466

Exemple 1

set.seed(1)
d = MASS::cats[sample(1:nrow(MASS::cats),10),]
m = lm(Hwt~Bwt, d)
# Student test
summary(m)
## 
## Call:
## lm(formula = Hwt ~ Bwt, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3417 -0.6850 -0.0854  0.6825  1.3884 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.2123     2.0167  -0.601 0.564390    
## Bwt           4.1746     0.7697   5.424 0.000628 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9641 on 8 degrees of freedom
## Multiple R-squared:  0.7862, Adjusted R-squared:  0.7595 
## F-statistic: 29.42 on 1 and 8 DF,  p-value: 0.000628
n = length(m$residuals)
xi = m$model$Bwt
s_beta1 = sqrt( sum(m$residuals^2)/(n-2) / sum((xi-mean(xi))^2))
s_beta1
## [1] 0.7696931
stat_t =  m$coefficient[2]/s_beta1
stat_t
##      Bwt 
## 5.423763
pval = pt(-stat_t, n-2) + 1-pt(stat_t, n-2)
pval
##          Bwt 
## 0.0006279806

Conditions d'applications du test

• Les observations sont ind?pendantes
• La variance des erreurs est constante \(\sigma^2\)
• La loi des erreurs est une loi normale \(\mathcal{N}(0, \sigma^2)\)

Pour tester la normalit? des r?sidus :

• Test de Shapiro-Wilk
• Test de Anderson-Darling
• Test de Jarque-Bera (package tseries, grands ?chantillons)
• Test de Kolmogorov-Smirnov (comparaison de distributions)
• QQ plot

Pour tester l'homoscedasticit? :

• Test de Goldfeld-Quandt
• Test de Hartley (??)
• Test de Bartlett (ANOVA)

shapiro.test(rnorm(100))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(100)
## W = 0.99612, p-value = 0.9942
shapiro.test(runif(100))
## 
##  Shapiro-Wilk normality test
## 
## data:  runif(100)
## W = 0.93637, p-value = 0.0001165
shapiro.test(m$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  m$residuals
## W = 0.96881, p-value = 0.8796
ks.test(rnorm(100), rnorm(100))
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(100) and rnorm(100)
## D = 0.08, p-value = 0.9062
## alternative hypothesis: two-sided
ks.test(rnorm(100), runif(100))
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(100) and runif(100)
## D = 0.58, p-value = 4.885e-15
## alternative hypothesis: two-sided
ks.test(rnorm(100), m$residuals)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(100) and m$residuals
## D = 0.13, p-value = 0.9941
## alternative hypothesis: two-sided
layout(matrix(1:6,2), respect=TRUE)
set.seed(1)
plot(sort(rnorm(10)  ), scale(sort(runif(10)  )))
abline(0, 1, col=2)
plot(sort(rnorm(10)  ), sort(rnorm(10)  ))
abline(0, 1, col=2)
plot(sort(rnorm(100) ), scale(sort(runif(100) )))
abline(0, 1, col=2)
plot(sort(rnorm(100) ), sort(rnorm(100) ))
abline(0, 1, col=2)
plot(sort(rnorm(1000)), scale(sort(runif(1000))))
abline(0, 1, col=2)
plot(sort(rnorm(1000)), sort(rnorm(1000)))
abline(0, 1, col=2)

layout(matrix(1:3,1), respect=TRUE)
qqnorm(scale(rnorm(100)))
abline(0, 1, col=2)
qqnorm(scale(runif(100)))
abline(0, 1, col=2)
qqnorm(scale(m$residuals))
abline(0, 1, col=2)

Travaux pratiques

• Construire le mod?le de r?gr?ssion lin?aire sur la totalit? des donn?es cats.
• Tracer la droite de r?gression correspondante.
• Tester l'hypoth?se \(\mathcal{H}_0 : {\beta_1 = 0}\) contre \(\mathcal{H}_1 : {\beta_1 \neq 0}\) sur ce mod?le.
• Comparer la p-valeur obtenue avec ce mod?le avec celle obtenue sur le mod?le du cours (30 observations).

layout(1, respect=TRUE)
plot(
  MASS::cats$Bwt, MASS::cats$Hwt, 
  xlab="Body weight", ylab="Heart weight", 
  main="Heart weight ~ Body weight"
)

---

1.4.2 Intervalle de confiance

On peut construire les intervalles de confiance suivants

\[IC_{1-\alpha}(\widehat\beta_1) = \Big]\widehat\beta_1 - t_{n-2;1-\alpha/2}*s_{\beta_1}; \widehat\beta_1 + t_{n-2;1-\alpha/2}*s_{\beta_1}\Big[\]

\[IC_{1-\alpha}(\widehat\beta_0) = \Big]\widehat\beta_0 - t_{n-2;1-\alpha/2}*s_{\beta_0}; \widehat\beta_0 + t_{n-2;1-\alpha/2}*s_{\beta_0}\Big[\]

confint(m)
##                 2.5 %   97.5 %
## (Intercept) -5.862783 3.438181
## Bwt          2.399718 5.949549

Ainsi que l'intervalle de pr?diction d'une valeur moyenne de \(Y\), sachant que \(x = x_0\).

L'estimation ponctuelle pour cette valeur de \(x_0\) est alors ?gale a \(\widehat Y_0 = \widehat\beta_0 + \widehat\beta_1x_0\)

\[IP_{1-\alpha}(Y) = \Big]\widehat Y_0 - t_{n-2;1-\alpha/2}\sqrt{\widehat\sigma^2_n\Big(\frac{1}{n}+\frac{(x_0-\overline{x}_n)^2}{(n-1)s^2_x}\Big)}; \widehat Y_0 + t_{n-2;1-\alpha/2}*\sqrt{\widehat\sigma^2_n\Big(\frac{1}{n}+\frac{(x_0-\overline{x}_n)^2}{(n-1)s^2_x}\Big)}\Big[\]

layout(1, respect=TRUE)
plot(
  d$Bwt, d$Hwt, 
  xlab="Body weight", ylab="Heart weight", 
  main="Heart weight ~ Body weight"
)
abline(m, col=2)
x = seq(min(m$model$Bwt), max(m$model$Bwt), length.out=100)
x = seq(2.5, 3, length.out=100)
pred =  predict(m, interval="predict", level=0.95)
## Warning in predict.lm(m, interval = "predict", level = 0.95): predictions on current data refer to _future_ responses
head(pred)
##           fit       lwr      upr
## 68   9.224283  6.887049 11.56152
## 129 12.981453 10.242097 15.72081
## 43  10.894136  8.498329 13.28994
## 14   7.971893  5.539546 10.40424
## 51   7.971893  5.539546 10.40424
## 85  10.059210  7.719281 12.39914
points(m$model$Bwt, pred[,2], col=2, pch=3)
points(m$model$Bwt, pred[,3], col=2, pch=3)
pred =  predict(m, interval="predict", level=0.99)
## Warning in predict.lm(m, interval = "predict", level = 0.99): predictions on current data refer to _future_ responses
points(m$model$Bwt, pred[,2], col=3, pch=3)
points(m$model$Bwt, pred[,3], col=3, pch=3)
pred =  predict(m, interval="predict", level=0.90)
## Warning in predict.lm(m, interval = "predict", level = 0.9): predictions on current data refer to _future_ responses
points(m$model$Bwt, pred[,2], col=4, pch=3)
points(m$model$Bwt, pred[,3], col=4, pch=3)
legend("topleft",c("95% ci", "99% ci", "90% ci"), col=2:4, pch=3)

Travaux pratiques

• Tracer l'intervalle de pr?diction du mod?le obtenu sur la totalit? des donn?es cats.
• Comparer la l'intervalle de pr?diction obtenu avec ce mod?le avec celui obtenu sur le mod?le du cours (30 observations).