Function Reference: multcompare
- statistics: C = multcompare (STATS)
- statistics: C = multcompare (STATS, "name", value)
- statistics: [C, M] = multcompare (...)
- statistics: [C, M, H] = multcompare (...)
- statistics: [C, M, H, GNAMES] = multcompare (...)
- statistics: padj = multcompare (p)
- statistics: padj = multcompare (p, "ctype", CTYPE)
Perform posthoc multiple comparison tests or p-value adjustments to control the family-wise error rate (FWER) or false discovery rate (FDR).
C = multcompare (STATS) performs a multiple comparison
using a STATS structure that is obtained as output from any of
the following functions: anova1, anova2, anovan, kruskalwallis, and friedman.
The return value C is a matrix with one row per comparison and six
columns. Columns 1-2 are the indices of the two samples being compared.
Columns 3-5 are a lower bound, estimate, and upper bound for their
difference, where the bounds are for 95% confidence intervals. Column 6-8 are
the multiplicity adjusted p-values for each individual comparison, the test
statistic and the degrees of freedom.
All tests by multcompare are two-tailed.
multcompare can take a number of optional parameters as name-value
pairs.
[…] = multcompare (STATS, "alpha", ALPHA)
- ALPHA sets the significance level of null hypothesis significance tests to ALPHA, and the central coverage of two-sided confidence intervals to 100*(1-ALPHA)%. (Default ALPHA is 0.05).
[…] = multcompare (STATS, "ControlGroup", REF)
- REF is the index of the control group to limit comparisons to. The index must be a positive integer scalar value. For each dimension (d) listed in DIM, multcompare uses STATS.grpnames{d}(idx) as the control group. (Default is empty, i.e. [], for full pairwise comparisons)
[…] = multcompare (STATS, "ctype", CTYPE)
-
CTYPE is the type of comparison test to use. In order of increasing
power, the choices are: "bonferroni", "scheffe", "mvt", "holm" (default),
"hochberg", "fdr", or "lsd". The first five methods control the family-wise
error rate. The "fdr" method controls false discovery rate (by the original
Benjamini-Hochberg step-up procedure). The final method, "lsd" (or "none"),
makes no attempt to control the Type 1 error rate of multiple comparisons.
The coverage of confidence intervals are only corrected for multiple
comparisons in the cases where CTYPE is "bonferroni", "scheffe" or
"mvt", which control the Type 1 error rate for simultaneous inference.
The "mvt" method uses the multivariate t distribution to assess the probability or critical value of the maximum statistic across the tests, thereby accounting for correlations among comparisons in the control of the family-wise error rate with simultaneous inference. In the case of pairwise comparisons, it simulates Tukey’s (or the Games-Howell) test, in the case of comparisons with a single control group, it simulates Dunnett’s test. CTYPE values "tukey-kramer" and "hsd" are recognised but set the value of CTYPE and REF to "mvt" and empty respectively. A CTYPE value "dunnett" is recognised but sets the value of CTYPE to "mvt", and if REF is empty, sets REF to 1. Since the algorithm uses a Monte Carlo method (of 1e+06 random samples), you can expect the results to fluctuate slightly with each call to multcompare and the calculations may be slow to complete for a large number of comparisons. If the parallel package is installed and loaded,
multcomparewill automatically accelerate computations by parallel processing. Note that p-values calculated by the "mvt" are truncated at 1e-06.
[…] = multcompare (STATS, "df", DF)
-
DF is an optional scalar value to set the number of degrees of freedom
in the calculation of p-values for the multiple comparison tests. By default,
this value is extracted from the STATS structure of the ANOVA test, but
setting DF maybe necessary to approximate Satterthwaite correction if
anovanwas performed using weights.
[…] = multcompare (STATS, "dim", DIM)
-
DIM is a vector specifying the dimension or dimensions over which the
estimated marginal means are to be calculated. Used only if STATS comes from
anovan. The value [1 3], for example, computes the estimated marginal mean
for each combination of the first and third predictor values. The default is
to compute over the first dimension (i.e. 1). If the specified dimension is,
or includes, a continuous factor then
multcomparewill return an error.
[…] = multcompare (STATS, "estimate", ESTIMATE)
- ESTIMATE is a string specifying the estimates to be compared when computing multiple comparisons after anova2; this argument is ignored by anovan and anova1. Accepted values for ESTIMATE are either "column" (default) to compare column means, or "row" to compare row means. If the model type in anova2 was "linear" or "nested" then only "column" is accepted for ESTIMATE since the row factor is assumed to be a random effect.
[…] = multcompare (STATS, "display", DISPLAY)
- DISPLAY is either "on" (the default): to display a table and graph of the comparisons (e.g. difference between means), their 100*(1-ALPHA)% intervals and multiplicity adjusted p-values in APA style; or "off": to omit the table and graph. On the graph, markers and error bars colored red have multiplicity adjusted p-values < ALPHA, otherwise the markers and error bars are blue.
[…] = multcompare (STATS, "seed", SEED)
- SEED is a scalar value used to initialize the random number generator so that CTYPE "mvt" produces reproducible results.
[C, M, H, GNAMES] = multcompare (…)
returns additional outputs. M is a matrix where columns 1-2 are the
estimated marginal means and their standard errors, and columns 3-4 are lower
and upper bounds of the confidence intervals for the means; the critical
value of the test statistic is scaled by a factor of 2^(-0.5) before
multiplying by the standard errors of the group means so that the intervals
overlap when the difference in means becomes significant at approximately
the level ALPHA. When ALPHA is 0.05, this corresponds to
confidence intervals with 83.4% central coverage. H is a handle to the
figure containing the graph. GNAMES is a cell array with one row for
each group, containing the names of the groups.
padj = multcompare (p) calculates and returns adjusted
p-values (padj) using the Holm-step down Bonferroni procedure to
control the family-wise error rate.
padj = multcompare (p, "ctype", CTYPE) calculates
and returns adjusted p-values (padj) computed using the method
CTYPE. In order of increasing power, CTYPE for p-value adjustment
can be either "bonferroni", "holm" (default), "hochberg", or "fdr". See
above for further information about the CTYPE methods.
See also: anova1, anova2, anovan, kruskalwallis, friedman, fitlm
Source Code: multcompare
Example: 1
## Demonstration using balanced one-way ANOVA from anova1
x = ones (50, 4) .* [-2, 0, 1, 5];
randn ("seed", 1); # for reproducibility
x = x + normrnd (0, 2, 50, 4);
groups = {"A", "B", "C", "D"};
[p, tbl, stats] = anova1 (x, groups, "off");
multcompare (stats);
HOLM Multiple Comparison (Post Hoc) Test for ANOVA1
Group ID Group ID LBoundDiff EstimatedDiff UBoundDiff p-value
----------------------------------------------------------------------
1 2 -2.789 -2.009 -1.228 <.001
1 3 -2.489 -1.709 -0.928 <.001
1 4 -7.365 -6.585 -5.804 <.001
2 3 -0.480 0.300 1.081 .449
2 4 -5.356 -4.576 -3.795 <.001
3 4 -5.656 -4.876 -4.095 <.001
|
Example: 2
## Demonstration using unbalanced one-way ANOVA example from anovan
dv = [ 8.706 10.362 11.552 6.941 10.983 10.092 6.421 14.943 15.931 ...
22.968 18.590 16.567 15.944 21.637 14.492 17.965 18.851 22.891 ...
22.028 16.884 17.252 18.325 25.435 19.141 21.238 22.196 18.038 ...
22.628 31.163 26.053 24.419 32.145 28.966 30.207 29.142 33.212 ...
25.694 ]';
g = [1 1 1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 3 3 3 ...
4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5]';
[P,ATAB, STATS] = anovan (dv, g, "varnames", "score", "display", "off");
[C, M, H, GNAMES] = multcompare (STATS, "dim", 1, "ctype", "holm", ...
"ControlGroup", 1, "display", "on")
HOLM Multiple Comparison (Post Hoc) Test for ANOVAN
Group ID Group ID LBoundDiff EstimatedDiff UBoundDiff p-value
----------------------------------------------------------------------
1 2 -11.343 -8.000 -4.657 <.001
1 3 -11.932 -9.000 -6.068 <.001
1 4 -14.035 -11.000 -7.965 <.001
1 5 -21.849 -19.000 -16.151 <.001
C =
Columns 1 through 6:
1.0000e+00 2.0000e+00 -1.1343e+01 -8.0000e+00 -4.6572e+00 2.8581e-05
1.0000e+00 3.0000e+00 -1.1932e+01 -9.0000e+00 -6.0682e+00 1.0459e-06
1.0000e+00 4.0000e+00 -1.4035e+01 -1.1000e+01 -7.9654e+00 6.3838e-08
1.0000e+00 5.0000e+00 -2.1849e+01 -1.9000e+01 -1.6151e+01 3.0198e-14
Columns 7 and 8:
-4.8748e+00 3.2000e+01
-6.2529e+00 3.2000e+01
-7.3834e+00 3.2000e+01
-1.3583e+01 3.2000e+01
M =
10.0000 1.0178 8.5341 11.4659
18.0000 1.2874 16.1458 19.8542
19.0000 1.0178 17.5341 20.4659
21.0001 1.0880 19.4330 22.5673
29.0001 0.9595 27.6180 30.3822
H = 1
GNAMES =
{
[1,1] = score=1
[2,1] = score=2
[3,1] = score=3
[4,1] = score=4
[5,1] = score=5
}
|
Example: 3
## Demonstration using factorial ANCOVA example from anovan
score = [95.6 82.2 97.2 96.4 81.4 83.6 89.4 83.8 83.3 85.7 ...
97.2 78.2 78.9 91.8 86.9 84.1 88.6 89.8 87.3 85.4 ...
81.8 65.8 68.1 70.0 69.9 75.1 72.3 70.9 71.5 72.5 ...
84.9 96.1 94.6 82.5 90.7 87.0 86.8 93.3 87.6 92.4 ...
100. 80.5 92.9 84.0 88.4 91.1 85.7 91.3 92.3 87.9 ...
91.7 88.6 75.8 75.7 75.3 82.4 80.1 86.0 81.8 82.5]';
treatment = {"yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" ...
"yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" ...
"yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" "yes" ...
"no" "no" "no" "no" "no" "no" "no" "no" "no" "no" ...
"no" "no" "no" "no" "no" "no" "no" "no" "no" "no" ...
"no" "no" "no" "no" "no" "no" "no" "no" "no" "no"}';
exercise = {"lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" ...
"mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" ...
"hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi" ...
"lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" "lo" ...
"mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" "mid" ...
"hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi" "hi"}';
age = [59 65 70 66 61 65 57 61 58 55 62 61 60 59 55 57 60 63 62 57 ...
58 56 57 59 59 60 55 53 55 58 68 62 61 54 59 63 60 67 60 67 ...
75 54 57 62 65 60 58 61 65 57 56 58 58 58 52 53 60 62 61 61]';
[P, ATAB, STATS] = anovan (score, {treatment, exercise, age}, "model", ...
[1 0 0; 0 1 0; 0 0 1; 1 1 0], "continuous", 3, ...
"sstype", "h", "display", "off", "contrasts", ...
{"simple","poly",""});
[C, M, H, GNAMES] = multcompare (STATS, "dim", [1 2], "ctype", "holm", ...
"display", "on")
HOLM Multiple Comparison (Post Hoc) Test for ANOVAN
Group ID Group ID LBoundDiff EstimatedDiff UBoundDiff p-value
----------------------------------------------------------------------
1 2 -4.544 -0.017 4.509 1.000
1 3 8.963 13.703 18.443 <.001
1 4 -6.002 -1.529 2.945 1.000
1 5 -6.174 -1.701 2.771 1.000
1 6 -0.692 3.957 8.605 .662
2 3 9.166 13.721 18.276 <.001
2 4 -6.060 -1.511 3.038 1.000
2 5 -6.195 -1.684 2.828 1.000
2 6 -0.533 3.974 8.481 .662
3 4 -20.018 -15.232 -10.446 <.001
3 5 -20.112 -15.404 -10.697 <.001
3 6 -14.228 -9.747 -5.265 <.001
4 5 -4.650 -0.172 4.305 1.000
4 6 0.798 5.485 10.172 .204
5 6 1.035 5.658 10.280 .174
C =
Columns 1 through 6:
1.0000e+00 2.0000e+00 -4.5437e+00 -1.7457e-02 4.5087e+00 1.0000e+00
1.0000e+00 3.0000e+00 8.9634e+00 1.3703e+01 1.8443e+01 4.5333e-06
1.0000e+00 4.0000e+00 -6.0019e+00 -1.5286e+00 2.9447e+00 1.0000e+00
1.0000e+00 5.0000e+00 -6.1735e+00 -1.7011e+00 2.7714e+00 1.0000e+00
1.0000e+00 6.0000e+00 -6.9212e-01 3.9565e+00 8.6051e+00 6.6201e-01
2.0000e+00 3.0000e+00 9.1656e+00 1.3721e+01 1.8276e+01 2.0212e-06
2.0000e+00 4.0000e+00 -6.0600e+00 -1.5111e+00 3.0378e+00 1.0000e+00
2.0000e+00 5.0000e+00 -6.1953e+00 -1.6836e+00 2.8281e+00 1.0000e+00
2.0000e+00 6.0000e+00 -5.3340e-01 3.9740e+00 8.4813e+00 6.6201e-01
3.0000e+00 4.0000e+00 -2.0018e+01 -1.5232e+01 -1.0446e+01 6.1827e-07
3.0000e+00 5.0000e+00 -2.0112e+01 -1.5404e+01 -1.0697e+01 3.4083e-07
3.0000e+00 6.0000e+00 -1.4228e+01 -9.7468e+00 -5.2654e+00 6.5704e-04
4.0000e+00 5.0000e+00 -4.6499e+00 -1.7249e-01 4.3050e+00 1.0000e+00
4.0000e+00 6.0000e+00 7.9808e-01 5.4851e+00 1.0172e+01 2.0411e-01
5.0000e+00 6.0000e+00 1.0354e+00 5.6576e+00 1.0280e+01 1.7402e-01
Columns 7 and 8:
-7.7360e-03 5.3000e+01
5.7988e+00 5.3000e+01
-6.8538e-01 5.3000e+01
-7.6287e-01 5.3000e+01
1.7071e+00 5.3000e+01
6.0416e+00 5.3000e+01
-6.6630e-01 5.3000e+01
-7.4847e-01 5.3000e+01
1.7684e+00 5.3000e+01
-6.3835e+00 5.3000e+01
-6.5634e+00 5.3000e+01
-4.3623e+00 5.3000e+01
-7.7268e-02 5.3000e+01
2.3473e+00 5.3000e+01
2.4550e+00 5.3000e+01
M =
86.9788 1.6031 84.7051 89.2525
86.9962 1.5774 84.7590 89.2334
73.2755 1.6514 70.9334 75.6176
88.5074 1.6166 86.2146 90.8002
88.6799 1.5948 86.4180 90.9417
83.0223 1.6130 80.7346 85.3100
H = 1
GNAMES =
{
[1,1] = X1=yes, X2=lo
[2,1] = X1=yes, X2=mid
[3,1] = X1=yes, X2=hi
[4,1] = X1=no, X2=lo
[5,1] = X1=no, X2=mid
[6,1] = X1=no, X2=hi
}
|
Example: 4
## Demonstration using one-way ANOVA from anovan, with fit by weighted least
## squares to account for heteroskedasticity.
g = [1, 1, 1, 1, 1, 1, 1, 1, ...
2, 2, 2, 2, 2, 2, 2, 2, ...
3, 3, 3, 3, 3, 3, 3, 3]';
y = [13, 16, 16, 7, 11, 5, 1, 9, ...
10, 25, 66, 43, 47, 56, 6, 39, ...
11, 39, 26, 35, 25, 14, 24, 17]';
[P,ATAB,STATS] = anovan(y, g, "display", "off");
fitted = STATS.X * STATS.coeffs(:,1); # fitted values
b = polyfit (fitted, abs (STATS.resid), 1);
v = polyval (b, fitted); # Variance as a function of the fitted values
[P,ATAB,STATS] = anovan (y, g, "weights", v.^-1, "display", "off");
[C, M] = multcompare (STATS, "display", "on", "ctype", "mvt")
MVT Multiple Comparison (Post Hoc) Test for ANOVAN
Group ID Group ID LBoundDiff EstimatedDiff UBoundDiff p-value
----------------------------------------------------------------------
1 2 -42.096 -26.750 -11.404 <.001
1 3 -26.805 -14.125 -1.445 .027
2 3 -5.239 12.625 30.489 .198
C =
Columns 1 through 6:
1.0000e+00 2.0000e+00 -4.2096e+01 -2.6750e+01 -1.1404e+01 7.0100e-04
1.0000e+00 3.0000e+00 -2.6805e+01 -1.4125e+01 -1.4448e+00 2.7414e-02
2.0000e+00 3.0000e+00 -5.2387e+00 1.2625e+01 3.0489e+01 1.9829e-01
Columns 7 and 8:
-4.3716e+00 2.1000e+01
-2.7937e+00 2.1000e+01
1.7725e+00 2.1000e+01
M =
9.7500 2.4770 5.3572 14.1428
36.5000 5.5953 26.5773 46.4227
23.8750 4.4076 16.0585 31.6915
|
Example: 5
## Demonstration of p-value adjustments to control the false discovery rate
## Data from Westfall (1997) JASA. 92(437):299-306
p = [.005708; .023544; .024193; .044895; ...
.048805; .221227; .395867; .693051; .775755];
padj = multcompare(p,'ctype','fdr')
padj =
0.051372
0.072579
0.072579
0.087849
0.087849
0.331840
0.508972
0.775755
0.775755
|