R Under development (unstable) (2023-06-12 r84534 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > library("multcomp") Loading required package: mvtnorm Loading required package: survival Loading required package: TH.data Loading required package: MASS Attaching package: 'TH.data' The following object is masked from 'package:MASS': geyser > set.seed(290875) > > testdata <- data.frame(y = rnorm(21), + f1 <- factor(c(rep(c("A", "B", "C"), 7))), + f2 <- factor(c(rep("D", 10), rep("E", 11))), + x <- rnorm(21)) > > # one-way ANOVA > coef(amod <- aov(y ~ f1, data = testdata)) (Intercept) f1B f1C -0.4394751 0.5151680 0.6886101 > glht(amod, linfct = mcp(f1 = "Dunnett")) General Linear Hypotheses Multiple Comparisons of Means: Dunnett Contrasts Linear Hypotheses: Estimate B - A == 0 0.5152 C - A == 0 0.6886 > > # and a continuous covariable: ANCOVA > coef(lmod <- lm(y ~ f1 + x, data = testdata)) (Intercept) f1B f1C x -0.434528566 0.509444592 0.686181780 -0.009491201 > glht(lmod, linfct = mcp(f1 = "Dunnett")) General Linear Hypotheses Multiple Comparisons of Means: Dunnett Contrasts Linear Hypotheses: Estimate B - A == 0 0.5094 C - A == 0 0.6862 > > # ANCOVA with an additional factor as covariable > coef(lmod <- lm(y ~ f1 + f2 + x, data = testdata)) (Intercept) f1B f1C f2E x -0.40849498 0.51296437 0.69200699 -0.05266965 -0.01613183 > glht(lmod, linfct = mcp(f1 = "Dunnett")) General Linear Hypotheses Multiple Comparisons of Means: Dunnett Contrasts Linear Hypotheses: Estimate B - A == 0 0.513 C - A == 0 0.692 > > # and with interaction terms > coef(lmod <- lm(y ~ f1 + f2 + f2:f1 + x, data = testdata)) (Intercept) f1B f1C f2E x f1B:f2E -0.44532319 0.70282663 0.65613337 0.05552324 -0.03443721 -0.37862471 f1C:f2E 0.02753451 > glht(lmod, linfct = mcp(f1 = "Dunnett")) General Linear Hypotheses Multiple Comparisons of Means: Dunnett Contrasts Linear Hypotheses: Estimate B - A == 0 0.7028 C - A == 0 0.6561 Warning message: In mcp2matrix(model, linfct = linfct) : covariate interactions found -- default contrast might be inappropriate > > # with contrasts as expressions > glht(lmod, linfct = mcp(f1 = c("B - A = 0", "C - A = 0"))) General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Linear Hypotheses: Estimate B - A == 0 0.7028 C - A == 0 0.6561 Warning message: In mcp2matrix(model, linfct = linfct) : covariate interactions found -- default contrast might be inappropriate > > tmp <- multcomp:::chrlinfct2matrix(c(l1 = "x1 - x2 = 2", + l2 = "x2 + 3 * x3 = 1"), + paste("x", 1:3, sep = "")) > > stopifnot(max(abs(tmp$K - rbind(c(1, -1, 0), c(0, 1, 3)))) < sqrt(.Machine$double.eps)) > stopifnot(max(abs(tmp$m - c(2, 1))) < sqrt(.Machine$double.eps)) > > ### coef.survreg and vcov.survreg need special tuning > ### thx to Z for pointing this out > if (require("survival")) { + smod <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, + data = ovarian, dist = 'weibull') + K <- diag(length(coef(smod))) + rownames(K) <- names(coef(smod)) + glht(smod, linfct = K) + } General Linear Hypotheses Linear Hypotheses: Estimate (Intercept) == 0 6.8967 ecog.ps == 0 -0.3850 rx == 0 0.5286 > > ### new `means' comparisons > amod <- aov(weight ~ dose + gesttime + number, data = litter) > confint(glht(amod, linfct = mcp(dose = "Means"))) Simultaneous Confidence Intervals Multiple Comparisons of Means: Mean Contrasts Fit: aov(formula = weight ~ dose + gesttime + number, data = litter) Quantile = 2.5558 95% family-wise confidence level Linear Hypotheses: Estimate lwr upr 0 == 0 32.3651 30.0805 34.6498 5 == 0 29.0127 26.6372 31.3883 50 == 0 30.0743 27.5239 32.6246 500 == 0 29.6899 27.1591 32.2207 > > proc.time() user system elapsed 1.28 0.18 1.45