############################################################################## # # Correction de la feuille de TD3 # ############################################################################### rm(list=ls()) # ---------- Exercice 1 : lien entre variables quantitatives ---------- # --- Partie 1 --- # 1. library(MASS) # 2. Animals help(Animals) # Animals contient pour plusieurs especes le poids moyen du corps en kilo et le poids moyen du cerveau en grammes. C'est un data.frame. # 3. plot(log(Animals$body), log(Animals$brain), main='Lien entre le poids du corps et celui du cerveau des animaux', xlab='log(Poids du corps)', ylab='log(Poids du cerveau)') # 4. On repere dans la table les indices des dinosaures : index= c(6,16,26) points(log(Animals$body)[index], log(Animals$brain)[index], col='red') legend('topleft', c('mammiferes','dinosaures'), pch=c(21,21), col=c('black','red')) # pch=21 permet d'afficher des cercles dans la legende. # Les dinosaures semblent correspondre a des points a part, les autres points donnes etant des mammiferes. # Il semble que le lien entre les log(corps) et log(cerveau) soit different selon que l'on est dinosaure ou mammifere. # 5. # On commence par exclure les dinosaures. X= log(Animals$body)[-index] Y= log(Animals$brain)[-index] plot(X, Y, main='Lien entre le poids du corps et celui du cerveau des mammiferes', xlab='log(Poids du corps)', ylab='log(Poids du cerveau)') # Il semble qu'il y a une relation lineaire du type y=a+bx. # a correspond a l'"intercept" et b a la pente le coefficient multiplicatif de x. # Mathematiquement cela veut dire que Y est proche de aX+b, c'est-a-dire Y=ax+b+bruit. # 6. reg= lm(Y~X) reg # Affiche les coeffcients a et b de la droite. # intercept est la valeur en 0 (c'est-a-dire b) et sous X il y a le coeffcient multiplicatif par rapport a X (c'est-a-dire a). # Pour verifier que cela correspond aux formules du cours, on recalcule ces coefficients à la main : n= length(X) hata= (sum(X*Y)-n*mean(X)*mean(Y))/(sum(X^2)-n*(mean(X))^2) hatb= mean(Y)-hata*mean(X) # 7. coefficients(reg) # Vecteur des coefficients. fitted.values(reg) # Vecteur des hatf(X), cad des predictions obtenues avec la droite. residuals(reg) # Vecteur des Y-hatf(X), cad des differences entre les observations et les predictions. # 8. points(X, fitted.values(reg), pch=4, col='red') legend('topleft', c('observation','estimation'), pch=c(21,4), col=c('black','red')) # NB : points permet de rajouter des points sur un graphique sans avoir a se soucier des xlim/ylim. # On garde l'echelle du premier plot quoiqu'il arrive. lines permet de rajouter des lignes. # Ici on aurait pu rajouter une ligne avec lines mais des que hatf est plus complique qu'une simple droite, le graphe va devenir affreux, points est la solution qui sera comprehensible quelquesoit le modele. # 9. # Calcul du coefficient de détermination (rapport entre variance expliquee et variance totale) rdeux= sum((fitted.values(reg)-mean(Y))^2)/sum((Y-mean(Y))^2) # Il est proche de 1, le modele est de bonne qualite. # Calcul du coefficient de correlation (covariance renormalisee par les ecarts-type) r= sum((X-mean(X))*(Y-mean(Y)))/sqrt(sum((X-mean(X))^2)*sum((Y-mean(Y))^2)) # Lui aussi est tres proche de 1 et on verifie que r^2=rdeux (attention formule vraie uniquement parce que c'est une regression lineaire !) # 10. par(mfrow=c(1,2)) hist(residuals(reg), freq=FALSE, main='histogramme des residus', xlab='residus', ylab='densite') boxplot(residuals(reg), main='boite a moustache des residus') # Sur l'histogramme on voit peu de choses car on a peu de donnees (25). # Si il y avait plus d'observations, on pourrait alors superposer la densite gaussienne correspondante et voir alors si les residus sont bien centres, symmetriques autour de 0 et de distribution en gros gaussienne. # Si on n'a pas assez d'observations, on peut au moins voir le caractere centre et symmetrique sur la boite a moustache. # Sur la boite a moustache on voit ici que les residus sont bien en gros centres autour de 0 et symmetriques. # Barre du milieu= mediane, fin du rectangle 1er et 3eme quartile, barre = 1*5 fois l'interquartile (representatif de l'etendue classique d'une gaussienne de son premier a son 99eme percentile). # Au dela= outliers - generalement des donnees "exceptionnelles" a eventuellement identifier plus avant. # Ici avec juste 25 points, les outliers ont un sens tres peu fiables, par contre on voit que le dernier pic de l'histogramme etait du en fait a seulement deux points, ces deux "outliers" qui rendait le graphe totalement asymmetrique alors que ce n'est pas vraiment le cas. # ------- Partie 2 ------- # 11. # On commence par creer les donnees. X= c(115.22, 135.98, 119.34, 114.96, 187.05, 243.92, 267.43, 238.71, 295.94, 317.78, 216, 240.35, 386.57, 261.53, 249.34, 309.87, 345.89, 165.54, 196.98, 395.26) Y= c(369, 439, 475, 603, 1247, 1298, 1420, 1476, 1532, 1639, 1735, 1777, 1793, 1843, 1855, 1880, 1881, 1904, 1922, 1993) # On reprend le modele precedent du type y=a+bx. # 12. plot(X, Y, xlab='Depense en nourriture (euros)', ylab='Revenu attribue (euros)', main='Etude du lien entre depense en nourriture et revenu attribue') # 13. summary(X) summary(Y) # Ces deux variables etant quantitatives continues, on etudie leur distribution en tracant les histogrammes de leurs frequences. par(mfrow=c(1,2)) hist(X, freq=FALSE, main='Histogramme des depenses en nourriture', xlab='Depense en nourriture (euros)') hist(Y, freq=FALSE, main='Histogramme des revenus attribues', xlab='Revenu attribue (euros)') # 14. r= sum((X-mean(X))*(Y-mean(Y)))/sqrt(sum((X-mean(X))^2)*sum((Y-mean(Y))^2)) # Le coefficient de correlation est ici moins proche de 1 que dans la 1ere partie. # On doit s'attendre a une regression lineaire de moins bonne qualite. # 15. reg2= lm(Y~X) reg2 # D'apres ce modele, le lien entre X et Y serait donc y = 328.626 + 4.686x plot(X, Y, xlab='Depense en nourriture (euros)', ylab='Revenu attribue (euros)', main='Etude du lien entre depense en nourriture et revenu attribue') points(X, fitted.values(reg2), pch=4, col='red') legend('topleft', c('observation','estimation'), pch=c(21,4), col=c('black','red')) # La droite suit grossierement la tendance des donnees, mais on constate que certaines predictions sont assez eloignees des observations. # Etudes des residus. par(mfrow=c(1,2)) hist(residuals(reg2), freq=FALSE, main='histogramme des residus', xlab='residus', ylab='densite') boxplot(residuals(reg2), main='boite a moustache des residus') # Les residus sont grossierement centrees autour de 0, de facon symetrique. # En revanche, on constate egalement un nombre important de residus eloignes de 0, ce qui indique la presence d'erreurs importantes. # ---------- Exercice 2 : lien entre variable qualitative et quantitative ---------- library(MASS) # 1. donnees1= cabbages par(mfrow=c(1,2)) boxplot(donnees1$HeadWt~donnees1$Cult, xlab='Variete de chou', ylab='Poids', main='Poids vs Variete du chou') boxplot(donnees1$VitC~donnees1$Cult, xlab='Variete de chou', ylab='Teneur de vitamine C', main='Teneur en vitamine C vs Variete du chou') # La variete c39 propose des choux plus lourds, mais contenant moins de vitamine que la variete c52. # 2. donnees2= anorexia # On fait la difference entre le poids apres et avant traitement. boxplot(donnees2$Postwt-donnees2$Prewt~donnees2$Treat, xlab='Traitement', ylab='Evolution du poids', main='Evolution du poids selon le traitement utilise') # FT semble donner en moyenne les meilleurs resultats. # 3. donnees3= Cars93 par(oma=c(0,0,3,0)) # Permet de regler les marges de la fenetre graphique. Sans cette option, le titre n'apparait qu'a moitie. par(mfrow=c(1,2)) boxplot(donnees3$Min.Price~donnees3$Origin, xlab='Origine', ylab='Prix minimum') boxplot(donnees3$Max.Price~donnees3$Origin, xlab='Origine', ylab='Prix maximum') title("Liens entre prix et origine du vehicule", outer=TRUE) # Rajoute un titre general a la fenetre # Les voitures etrangeres ne sont pas en moyenne plus abordables que les voitures americaines. # Le marche etranger propose une plus grande gamme de prix (plus grandes etendues des boites a moustache). # 4. # On multiplie la capacite du reservoir et la consommation du vehicule sur voie rapide pour obtenir son autonomie. boxplot(donnees3$Fuel.tank.capacity*donnees3$MPG.highway~donnees3$Type, xlab='Type de voiture', ylab='Autonomie sur voie rapide', main="Comparaison de l'autonomie des vehicules sur voie rapide") # Les petites voitures ont en moyenne la plus faible autonomie.