FONCTION lm() et PREVISION COMMANDES S-Plus et quelques sorties ==================================== Creer un objet Splus a partir d'un fichier Unix >d652446.dat_scan("d652446.dat") Prendre les 78 premières données : >donnees<-d652446.dat[1:78] Generer la variable temps > temps<-1:78 ou > temps_seq(1,78,1) > temps [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 [76] 76 77 78 Graphique de la serie >plot(temps, donnees, type = "b") >tsplot(donnees,type = "b", pch = "*") - type ="b" pour tracer des lignes et des points - pch = "*" pour que le symbole des points soient * >tsplot(donnees,type = "b", pch = "*", ylab = "Dollars", xlab = "temps en mois") - Noter que la variable temps n'a pas été utilisée - tsplot: "time series plot" - pour nommer les coordonnées. >title(" Prix moyen par unite, voitures part. Can & EU") Estimation par moindres carres avec la fonction "lm" =================================================== Estimation de la droite de regression avec la variable explicative "temps" Les resultats sont dans l'objet "lm.prix.out" >lm.prix.out <- lm(donnees ~ temps) > lm.prix.out Call: lm(formula = donnees ~ temps) Coefficients: (Intercept) temps 2423.92 10.9595 Degrees of freedom: 78 total; 76 residual Residual standard error: 78.13619 Certaines fonctions Splus operent sur des sorties "lm" : summary, anova, coefficients, formula, fitted, residuals, etc. > summary(lm.prix.out) Call: lm(formula = donnees ~ temps) Residuals: Min 1Q Median 3Q Max -143.6 -63.36 6.545 63.27 150.6 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) 2423.9204 17.8659 135.6732 0.0000 temps 10.9595 0.3929 27.8903 0.0000 Residual standard error: 78.14 on 76 degrees of freedom Multiple R-Squared: 0.911 F-statistic: 777.9 on 1 and 76 degrees of freedom, the p-value is 0 Correlation of Coefficients: (Intercept) temps -0.8688 > anova(lm.prix.out) Analysis of Variance Table Response: donnees Terms added sequentially (first to last) Df Sum of Sq Mean Sq F Value Pr(F) temps 1 4749111 4749111 777.8715 0 Residuals 76 464000 6105 > formula(lm.prix.out) donnees ~ temps > coef(lm.prix.out) (Intercept) temps 2423.92 10.9595 > resid(lm.prix.out) 1 2 3 4 5 6 7 8 73.12009 104.1606 99.2011 112.2416 102.2821 75.32261 35.36311 59.40362 9 10 11 12 13 14 15 16 6.444121 -29.51537 4.525129 -12.43437 -11.39386 6.64664 54.68714 16.72765 17 18 19 20 21 22 23 24 23.76815 24.80866 -20.15084 -38.11034 -92.06983 -131.0293 -67.98883 32.05168 25 26 27 28 29 30 31 32 4.092182 -33.86731 -137.8268 -54.78631 -97.7458 -104.7053 -124.6648 -143.6243 33 34 35 36 37 38 39 40 -140.5838 40.45672 113.4972 116.5377 150.5782 148.6187 56.65924 -17.30026 41 42 43 44 45 46 47 -42.25976 -66.21925 -109.1787 -102.1382 -93.09774 -133.0572 -131.0167 48 49 50 51 52 53 54 55 -42.97623 67.06427 67.10478 20.14528 -6.814215 -41.77371 -73.73321 -86.6927 56 57 58 59 60 61 62 63 -87.6522 -99.6117 -75.57119 -15.53069 63.50982 62.55032 34.59082 100.6313 64 65 66 67 68 69 70 71 72.67183 53.71233 27.75284 13.79334 -5.166155 13.87435 75.91485 78.95536 72 73 74 75 76 77 78 112.9959 79.03636 65.07687 44.11737 16.15788 -23.80162 -36.76112 > summary(resid(lm.prix.out)) Min. 1st Qu. Median Mean 3rd Qu. -1.436243e+02 -6.336102e+01 6.545381e+00 6.376666e-16 6.326994e+01 Max. 1.505782e+02 > La commande suivante permet de visualiser les 6 graphiques suivants: > plot(lm.prix.out,ask=T) Make a plot selection (or 0 to exit): 1: plot: All 2: plot: Residuals vs Fitted Values 3: plot: Sqrt of abs(Residuals) vs Fitted Values 4: plot: Response vs Fitted Values 5: plot: Normal QQplot of Residuals 6: plot: r-f spread plot 7: plot: Cook's Distances Selection: On peut visualiser les 6 graphiques simultanement en utilisant la fonction par(): > par(mfrow=c(3,2)) - divise le graphe en un graphe de dimension 3 x 2. > plot(lm.prix.out) REMARQUE Ce n'est pas pertinent ici mais si on veut enlever la constante dans l'equation de regression, il suffit d'ajouter "-1" dans l'argument de la fonction "lm". > lm1.prix.out<-lm(donnees ~ temps -1) > summary(lm1.prix.out) Call: lm(formula = donnees ~ temps - 1) Residuals: Min 1Q Median 3Q Max -1226 -322.2 614.1 1434 2451 Coefficients: Value Std. Error t value Pr(>|t|) temps 57.2764 3.0148 18.9983 0.0000 Residual standard error: 1211 on 77 degrees of freedom Multiple R-Squared: 0.8242 F-statistic: 360.9 on 1 and 77 degrees of freedom, the p-value is 0 Autocorrelations des residus: ============================ On obtient les 24 premieres autocorrelations des residus avec la commande >acf(resid(lm.prix.out),lag.max=24) L'ajout d'un terme au carre ameliore-t-il le modele? ==================================================== > lm2.prix.out <- lm(donnees ~ temps + temps^2) Les variables sont incluses dans le modele de facon sequentielle: > anova(lm2.prix.out) Analysis of Variance Table Response: donnees Terms added sequentially (first to last) Df Sum of Sq Mean Sq F Value Pr(F) temps 1 4749111 4749111 960.4362 0.000000e+00 I(temps^2) 1 93144 93144 18.8370 4.387829e-05 Residuals 75 370856 4945 > coef(lm2.prix.out) (Intercept) temps I(temps^2) 2504.222 4.936904 0.07623535 Bien que le coefficient de temps^2 soit significatif a 5%, le terme quadratique ne ressort pas dans le graphique de "valeurs ajustees". Effectuer des previsions ======================== NOTE: La fonction predict() utilisee ici est une fonction "locale" et n'est donc pas dans la librairie de base de S-Plus. Avec le premier modele, nous pouvons obtenir des previsions pour les 6 derniers mois de l'annee 98: > length(donnees) > [1] 78 > temps.futur_cbind(1, 79:84) > temps.futur [,1] [,2] [1,] 1 79 [2,] 1 80 [3,] 1 81 [4,] 1 82 [5,] 1 83 [6,] 1 84 > predict(lm.prix.out,temps.futur) [1] 3289.721 3300.680 3311.640 3322.599 3333.559 3344.518 Evaluation de la performance previsionnelle =========================================== Laissons de côté les 6 observations; reestimons le modele avec les 72 premieres observations et ensuite, nous comparerons les previsions aux observations pour les 6 premiers mois de la prochaine année. > lm3.prix.out_lm(donnees[1:72] ~ temps[1:72]) > coef(lm3.prix.out) (Intercept) temps[1:72] 2427.994 10.79317 > anova(lm3.prix.out) Analysis of Variance Table Response: donnees[1:72] Terms added sequentially (first to last) Df Sum of Sq Mean Sq F Value Pr(F) temps[1:72] 1 3622684 3622684 565.7345 0 Residuals 70 448245 6404 > temps.futur1_cbind(1, 73:78) > temps.futur1 [,1] [,2] [1,] 1 73 [2,] 1 74 [3,] 1 75 [4,] 1 76 [5,] 1 77 [6,] 1 78 > prev_predict(lm3.prix.out, temps.futur1) > prev [1] 3215.895 3226.688 3237.481 3248.275 3259.068 3269.861 > prev.obs_cbind(73:78,d652446.dat[73:78],prev) > prev.obs prev [1,] 73 3303 3215.895 [2,] 74 3300 3226.688 [3,] 75 3290 3237.481 [4,] 76 3273 3248.275 [5,] 77 3244 3259.068 [6,] 78 3242 3269.861 On peut alors calculer l'erreur relative de prevision en % pour chacun des 6 mois: > 100*(prev.obs[,2] - prev.obs[,3])/prev.obs[,2] [1] 2.6371435 2.2215661 1.5963073 0.7554336 -0.4644830 -0.8593769 En conclusion, bien que le modele ne semble pas tres adéquat ( il faut faire une analyse residuelle plus poussee), les previsions ne sont pas si pires!