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().

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

Arguments

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.

Value

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

See also

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

Examples

library(doBy) 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
# Single factor. m0 <- lm(weight ~ feed, data = chickwts) anova(m0)
#> Analysis of Variance Table #> #> Response: weight #> Df Sum Sq Mean Sq F value Pr(>F) #> feed 5 231129 46226 15.365 5.936e-10 *** #> Residuals 65 195556 3009 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 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")
#> feed fit lwr upr cld #> 1 casein 323.5833 291.9608 355.2058 ab #> 2 horsebean 160.2000 125.5593 194.8407 e #> 3 linseed 218.7500 187.1275 250.3725 d #> 4 meatmeal 276.9091 243.8805 309.9377 bc #> 5 soybean 246.4286 217.1518 275.7053 cd #> 6 sunflower 328.9167 297.2942 360.5392 a
data(warpbreaks) # Two factors (complete factorial). m1 <- lm(breaks ~ wool * tension, data = warpbreaks) anova(m1)
#> Analysis of Variance Table #> #> Response: breaks #> Df Sum Sq Mean Sq F value Pr(>F) #> wool 1 450.7 450.67 3.7653 0.0582130 . #> tension 2 2034.3 1017.13 8.4980 0.0006926 *** #> wool:tension 2 1002.8 501.39 4.1891 0.0210442 * #> Residuals 48 5745.1 119.69 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L <- LE_matrix(m1, effect = c("wool", "tension")) attributes(L)
#> $dim #> [1] 6 6 #> #> $dimnames #> $dimnames[[1]] #> NULL #> #> $dimnames[[2]] #> [1] "(Intercept)" "woolB" "tensionM" "tensionH" #> [5] "woolB:tensionM" "woolB:tensionH" #> #> #> $grid #> wool tension #> 1 A L #> 2 B L #> 3 A M #> 4 B M #> 5 A H #> 6 B H #> #> $class #> [1] "LE_matrix_class" "matrix" #>
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)
#> $H #> wool fit lwr upr cld #> 1 A 24.55556 18.43912 30.67199 a #> 2 B 18.77778 12.66134 24.89421 a #> #> $L #> wool fit lwr upr cld #> 1 A 44.55556 38.43912 50.67199 a #> 2 B 28.22222 22.10579 34.33866 b #> #> $M #> wool fit lwr upr cld #> 1 A 24.00000 17.88356 30.11644 a #> 2 B 28.77778 22.66134 34.89421 a #>
# Two factors (incomplete factorial). warpbreaks <- subset(warpbreaks, !(tension == "H" & wool == "A")) xtabs(~tension + wool, data = warpbreaks)
#> wool #> tension A B #> L 9 9 #> M 9 9 #> H 0 9
# There is NA in the estimated parameters. m2 <- lm(breaks ~ wool * tension, data = warpbreaks) coef(m2)
#> (Intercept) woolB tensionM tensionH woolB:tensionM #> 44.555556 -16.333333 -20.555556 -9.444444 21.111111 #> woolB:tensionH #> NA
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)
#> Analysis of Variance Table #> #> Model 1: breaks ~ wool * tension #> Model 2: breaks ~ X - 1 #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 40 4900.9 #> 2 40 4900.9 0 0
# 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)
#> $A #> tension fit lwr upr cld #> 1 L 44.55556 38.34272 50.76839 b #> 2 M 24.00000 17.78716 30.21284 a #> #> $B #> tension fit lwr upr cld #> 1 L 28.22222 22.00939 34.43506 a #> 2 M 28.77778 22.56494 34.99061 a #> 3 H 18.77778 12.56494 24.99061 a #>
# Sample means. aggregate(breaks ~ tension + wool, data = warpbreaks, FUN = mean)
#> tension wool breaks #> 1 L A 44.55556 #> 2 M A 24.00000 #> 3 L B 28.22222 #> 4 M B 28.77778 #> 5 H B 18.77778