TP1 : Régression linéaire avec R
2023-11-09
Berté Jean-Luc
L2 MGA
#TP1 : Régression Linéaire avec R
#1 Calcul des coefficients de la droite de régression
#On suppose qu’on dispose de n données (xi, yi) contenues dans un fichier .csv. On commence par #lire ce fichier en affectant
#les valeurs dans une data frame nommées ici ”données” puis on trace le nuage de points pour #vérifier visuellement la potentielle linéarité.
# Lecture des données dans un fichier csv (séparateur point virgule)
Droite <- read.csv("/cloud/project/Droite.csv", sep=";")
print(t(Droite))# x=variable explicative,y=variable expliquée## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## Taille 1.56 1.63 1.67 1.67 1.69 1.64 1.63 1.56 1.67 1.72 1.67 1.56
## Poids 47.00 44.00 52.00 61.00 54.00 47.00 44.00 47.00 52.00 57.00 52.00 47.00
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## Taille 1.69 1.7 1.63 1.67 1.55 1.63 1.64 1.56 1.75 1.72 1.67 1.68
## Poids 54.00 52.0 53.00 52.00 46.00 44.00 47.00 47.00 54.00 56.00 52.00 55.00
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## Taille 1.72 1.62 1.69 1.6 1.56 1.65 1.79 1.69 1.73 1.64 1.63 1.63
## Poids 62.00 48.00 57.00 45.0 45.00 58.00 64.00 51.00 54.00 47.00 52.00 53.00
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## Taille 1.7 1.62 1.67 1.69 1.63 1.6 1.73 1.61 1.56 1.48 1.6 1.69
## Poids 52.0 50.00 47.00 54.00 52.00 44.0 57.00 43.00 47.00 38.00 44.0 54.00
## [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
## Taille 1.64 1.56 1.63 1.55 1.67 1.73 1.67 1.69 1.63 1.67 1.56 1.68
## Poids 47.00 45.00 44.00 47.00 47.00 57.00 61.00 51.00 52.00 52.00 47.00 50.00
## [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
## Taille 1.63 1.6 1.75 1.73 1.75 1.61 1.67 1.62 1.69 1.64 1.62 1.6
## Poids 44.00 48.0 54.00 66.00 53.00 50.00 47.00 56.00 57.00 52.00 52.00 50.0
## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84]
## Taille 1.72 1.64 1.7 1.7 1.65 1.67 1.65 1.62 1.57 1.6 1.6 1.6
## Poids 56.00 65.00 55.0 59.0 55.00 50.00 55.00 47.00 48.00 48.0 48.0 44.0
## [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96]
## Taille 1.71 1.67 1.6 1.79 1.64 1.7 1.61 1.68 1.68 1.73 1.68 1.74
## Poids 53.00 52.00 49.0 64.00 56.00 51.0 51.00 50.00 56.00 54.00 55.00 70.00
## [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] [,105] [,106]
## Taille 1.6 1.48 1.65 1.67 1.67 1.65 1.7 1.67 1.62 1.7
## Poids 43.0 38.00 54.00 55.00 55.00 48.00 52.0 54.00 48.00 55.0
## [,107] [,108] [,109] [,110] [,111] [,112] [,113] [,114] [,115] [,116]
## Taille 1.6 1.51 1.67 1.64 1.68 1.55 1.65 1.63 1.66 1.68
## Poids 45.0 54.00 48.00 50.00 63.00 53.00 48.00 59.00 53.00 51.00
## [,117] [,118] [,119] [,120] [,121] [,122] [,123] [,124] [,125] [,126]
## Taille 1.65 1.64 1.7 1.71 1.7 1.62 1.5 1.6 1.65 1.69
## Poids 60.00 59.00 60.0 57.00 53.0 60.00 40.0 50.0 53.00 58.00
## [,127] [,128] [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136]
## Taille 1.6 1.65 1.67 1.65 1.55 1.65 1.56 1.63 1.65 1.72
## Poids 42.0 44.00 60.00 63.00 43.00 61.00 45.00 52.00 56.00 56.00
## [,137] [,138] [,139] [,140] [,141] [,142] [,143] [,144]
## Taille 1.58 1.62 1.61 1.55 1.61 1.72 1.64 1.76
## Poids 47.00 48.00 56.00 49.00 43.00 56.00 63.00 58.00# t pour transposer le tableau
# pour accéder aux variables de donnees par leur nom (seul) (simplifier leur nom)
attach(Droite)# par ex donnees$x devient x
# tracé: nuage de points
eq = paste0("Nuage de points: observations")
plot(Taille, Poids,type="p",pch=4,main=eq,col="blue")
#y = a1 + a2x #dont il faut déterminer les coefficients. #1.1 Calcul matriciel #La méthode des moindres carrés consiste à déterminer les valeurs a1 et a2 qui minimisent les #écarts au carré entre les #valeurs expérimentales et les valeurs prédites par le modèle (écarts aussi appelés résidus). #C’est un problème de recherche #du minimum d’une fonction de deux variables a1 et a2 (cf cours): #S(a1, a2) = Xn #i=1 #(a1 + a2xi − yi) #2 #Cela conduit à un système linéaire qu’il suffit d’inverser
y=a1+a2*x détermination de a1 et a2 par la méthode des moindres carrés
Construction de X et Y (cf cours)
c1=rep(1,144)
X<-matrix(c(c1,Taille),ncol = 2)
Y<-matrix(Poids,ncol=1)
M=t(X)%*%X # transposée(X) *X
b<-t(X)%*%Y # transposée(X) *Y
a=solve(M,b) # inversion du système donne le vecteur a=les coefficients a1 et a2
plot(Taille, Poids,type="p",pch=4,main=eq,col="blue")# tracé du nuage de points (x,y)
abline(a, col="red",main=eq)# tracé droite de regression (courbe de tendance linéaire)
# equation de la droite sous forme de chaine de caractères
eq = paste0("Droite regression: y = ", round(a[2],3), "*x +",round(a[1],3))
title(sub=eq)# ajout d'un sous titre contenant l'équation de la droite
#Nous avons ainsi déterminé la droite de régression par la méthode des moindres carrés c’est à-#dire en minimisant les #écarts au carré entre les points expérimentaux et la droite. #1.2 Coefficients de la droite de régression avec lm #Il existe sous R une fonction qui permet le calcul des coefficients de la droite de régression #sans avoir besoin de programmer #les matrices X et Y (la fonction le fait elle mˆeme). C’est la fonction lm (pour Linear #Models).
droite <- lm(Poids~Taille, data=Droite)
# y~x signifie chercher une relation entre y et x du type y=a1 + a2*x sur les données de la #data frame donnees
print(droite) # affichage du résultat nommé ici droite qui est une liste##
## Call:
## lm(formula = Poids ~ Taille, data = Droite)
##
## Coefficients:
## (Intercept) Taille
## -56.15 65.59a1 <- coef(droite)[1] # ordonnée à l'origine appelé intercept
a2 <- coef(droite)[2] # la pente
# tracé: nuage de points avec un titre principal et en sous titre l'équation
eq = paste0("y = ", round(a2,3), "*x +",round(a1,3))
plot(Taille, Poids,main="Droite de régression", sub=eq)
# tracé: droite de regression (courbe de tendance linéaire)
abline(droite, col="red")
#On retrouve bien les mêmes coefficients qu’avec le calcul matriciel (évidemment, c’est le même #calcul !). Mais on obtient en plus d’autres résultats. #2 Approche probabiliste: qualité de la prédiction #Si on veut estimer la qualité de la prédiction, on peut avoir une approche probabiliste: elle #consiste à faire l’hypothèse #que dans les observations notées (xi, yi), les xi sont des valeurs exactes et certaines d’une #variable non aléatoire x et les #yi sont des réalisations d’une variable aléatoire notée Y . #On fait la supposition que cette V.A. vérifie #Y = α1 + α2x + #où l’erreur est une variable aléatoire vérifiant ∼ N (0, σ2). On peut montrer que les # paramètres, a1 et a2, déterminés #par les moindres carrés, sont alors les meilleurs estimateurs sans biais de α1 et α2 #1 de la méthode de maximum de vraissemblance. On peut utiliser cette régression pour répondre #à certaines questions statistiques (signification, comparaison, #prédiction, …) que nous allons détailler. Mais on doit vérifier que ∼ N (0, σ) en étudiant #la répartition des réalisations #de la variable aléatoire, les résidus ri (où ri = a1 + a2xi − yi). #2.1 Qualité de la prédiction: R2, résidus, écart-type, tests… #2.1.1 Analyse de la variance #Lors de l’appel de lm, l’analyse de la variance a été faite
anova(droite)# table d'analyse de la variance## Analysis of Variance Table
##
## Response: Poids
## Df Sum Sq Mean Sq F value Pr(>F)
## Taille 1 2129.5 2129.46 103.81 < 2.2e-16 ***
## Residuals 142 2912.8 20.51
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # Response: y
# Df Sum Sq Mean Sq F value Pr(>F)
# x Df1 SSR MS1=SSR/Df1 MS1/MS2 pvalue prob(F(alpha,Df1,Df2)>F)
# Residuals Df2 SSE MS2=SSE/Df2
# H0 "pente nulle" est rejetée si pvalue<risque_alpha
# SSR=sum square regression (variance expliquée par le modèle)
# SSE= sum square error (variance non expliquée par le modèle)
# SST=SSE+SSR (variance totale)
# Dans la table, on obtient en colonnes
# les degrés de liberté ”Df” (pour les lois des tests sur les paramètres)
# les sommes des carrés (Sum Sq): Sum Square Regression: variance expliquée par le modèle et #Sum Square Erreur:
#variance résiduelle non expliquée par le modèle.
# les moyennes de sommes des carrés (Mean Sq): SSR/Df1 et SSE/Df2 (estimation de la variance σ
#2 de l’erreur )
# la statistique (F value) (quotient de deux valeurs précédentes-suit une distribution de #Fisher F1,n−2).
# la p-value (Pr(>F)) associée au test d’hypothèse H0 : ”la pente a2 est nulle” contre H1 : #”la pente n’est pas nulle”.
#On rejette H0 au risque α, si la p-value est inférieure à α. Dans notre cas, par exemple, on #rejette l’hypothèse ”la pente est nulle”.
#2.1.2 Commande Summary
#On peut obtenir d’autres informations statistiques, en tapant la commande summary().
summary(droite)# résumé des calculs effectués par l'instruction lm##
## Call:
## lm(formula = Poids ~ Taille, data = Droite)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0756 -3.3553 -0.6672 1.9002 13.5803
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -56.154 10.611 -5.292 4.48e-07 ***
## Taille 65.594 6.438 10.189 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.529 on 142 degrees of freedom
## Multiple R-squared: 0.4223, Adjusted R-squared: 0.4183
## F-statistic: 103.8 on 1 and 142 DF, p-value: < 2.2e-16 # Coefficients:
# Estimation a1, sa1=ecart-type de a1, t=a1/sa1, probabilité critique PCa1
# H0: la droite passe par l'origine"vs H1 la droite ne passe pas pas l'origine
# H0 est rejetée si t<student(n-2,1-risque/2) ie si la pvalue PCa1<risque.
# Estimation a2, sa2=ecart-type de a, t=a2/sa2, probabilité critique PCa2
# H0: la pente de la droite est nulle" rejetée si t<student(n-2,1-risque/2)
# H0: la pente de la droite est nulle" rejetée si la pvalue PCa2<risque.
# Residual standard error: SQRT(MS2)
# doit être faible pour que le modèle soit considéré comme bon (=prédictif)
# Multiple R-squared: SSR/SST,
# Plus cette valeure sera proche de 1 meilleur sera l'ajustement.
# Adjusted R-squared: ajustement du R^2 au nombre p de variables explicatives
#F-statistic: MS1/MS2 on Df1 and Df2, p-value
# Dans la 1 ère table ”Residuals”, on trouve les minimum et maximum ainsi que les quartiles des #résidus. On peut aussi
#afficher les valeurs des résidus avec la commande
resid(droite)# les résidus## 1 2 3 4 5 6 7
## 0.8277810 -6.7637754 -1.3875219 7.6124781 -0.6993951 -4.4197120 -6.7637754
## 8 9 10 11 12 13 14
## 0.8277810 -1.3875219 0.3327950 -1.3875219 0.8277810 -0.6993951 -3.3553317
## 15 16 17 18 19 20 21
## 2.2362246 -1.3875219 0.4837176 -6.7637754 -4.4197120 0.8277810 -4.6350149
## 22 23 24 25 26 27 28
## -0.6672050 -1.3875219 0.9565415 5.3327950 -2.1078387 2.3006049 -3.7959655
## 29 30 31 32 33 34 35
## -1.1722190 5.9243514 2.7412386 -3.6993951 -3.3231416 -4.4197120 1.2362246
## 36 37 38 39 40 41 42
## 2.2362246 -3.3553317 -0.1078387 -6.3875219 -0.6993951 1.2362246 -4.7959655
## 43 44 45 46 47 48 49
## -0.3231416 -6.4519021 0.8277810 -2.9247260 -4.7959655 -0.6993951 -4.4197120
## 50 51 52 53 54 55 56
## -1.1722190 -6.7637754 1.4837176 -6.3875219 -0.3231416 7.6124781 -3.6993951
## 57 58 59 60 61 62 63
## 1.2362246 -1.3875219 0.8277810 -4.0434585 -6.7637754 -0.7959655 -4.6350149
## 64 65 66 67 68 69 70
## 8.6768584 -5.6350149 0.5480979 -6.3875219 5.8921613 2.3006049 0.5802880
## 71 72 73 74 75 76 77
## 1.8921613 1.2040345 -0.6672050 13.5802880 -0.3553317 3.6446683 2.9243514
## 78 79 80 81 82 83 84
## -3.3875219 2.9243514 -3.1078387 1.1718444 -0.7959655 -0.7959655 -4.7959655
## 85 86 87 88 89 90 91
## -3.0112684 -1.3875219 0.2040345 2.7412386 4.5802880 -4.3553317 1.5480979
## 92 93 94 95 96 97 98
## -4.0434585 1.9565415 -3.3231416 0.9565415 12.0209218 -5.7959655 -2.9247260
## 99 100 101 102 103 104 105
## 1.9243514 1.6124781 1.6124781 -4.0756486 -3.3553317 0.6124781 -2.1078387
## 106 107 108 109 110 111 112
## -0.3553317 -3.7959655 11.1074641 -5.3875219 -1.4197120 8.9565415 7.4837176
## 113 114 115 116 117 118 119
## -4.0756486 8.2362246 0.2684148 -3.0434585 7.9243514 7.5802880 4.6446683
## 120 121 122 123 124 125 126
## 0.9887316 -2.3553317 9.8921613 -2.2365993 1.2040345 0.9243514 3.3006049
## 127 128 129 130 131 132 133
## -6.7959655 -8.0756486 6.6124781 10.9243514 -2.5162824 8.9243514 -1.1722190
## 134 135 136 137 138 139 140
## 1.2362246 3.9243514 -0.6672050 -0.4840923 -2.1078387 6.5480979 3.4837176
## 141 142 143 144
## -6.4519021 -0.6672050 11.5802880 -1.2909515#Dans la seconde table ”Coefficients”
# la colonne ”estimate” contient les estimations de l’ordonnée à l’origine, a1 #appelé(intercept) puis de la pente, a2,
# leurs écarts-types (Std. Error) respectifs,
# pour chaque paramètre, la statistique observée (t value) ainsi que la p-value (Pr(>|t|)) #associée au test d’hypothèse
#H0 : ”le paramètre est nul” contre H1 : ”le paramètre n’est pas nul”. On rejette H0 au risque #α, si la p-value est
#inférieure à α. Dans notre cas, par exemple, on rejette les hypothèses ”la pente est nulle” et #”la droite passe par
#l’origine”.
#Ensuite, on trouve
# ”Residual standard error” l’estimation sr de l’écart-type σ de la VA o`u sr =
#qSSE
#Df2
# ”Multiple R-squared” le coefficient de détermination R2 =
#SSR
#SST qui exprime le rapport entre la variance expliquée
#par le modèle et la variance totale. Ce coefficient varie entre 0 et 1. S’il est égal à 1, #cela signifie que le modèle
#explique parfaitement les données.
# ”Adjusted R-squared” le coefficient de détermination ajusté (prend en compte le nombre de #variables)
# la statistique F déjà définie dans l’anova.
#2.1.3 Analyse des résidus
#Une étape fondamentale dans la démarche de la régression est l’étude des résidus. En effet, #les résultats statistiques donnés
#précédemment sont basés sur l’hypotèse que les érreurs (les résidus) sont indépendantes, #normalement distribuées, de
#moyenne nulle et de variance constante. La fonction lm fournit quatre graphiques qui #permettent d’en juger.
#layout(matrix(1:4,2,2))# fenètre graphique coupée en 4
#plot(droite) # 4 graphiques
#Le premier graphe ”Residual vs fitted” donne les résidus en fonction des valeurs prédites. Les #points doivent être
#répartis aléatoirement autour de l’axe horizontal y = 0 et ne pas montrer de tendance. Le #graphe ”Scale-Location”
#donne la racine des résidus standardisés en fonction des valeurs prédites et ne doit pas non #plus montrer de tendance.
#Le graphe ”Normal Q-Q” permet de vérifier la normalité des résidus en comparant les quantiles #de la population avec
#ceux de la loi normale. Le dernier graphique ”Residuals vs Leverage” met en valeur #l’importance de chaque point dans la regression. On se questionnera particulièrement sur les #points ayant une distance de Cook supérieure à 1 (rend la donnée
#suspecte=>point aberrant ?).
#On peut aussi afficher les résidus studentisés qui en pratique devront ˆetre compris entre les #bornes −2 et 2.
# résidus studentisés
res <- rstudent( droite)
# tracé des résidus studentisés (ylim échelle pour y)
plot(res,ylab="Résidus studentisés",xlab="",main="Résidus de student",ylim=c(-2.5,2.5))
# tracés des droites y=0 (trait plein bleu), y=+/- 2 (trait espacé rouge)
# h= contient équation y=(-2,0,2) et lty le type de lignes
abline(h=c(-2,0,2),lty=c(2,1,2),col=c("red","blue","red"))
#En conclusion, les points suspects sont les points dont le résidu studentisé est supérieur à 2 #en valeur absolue et/ou la
#distance de Cook est supérieure à 1. Dans ce dernier cas, le point contribue très/trop #fortement à la détermination des
#coefficients du modèle comparativement aux autres. Il n’y a pas de méthode universelle pour #traiter ce type de points.
#2.1.4 Les prédictions: intervalle de confiance et de prédiction
#On peut obtenir les valeurs prédites par la commande fitted
fitted(droite)# les valeurs prédites## 1 2 3 4 5 6 7 8
## 46.17222 50.76378 53.38752 53.38752 54.69940 51.41971 50.76378 46.17222
## 9 10 11 12 13 14 15 16
## 53.38752 56.66720 53.38752 46.17222 54.69940 55.35533 50.76378 53.38752
## 17 18 19 20 21 22 23 24
## 45.51628 50.76378 51.41971 46.17222 58.63501 56.66720 53.38752 54.04346
## 25 26 27 28 29 30 31 32
## 56.66720 50.10784 54.69940 48.79597 46.17222 52.07565 61.25876 54.69940
## 33 34 35 36 37 38 39 40
## 57.32314 51.41971 50.76378 50.76378 55.35533 50.10784 53.38752 54.69940
## 41 42 43 44 45 46 47 48
## 50.76378 48.79597 57.32314 49.45190 46.17222 40.92473 48.79597 54.69940
## 49 50 51 52 53 54 55 56
## 51.41971 46.17222 50.76378 45.51628 53.38752 57.32314 53.38752 54.69940
## 57 58 59 60 61 62 63 64
## 50.76378 53.38752 46.17222 54.04346 50.76378 48.79597 58.63501 57.32314
## 65 66 67 68 69 70 71 72
## 58.63501 49.45190 53.38752 50.10784 54.69940 51.41971 50.10784 48.79597
## 73 74 75 76 77 78 79 80
## 56.66720 51.41971 55.35533 55.35533 52.07565 53.38752 52.07565 50.10784
## 81 82 83 84 85 86 87 88
## 46.82816 48.79597 48.79597 48.79597 56.01127 53.38752 48.79597 61.25876
## 89 90 91 92 93 94 95 96
## 51.41971 55.35533 49.45190 54.04346 54.04346 57.32314 54.04346 57.97908
## 97 98 99 100 101 102 103 104
## 48.79597 40.92473 52.07565 53.38752 53.38752 52.07565 55.35533 53.38752
## 105 106 107 108 109 110 111 112
## 50.10784 55.35533 48.79597 42.89254 53.38752 51.41971 54.04346 45.51628
## 113 114 115 116 117 118 119 120
## 52.07565 50.76378 52.73159 54.04346 52.07565 51.41971 55.35533 56.01127
## 121 122 123 124 125 126 127 128
## 55.35533 50.10784 42.23660 48.79597 52.07565 54.69940 48.79597 52.07565
## 129 130 131 132 133 134 135 136
## 53.38752 52.07565 45.51628 52.07565 46.17222 50.76378 52.07565 56.66720
## 137 138 139 140 141 142 143 144
## 47.48409 50.10784 49.45190 45.51628 49.45190 56.66720 51.41971 59.29095#On peut obtenir de nouvelles prédictions avec intervalle de confiance.
newx=seq(1.5,10,1);
pc=predict(droite,data.frame(x= newx), level = 0.95, interval = "confidence")## Warning: 'newdata' had 9 rows but variables found have 144 rowsprint(pc)## fit lwr upr
## 1 46.17222 44.83550 47.50894
## 2 50.76378 49.98641 51.54114
## 3 53.38752 52.58678 54.18827
## 4 53.38752 52.58678 54.18827
## 5 54.69940 53.77528 55.62351
## 6 51.41971 50.66809 52.17134
## 7 50.76378 49.98641 51.54114
## 8 46.17222 44.83550 47.50894
## 9 53.38752 52.58678 54.18827
## 10 56.66720 55.47720 57.85721
## 11 53.38752 52.58678 54.18827
## 12 46.17222 44.83550 47.50894
## 13 54.69940 53.77528 55.62351
## 14 55.35533 54.35086 56.35981
## 15 50.76378 49.98641 51.54114
## 16 53.38752 52.58678 54.18827
## 17 45.51628 44.07222 46.96034
## 18 50.76378 49.98641 51.54114
## 19 51.41971 50.66809 52.17134
## 20 46.17222 44.83550 47.50894
## 21 58.63501 57.12844 60.14159
## 22 56.66720 55.47720 57.85721
## 23 53.38752 52.58678 54.18827
## 24 54.04346 53.18824 54.89867
## 25 56.66720 55.47720 57.85721
## 26 50.10784 49.28561 50.93007
## 27 54.69940 53.77528 55.62351
## 28 48.79597 47.83850 49.75343
## 29 46.17222 44.83550 47.50894
## 30 52.07565 51.32868 52.82262
## 31 61.25876 59.29371 63.22382
## 32 54.69940 53.77528 55.62351
## 33 57.32314 56.03153 58.61476
## 34 51.41971 50.66809 52.17134
## 35 50.76378 49.98641 51.54114
## 36 50.76378 49.98641 51.54114
## 37 55.35533 54.35086 56.35981
## 38 50.10784 49.28561 50.93007
## 39 53.38752 52.58678 54.18827
## 40 54.69940 53.77528 55.62351
## 41 50.76378 49.98641 51.54114
## 42 48.79597 47.83850 49.75343
## 43 57.32314 56.03153 58.61476
## 44 49.45190 48.56861 50.33520
## 45 46.17222 44.83550 47.50894
## 46 40.92473 38.67045 43.17900
## 47 48.79597 47.83850 49.75343
## 48 54.69940 53.77528 55.62351
## 49 51.41971 50.66809 52.17134
## 50 46.17222 44.83550 47.50894
## 51 50.76378 49.98641 51.54114
## 52 45.51628 44.07222 46.96034
## 53 53.38752 52.58678 54.18827
## 54 57.32314 56.03153 58.61476
## 55 53.38752 52.58678 54.18827
## 56 54.69940 53.77528 55.62351
## 57 50.76378 49.98641 51.54114
## 58 53.38752 52.58678 54.18827
## 59 46.17222 44.83550 47.50894
## 60 54.04346 53.18824 54.89867
## 61 50.76378 49.98641 51.54114
## 62 48.79597 47.83850 49.75343
## 63 58.63501 57.12844 60.14159
## 64 57.32314 56.03153 58.61476
## 65 58.63501 57.12844 60.14159
## 66 49.45190 48.56861 50.33520
## 67 53.38752 52.58678 54.18827
## 68 50.10784 49.28561 50.93007
## 69 54.69940 53.77528 55.62351
## 70 51.41971 50.66809 52.17134
## 71 50.10784 49.28561 50.93007
## 72 48.79597 47.83850 49.75343
## 73 56.66720 55.47720 57.85721
## 74 51.41971 50.66809 52.17134
## 75 55.35533 54.35086 56.35981
## 76 55.35533 54.35086 56.35981
## 77 52.07565 51.32868 52.82262
## 78 53.38752 52.58678 54.18827
## 79 52.07565 51.32868 52.82262
## 80 50.10784 49.28561 50.93007
## 81 46.82816 45.59498 48.06133
## 82 48.79597 47.83850 49.75343
## 83 48.79597 47.83850 49.75343
## 84 48.79597 47.83850 49.75343
## 85 56.01127 54.91749 57.10505
## 86 53.38752 52.58678 54.18827
## 87 48.79597 47.83850 49.75343
## 88 61.25876 59.29371 63.22382
## 89 51.41971 50.66809 52.17134
## 90 55.35533 54.35086 56.35981
## 91 49.45190 48.56861 50.33520
## 92 54.04346 53.18824 54.89867
## 93 54.04346 53.18824 54.89867
## 94 57.32314 56.03153 58.61476
## 95 54.04346 53.18824 54.89867
## 96 57.97908 56.58164 59.37651
## 97 48.79597 47.83850 49.75343
## 98 40.92473 38.67045 43.17900
## 99 52.07565 51.32868 52.82262
## 100 53.38752 52.58678 54.18827
## 101 53.38752 52.58678 54.18827
## 102 52.07565 51.32868 52.82262
## 103 55.35533 54.35086 56.35981
## 104 53.38752 52.58678 54.18827
## 105 50.10784 49.28561 50.93007
## 106 55.35533 54.35086 56.35981
## 107 48.79597 47.83850 49.75343
## 108 42.89254 40.99432 44.79075
## 109 53.38752 52.58678 54.18827
## 110 51.41971 50.66809 52.17134
## 111 54.04346 53.18824 54.89867
## 112 45.51628 44.07222 46.96034
## 113 52.07565 51.32868 52.82262
## 114 50.76378 49.98641 51.54114
## 115 52.73159 51.96779 53.49538
## 116 54.04346 53.18824 54.89867
## 117 52.07565 51.32868 52.82262
## 118 51.41971 50.66809 52.17134
## 119 55.35533 54.35086 56.35981
## 120 56.01127 54.91749 57.10505
## 121 55.35533 54.35086 56.35981
## 122 50.10784 49.28561 50.93007
## 123 42.23660 40.22075 44.25245
## 124 48.79597 47.83850 49.75343
## 125 52.07565 51.32868 52.82262
## 126 54.69940 53.77528 55.62351
## 127 48.79597 47.83850 49.75343
## 128 52.07565 51.32868 52.82262
## 129 53.38752 52.58678 54.18827
## 130 52.07565 51.32868 52.82262
## 131 45.51628 44.07222 46.96034
## 132 52.07565 51.32868 52.82262
## 133 46.17222 44.83550 47.50894
## 134 50.76378 49.98641 51.54114
## 135 52.07565 51.32868 52.82262
## 136 56.66720 55.47720 57.85721
## 137 47.48409 46.34963 48.61855
## 138 50.10784 49.28561 50.93007
## 139 49.45190 48.56861 50.33520
## 140 45.51628 44.07222 46.96034
## 141 49.45190 48.56861 50.33520
## 142 56.66720 55.47720 57.85721
## 143 51.41971 50.66809 52.17134
## 144 59.29095 57.67259 60.90931#Attention, un intervalle de confiance n’est pas un intervalle dans lequel la valeur a une #probabilité de 95 % de se trouver.
#L’IC a 95 % de chance de contenir la vraie valeur si on répète les estimations un grand nombre #de fois. Autrement dit,
#un intervalle de confiance à 95 % donnera un encadrement correct 95 fois sur 100.
#On peut aussi obtenir un intervalle de prédiction
newx=seq(1.5,10,1);
pp=predict(droite,data.frame(x= newx), level = 0.95, interval = "prediction")## Warning: 'newdata' had 9 rows but variables found have 144 rowsprint(pp)## fit lwr upr
## 1 46.17222 37.11988 55.22456
## 2 50.76378 41.77699 59.75056
## 3 53.38752 44.39869 62.37636
## 4 53.38752 44.39869 62.37636
## 5 54.69940 45.69873 63.70006
## 6 51.41971 42.43512 60.40430
## 7 50.76378 41.77699 59.75056
## 8 46.17222 37.11988 55.22456
## 9 53.38752 44.39869 62.37636
## 10 56.66720 47.63537 65.69904
## 11 53.38752 44.39869 62.37636
## 12 46.17222 37.11988 55.22456
## 13 54.69940 45.69873 63.70006
## 14 55.35533 46.34606 64.36460
## 15 50.76378 41.77699 59.75056
## 16 53.38752 44.39869 62.37636
## 17 45.51628 36.44747 54.58509
## 18 50.76378 41.77699 59.75056
## 19 51.41971 42.43512 60.40430
## 20 46.17222 37.11988 55.22456
## 21 58.63501 49.55604 67.71399
## 22 56.66720 47.63537 65.69904
## 23 53.38752 44.39869 62.37636
## 24 54.04346 45.04961 63.03731
## 25 56.66720 47.63537 65.69904
## 26 50.10784 41.11706 59.09861
## 27 54.69940 45.69873 63.70006
## 28 48.79597 39.79182 57.80012
## 29 46.17222 37.11988 55.22456
## 30 52.07565 43.09144 61.05985
## 31 61.25876 52.09255 70.42497
## 32 54.69940 45.69873 63.70006
## 33 57.32314 48.27736 66.36893
## 34 51.41971 42.43512 60.40430
## 35 50.76378 41.77699 59.75056
## 36 50.76378 41.77699 59.75056
## 37 55.35533 46.34606 64.36460
## 38 50.10784 41.11706 59.09861
## 39 53.38752 44.39869 62.37636
## 40 54.69940 45.69873 63.70006
## 41 50.76378 41.77699 59.75056
## 42 48.79597 39.79182 57.80012
## 43 57.32314 48.27736 66.36893
## 44 49.45190 40.45534 58.44847
## 45 46.17222 37.11988 55.22456
## 46 40.92473 31.69219 50.15726
## 47 48.79597 39.79182 57.80012
## 48 54.69940 45.69873 63.70006
## 49 51.41971 42.43512 60.40430
## 50 46.17222 37.11988 55.22456
## 51 50.76378 41.77699 59.75056
## 52 45.51628 36.44747 54.58509
## 53 53.38752 44.39869 62.37636
## 54 57.32314 48.27736 66.36893
## 55 53.38752 44.39869 62.37636
## 56 54.69940 45.69873 63.70006
## 57 50.76378 41.77699 59.75056
## 58 53.38752 44.39869 62.37636
## 59 46.17222 37.11988 55.22456
## 60 54.04346 45.04961 63.03731
## 61 50.76378 41.77699 59.75056
## 62 48.79597 39.79182 57.80012
## 63 58.63501 49.55604 67.71399
## 64 57.32314 48.27736 66.36893
## 65 58.63501 49.55604 67.71399
## 66 49.45190 40.45534 58.44847
## 67 53.38752 44.39869 62.37636
## 68 50.10784 41.11706 59.09861
## 69 54.69940 45.69873 63.70006
## 70 51.41971 42.43512 60.40430
## 71 50.10784 41.11706 59.09861
## 72 48.79597 39.79182 57.80012
## 73 56.66720 47.63537 65.69904
## 74 51.41971 42.43512 60.40430
## 75 55.35533 46.34606 64.36460
## 76 55.35533 46.34606 64.36460
## 77 52.07565 43.09144 61.05985
## 78 53.38752 44.39869 62.37636
## 79 52.07565 43.09144 61.05985
## 80 50.10784 41.11706 59.09861
## 81 46.82816 37.79053 55.86578
## 82 48.79597 39.79182 57.80012
## 83 48.79597 39.79182 57.80012
## 84 48.79597 39.79182 57.80012
## 85 56.01127 46.99161 65.03093
## 86 53.38752 44.39869 62.37636
## 87 48.79597 39.79182 57.80012
## 88 61.25876 52.09255 70.42497
## 89 51.41971 42.43512 60.40430
## 90 55.35533 46.34606 64.36460
## 91 49.45190 40.45534 58.44847
## 92 54.04346 45.04961 63.03731
## 93 54.04346 45.04961 63.03731
## 94 57.32314 48.27736 66.36893
## 95 54.04346 45.04961 63.03731
## 96 57.97908 48.91758 67.04058
## 97 48.79597 39.79182 57.80012
## 98 40.92473 31.69219 50.15726
## 99 52.07565 43.09144 61.05985
## 100 53.38752 44.39869 62.37636
## 101 53.38752 44.39869 62.37636
## 102 52.07565 43.09144 61.05985
## 103 55.35533 46.34606 64.36460
## 104 53.38752 44.39869 62.37636
## 105 50.10784 41.11706 59.09861
## 106 55.35533 46.34606 64.36460
## 107 48.79597 39.79182 57.80012
## 108 42.89254 33.74042 52.04465
## 109 53.38752 44.39869 62.37636
## 110 51.41971 42.43512 60.40430
## 111 54.04346 45.04961 63.03731
## 112 45.51628 36.44747 54.58509
## 113 52.07565 43.09144 61.05985
## 114 50.76378 41.77699 59.75056
## 115 52.73159 43.74597 61.71720
## 116 54.04346 45.04961 63.03731
## 117 52.07565 43.09144 61.05985
## 118 51.41971 42.43512 60.40430
## 119 55.35533 46.34606 64.36460
## 120 56.01127 46.99161 65.03093
## 121 55.35533 46.34606 64.36460
## 122 50.10784 41.11706 59.09861
## 123 42.23660 33.05936 51.41383
## 124 48.79597 39.79182 57.80012
## 125 52.07565 43.09144 61.05985
## 126 54.69940 45.69873 63.70006
## 127 48.79597 39.79182 57.80012
## 128 52.07565 43.09144 61.05985
## 129 53.38752 44.39869 62.37636
## 130 52.07565 43.09144 61.05985
## 131 45.51628 36.44747 54.58509
## 132 52.07565 43.09144 61.05985
## 133 46.17222 37.11988 55.22456
## 134 50.76378 41.77699 59.75056
## 135 52.07565 43.09144 61.05985
## 136 56.66720 47.63537 65.69904
## 137 47.48409 38.45941 56.50878
## 138 50.10784 41.11706 59.09861
## 139 49.45190 40.45534 58.44847
## 140 45.51628 36.44747 54.58509
## 141 49.45190 40.45534 58.44847
## 142 56.66720 47.63537 65.69904
## 143 51.41971 42.43512 60.40430
## 144 59.29095 50.19276 68.38914#Un intervalle de prédiction est un intervalle qui est susceptible de contenir une observation #individuelle future.
#L’intervalle de prédiction est toujours plus large que l’intervalle de confiance.#On peut afficher les prédictions, les bandes de confiance et de prédiction
#plot( pp[,1] ~ newx, type='p',pch=3,ylab="prédictions",xlab="x" )
#points( pc[,2] ~ newx, type='l', col="green" )
#points( pc[,3] ~ newx, type='l', col="green" )
#points( pp[,2] ~ newx, type='l', col="red" )
#points( pp[,3] ~ newx, type='l', col="red" )
#title(main="Bandes de confiance et de prédiction")
#legend("topleft", c("Bande de confiance", "Bande de prédiction"),lwd=1, lty=1, col=c("green", "red") )