The margins command introduced in Stata 11 provides an easy way to compute contrasts in analysis of variance designs. We will begin by loading the hsbanova dataset and running a factorial anova with read as the response variable.
use https://stats.idre.ucla.edu/stat/data/hsbanova, clear anova read grp##female Number of obs = 200 R-squared = 0.2104 Root MSE = 9.27502 Adj R-squared = 0.1817 Source | Partial SS df MS F Prob > F -----------+---------------------------------------------------- Model | 4402.42 7 628.917143 7.31 0.0000 | grp | 3909.6042 3 1303.2014 15.15 0.0000 female | 117.305356 1 117.305356 1.36 0.2444 grp#female | 427.282493 3 142.427498 1.66 0.1780 | Residual | 16517 192 86.0260417 -----------+---------------------------------------------------- Total | 20919.42 199 105.122714
Neither the grp by female interaction nor the female main effect were significant. Prior to the experiment we decided that four comparisons among the grp means would be theoretically interesting and that we would test these if grp was statistically significant. It was and so we will test the following contrasts:
-
- grp1 versus grp2
-
- grp3 versus grp4
-
- grp1 & grp2 versus grp3 & grp4
- grp1, grp2 & grp3 versus grp4
It can also be useful to look at these contrasts in terms of the weights applied to each of the means. The four contrasts above would look something like this:
grp1 grp2 grp3 grp4 1 -1 0 0 contrast 1 0 0 1 -1 contrast 2 1 1 -1 -1 contrast 3 1/3 1/3 1/3 -1 contrast 4
Note that for each of the above contrasts the sum of the weights is zero.
We can now begin our multiple comparisons by running the margins command with the asbalanced and post options.
margins grp, asbalanced post Adjusted predictions Number of obs = 200 Expression : Linear prediction, predict() at : grp (asbalanced) female (asbalanced) ------------------------------------------------------------------------------ | Delta-method | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- grp | 1 | 46.15942 1.315904 35.08 0.000 43.5803 48.73854 2 | 49.95536 1.385722 36.05 0.000 47.23939 52.67132 3 | 54.99107 1.20007 45.82 0.000 52.63898 57.34317 4 | 57.86032 1.399677 41.34 0.000 55.11701 60.60364 ------------------------------------------------------------------------------
The margins command displays the adjusted marginal means (aka lsmeans, aka estimated marginal means) for each of the four levels of grp. We will now proceed to compute each of the one degree of freedom contrasts two ways; first using the test command and then using lincom. The test command produces a chi-square while the lincom command produces a z-test. Squaring the z yields the same value as the chi-square.
/* contrast 1: grp1 vs grp2 */ test 1.grp - 2.grp = 0 ( 1) 1bn.grp - 2.grp = 0 chi2( 1) = 3.95 Prob > chi2 = 0.0470 lincom 1.grp - 2.grp ( 1) 1bn.grp - 2.grp = 0 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -3.795937 1.910975 -1.99 0.047 -7.54138 -.0504938 ------------------------------------------------------------------------------ /* contrast 2: grp3 vs grp4 */ test 3.grp - 4.grp = 0 ( 1) 3.grp - 4.grp = 0 chi2( 1) = 2.42 Prob > chi2 = 0.1197 lincom 3.grp - 4.grp ( 1) 3.grp - 4.grp = 0 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -2.869252 1.843709 -1.56 0.120 -6.482856 .7443509 ------------------------------------------------------------------------------ /* contrast 3: (grp1 & grp2) vs (grp3 & grp4) */ test 1.grp + 2.grp - 3.grp - 4.grp = 0 ( 1) 1bn.grp + 2.grp - 3.grp - 4.grp = 0 chi2( 1) = 39.73 Prob > chi2 = 0.0000 lincom 1.grp + 2.grp - 3.grp - 4.grp ( 1) 1bn.grp + 2.grp - 3.grp - 4.grp = 0 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -16.73662 2.655389 -6.30 0.000 -21.94108 -11.53215 ------------------------------------------------------------------------------ /* contrast 4: (grp1, grp2 &grp3) vs (grp4) */ test 1/3*1.grp + 1/3*2.grp + 1/3*3.grp - 4.grp = 0 ( 1) .3333333*1bn.grp + .3333333*2.grp + .3333333*3.grp - 4.grp = 0 chi2( 1) = 22.23 Prob > chi2 = 0.0000 lincom 1/3*1.grp + 1/3*2.grp + 1/3*3.grp - 4.grp ( 1) .3333333*1bn.grp + .3333333*2.grp + .3333333*3.grp - 4.grp = 0 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -7.491708 1.588985 -4.71 0.000 -10.60606 -4.377355 ------------------------------------------------------------------------------
We did the last comparison using the weights (1/3 1/3 1/3 -1). We could get the same chi-square and z using the weights (1 1 1 -3). However, you will note that the coefficient in the lincom will have a different numerical value.
test 1.grp + 2.grp + 3.grp - 3*4.grp = 0 ( 1) 1bn.grp + 2.grp + 3.grp - 3*4.grp = 0 chi2( 1) = 22.23 Prob > chi2 = 0.0000 lincom 1.grp + 2.grp + 3.grp - 3*4.grp ( 1) 1bn.grp + 2.grp + 3.grp - 3*4.grp = 0 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | -22.47512 4.766954 -4.71 0.000 -31.81818 -13.13206 ------------------------------------------------------------------------------
At this point we should consider some kind of adjustment for having done four contrasts. The easiest approach would be a Bonferroni adjustment made by dividing our alpha level by the number of comparisons (.05/4 = .0125). Thus, any contrast less then .0125 would be significant, namely contrasts 3 and 4.
Alternatively, we might want to consider using a Scheffé adjustment since two of the contrasts (3 and 4) were non-pairwise contrasts. The critical value for Scheffé can be computed from the following:
(p-1)*critical_value_F(p-1, dfe) = (4-1)*critical_value_F(4-1, 192) = (4-1)*2.65 = 3*2.65 = 7.95
Although, the critical value for Scheffé is an F-ratio we can compare it to the chi-square with one degree of freedom from the test commands. Once again only contrasts 3 and 4 were statistically significant.