> library(MASS) > data(hills) > (hills.lm <- lm(time ~ dist + climb, data = hills)) Call: lm(formula = time ~ dist + climb, data = hills) Coefficients: (Intercept) dist climb -8.99204 6.21796 0.01105 > summary(hills.lm) Call: lm(formula = time ~ dist + climb, data = hills) Residuals: Min 1Q Median 3Q Max -16.215 -7.129 -1.186 2.371 65.121 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -8.992039 4.302734 -2.090 0.0447 * dist 6.217956 0.601148 10.343 9.86e-12 *** climb 0.011048 0.002051 5.387 6.45e-06 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 14.68 on 32 degrees of freedom Multiple R-Squared: 0.9191, Adjusted R-squared: 0.914 F-statistic: 181.7 on 2 and 32 DF, p-value: < 2.2e-16 > plot(fitted(hills.lm), studres(hills.lm)) > title("externally studentized residuals vs. fitted values") > plot(fitted(hills.lm), stdres(hills.lm)) > title("internally studentized residuals vs. fitted values") > plot(fitted(hills.lm), studres(hills.lm)) > points(fitted(hills.lm), stdres(hills.lm), col = "green") #### robust regression, section 6.5 in MASS > data(phones) > attach(phones) > out <- lm(calls ~ year) > summary(out) Call: lm(formula = calls ~ year) Residuals: Min 1Q Median 3Q Max -78.97 -33.52 -12.04 23.38 124.20 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -260.059 102.607 -2.535 0.0189 * year 5.041 1.658 3.041 0.0060 ** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Residual standard error: 56.22 on 22 degrees of freedom Multiple R-Squared: 0.2959, Adjusted R-squared: 0.2639 F-statistic: 9.247 on 1 and 22 DF, p-value: 0.005998 > plot(year, calls) > abline(out) > library(lqs) > out.lqs <- lqs(calls ~ year) > summary(out.lqs) Length Class Mode crit 1 -none- numeric sing 1 -none- character coefficients 2 -none- numeric bestone 2 -none- numeric fitted.values 24 -none- numeric residuals 24 -none- numeric scale 2 -none- numeric terms 3 terms call call 2 -none- call xlevels 0 -none- list model 2 data.frame list > out.lqs Call: lqs.formula(formula = calls ~ year) Coefficients: (Intercept) year -56.162 1.159 Scale estimates 1.249 1.131 > out Call: lm(formula = calls ~ year) Coefficients: (Intercept) year -260.059 5.041 > abline(out.lqs, col = "red") > out.rlm <- rlm(calls ~ year, method = "MM") > summary(out.rlm) Call: rlm.formula(formula = calls ~ year, method = "MM") Residuals: Min 1Q Median 3Q Max -1.7442 -0.4543 0.2148 39.0082 188.4577 Coefficients: Value Std. Error t value (Intercept) -52.4230 2.9159 -17.9783 year 1.1009 0.0471 23.3669 Residual standard error: 2.129 on 22 degrees of freedom Correlation of Coefficients: (Intercept) year -0.9937 > abline(out.rlm, col = "seagreen") # ANOVA, section 6.7 in MASS > aov(yield ~ block + N * P * K, data = npk) Call: aov(formula = yield ~ block + N * P * K, data = npk) Terms: block N P K N:P N:K P:K Sum of Squares 343.2950 189.2817 8.4017 95.2017 21.2817 33.1350 0.4817 Deg. of Freedom 5 1 1 1 1 1 1 Residuals Sum of Squares 185.2867 Deg. of Freedom 12 Residual standard error: 3.929447 1 out of 13 effects not estimable Estimated effects may be unbalanced > out <- aov(yield ~ block + N * P * K, data = npk) > out2 <- aov(yield ~ block + N + P + K, data = npk) > anova(out2, out) Analysis of Variance Table Model 1: yield ~ block + N + P + K Model 2: yield ~ block + N * P * K Res.Df RSS Df Sum of Sq F Pr(>F) 1 15 240.185 2 12 185.287 3 54.898 1.1852 0.3566 > out2 Call: aov(formula = yield ~ block + N + P + K, data = npk) Terms: block N P K Residuals Sum of Squares 343.2950 189.2817 8.4017 95.2017 240.1850 Deg. of Freedom 5 1 1 1 15 Residual standard error: 4.001541 Estimated effects may be unbalanced > plot(fitted(out2), studres(out2)) > abline(h = 0)