A Wraper of glht to Get All Pairwise Mean Comparisons

Source: R/pairwise.R

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