Cours CMI / M1 Statistiques - séance 2   Résumé de la séance 2 Dans cette séance, on introduit le vocabulaire classique de la régression et on explique comment on résoud mathématiquement le problème linéaire à l'aide d'un exemple simple. On y présente aussi tout ce qui est lié à la régression : le choix des régresseurs, l'analyse de la qualité de la régression... Des exemples sur données réelles complètent cette séance.   Liens utiles Sur la notion de régression en général, on peut lire notre page EDA2 CRS1 et le document siab5 dont une copie est ici. Comme premier exemple de régression linéaire simple, on utilisera EDA2 CRS2 et le document modStat_C1 dont une copie est ici. Pour la régression linéaire multiple, on consultera la fin du document siab5 (à partir de la page 51), notre cours EDA2 CRS3, et le document modStat_C2 dont une copie est ici. Exercices de régression linéaire simple Si ml est le résultat d'une régression linéaire simple, quelle est la différence entre les trois instructions ci-dessous : # remarque : ml est ici issu d'un calcul de régression linéaire, # comme par exemple ml <- lm( y ~ . , data=...) print(ml) print(anova(ml)) print(summary(ml)) Effectuer, par copier/coller à partir de EDA2 CRS2, tous les calculs et graphiques pour la modélisation linéaire essence en fonction de distance. Faire une fonction nommée reglin() qui permet d'automatiser les traitements usuels d'une régression linéaire. Essayer de modéliser linéairement Sepal.Length à l'aide la variable Sepal.Width pour les données iris. Est-ce un "bon" modèle  ? Est-ce graphiquement satisfaisant ? Essayer de modéliser linéairement Sepal.Length à l'aide la variable Petal.Width pour les données iris. Est-ce un "bon" modèle  ? Est-ce graphiquement satisfaisant ? Pour les données LEA, essayer de modéliser linéairement MW à l'aide la variable LENGTH. Est-ce un "bon" modèle  ? Y a-t-il causalité ? Est-ce graphiquement satisfaisant ? Peut-on fournir une explication biologique aux coefficients ? Ou était le piège ? Toujours pour ces données LEA, avec la modélisation linéaire de MW à l'aide la variable LENGTH, quels sont les points-levier les plus importants ? Et les points influents les plus importants ?  Solutions  print(ml) n'affiche que les coefficients du modèle. print(anova(ml)) affiche uniquement la décomposition de la variance, mais avec une p-value qui indique si le modèle linéaire est approprié. print(summary(ml)) fournit tous les résultats utiles de la modélisation et en particulier le résultat du test de Fisher pour savoir si le modèle linéaire est approprié. Voici le code de la fonction demandée et son utilisation, fichier cmi02reglin.r : ## automatisation des traitements usuels d'une régression linéaire simple ## on utilise des fonctions de statgh.r pour les titres if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } ############################################################################# reglin <- function(x,y,titre) { ############################################################################# # il serait bon de rajouter les unités, # de commenter le résultat du test de Fisher et de l'anova catss(titre) cats("Résumé des données") print(summary(data.frame(x,y))) # calcul du modèle linéaire ml <- lm( y ~ x ) cats("Coefficients") print(ml) cats("Modèle complet") print( summary(ml) ) cats("ANOVA du modèle") print( anova(ml) ) idx <- order(x) plot(y[idx]~x[idx],type="b",main=titre) abline( ml, col="red",lwd=2) } # fin fonction reglin ## exemples d'utilisation : # km km <- lit.dar("http://forge.info.univ-angers.fr/~gh/wstat/Eda/km.dar") reglin(km$distance,km$essence,"Essence en fonction de distance") # iris (1) data(iris) reglin(iris$Sepal.Width,iris$Sepal.Length ,"Length en fonction Width pour les sépales d'iris") # iris (2) reglin(iris$Petal.Width,iris$Sepal.Length ,"Sepal.Length en fonction de Petal.Width pour les données d'iris") Il est possible de la télécharger sur internet, via l'instruction : source("http://forge.info.univ-angers.fr/~gh/Cmi/reglin.r",encoding="latin1") A l'aide de la commande reglin(iris$Sepal.Width,iris$Sepal.Length ,"Length en fonction Width pour les sépales d'iris") on voit que le modèle linéaire n'est pas approprié ici (p-value de Fisher > 0,15). + ================================================ + + + + Length en fonction Width pour les sépales d'iris + + + + ================================================ + Résumé des données ================== x y Min. :2.000 Min. :4.300 1st Qu.:2.800 1st Qu.:5.100 Median :3.000 Median :5.800 Mean :3.057 Mean :5.843 3rd Qu.:3.300 3rd Qu.:6.400 Max. :4.400 Max. :7.900 Coefficients ============ Call: lm(formula = y ~ x) Coefficients: (Intercept) x 6.5262 -0.2234 Modèle complet ============== Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -1.5561 -0.6333 -0.1120 0.5579 2.2226 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5262 0.4789 13.63 <2e-16 *** x -0.2234 0.1551 -1.44 0.152 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.8251 on 148 degrees of freedom Multiple R-squared: 0.01382, Adjusted R-squared: 0.007159 F-statistic: 2.074 on 1 and 148 DF, p-value: 0.1519 ANOVA du modèle =============== Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x 1 1.412 1.41224 2.0744 0.1519 Residuals 148 100.756 0.68078 On utilise ici la commande reglin(iris$Petal.Width,iris$Sepal.Length ,"Sepal.Length en fonction de Petal.Width pour les données d'iris"). Le modèle linéaire est approprié (p-value de Fisher < 2.2e-16), la régression est de moyenne qualité (R2a=0.6668) et le résultat graphique moyennement satisfaisant. + =============================================================== + + + + Sepal.Length en fonction de Petal.Width pour les données d'iris + + + + =============================================================== + Résumé des données ================== x y Min. :0.100 Min. :4.300 1st Qu.:0.300 1st Qu.:5.100 Median :1.300 Median :5.800 Mean :1.199 Mean :5.843 3rd Qu.:1.800 3rd Qu.:6.400 Max. :2.500 Max. :7.900 Coefficients ============ Call: lm(formula = y ~ x) Coefficients: (Intercept) x 4.7776 0.8886 Modèle complet ============== Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -1.38822 -0.29358 -0.04393 0.26429 1.34521 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.77763 0.07293 65.51 <2e-16 *** x 0.88858 0.05137 17.30 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.478 on 148 degrees of freedom Multiple R-squared: 0.669, Adjusted R-squared: 0.6668 F-statistic: 299.2 on 1 and 148 DF, p-value: < 2.2e-16 ANOVA du modèle =============== Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x 1 68.353 68.353 299.17 < 2.2e-16 *** Residuals 148 33.815 0.228 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Il est prudent de lire toutes les sorties des fonctions ! Nous avons pris la précaution, dans la fonction reglin(), d'afficher un résumé des données à traiter. Une lecture attentive de ces résumés montre qu'il y a des poids moléculaires égaux à -9, ce qui n'est bien sûr pas possible. Il s'agit en fait de données incomplètes : la récupération des données a été faite sur Internet, et les données manquantes sont repérées, pour les variables quantitatives, par -9. Il faut donc retirer ces données, heureusement en faible quantité (il y en a 5) avant de continuer. Voici le code associé : ## régression linéaire simple de MW à l'aide de LENGTH dans les données LEA # lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } # lecture des données lea <- lit.dar("http://forge.info.univ-angers.fr/~gh/Datasets/lea.dar") cats("Données de départ") cols <- c("length","mw") leaLM <- lea[ , cols ] print(summary(leaLM)) # suppression des données incorrectes (après une lecture rapide de summary) filtre <- (leaLM[ ,"mw"] > 0 ) cat("\n il y a ",sum(!filtre)," données incorrectes (MW<0) \n") # affichage cats("Données correctes") leaLM <- leaLM[ filtre, ] print(summary(leaLM)) # sauvegarde save(leaLM,file="leaLM.rdata") # sauvegarde dans un format adapté à R # analyse via la fonction de l'exercice précédent source("reglin.r",encoding="latin1") # lecture du code de la fonction reglin() reglin(leaLM$length,leaLM$mw,"MW en fonction de LENGTH, données LEA") Et ses résultats : Données de départ ================= length mw Min. : 68.0 Min. : -9 1st Qu.: 130.0 1st Qu.: 13961 Median : 168.0 Median : 17447 Mean : 205.7 Mean : 21839 3rd Qu.: 236.0 3rd Qu.: 25100 Max. :1864.0 Max. :197129 il y a 5 données incorrectes (MW<0) Données correctes ================= length mw Min. : 68.0 Min. : 7521 1st Qu.: 130.0 1st Qu.: 14011 Median : 168.0 Median : 17473 Mean : 205.9 Mean : 21981 3rd Qu.: 235.2 3rd Qu.: 25149 Max. :1864.0 Max. :197129 + ===================================== + + + + MW en fonction de LENGTH, données LEA + + + + ===================================== + Résumé des données ================== x y Min. : 68.0 Min. : 7521 1st Qu.: 130.0 1st Qu.: 14011 Median : 168.0 Median : 17473 Mean : 205.9 Mean : 21981 3rd Qu.: 235.2 3rd Qu.: 25149 Max. :1864.0 Max. :197129 Coefficients ============ Call: lm(formula = y ~ x) Coefficients: (Intercept) x -71.01 107.12 ANOVA du modèle =============== Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x 1 1.9554e+11 1.9554e+11 152317 < 2.2e-16 *** Residuals 766 9.8338e+08 1.2838e+06 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Modèle complet ============== Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -6359.3 -590.0 77.5 518.2 6203.6 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -71.0072 69.7444 -1.018 0.309 x 107.1227 0.2745 390.278 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1133 on 766 degrees of freedom Multiple R-squared: 0.995, Adjusted R-squared: 0.995 F-statistic: 1.523e+05 on 1 and 766 DF, p-value: < 2.2e-16 Pour disposer des indicateurs de qualité comme les valeurs des leviers, des points influents, il faudrait passer par plusieurs fonctions, dont hatvalues(), dfits()... Nous fournissons, avec statgh.r une fonction analyseRegression() qui regroupe tous ces appels de fonctions : Code associé : ## compléments sur la régression linéaire simple de MW à l'aide de LENGTH dans les données LEA # lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } # reprise des données load("leaLM.rdata") # df = leaLM cats("Rappel des données") print(summary(leaLM)) analyseRegression("MW en fonction de LENGTH",xi = leaLM$length,yi=leaLM$mw,interactif = FALSE) Résultats : Rappel des données de départ ============================ length mw Min. : 68.0 Min. : 7521 1st Qu.: 130.0 1st Qu.: 14011 Median : 168.0 Median : 17473 Mean : 205.9 Mean : 21981 3rd Qu.: 235.2 3rd Qu.: 25149 Max. :1864.0 Max. :197129 Analyse de la régression "MW en fonction de LENGTH" sur 768 valeurs =================================================================== Affichage de summary(lm(yi~xi) : -------------------------------- Call: lm(formula = yi ~ xi) Residuals: Min 1Q Median 3Q Max -6359.3 -590.0 77.5 518.2 6203.6 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -71.0072 69.7444 -1.018 0.309 xi 107.1227 0.2745 390.278 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1133 on 766 degrees of freedom Multiple R-squared: 0.995, Adjusted R-squared: 0.995 F-statistic: 1.523e+05 on 1 and 766 DF, p-value: < 2.2e-16 Résultats de la régression : ---------------------------- valeur de R : rho = 0.997 donc R2= 0.995 ; p-value : 0.000 *** coefficients de la régression et intervalles de confiance à 95 % Estimate Std. Error t value Pr(>|t|) 2.5 % 97.5 % (Intercept) -71.007 69.744 -1.018 0.309 -207.920 65.906 xi 107.123 0.274 390.278 0.000 106.584 107.662 Résidus et points remarquables ------------------------------ Début de ces résultats, 10 lignes au maximum seront affichées obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 1 200 20295.80 21353.53 -1057.73103 -0.934142499 0.001304098 -0.93406474 -0.0337532407 5.697355e-04 2 153 15909.32 16318.77 -409.44485 -0.361633387 0.001466055 -0.36142811 -0.0138489104 9.600512e-05 3 132 14103.76 14069.19 34.57125 0.030536701 0.001622220 0.03051678 0.0012301162 7.575810e-07 4 143 15209.61 15247.54 -37.92780 -0.033500056 0.001533964 -0.03347821 -0.0013122092 8.620707e-07 5 144 15422.86 15354.66 68.19950 0.060237577 0.001526645 0.06019839 0.0023538870 2.774000e-06 6 143 15237.67 15247.54 -9.86830 -0.008716261 0.001533964 -0.00871057 -0.0003414189 5.835960e-08 768 résidu(s) fort(s) >= seuil = 0 on n'affiche que les 10 premiers résidus obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 251 507 47880.90 54240.20 -6359.301 -5.631 0.007 -5.748 -0.469 0.106 652 502 47658.48 53704.59 -6046.105 -5.353 0.006 -5.453 -0.439 0.093 661 486 46163.72 51990.62 -5826.904 -5.158 0.006 -5.247 -0.404 0.079 328 647 64202.23 69237.38 -5035.146 -4.472 0.013 -4.529 -0.514 0.129 327 618 61093.58 66130.82 -5037.240 -4.471 0.011 -4.528 -0.483 0.114 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 172 535 61524.10 57239.64 4284.467 3.796 0.008 3.830 0.336 0.056 667 385 37217.05 41171.23 -3954.180 -3.495 0.003 -3.521 -0.199 0.020 51 levier(s) important(s) >= seuil 3/n = 0.00390625 on n'affiche que les 10 premiers leviers obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 759 1864 197129.31 199605.70 -2476.384 -2.388 0.163 -2.396 -1.056 0.554 122 1509 159774.88 161577.14 -1802.265 -1.678 0.101 -1.680 -0.563 0.158 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 745 1236 132565.45 132332.64 232.809 0.212 0.064 0.212 0.055 0.002 397 847 94075.55 90661.92 3413.631 3.052 0.025 3.069 0.496 0.121 76 742 83297.49 79414.03 3883.460 3.459 0.018 3.484 0.474 0.111 722 733 77051.16 78449.93 -1398.772 -1.246 0.018 -1.246 -0.167 0.014 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 744 684 75708.18 73200.92 2507.263 2.229 0.015 2.235 0.273 0.037 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 autres seuils possibles : 2/n = 0.003 Hoaglin et Welsch (1978 ) ; 0.5 Huber (1987) 39 résidu(s) studentisé(s) fort(s) >= seuil = 2 on n'affiche que les 10 premiers résidus studentisés obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 251 507 47880.90 54240.20 -6359.301 -5.631 0.007 -5.748 -0.469 0.106 652 502 47658.48 53704.59 -6046.105 -5.353 0.006 -5.453 -0.439 0.093 661 486 46163.72 51990.62 -5826.904 -5.158 0.006 -5.247 -0.404 0.079 328 647 64202.23 69237.38 -5035.146 -4.472 0.013 -4.529 -0.514 0.129 327 618 61093.58 66130.82 -5037.240 -4.471 0.011 -4.528 -0.483 0.114 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 172 535 61524.10 57239.64 4284.467 3.796 0.008 3.830 0.336 0.056 667 385 37217.05 41171.23 -3954.180 -3.495 0.003 -3.521 -0.199 0.020 39 observation(s) mal reconstituée(s), rstudent >= seuil = 1.96307 on n'affiche que les 10 premières observation(s) mal reconstituée(s) obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 251 507 47880.90 54240.20 -6359.301 -5.631 0.007 -5.748 -0.469 0.106 652 502 47658.48 53704.59 -6046.105 -5.353 0.006 -5.453 -0.439 0.093 661 486 46163.72 51990.62 -5826.904 -5.158 0.006 -5.247 -0.404 0.079 328 647 64202.23 69237.38 -5035.146 -4.472 0.013 -4.529 -0.514 0.129 327 618 61093.58 66130.82 -5037.240 -4.471 0.011 -4.528 -0.483 0.114 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 172 535 61524.10 57239.64 4284.467 3.796 0.008 3.830 0.336 0.056 667 385 37217.05 41171.23 -3954.180 -3.495 0.003 -3.521 -0.199 0.020 38 dfits significatif(s) >= seuil 2*racine(2/n) = 0.1020621 on n'affiche que les 10 premiers dfits significatifs obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 759 1864 197129.31 199605.70 -2476.384 -2.388 0.163 -2.396 -1.056 0.554 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 122 1509 159774.88 161577.14 -1802.265 -1.678 0.101 -1.680 -0.563 0.158 328 647 64202.23 69237.38 -5035.146 -4.472 0.013 -4.529 -0.514 0.129 397 847 94075.55 90661.92 3413.631 3.052 0.025 3.069 0.496 0.121 327 618 61093.58 66130.82 -5037.240 -4.471 0.011 -4.528 -0.483 0.114 76 742 83297.49 79414.03 3883.460 3.459 0.018 3.484 0.474 0.111 251 507 47880.90 54240.20 -6359.301 -5.631 0.007 -5.748 -0.469 0.106 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 38 distance(s) de cook >= seuil = 0.005208333 on n'affiche que les 10 premières distance(s) de cook >= seuil obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 487 1429 159210.94 153007.32 6203.613 5.737 0.089 5.860 1.833 1.610 759 1864 197129.31 199605.70 -2476.384 -2.388 0.163 -2.396 -1.056 0.554 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 122 1509 159774.88 161577.14 -1802.265 -1.678 0.101 -1.680 -0.563 0.158 328 647 64202.23 69237.38 -5035.146 -4.472 0.013 -4.529 -0.514 0.129 397 847 94075.55 90661.92 3413.631 3.052 0.025 3.069 0.496 0.121 327 618 61093.58 66130.82 -5037.240 -4.471 0.011 -4.528 -0.483 0.114 76 742 83297.49 79414.03 3883.460 3.459 0.018 3.484 0.474 0.111 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 251 507 47880.90 54240.20 -6359.301 -5.631 0.007 -5.748 -0.469 0.106 autres seuils possibles : souhaitable f(1,n-1,0.1) = 0.016 ; précoccupant f(1,n-1,0.5) = 0.455 Cook (1977) 46 point(s) remarquable(s) vus par influence.measures (sans tri) on n'affiche que les 10 premiers points remarquables (sans tri) obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 16 484 54509.96 51776.38 2733.584 2.420 0.006 2.427 0.186 0.017 76 742 83297.49 79414.03 3883.460 3.459 0.018 3.484 0.474 0.111 80 643 67887.34 68808.89 -921.542 -0.818 0.013 -0.818 -0.092 0.004 118 417 47919.19 44599.16 3320.035 2.936 0.004 2.951 0.185 0.017 122 1509 159774.88 161577.14 -1802.265 -1.678 0.101 -1.680 -0.563 0.158 127 723 82204.35 77378.70 4825.650 4.296 0.017 4.346 0.571 0.160 128 678 76856.36 72558.18 4298.179 3.821 0.014 3.855 0.466 0.107 172 535 61524.10 57239.64 4284.467 3.796 0.008 3.830 0.336 0.056 201 391 38903.89 41813.97 -2910.076 -2.573 0.003 -2.582 -0.149 0.011 209 575 58560.91 61524.54 -2963.632 -2.628 0.009 -2.638 -0.256 0.032 DESCRIPTION STATISTIQUE DE LA VARIABLE Résidus Taille 768 individus Moyenne 0.0000 Ecart-type 1132.3018 Coef. de variation -1 % 1er Quartile -590.0171 Mediane 77.5062 3eme Quartile 518.1828 iqr absolu 1108.1999 iqr relatif 1430.0000 % Minimum -6359.301 Maximum 6203.613 Tracé tige et feuilles The decimal point is 3 digit(s) to the right of the | -6 | 40 -5 | 800 -4 | 0 -3 | 64443210000 -2 | 9855440 -1 | 887766655555544444444444444333333333333222222222222111111111111111111111000000000000 -0 | 9999999999999999988888888888888888888888888888888888888888887777777777777777766666666666+132 0 | 0000000000000000000000000000111111111111111111111111111111111111111111111122222222222222+227 1 | 0000000000011111122222223333333344444444444444455555555555566666667777778888899 2 | 001111345777 3 | 013449 4 | 338 5 | 6 | 2 Test de normalité de Shapiro-Wilk Shapiro-Wilk normality test data: ei W = 0.9079, p-value < 2.2e-16 Test de l'homoscédasticité des résidus ====================================== studentized Breusch-Pagan test data: lm(yi ~ xi) BP = 209.14, df = 1, p-value < 2.2e-16 $modele Call: lm(formula = yi ~ xi) Coefficients: (Intercept) xi -71.01 107.12 $coefs Estimate Std. Error t value Pr(>|t|) 2.5 % 97.5 % (Intercept) -71.00722 69.744370 -1.018107 0.3089486 -207.9200 65.90556 xi 107.12270 0.274478 390.277892 0.0000000 106.5839 107.66151 $indicateurs obs xi yi yi_chapeau résidu résidu_studentisé levier rstudent dfits cook_dist 1 200 20295.80 21353.532 -1057.731031 -0.934142499 0.001304098 -0.934064738 -0.0337532407 5.697355e-04 2 153 15909.32 16318.765 -409.444853 -0.361633387 0.001466055 -0.361428111 -0.0138489104 9.600512e-05 3 132 14103.76 14069.189 34.571248 0.030536701 0.001622220 0.030516781 0.0012301162 7.575810e-07 4 143 15209.61 15247.538 -37.927800 -0.033500056 0.001533964 -0.033478206 -0.0013122092 8.620707e-07 [...] Exercices de régression linéaire multiple Comment calculer et analyser la matrice des coefficients de corrélation linéaire pour les données iris ? Comment utiliser la fonction pairsi() de statgh.r au lieu de la simple fonction pairs() de R pour les données iris ? Ecrire une fonction nommée mrls() qui, étant donné un data.frame qui contient en colonne 1 une variable à expliquer et ensuite les variables explicatives, fournit un tableau résumé de toutes les régressions linéaires simples associées. Quel est le meilleur critère de tri pour les lignes du tableau ? Pour mettre au point la fonction, on utilisera les données chenilles de EDA2CR2 avec LN_NID comme cible car on en connait les résultats. On retirera la variable NIDS. On pourra consulter, si besoin est, la page Linear Regression with R. Quelle est la meilleure régression linéaire simple possible pour expliquer la variable PI parmi toutes les régressions linéaires simples possibles dans les données LEA ? On pourra utiliser la fonction mrls() de la question précédente. On commencera par ne retenir que les données "propres" pour les variables quantitatives. Effectuer ensuite une régression linéaire multiple de PI en fonction de toutes les variables. Que faut-il en conclure ? Les résultats sont-ils meilleurs si on se restreint aux protéines dont la longueur est inférieure à 900 acides aminés ? Quelle est la meilleure régression linéaire multiple pour expliquer la variable MW parmi toutes les régressions linéaires simples possibles dans les données LEA ? Comment tester si le modèle qui utilise pi length foldindex gravy est significativement meilleur que celui qui utilise length foldindex gravy seulement ?  Solutions  On peut utiliser la fonction cor() de R pour la calculer, mais il n'y a aucun outil ni aucune fonction pour l'analyser. La fonction mdc() de statgh.r permet de le faire. Avant d'utiliser ces fonctions, il faut retirer la colonne 5 qui n'est pas quantitative. Code associé : ## matrice des coefficients de corrélation linéaire pour les données iris # lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } # préparation des données, on retire la colonne 5 data(iris) iris4 <- iris[, -5 ] # solution R classique print(round(cor(iris4),3)) # solution gh mdc(iris4,meilcor=TRUE) Résultats : Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1.000 -0.118 0.872 0.818 Sepal.Width -0.118 1.000 -0.428 -0.366 Petal.Length 0.872 -0.428 1.000 0.963 Petal.Width 0.818 -0.366 0.963 1.000 Matrice des corrélations au sens de Pearson pour 150 lignes et 4 colonnes Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1.000 Sepal.Width -0.118 1.000 Petal.Length 0.872 -0.428 1.000 Petal.Width 0.818 -0.366 0.963 1.000 Meilleure corrélation 0.9628654 pour Petal.Width et Petal.Length p-value 4.68e-86 Formules Petal.Length = 2.230 * Petal.Width + 1.084 et Petal.Width = 0.416 * Petal.Length -0.363 Coefficients de corrélation par ordre décroissant 0.963 p-value 0.0000 pour Petal.Width et Petal.Length : Petal.Width = 0.416 x Petal.Length - 0.363 0.872 p-value 0.0000 pour Petal.Length et Sepal.Length : Petal.Length = 1.858 x Sepal.Length - 7.101 0.818 p-value 0.0000 pour Petal.Width et Sepal.Length : Petal.Width = 0.753 x Sepal.Length - 3.200 -0.428 p-value 0.0000 pour Petal.Length et Sepal.Width : Petal.Length = -1.735 x Sepal.Width + 9.063 -0.366 p-value 0.0000 pour Petal.Width et Sepal.Width : Petal.Width = -0.640 x Sepal.Width + 3.157 -0.118 p-value 0.1519 pour Sepal.Width et Sepal.Length : Sepal.Width = -0.062 x Sepal.Length + 3.419 Il faut utiliser les options de pairsi() pour avoir soit les valeurs de corrélation linéaire, soit la couleur. ## tracé en paires pour les données iris # lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } # préparation des données, on retire la colonne 5 data(iris) iris4 <- iris[, -5 ] # solution R classique pairs(iris4) # trois extensions gh pairsi(iris4) # équivalent à pairsi(iris4,opt=2) couleursIris <- c("red","green","blue")[ as.numeric(iris$Species) ] pairsi(iris4,opt=1,col=couleursIris,pch=19) print(table(iris$Species,couleursIris)) blue green red setosa 0 0 50 versicolor 0 50 0 virginica 50 0 0 Une seule boucle R suffit pour passer en revue toutes les régressions linéaires simples : ## tableau résumé de toutes les régressions linéaires simples ## pour un dataframe de variables quantitatives dont la colonne 1 est la cible # lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } ####################################################################### mrls <- function(df) { ####################################################################### nbvar <- ncol(df) # rappel rapide des données print(summary(df)) mdc(df) # préparation du tableau des résultats cols <- c("pvalue","R2a","cor","AIC","BIC") nbc <- length(cols) tabRes <- data.frame(matrix(NA,nrow=nbvar-1,ncol=nbc)) names(tabRes) <- cols row.names(tabRes) <- names(df)[2:nbvar] # remplissage for (idv in (2:nbvar)) { ml <- lm( df[,1] ~ df[ ,idv] ) res <- summary(ml) fstat <- res$fstatistic tabRes$pvalue[ idv-1 ] <- pf(q=fstat[1], df1=fstat[2], df2=fstat[3], lower.tail=FALSE) tabRes$R2a[ idv-1 ] <- round(abs(res$adj.r.squared) , 4) tabRes$cor[ idv-1 ] <- round(cor(df[,1],df[,idv]),4) tabRes$AIC[ idv-1 ] <- round(AIC(ml),4) tabRes$BIC[ idv-1 ] <- round(BIC(ml),4) } # fin pour idv # affichage cat("\nTableau de toutes les régressions linéaires simples\n\n") idx <- order(tabRes$R2a) print(tabRes[idx,]) cat("\n") } # fin de fonction mrls ####################################################################### ## exemples d'applications # iris data(iris) mrls( iris[,-5] ) # chenilles url <- "http://forge.info.univ-angers.fr/~gh/wstat/Eda/chenilles.dar" chen <- lit.dar(url) chen$NIDS <- NULL nbc <- ncol(chen) chen <- chen[ , names(chen)[c(nbc,2:(nbc-1)) ] ] mrls( chen ) Sepal.Length Sepal.Width Petal.Length Petal.Width Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 Median :5.800 Median :3.000 Median :4.350 Median :1.300 Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800 Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500 Matrice des corrélations au sens de Pearson pour 150 lignes et 4 colonnes Sepal.Length Sepal.Width Petal.Length Petal.Width Sepal.Length 1.000 Sepal.Width -0.118 1.000 Petal.Length 0.872 -0.428 1.000 Petal.Width 0.818 -0.366 0.963 1.000 Tableau de toutes les régressions linéaires simples pvalue R2a cor AIC BIC Sepal.Width 1.518983e-01 0.0072 -0.1176 371.9917 381.0236 Petal.Width 2.325498e-37 0.6668 0.8179 208.2215 217.2534 Petal.Length 1.038667e-47 0.7583 0.8718 160.0404 169.0723 LN_NID PENTE PINS HAUT DIAM DENS ORIEN HAUTDOM STRAT MELAN Min. :-3.5066 Min. :15.00 Min. : 0.00 Min. :2.400 Min. : 5.80 Min. :1.000 Min. :1.100 Min. : 3.600 Min. :1.100 Min. :1.300 1st Qu.:-1.7148 1st Qu.:24.00 1st Qu.: 4.00 1st Qu.:3.700 1st Qu.:11.50 1st Qu.:1.200 1st Qu.:1.600 1st Qu.: 5.900 1st Qu.:1.500 1st Qu.:1.600 Median :-0.4005 Median :28.00 Median : 8.00 Median :4.400 Median :15.70 Median :1.500 Median :1.700 Median : 7.200 Median :2.000 Median :1.800 Mean :-0.8133 Mean :29.03 Mean :11.45 Mean :4.452 Mean :15.25 Mean :1.791 Mean :1.658 Mean : 7.539 Mean :1.982 Mean :1.752 3rd Qu.: 0.1222 3rd Qu.:34.00 3rd Qu.:18.00 3rd Qu.:5.300 3rd Qu.:18.30 3rd Qu.:2.400 3rd Qu.:1.800 3rd Qu.: 9.100 3rd Qu.:2.500 3rd Qu.:2.000 Max. : 1.0986 Max. :46.00 Max. :32.00 Max. :6.500 Max. :21.80 Max. :3.300 Max. :1.900 Max. :13.700 Max. :2.900 Max. :2.000 Matrice des corrélations au sens de Pearson pour 33 lignes et 10 colonnes LN_NID PENTE PINS HAUT DIAM DENS ORIEN HAUTDOM STRAT MELAN LN_NID 1.000 PENTE -0.429 1.000 PINS -0.518 0.322 1.000 HAUT -0.425 0.137 0.414 1.000 DIAM -0.201 0.113 0.295 0.905 1.000 DENS -0.528 0.301 0.980 0.439 0.306 1.000 ORIEN -0.230 -0.152 0.128 0.058 -0.079 0.151 1.000 HAUTDOM -0.541 0.262 0.759 0.772 0.596 0.810 0.060 1.000 STRAT -0.594 0.326 0.877 0.460 0.267 0.909 0.063 0.854 1.000 MELAN -0.063 0.129 0.206 -0.045 -0.025 0.130 0.138 0.054 0.175 1.000 Tableau de toutes les régressions linéaires simples pvalue R2a cor AIC BIC DIAM 0.2621705015 0.0094 -0.2009 111.7342 116.2237 ORIEN 0.1985066960 0.0222 -0.2297 111.3058 115.7953 MELAN 0.7282888025 0.0282 -0.0628 112.9637 117.4532 HAUT 0.0136117880 0.1545 -0.4253 106.5101 110.9997 PENTE 0.0126307606 0.1581 -0.4294 106.3671 110.8566 PINS 0.0020270163 0.2445 -0.5178 102.7944 107.2839 DENS 0.0015772396 0.2558 -0.5283 102.2959 106.7854 HAUTDOM 0.0011400381 0.2703 -0.5414 101.6485 106.1380 STRAT 0.0002680646 0.3319 -0.5940 98.7348 103.2243 Il est possible de télécharger la fonction mrls() sur internet, via l'instruction : source("http://forge.info.univ-angers.fr/~gh/Cmi/mrls.r",encoding="latin1") Une fois la préparation des données terminée, il suffit d'un seul appel de la fonction de l'exercice précédent pour obtenir la réponse sur les régressions linéaires simples. ### recherche de la meilleure régression linéaire simple pour PI dans les données LEA ## lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } catss("Régressions linéaires simples") ## lecture des données, on se restreient au QT, suppression des données incorrectes lea <- lit.dar("http://forge.info.univ-angers.fr/~gh/Datasets/lea.dar") cats("Données de départ") print(summary(lea)) filtre <- (lea[ ,"mw"] > 0 ) cat("\n il y a ",sum(!filtre)," données incorrectes (MW<0) \n") leaOK <- lea[ filtre, ] # on ne retient que les quantitatives et on met pi en premier leaOK <- leaOK[ ,c("pi","length","foldindex","mw","gravy") ] cats("Données correctes") ## analyse via la fonction de l'exercice précédent source("mrls.r",encoding="latin1") # lecture du code de la fonction mrls() mrls(leaOK) # tracé pairsi(leaOK) catss("Régressions linéaires multiples") # esai de régression linéaire multiple avec toutes les données cats("Avec toutes les données") mlm <- lm( pi ~ . , data=leaOK) print(mlm) print(anova(mlm)) print(summary(mlm)) # régression linéaire multiple pour les protéines dont la longueur est inférieure à 900 acides aminés cats("On se limite aux petites LEA") petitesLEA <- leaOK[ (leaOK$length<900) , ] mlm2 <- lm( pi ~ . , data=petitesLEA) print(mlm2) print(anova(mlm2)) print(summary(mlm2)) + ============================= + + + + Régressions linéaires simples + + + + ============================= + Données de départ ================= length reign pfam cdd foldindex pi mw gravy genre Min. : 68.0 Length:773 Length:773 Length:773 Min. :-0.45034 Min. :-9.000 Min. : -9 Min. :-9.0000 Length:773 1st Qu.: 130.0 Class :character Class :character Class :character 1st Qu.:-0.19701 1st Qu.: 5.590 1st Qu.: 13961 1st Qu.:-1.3322 Class :character Median : 168.0 Mode :character Mode :character Mode :character Median :-0.11690 Median : 7.100 Median : 17447 Median :-1.1003 Mode :character Mean : 205.7 Mean :-0.09963 Mean : 7.377 Mean : 21839 Mean :-1.0610 3rd Qu.: 236.0 3rd Qu.:-0.03098 3rd Qu.: 9.560 3rd Qu.: 25100 3rd Qu.:-0.8095 Max. :1864.0 Max. : 0.41210 Max. :12.610 Max. :197129 Max. : 0.7991 espece Length:773 Class :character Mode :character iil y a 5 données incorrectes (MW<0) Données correctes ================= pi length foldindex mw gravy Min. : 3.810 Min. : 68.0 Min. :-0.45034 Min. : 7521 Min. :-2.2056 1st Qu.: 5.620 1st Qu.: 130.0 1st Qu.:-0.19496 1st Qu.: 14011 1st Qu.:-1.3210 Median : 7.130 Median : 168.0 Median :-0.11665 Median : 17473 Median :-1.0985 Mean : 7.484 Mean : 205.9 Mean :-0.09975 Mean : 21981 Mean :-1.0093 3rd Qu.: 9.562 3rd Qu.: 235.2 3rd Qu.:-0.03051 3rd Qu.: 25149 3rd Qu.:-0.8059 Max. :12.610 Max. :1864.0 Max. : 0.41210 Max. :197129 Max. : 0.7991 Matrice des corrélations au sens de Pearson pour 768 lignes et 5 colonnes pi length foldindex mw gravy pi 1.000 length -0.241 1.000 foldindex 0.038 0.174 1.000 mw -0.252 0.997 0.165 1.000 gravy 0.017 0.179 0.985 0.173 1.000 Tableau de toutes les régressions linéaires simples pvalue R2a cor AIC BIC foldindex 2.919533e-01 0.0001 0.0381 3291.788 3305.719 gravy 6.353834e-01 0.0010 0.0171 3292.676 3306.608 length 1.336249e-11 0.0568 -0.2409 3246.992 3260.923 mw 1.513903e-12 0.0620 -0.2515 3242.708 3256.640 + =============================== + + + + Régressions linéaires multiples + + + + =============================== + Avec toutes les données ======================= Call: lm(formula = pi ~ ., data = leaOK) Coefficients: (Intercept) length foldindex mw gravy 6.9163040 0.0223285 6.9991193 -0.0002405 -1.9372357 Analysis of Variance Table Response: pi Df Sum Sq Mean Sq F value Pr(>F) length 1 188.46 188.459 48.7938 6.197e-12 *** foldindex 1 21.48 21.480 5.5614 0.01861 * mw 1 72.47 72.467 18.7622 1.679e-05 *** gravy 1 18.41 18.406 4.7654 0.02934 * Residuals 763 2946.99 3.862 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Call: lm(formula = pi ~ ., data = leaOK) Residuals: Min 1Q Median 3Q Max -3.6110 -1.6437 -0.3399 1.8622 6.4797 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.916e+00 6.328e-01 10.930 < 2e-16 *** length 2.233e-02 6.994e-03 3.192 0.001469 ** foldindex 6.999e+00 2.859e+00 2.448 0.014584 * mw -2.405e-04 6.508e-05 -3.695 0.000235 *** gravy -1.937e+00 8.874e-01 -2.183 0.029341 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.965 on 763 degrees of freedom Multiple R-squared: 0.09262, Adjusted R-squared: 0.08786 F-statistic: 19.47 on 4 and 763 DF, p-value: 2.863e-15 On se limite aux petites LEA ============================ Call: lm(formula = pi ~ ., data = petitesLEA) Coefficients: (Intercept) length foldindex mw gravy 7.0463796 0.0242442 7.7530926 -0.0002717 -2.1490011 Analysis of Variance Table Response: pi Df Sum Sq Mean Sq F value Pr(>F) length 1 197.63 197.630 52.1368 1.260e-12 *** foldindex 1 25.63 25.630 6.7615 0.009496 ** mw 1 90.64 90.644 23.9129 1.231e-06 *** gravy 1 21.25 21.249 5.6058 0.018151 * Residuals 759 2877.07 3.791 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Call: lm(formula = pi ~ ., data = petitesLEA) Residuals: Min 1Q Median 3Q Max -3.717 -1.540 -0.275 1.788 5.966 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.046e+00 6.354e-01 11.090 < 2e-16 *** length 2.424e-02 7.224e-03 3.356 0.00083 *** foldindex 7.753e+00 2.940e+00 2.637 0.00854 ** mw -2.717e-04 6.694e-05 -4.059 5.43e-05 *** gravy -2.149e+00 9.077e-01 -2.368 0.01815 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.947 on 759 degrees of freedom Multiple R-squared: 0.1043, Adjusted R-squared: 0.09962 F-statistic: 22.1 on 4 and 759 DF, p-value: < 2.2e-16 Les résultats montrent que PI ne peut pas se modéliser linéairement à partir des autres variables. Se restreindre aux petites LEA ne change rien, il n'y a que 4 protéines plus grandes que 900 aa. Suite à l'exécution du code R ci-dessous, la meilleure régression linéaire multiple pour MW dans les données LEA est, celle qui utilise toutes les variables. ### recherche de la meilleure régression linéaire simple pour MW dans les données LEA ## lecture des fonctions gh if (!exists("cats")) { void <- capture.output( source("http://forge.info.univ-angers.fr/~gh/statgh.r",encoding="latin1") ) } catss("Régressions linéaires simples") ## lecture des données, on se restreient au QT, suppression des données incorrectes lea <- lit.dar("http://forge.info.univ-angers.fr/~gh/Datasets/lea.dar") cats("Données de départ") print(summary(lea)) filtre <- (lea[ ,"mw"] > 0 ) cat("\n il y a ",sum(!filtre)," données incorrectes (MW<0) \n") leaOK <- lea[ filtre, ] # on ne retient que les quantitatives et on met mw en premier leaOK <- leaOK[ ,c("mw","pi","length","foldindex","gravy") ] cats("Données correctes") ## analyse via la fonction de l'exercice précédent source("mrls.r",encoding="latin1") # lecture du code de la fonction mrls() mrls(leaOK) # tracé pairsi(leaOK) catss("Régressions linéaires multiples") # esai de régression linéaire multiple avec toutes les données cats("Avec toutes les données") mlm <- lm( mw ~ . , data=leaOK) print(mlm) print(anova(mlm)) print(summary(mlm)) # utilisation de stepAIC pour trouver la meilleure régression linéaire multiple cats("Utilisation du package MASS") library(MASS) msat <- lm( mw ~ . ,data=leaOK) msel <- stepAIC(msat,direction="both",trace=TRUE) print(summary(msel)) # utilisation de olsrr pour trouver la meilleure régression linéaire multiple cats("Utilisation du package olsrr") library(olsrr) ols_step_all_possible(msat) ols_step_best_subset(msat) # comparaison de m1 avec les variables pi length foldindex gravy # vs m2 avec les variables length foldindex gravy cats("Comparaison des modèles") m1 <- lm( mw ~ pi + length + foldindex + gravy, data=leaOK) m2 <- lm( mw ~ length + foldindex + gravy, data=leaOK) print(anova(m1,m2)) + ============================= + + + + Régressions linéaires simples + + + + ============================= + Matrice des corrélations au sens de Pearson pour 768 lignes et 5 colonnes mw pi length foldindex gravy mw 1.000 pi -0.252 1.000 length 0.997 -0.241 1.000 foldindex 0.165 0.038 0.174 1.000 gravy 0.173 0.017 0.179 0.985 1.000 Tableau de toutes les régressions linéaires simples pvalue R2a cor AIC BIC foldindex 4.547178e-06 0.0258 0.1646 17033.09 17047.02 gravy 1.516555e-06 0.0285 0.1725 17030.97 17044.90 pi 1.513903e-12 0.0620 -0.2515 17003.98 17017.91 length 0.000000e+00 0.9950 0.9975 12985.65 12999.58 + =============================== + + + + Régressions linéaires multiples + + + + =============================== + Avec toutes les données ======================= Call: lm(formula = mw ~ ., data = leaOK) Coefficients: (Intercept) pi length foldindex gravy 2506.69 -73.11 106.97 -10462.83 3014.63 Analysis of Variance Table Response: mw Df Sum Sq Mean Sq F value Pr(>F) pi 1 1.2433e+10 1.2433e+10 10587.944 < 2.2e-16 *** length 1 1.8313e+11 1.8313e+11 155952.594 < 2.2e-16 *** foldindex 1 1.4415e+07 1.4415e+07 12.276 0.0004856 *** gravy 1 4.6599e+07 4.6599e+07 39.682 5.046e-10 *** Residuals 763 8.9599e+08 1.1743e+06 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Call: lm(formula = mw ~ ., data = leaOK) Residuals: Min 1Q Median 3Q Max -5882.2 -534.5 62.2 517.3 6888.5 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.507e+03 3.641e+02 6.885 1.21e-11 *** pi -7.311e+01 1.979e+01 -3.695 0.000235 *** length 1.070e+02 2.755e-01 388.292 < 2e-16 *** foldindex -1.046e+04 1.537e+03 -6.809 1.99e-11 *** gravy 3.015e+03 4.786e+02 6.299 5.05e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1084 on 763 degrees of freedom Multiple R-squared: 0.9954, Adjusted R-squared: 0.9954 F-statistic: 4.165e+04 on 4 and 763 DF, p-value: < 2.2e-16 Utilisation du package MASS =========================== Start: AIC=10738.69 mw ~ pi + length + foldindex + gravy Df Sum of Sq RSS AIC <none> 8.9599e+08 10739 - pi 1 1.6035e+07 9.1202e+08 10750 - gravy 1 4.6599e+07 9.4259e+08 10776 - foldindex 1 5.4445e+07 9.5043e+08 10782 - length 1 1.7705e+11 1.7795e+11 14800 Call: lm(formula = mw ~ pi + length + foldindex + gravy, data = leaOK) Residuals: Min 1Q Median 3Q Max -5882.2 -534.5 62.2 517.3 6888.5 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.507e+03 3.641e+02 6.885 1.21e-11 *** pi -7.311e+01 1.979e+01 -3.695 0.000235 *** length 1.070e+02 2.755e-01 388.292 < 2e-16 *** foldindex -1.046e+04 1.537e+03 -6.809 1.99e-11 *** gravy 3.015e+03 4.786e+02 6.299 5.05e-10 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1084 on 763 degrees of freedom Multiple R-squared: 0.9954, Adjusted R-squared: 0.9954 F-statistic: 4.165e+04 on 4 and 763 DF, p-value: < 2.2e-16 Utilisation du package olsrr ============================ # A tibble: 15 x 6 Index N Predictors `R-Square` `Adj. R-Square` `Mallow's Cp` <int> <int> <chr> <dbl> <dbl> <dbl> 1 1 1 length 0.995 0.995 73.4 2 2 1 pi 0.0633 0.0620 156004. 3 3 1 gravy 0.0298 0.0285 161611. 4 4 1 foldindex 0.0271 0.0258 162059. 5 5 2 pi length 0.995 0.995 53.0 6 6 2 length foldindex 0.995 0.995 60.3 7 7 2 length gravy 0.995 0.995 68.3 8 8 2 pi gravy 0.0945 0.0922 150772. 9 9 2 pi foldindex 0.0936 0.0913 150923. 10 10 2 foldindex gravy 0.0308 0.0282 161445. 11 11 3 length foldindex gravy 0.995 0.995 16.7 12 12 3 pi length foldindex 0.995 0.995 42.7 13 13 3 pi length gravy 0.995 0.995 49.4 14 14 3 pi foldindex gravy 0.0945 0.0910 150774. 15 15 4 pi length foldindex gravy 0.995 0.995 5.00 Best Subsets Regression ---------------------------------------- Model Index Predictors ---------------------------------------- 1 length 2 pi length 3 length foldindex gravy 4 pi length foldindex gravy Subsets Regression Summary ----------------------------------------------------------------------------------------------------------------------------------------------------- Adj. Pred Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC ----------------------------------------------------------------------------------------------------------------------------------------------------- 1 0.9950 0.9950 0.9949 73.4179 12985.6529 10805.8180 12999.5843 1287132.9731 1287124.2327 1678.1452 0.0050 2 0.9951 0.9951 0.995 52.9580 12966.7735 10786.9316 12985.3487 1255890.5853 1255869.2649 1637.4118 0.0049 3 0.9954 0.9953 0.9952 16.6553 12931.8020 10752.2231 12955.0210 1200001.4536 1199964.7846 1564.5444 0.0047 4 0.9954 0.9954 0.9953 5.0000 12920.1788 10740.7547 12948.0416 1181997.0511 1181940.8662 1541.0705 0.0046 ----------------------------------------------------------------------------------------------------------------------------------------------------- AIC: Akaike Information Criteria SBIC: Sawa's Bayesian Information Criteria SBC: Schwarz Bayesian Criteria MSEP: Estimated error of prediction, assuming multivariate normality FPE: Final Prediction Error HSP: Hocking's Sp APC: Amemiya Prediction Criteria Comparaison des modèles ======================= Analysis of Variance Table Model 1: mw ~ pi + length + foldindex + gravy Model 2: mw ~ length + foldindex + gravy Res.Df RSS Df Sum of Sq F Pr(>F) 1 763 895987628 2 764 912022976 -1 -16035348 13.655 0.0002353 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1   Cliquer  ici  pour revenir à la page de départ des cours CMI / M1 Statistiques.

 

 

   Retour à la page principale de   (gH)