options(digits = 4) library("multcomp") RNGversion("3.5.2") set.seed(290875) ### mcp didn't accept objects of class `contrMat' ### spotted by Yves Brostaux amod <- aov(response ~ trt, data = cholesterol) cht1 <- glht(amod, linfct = mcp(trt = "Tukey")) K <- contrMat(table(cholesterol$trt), type = "Tukey") cht2 <- glht(amod, linfct = mcp(trt = K)) stopifnot(all.equal(coef(cht1), coef(cht2))) ### several inconsistencies spotted by ### Rich Heiberger 2006-11-28 ### need to be identical stopifnot(identical(cht1, print(cht1))) ### was: error summary(cht1)$test$coefficients ### NAs in coefficients tmp.data <- data.frame(EE=gl(2, 1, 24, letters[1:2]), FF=gl(3, 2, 24, LETTERS[3:5]), GG=gl(4, 6, 24, letters[6:9])) tmp.data$x <- rep(12, 24) tmp.data$y <- rep(7, 24) tmp.data$z <- c(9, 14, 3, 4, 15, 1, 11, 13, 24, 10, 22, 18, 20, 21, 6, 7, 16, 2, 19, 12, 17, 8, 23, 5) tmp.data$w <- c(15, 9, 18, 21, 17, 11, 23, 12, 1, 10, 2, 14, 24, 7, 13, 4, 5, 19, 16, 20, 3, 8, 22, 6) tmp.aov <- aov(z ~ EE+FF*GG + x*y +x*EE + y*FF, data=tmp.data) try(glht(tmp.aov, linfct=mcp(EE="Tukey"))) try(glht(tmp.aov, linfct=mcp(FF="Tukey"))) glht(tmp.aov, linfct=mcp(GG="Tukey")) ### covariate interactions: fire a warning tmp.aov <- aov(z ~ w*GG , data=tmp.data) glht(tmp.aov, linfct = mcp(GG = "Tukey")) ### stop with informative error message amod <- aov(breaks ~ tension + Error(wool), data = warpbreaks) try(glht(amod, linfct = mcp(tension = "Tukey"))) ### print error, spotted by Rich amod <- aov(breaks ~ wool * tension, data = warpbreaks) wht <- glht(amod, linfct = mcp(tension = "Tukey")) tmp <- confint(wht, calpha=2) print(tmp) ### coef. and vcov. didn't pass through ### bug report by John Deke lmod <- lm(Fertility ~ ., data = swiss) my.model <- list(coef(lmod),vcov(lmod)) coef2 <- function(model) return(model[[1]]) vcov2 <- function(model) return(model[[2]]) a <- glht(model = my.model, linfct = c("Agriculture=0","Catholic=0"), coef. = coef2, vcov. = vcov2, df = 100) b <- glht(model = lmod, linfct = c("Agriculture=0","Catholic=0"), df = 100) stopifnot(all.equal(coef(a), coef(b))) ### checks in mcp (spotted by Rich) amod <- aov(breaks ~ tension, data = warpbreaks) try(glht(amod, linfct = mcp(group = "Tukey"))) tmp <- warpbreaks class(tmp$tension) <- "numeric" amod <- aov(breaks ~ tension, data = tmp) try(glht(amod, linfct = mcp(tension = "Tukey"))) ### symbolic description and interactions ### spotted by Antonio Fabio Di Narzo dat <- data.frame(y = rnorm(6), x = seq_len(6), f = gl(2, 3)) lf <- glht(lm(y ~ x * f, data = dat), 'x + x:f2 = 0')$linfct stopifnot(all.equal(max(abs(lf - c(0, 1, 0, 1))), 0)) lf <- glht(lm(y ~ x * f, data = dat), 'x + 2.5 * x:f2 = 0')$linfct stopifnot(all.equal(max(abs(lf - c(0, 1, 0, 2.5))), 0)) ### example from Bretz 2001 JSCS `tmp` <- structure(list(gr = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", "2", "3"), class = "factor"), age = c(39L, 40L, 41L, 41L, 45L, 49L, 52L, 47L, 61L, 65L, 58L, 59L, 29L, 29L, 33L, 32L, 31L, 29L, 29L, 30L, 21L, 28L, 23L, 35L, 38L, 38L, 43L, 39L, 38L, 42L, 43L, 43L, 37L, 50L, 50L, 45L, 48L, 51L, 46L, 58L, 27L, 25L, 24L, 32L, 23L, 25L, 32L, 18L, 19L, 26L, 33L, 27L, 33L, 25L, 42L, 35L, 35L, 41L, 38L, 41L, 36L, 36L, 41L, 41L, 37L, 42L, 39L, 41L, 43L, 41L, 48L, 47L, 53L, 49L, 54L, 48L, 49L, 47L, 52L, 58L, 62L, 65L, 62L, 59L), y = c(4.62, 5.29, 5.52, 3.71, 4.02, 5.09, 2.7, 4.31, 2.7, 3.03, 2.73, 3.67, 5.21, 5.17, 4.88, 4.5, 4.47, 5.12, 4.51, 4.85, 5.22, 4.62, 5.07, 3.64, 3.64, 5.09, 4.61, 4.73, 4.58, 5.12, 3.89, 4.62, 4.3, 2.7, 3.5, 5.06, 4.06, 4.51, 4.66, 2.88, 5.29, 3.67, 5.82, 4.77, 5.71, 4.47, 4.55, 4.61, 5.86, 5.2, 4.44, 5.52, 4.97, 4.99, 4.89, 4.09, 4.24, 3.88, 4.85, 4.79, 4.36, 4.02, 3.77, 4.22, 4.94, 4.04, 4.51, 4.06, 4.02, 4.99, 3.86, 4.68, 4.74, 3.76, 3.98, 5, 3.31, 3.11, 4.76, 3.95, 4.6, 4.83, 3.18, 3.03)), .Names = c("gr", "age", "y"), row.names = c(NA, -84L), class = "data.frame") amod <- aov(y ~ gr + age, data = tmp) glht(amod, linfct = mcp(gr = "Tukey")) ### better error message ### suggested by Rich amod <- aov(breaks ~ tension, data = warpbreaks) try(glht(amod, linfct = mcp(tension = "Warp"))) ### cld did not find a terms component ### spotted by Peter B. Mandeville if (require("nlme")) { data("Orthodont") fm1 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1) hsd1 <- glht(fm1, linfct = mcp(Sex = "Tukey")) cld(hsd1) } ### spotted by ### example code by Achim Zeileis ## various models with and without intercept m1a <- lm(breaks ~ tension, data = warpbreaks) m1b <- lm(breaks ~ 0 + tension, data = warpbreaks) m2a <- lm(breaks ~ wool + tension, data = warpbreaks) m2b <- lm(breaks ~ 0 + wool + tension, data = warpbreaks) ## these two are equivalent: one factor with/without intercept stopifnot(all.equal( coef(glht(m1a, linfct = mcp(tension = "Tukey"))), coef(glht(m1b, linfct = mcp(tension = "Tukey"))))) ## these two should be equivalent: two factors with/without intercept ## but the latter fails stopifnot(all.equal( coef(glht(m2a, linfct = mcp(tension = "Tukey"))), coef(glht(m2b, linfct = mcp(tension = "Tukey"))))) library("MASS") xdf <- data.frame(y = gl(3, 10, ordered = TRUE), grp = sample(gl(3, 10))) glht(polr(y ~ grp, data = xdf), mcp(grp = "Dunnett")) ### interactions of two factors dat <- expand.grid(f = gl(2, 3), f2 = gl(3, 2)) dat$y <- rnorm(nrow(dat)) lf <- glht(lm(y ~ f : f2 - 1, data = dat), 'f1:f21 - f2:f22 = 0')$linfct stopifnot(all.equal(max(abs(lf - c(1, 0, 0, -1, 0, 0))), 0)) ### plotting one-sided confidence intervals amod <- aov(breaks ~ wool + tension, data = warpbreaks) wht <- glht(amod, linfct = mcp(tension = "Tukey"), alternative="greater") plot(wht, xlim=c(-30, 30), main="right side was missing") wht <- glht(amod, linfct = mcp(tension = "Tukey"), alternative="less") plot(wht, xlim=c(-40, 20), main="left side was missing") ### reported by Christian Ritz summary(glht(parm(1:4,matrix(c(1,0.97,0.89,0.74, 0.97,1,0.97,0.89, 0.89,0.97,1,0.97, 0.74,0.89,0.97,1), 4, 4)))) ### reported by Melissa Chester Key (Apr 22, 2016) set.seed(2343) X <- data.frame(X1 = rep(c(1,0),c(20,30)), X2 = rep(rep(c(1,0),3),c(rep(10,4),0,10)), X3 = rep(rep(c(1,0),5),each=5)) Y <- rnorm(50,4 + 4*X[,1] + 4*X[,2] + X[,3] + .5*X[,1]*X[,3] + .4*X[,2]*X[,3],.25) model <- lm(Y ~ (X1 + X2) * X3,data=X) coef(model) my.contrasts<- c( "X1 - X2 + .5*X1:X3 - .5*X2:X3 = 0", # previously wrong answer (actually got X1 + X2 + 0.5* X1) "X1 + .5*X1:X3 - X2 - .5*X2:X3 = 0", # previously wrong answer "X1 + .5*X1:X3 - .5*X2:X3 - X2 = 0") # right answer (contrast.result <- glht(model,lin = my.contrasts)) # right calculation (ok <- sum(coef(model) * c(0,1,-1,0,.5,-.5))) stopifnot(all.equal(as.numeric(coef(contrast.result)), rep(sum(coef(model) * c(0,1,-1,0,.5,-.5)),3))) # actual calculation - note that -1 has changed to 1 #sum(coef(model) * c(0, 1, 1, 0, .5, -.5)) (mc <- multcomp:::chrlinfct2matrix(my.contrasts, names(coef(model)))) stopifnot(all.equal(as.numeric(mc$K[1,c('(Intercept)', 'X1','X2', 'X3','X1:X3','X2:X3')]),c( 0,1,-1 ,0, 0.5,-0.5))) stopifnot(all.equal(as.numeric(mc$K[2,c('(Intercept)', 'X1','X2', 'X3','X1:X3','X2:X3')]),c( 0,1,-1 ,0, 0.5,-0.5))) stopifnot(all.equal(as.numeric(mc$K[3,c('(Intercept)', 'X1','X2', 'X3','X1:X3','X2:X3')]),c( 0,1,-1 ,0, 0.5,-0.5))) ### "(Intercept)" in char exprs for linfct x <- runif(100) y <- rnorm(length(x)) m <- lm(y ~ x) stopifnot(all.equal(coef(glht(m, linfct = "(Intercept) = 0")), coef(m)["(Intercept)"], check.attrributes = FALSE)) stopifnot(all.equal(coef(glht(m, linfct = "(Intercept) + x = 0")), sum(coef(m)), check.attributes = FALSE))