## ----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) ## ------------------------------------------------------------------------ mInter <- lm(consumption ~ fdrive + lweight + fdrive:lweight, data = CarsNow) summary(mInter) ## ------------------------------------------------------------------------ anova(mInter) ## ------------------------------------------------------------------------ ldrive <- levels(CarsNow[, "fdrive"]) CZ <- contr.treatment(3) rownames(CZ) <- ldrive print(CZ) L.slope <- cbind(0, 0, 0, 1, CZ) colnames(L.slope) <- names(coef(mInter)) rownames(L.slope) <- ldrive print(L.slope) library("mffSM") LSest(mInter, L = L.slope) print(LSest(mInter, L = L.slope), digits = 3) ## ------------------------------------------------------------------------ L.diff.slope <- cbind(0, 0, 0, 0, rbind(CZ[2, ] - CZ[1, ], CZ[3, ] - CZ[1, ], CZ[3, ] - CZ[2, ])) colnames(L.diff.slope) <- names(coef(mInter)) rownames(L.diff.slope) <- paste(ldrive[c(2, 3, 3)], "-", ldrive[c(1, 1, 2)], sep = "") print(L.diff.slope) LSunadj <- LSest(mInter, L = L.diff.slope) print(LSunadj, digits = 3) ## ------------------------------------------------------------------------ LSBonf <- LSest(mInter, L = L.diff.slope, conf.level = 1 - 0.05/3) LSBonf[["P value"]] <- 3 * LSBonf[["P value"]] print(LSBonf, digits = 3) ## ------------------------------------------------------------------------ library("multcomp") ## ------------------------------------------------------------------------ mcp <- glht(mInter, linfct = L.diff.slope, rhs = 0) print(mcp) ## ------------------------------------------------------------------------ set.seed(20010911) smcp <- summary(mcp) print(smcp) ## ------------------------------------------------------------------------ set.seed(20010911) cimcp <- confint(mcp, level = 0.95) print(cimcp) ## ------------------------------------------------------------------------ MCP <- data.frame(Estimate = LSunadj[["Estimate"]], PvalUnadj = format.pval(LSunadj[["P value"]], digits = 2, eps = 1e-4), PvalBonf = format.pval(LSBonf[["P value"]], digits = 2, eps = 1e-4), PvalHBW = format.pval(smcp[["test"]][["pvalues"]], digits = 2, eps = 1e-4), LowerUnadj = LSunadj[["Lower"]], UpperUnadj = LSunadj[["Upper"]], LowerBonf = LSBonf[["Lower"]], UpperBonf = LSBonf[["Upper"]], LowerHBW = cimcp[["confint"]][, "lwr"], UpperHBW = cimcp[["confint"]][, "upr"], stringsAsFactors = FALSE) print(MCP[, 1:4], digits = 3) print(MCP[, 5:10], digits = 3)