## ----echo=FALSE, results='hide'------------------------------------------ source(".Rprofile") ## ------------------------------------------------------------------------ data(Cars2004nh, package = "mffSM") head(Cars2004nh) dim(Cars2004nh) summary(Cars2004nh) ## ------------------------------------------------------------------------ isComplete <- complete.cases(Cars2004nh[, c("consumption", "lweight", "engine.size")]) sum(!isComplete) CarsNow <- subset(Cars2004nh, isComplete, select = c("consumption", "drive", "fdrive", "weight", "lweight", "engine.size", "horsepower")) dim(CarsNow) summary(CarsNow) ## ------------------------------------------------------------------------ m <- lm(consumption ~ lweight, data = CarsNow) summary(m) ## ----fig-CheckModelAssumpt-03-01, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1) plot(consumption ~ lweight, data = CarsNow, pch = PCH, col = COL2, bg = BGC2, xlab = "Log(weight) [log(kg)]", ylab = "Consumption [l/100 km]") abline(m, col = "red2", lwd = 2) #lines(lowess(CarsNow[, "lweight"], CarsNow[, "consumption"]), col = "blue", lwd = 2) par(mar = c(4, 4, 2, 1) + 0.1) ## ----fig-CheckModelAssumpt-03-02, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- plot(m, which = 1, pch = 21, col = "blue4", bg = "skyblue") ## ----fig-CheckModelAssumpt-03-03, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- plot(m, which = 2, pch = 21, col = "blue4", bg = "skyblue") ## ----fig-CheckModelAssumpt-03-04, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- plot(m, which = 3, pch = 21, col = "blue4", bg = "skyblue") ## ----fig-CheckModelAssumpt-03-05, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- library("mffSM") plotLM(m) ## ------------------------------------------------------------------------ m <- lm(consumption ~ lweight, data = CarsNow, x = TRUE) ## ----fig-CheckModelAssumpt-03-06, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(1, 3), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) ### Scatterplot + fitted line plot(consumption ~ lweight, data = CarsNow, pch = PCH, col = COL2, bg = BGC2, xlab = "Log(weight) [log(kg)]", ylab = "Consumption [l/100 km]") abline(m, col = "red2", lwd = 2) ### Residuals versus fitted values (+ lowess) plot(fitted(m), residuals(m), pch = 21, col = "blue4", bg = "skyblue", xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs Fitted") lines(lowess(fitted(m), residuals(m)), col = "red") ### Residuals versus the covariate (+ lowess) plot(m[["x"]][, "lweight"], residuals(m), pch = 21, col = "blue4", bg = "skyblue", xlab = "Log(weight) [log(kg)]", ylab = "Residuals", main = "Residuals vs Covariate") lines(lowess(m[["x"]][, "lweight"], residuals(m)), col = "red") ## ----fig-CheckModelAssumpt-03-07, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(1, 2), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) ### Residuals versus engine.size (+ lowess) plot(CarsNow[, "engine.size"], residuals(m), pch = 21, col = COL2, bg = BGC2, xlab = "Engine size [l]", ylab = "Residuals", main = "Residuals vs Omitted covariate 1") lines(lowess(CarsNow[, "engine.size"], residuals(m)), col = "red") ### Residuals versus horsepower (+ lowess) plot(CarsNow[, "horsepower"], residuals(m), pch = 21, col = COL2, bg = BGC2, xlab = "Horsepower", ylab = "Residuals", main = "Residuals vs Omitted covariate 2") lines(lowess(CarsNow[, "horsepower"], residuals(m)), col = "red") ## ------------------------------------------------------------------------ m2a <- lm(consumption ~ lweight + engine.size, data = CarsNow, x = TRUE) summary(m2a) ## ----fig-CheckModelAssumpt-03-08, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(2, 2), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) ### Residuals versus fitted values (+ lowess) plot(fitted(m2a), residuals(m2a), pch = 21, col = "blue4", bg = "skyblue", xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs Fitted") lines(lowess(fitted(m2a), residuals(m2a)), col = "red") ### Residuals versus lweight (included covariate) plot(m2a[["x"]][, "lweight"], residuals(m2a), pch = 21, col = "blue4", bg = "skyblue", xlab = "Log(weight) [log(kg)]", ylab = "Residuals", main = "Residuals vs Covariate 1") lines(lowess(m2a[["x"]][, "lweight"], residuals(m2a)), col = "red") ### Residuals versus engine.size (included covariate) plot(m2a[["x"]][, "engine.size"], residuals(m2a), pch = 21, col = "blue4", bg = "skyblue", xlab = "Engine size [l]", ylab = "Residuals", main = "Residuals vs Covariate 2") lines(lowess(m2a[["x"]][, "engine.size"], residuals(m2a)), col = "red") ### Residuals versus horsepower (potentially omitted covariate) plot(CarsNow[, "horsepower"], residuals(m2a), pch = 21, col = COL2, bg = BGC2, xlab = "Horsepower", ylab = "Residuals", main = "Residuals vs Omitted covariate") lines(lowess(CarsNow[, "horsepower"], residuals(m2a)), col = "red") ## ------------------------------------------------------------------------ m2b <- lm(consumption ~ lweight + horsepower, data = CarsNow, x = TRUE) summary(m2b) ## ----fig-CheckModelAssumpt-03-09, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(2, 2), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) ### Residuals versus fitted values (+ lowess) plot(fitted(m2b), residuals(m2b), pch = 21, col = "blue4", bg = "skyblue", xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs Fitted") lines(lowess(fitted(m2b), residuals(m2b)), col = "red") ### Residuals versus lweight (included covariate) plot(m2b[["x"]][, "lweight"], residuals(m2b), pch = 21, col = "blue4", bg = "skyblue", xlab = "Log(weight) [log(kg)]", ylab = "Residuals", main = "Residuals vs Covariate 1") lines(lowess(m2b[["x"]][, "lweight"], residuals(m2b)), col = "red") ### Residuals versus horsepower (included covariate) plot(m2b[["x"]][, "horsepower"], residuals(m2b), pch = 21, col = "blue4", bg = "skyblue", xlab = "Horsepower", ylab = "Residuals", main = "Residuals vs Covariate 2") lines(lowess(m2b[["x"]][, "horsepower"], residuals(m2b)), col = "red") ### Residuals versus engine.size (potentially omitted covariate) plot(CarsNow[, "engine.size"], residuals(m2b), pch = 21, col = COL2, bg = BGC2, xlab = "Engine size [l]", ylab = "Residuals", main = "Residuals vs Omitted covariate") lines(lowess(CarsNow[, "engine.size"], residuals(m2b)), col = "red") ## ------------------------------------------------------------------------ m3 <- lm(consumption ~ lweight + engine.size + horsepower, data = CarsNow, x = TRUE) summary(m3) ## ----fig-CheckModelAssumpt-03-10, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(2, 2), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) ### Residuals versus fitted values (+ lowess) plot(fitted(m3), residuals(m3), pch = 21, col = "blue4", bg = "skyblue", xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs Fitted") lines(lowess(fitted(m3), residuals(m3)), col = "red") ### Residuals versus lweight (included covariate) plot(m3[["x"]][, "lweight"], residuals(m3), pch = 21, col = "blue4", bg = "skyblue", xlab = "Log(weight) [log(kg)]", ylab = "Residuals", main = "Residuals vs Covariate 1") lines(lowess(m3[["x"]][, "lweight"], residuals(m3)), col = "red") ### Residuals versus engine.size (included covariate) plot(m3[["x"]][, "engine.size"], residuals(m3), pch = 21, col = "blue4", bg = "skyblue", xlab = "Engine size [l]", ylab = "Residuals", main = "Residuals vs Covariate 2") lines(lowess(m3[["x"]][, "engine.size"], residuals(m3)), col = "red") ### Residuals versus horsepower (included covariate) plot(m3[["x"]][, "horsepower"], residuals(m3), pch = 21, col = "blue4", bg = "skyblue", xlab = "Horsepower", ylab = "Residuals", main = "Residuals vs Covariate 3") lines(lowess(m3[["x"]][, "horsepower"], residuals(m3)), col = "red") ## ------------------------------------------------------------------------ ybar <- with(CarsNow, mean(consumption)) prY <- residuals(m3, type = "partial") + ybar YLIM <- range(c(CarsNow[, "consumption"], prY)) ## ----fig-CheckModelAssumpt-03-11, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(1, 3), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) plot(CarsNow[, "lweight"], prY[, "lweight"], ylim = YLIM, main = "Log(weight)", xlab = "Log(weight) [log(kg)]", ylab = "Response-mean partial residuals [l/100 km]", pch = PCH3, col = COL3, bg = BGC3) lines(lowess(CarsNow[, "lweight"], prY[, "lweight"]), col = "red3") # plot(CarsNow[, "engine.size"], prY[, "engine.size"], ylim = YLIM, main = "Engine size", xlab = "Engine size", ylab = "Response-mean partial residuals [l/100 km]", pch = PCH3, col = COL3, bg = BGC3) lines(lowess(CarsNow[, "engine.size"], prY[, "engine.size"]), col = "red3") # plot(CarsNow[, "horsepower"], prY[, "horsepower"], ylim = YLIM, main = "Horsepower", xlab = "Horsepower", ylab = "Response-mean partial residuals [l/100 km]", pch = PCH3, col = COL3, bg = BGC3) lines(lowess(CarsNow[, "horsepower"], prY[, "horsepower"]), col = "red3") ## ------------------------------------------------------------------------ mW <- lm(consumption ~ lweight, data = CarsNow) mWE <- lm(consumption ~ lweight + engine.size, data = CarsNow) mWH <- lm(consumption ~ lweight + horsepower, data = CarsNow) mWEH <- lm(consumption ~ lweight + engine.size + horsepower, data = CarsNow) ## ------------------------------------------------------------------------ anova(mW, mWE) ## ------------------------------------------------------------------------ anova(mW, mWH) ## ------------------------------------------------------------------------ anova(mW, mWEH) ## ----fig-CheckModelAssumpt-03-12, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(2, 2), bty = BTY, mar = c(5, 4, 3, 1) + 0.1) # plot(m3, which = 3, pch = 21, col = "blue4", bg = "skyblue", lwd = 2) # XLAB <- c("Weight [kg]", "Engine size [l]", "Horsepower") for (j in 1:3){ plot(m3[["x"]][, j + 1], sqrt(abs(rstandard(m3))), pch = 21, col = "blue4", bg = "skyblue", xlab = XLAB[j], ylab = as.expression(substitute(sqrt(abs(yL)), list(yL = as.name("Standardized residuals")))), main = paste("Scale-Location (", colnames(m3[["x"]])[j + 1], ")", sep = "")) lines(lowess(m3[["x"]][, j + 1], sqrt(abs(rstandard(m3)))), col = "red3", lwd = 2) } ## ------------------------------------------------------------------------ library("lmtest") bptest(m3, varformula = ~fitted(m3)) bptest(m3, varformula = ~lweight, data = CarsNow) bptest(m3, varformula = ~engine.size, data = CarsNow) bptest(m3, varformula = ~horsepower, data = CarsNow) ## ----fig-CheckModelAssumpt-03-13, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 3, 1) + 0.1) plot(m3, which = 2, pch = 21, col = "blue4", bg = "skyblue") ## ------------------------------------------------------------------------ shapiro.test(residuals(m3)) ## ------------------------------------------------------------------------ shapiro.test(rstandard(m3)) ## ------------------------------------------------------------------------ library("nortest") lillie.test(residuals(m3)) ## ------------------------------------------------------------------------ ad.test(residuals(m3))