## ----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")) dim(CarsNow) summary(CarsNow) with(CarsNow, tapply(consumption, fdrive, mean)) ## ----fig-ParamCovar-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 ~ drive, data = CarsNow, xaxt = "n", col = COL, bg = BGC, pch = PCH, xlab = "Drive", ylab = "Consumption [l/100 km]", cex.lab = 1.2, cex.axis = 1.2) axis(1, at = 1:3, labels = levels(CarsNow[, "fdrive"]), cex.axis = 1.2) ## ----fig-ParamCovar-03-02, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- set.seed(20010911) par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1) plot(CarsNow[, "drive"] + runif(nrow(CarsNow), -0.2, 0.2), CarsNow[, "consumption"], xaxt = "n", col = COL, bg = BGC, pch = PCH, xlab = "Drive", ylab = "Consumption [l/100 km]", cex.lab = 1.2, cex.axis = 1.2) axis(1, at = 1:3, labels = levels(CarsNow[, "fdrive"]), cex.axis = 1.2) ## ------------------------------------------------------------------------ (ybar <- with(CarsNow, mean(consumption))) (ybargr <- with(CarsNow, tapply(consumption, fdrive, mean))) ## ----fig-ParamCovar-03-03, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- layout(matrix(c(1, 2,2,2,2), nrow = 1)) par(bty = BTY, mar = c(4, 4, 1, 1) + 0.1) set.seed(20010911) plot(runif(nrow(CarsNow), -0.1, 0.1), CarsNow[, "consumption"], xaxt = "n", col = "darkblue", bg = "skyblue", pch = PCH, xlab = "", ylab = "Consumption [l/100 km]", xlim = c(-0.2, 0.2), cex.lab = 1.2, cex.axis = 1.2) axis(1, at = 0, labels = "All") points(0, ybar, pch = 22, col = "red", bg = "yellow", cex = 3) set.seed(20010911) plot(CarsNow[, "drive"] + runif(nrow(CarsNow), -0.1, 0.1), CarsNow[, "consumption"], xaxt = "n", col = COL, bg = BGC, pch = PCH, xlab = "Drive", ylab = "Consumption [l/100 km]", cex.lab = 1.2, cex.axis = 1.2) points(1:3, ybargr, pch = 22, col = "darkgreen", bg = "seagreen", cex = 3) axis(1, at = 1:3, labels = levels(CarsNow[, "fdrive"]), cex.axis = 1.2) ## ----fig-ParamCovar-03-04, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- FCOL <- rainbow_hcl(3) names(FCOL) <- levels(CarsNow[, "fdrive"]) par(mfrow = c(1, 1), bty = BTY, mar = c(4, 4, 1, 1) + 0.1) plot(consumption ~ fdrive, data = CarsNow, col = FCOL, xlab = "Drive", ylab = "Consumption [l/100 km]", cex.lab = 1.2, cex.axis = 1.2) ## ----fig-ParamCovar-03-05, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- layout(matrix(c(1, 2,2,2,2), nrow = 1)) par(bty = BTY, mar = c(4, 4, 1, 1) + 0.1) boxplot(CarsNow[, "consumption"], col = "sandybrown", ylab = "Consumption [l/100 km]", cexc.lab = 1.2, cex.axis = 1.2) plot(consumption ~ fdrive, data = CarsNow, col = FCOL, xlab = "Drive", ylab = "Consumption [l/100 km]", cex.lab = 1.2, cex.axis = 1.2) ## ----fig-ParamCovar-03-06, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- XLIM <- with(CarsNow, range(consumption)) layout(matrix(c(0,1,1,0, 2,2,3,3), nrow = 2, byrow = TRUE)) par(bty = BTY, mar = c(4, 4, 4, 1) + 0.1) for (fl in levels(CarsNow[, "fdrive"])){ hist(subset(CarsNow, fdrive == fl)[, "consumption"], prob = TRUE, col = FCOL[fl], xlab = "Consumption [l/100 km]", xlim = XLIM, ylim = c(0, 0.42), main = fl, cex.main = 1.6, cex.lab = 1.4, cex.axis = 1.4) } par(mfrow = c(1, 1)) ## ----fig-ParamCovar-03-07, fig.ext='png', fig.path=FIGDIRPNG, fig.width=9, fig.height=6---- layout(matrix(c(0,1,1,0, 2,2,3,3), nrow = 2, byrow = TRUE)) par(bty = BTY, mar = c(4, 4, 4, 2) + 0.1) for (fl in levels(CarsNow[, "fdrive"])){ qqnorm(subset(CarsNow, fdrive == fl)[, "consumption"], main = fl, pch = PCH, col = COL, bg = BGC, cex.main = 1.6, cex.lab = 1.4, cex.axis = 1.4) qqline(subset(CarsNow, fdrive == fl)[, "consumption"], lwd = 2, col = "blue") } par(mfrow = c(1, 1)) ## ------------------------------------------------------------------------ mF <- lm(consumption ~ -1 + fdrive, data = CarsNow) summary(mF) ## ------------------------------------------------------------------------ deviance(mF) ## SS_e ## ------------------------------------------------------------------------ m0 <- lm(consumption ~ 1, data = CarsNow) summary(m0) ## ------------------------------------------------------------------------ deviance(m0) ## SS_T ## ------------------------------------------------------------------------ deviance(m0) - deviance(mF) ## SS_R ## ------------------------------------------------------------------------ a1 <- aov(consumption ~ fdrive, data = CarsNow) summary(a1) ## ------------------------------------------------------------------------ (ybar <- with(CarsNow, mean(consumption))) ## ------------------------------------------------------------------------ (ybargr <- with(CarsNow, tapply(consumption, fdrive, mean))) ## ------------------------------------------------------------------------ (meanybargr <- mean(ybargr)) ## ------------------------------------------------------------------------ (TAB <- with(CarsNow, table(fdrive))) (PTAB <- prop.table(TAB)) (wmeanybargr <- sum(ybargr * PTAB)) ## ------------------------------------------------------------------------ levels(CarsNow[["fdrive"]]) ## Which group is the first one? (G <- length(levels(CarsNow[, "fdrive"]))) ## ------------------------------------------------------------------------ mTrt <- lm(consumption ~ fdrive, data = CarsNow) summary(mTrt) ## ------------------------------------------------------------------------ (CTrt <- contr.treatment(G)) ## ------------------------------------------------------------------------ coef(mTrt)["(Intercept)"] + CTrt %*% coef(mTrt)[-1] muhat <- as.numeric(coef(mTrt)["(Intercept)"] + CTrt %*% coef(mTrt)[-1]) names(muhat) <- levels(CarsNow[, "fdrive"]) print(muhat) ## ------------------------------------------------------------------------ mSAS <- lm(consumption ~ fdrive, data = CarsNow, contrasts = list(fdrive = contr.SAS)) summary(mSAS) ## ------------------------------------------------------------------------ (CSAS <- contr.SAS(G)) ## ------------------------------------------------------------------------ coef(mSAS)["(Intercept)"] + CSAS %*% coef(mSAS)[-1] muhat <- as.numeric(coef(mSAS)["(Intercept)"] + CSAS %*% coef(mSAS)[-1]) names(muhat) <- levels(CarsNow[, "fdrive"]) print(muhat) ## ------------------------------------------------------------------------ mSum <- lm(consumption ~ fdrive, data = CarsNow, contrasts = list(fdrive = contr.sum)) summary(mSum) ## ------------------------------------------------------------------------ (CSum <- contr.sum(G)) ## ------------------------------------------------------------------------ coef(mSum)["(Intercept)"] + CSum %*% coef(mSum)[-1] muhat <- as.numeric(coef(mSum)["(Intercept)"] + CSum %*% coef(mSum)[-1]) names(muhat) <- levels(CarsNow[, "fdrive"]) print(muhat) ## ------------------------------------------------------------------------ alphaSum <- as.numeric(CSum %*% coef(mSum)[-1]) names(alphaSum) <- levels(CarsNow[, "fdrive"]) print(alphaSum) ## ------------------------------------------------------------------------ library("mffSM") LSum <- cbind(0, CSum) rownames(LSum) <- levels(CarsNow[["fdrive"]]) colnames(LSum) <- names(coef(mSum)) print(LSum) LSest(mSum, L = LSum) ## ------------------------------------------------------------------------ (ng <- with(CarsNow, table(fdrive))) CwSum <- rbind(diag(G - 1), - ng[-G] / ng[G]) rownames(CwSum) <- levels(CarsNow[["fdrive"]]) print(CwSum) ## ------------------------------------------------------------------------ mwSum <- lm(consumption ~ fdrive, data = CarsNow, contrasts = list(fdrive = CwSum)) summary(mwSum) ## ------------------------------------------------------------------------ coef(mwSum)["(Intercept)"] + CwSum %*% coef(mwSum)[-1] muhat <- as.numeric(coef(mwSum)["(Intercept)"] + CwSum %*% coef(mwSum)[-1]) names(muhat) <- levels(CarsNow[, "fdrive"]) print(muhat) ## ------------------------------------------------------------------------ alphawSum <- as.numeric(CwSum %*% coef(mwSum)[-1]) names(alphawSum) <- levels(CarsNow[, "fdrive"]) print(alphawSum) ## ------------------------------------------------------------------------ LwSum <- cbind(0, CwSum) rownames(LwSum) <- levels(CarsNow[["fdrive"]]) colnames(LwSum) <- names(coef(mwSum)) print(LwSum) LSest(mwSum, L = LwSum) ## ------------------------------------------------------------------------ mHelmert <- lm(consumption ~ fdrive, data = CarsNow, contrasts = list(fdrive = contr.helmert)) summary(mHelmert) ## ------------------------------------------------------------------------ (CHelmert <- contr.helmert(G)) ## ------------------------------------------------------------------------ coef(mHelmert)["(Intercept)"] + CHelmert %*% coef(mHelmert)[-1] muhat <- as.numeric(coef(mHelmert)["(Intercept)"] + CHelmert %*% coef(mHelmert)[-1]) names(muhat) <- levels(CarsNow[, "fdrive"]) print(muhat)