apmc: A Wraper of glht to Get All Pairwise Mean Comparisons
In walmes/wzRfun: Walmes Zeviani's collection of functions

Description Usage Arguments Value Author(s) See Also Examples

This function performs all pairwise compararisons among means returning pontual and intervalar estimates followed by letters to easy discriminate values. It is in fact a wraper of glht().

1	apmc(X, model, focus, test = "single-step", level = 0.05, cld2 = FALSE)

`X`	a matrix where each line is a linear function of the model parameters to estimate a least squares mean. In most pratical cases, it is an object from the `LE_matrix()`.
`model`	a model with class recognized by `glht()`.
`focus`	a string with the name of the factor which levels are being compared.
`test`	a p-value correction method. See `p.adjust.methods()`.
`level`	the experimentwise significance level for the multiple comparisons. The individual coverage of the confidence interval is `1-level`. Default is 0.05.
`cld2`	Logical, if `TRUE` uses the `cld2()` functions, otherwise uses the `cld()` function.

a data.frame with interval estimates and compact letter display for the means comparisons.

Walmes Zeviani, walmes@ufpr.br.

apc(), LE_matrix(), glht().

library(doBy)
library(multcomp)

# Single factor.
m0 <- lm(weight ~ feed, data = chickwts)
anova(m0)

# Prepare the matrix to estimate lsmeans.
L <- LE_matrix(m0, effect = "feed")
rownames(L) <- levels(chickwts$feed)
apmc(L, model = m0, focus = "feed", test = "fdr")

data(warpbreaks)

# Two factors (complete factorial).
m1 <- lm(breaks ~ wool * tension, data = warpbreaks)
anova(m1)

L <- LE_matrix(m1, effect = c("wool", "tension"))
attributes(L)

Ls <- by(L, INDICES = attr(L, "grid")$tension, FUN = as.matrix)
Ls <- lapply(Ls, "rownames<-", levels(warpbreaks$wool))

# Comparing means of wool in each tension.
lapply(Ls, apmc, model = m1, focus = "wool",
       test = "single-step", level = 0.1)

# Two factors (incomplete factorial).
warpbreaks <- subset(warpbreaks, !(tension == "H" & wool == "A"))
xtabs(~tension + wool, data = warpbreaks)

# There is NA in the estimated parameters.
m2 <- lm(breaks ~ wool * tension, data = warpbreaks)
coef(m2)

X <- model.matrix(m2)
b <- coef(m2)

X <- X[, !is.na(b)]

# unique(X)

# Uses the full estimable model matriz.
m3 <- update(m2, . ~ 0 + X)

# These models are in fact the same.
anova(m2, m3)

# LS matrix has all cells.
L <- LE_matrix(m2, effect = c("wool", "tension"))
g <- attr(L, "grid")
L <- L[, !is.na(b)]
i <- 5
L <- L[-i, ]
g <- g[-i, ]

rownames(L) <- make.names(g$tension, unique = FALSE)
Ls <- split.data.frame(L, g$wool)

# LSmeans with MCP test.
lapply(Ls, apmc, model = m3, focus = "tension",
       test = "single-step", level = 0.1, cld2 = TRUE)

# Sample means.
aggregate(breaks ~ tension + wool, data = warpbreaks, FUN = mean)