1
$\begingroup$

I would like to understand how to specify a contrast in a follow-up analysis using the afex package. Below I included code from the documentation and explain why I do not follow all steps here.

Edit: I was doubting whether to post here or on stackoverflow. I decided to post here because I am also interested in understanding differences in output for different tests (aside from understanding how to use R-code to specify contrasts in afex).

Code from the afex documentation:

    ########################### ## 2: Follow-up Analysis ## ###########################
    # use data as above 
    data(obk.long, package = "afex")
    # 1. obtain afex_aov object: 
    a1 <- aov_ez("id", "value", 
                 obk.long, 
                 between = c("treatment", "gender"), 
                 within = c("phase", "hour"), 
                 observed = "gender")
    # 2. obtain reference grid object (default is uses univariate model): 
    r1 <- emmeans(a1, ~treatment +phase, model = "multivariate") 
    r1
    # treatment phase   emmean        SE df lower.CL upper.CL
    # control   fup   4.333333 0.5511982 10 3.105187 5.561479
    # A         fup   7.250000 0.6038074 10 5.904633 8.595367
    # B         fup   7.291667 0.4611655 10 6.264126 8.319207
    # control   post  4.083333 0.6277716 10 2.684571 5.482096
    # A         post  6.500000 0.6876894 10 4.967733 8.032267
    # B         post  6.625000 0.5252314 10 5.454711 7.795289
    # control   pre   4.250000 0.7660323 10 2.543174 5.956826
    # A         pre   5.000000 0.8391464 10 3.130265 6.869735
    # B         pre   4.166667 0.6409086 10 2.738633 5.594700
    # 
    # Results are averaged over the levels of: gender, hour 
    # Confidence level used: 0.95 

Based on this reference grid, I would specify the contrast A_B_pre (A vs. B, pretest) as follows: c(rep(0, 6), 0, -1, 1). However, the documentation specifies the same contrast as A_B_pre = c(0, -1, 1, rep(0, 6)) (which to me seems to compare A vs. B in the follow-up treatment):

    # 3. create list of contrasts on the reference grid: 
    c1 <- list( A_B_pre = c(0, -1, 1, rep(0, 6)), 
                # A versus B for pretest 
                A_B_comb = c(0, 0, 0, 0, -0.5, 0.5, 0, -0.5, 0.5), 
                # A vs. B for post and follow-up combined 
                effect_post = c(0, 0, 0, -1, 0.5, 0.5, 0, 0, 0), 
                # control versus A&B post 
                effect_fup = c(0, 0, 0, 0, 0, 0, -1, 0.5, 0.5), 
                # control versus A&B follow-up 
                effect_comb = c(0, 0, 0, -0.5, 0.25, 0.25, -0.5, 0.25, 0.25) 
                # control versus A&B combined
    )
    # 4. test contrasts on reference grid: 
    contrast(r1, c1)
    # contrast       estimate        SE df t.ratio p.value
    # A_B_pre      0.04166667 0.7597743 10   0.055  0.9573
    # A_B_comb    -0.35416667 0.8327863 10  -0.425  0.6796
    # effect_post  2.47916667 0.7624260 10   3.252  0.0087
    # effect_fup   0.33333333 0.9303431 10   0.358  0.7276
    # effect_comb  1.40625000 0.7337579 10   1.917  0.0843
    # 
    # Results are averaged over the levels of: gender, hour 

For the comparison of treatment A vs. B in the phase pre, the estimate of the difference (based on the contrast from the documentation) is 0.04166667. However, when I calculate that difference, I obtain 0.8571429: 1)

    mean(obk.long[obk.long$treatment=="A" & obk.long$phase=="pre",]$value, na.rm=TRUE) -
    mean(obk.long[obk.long$treatment=="B" & obk.long$phase=="pre",]$value, na.rm=TRUE)
    # [1] 0.8571429

As I mentioned above, to me it seems that the contrast from the documentation compares A vs. B in the follow-up treatment. Indeed, when I manually calculate this comparison in the follow-up treatment, I come closer to the output from afex, although I am still far away from an exact match:

    mean(obk.long[obk.long$treatment=="A" & obk.long$phase=="fup",]$value, na.rm=TRUE) -
            mean(obk.long[obk.long$treatment=="B" & obk.long$phase=="fup",]$value, na.rm=TRUE)
    #[1] -0.03571429

    t.test(obk.long[obk.long$treatment=="A" & obk.long$phase=="fup",]$value,
           obk.long[obk.long$treatment=="B" & obk.long$phase=="fup",]$value)
    #p-value = 0.9484

1)

Edit 2: When obtaining all pairwise comparisons (i.e., without using the contrasts specified in the documentation, I obtain 0.83333333, which comes relatively close to 0.8571429 (differences likely arise because I do not control for gender and hour here, see also below):

emmeans(a1, c("treatment", "phase"), contr = "pairwise", model = "multivariate")
#...
#$contrasts
#contrast                    estimate        SE   df t.ratio p.value
#...
#A,pre - B,pre               0.83333333 1.0559027 10   0.789  0.9946

Additional thoughts / bonus question:

Perhaps the difference arises because in the contrast for the afex package, hour and gender are controlled for? Usually I would regress the dependent variable (value) on these variables (hour and gender) and then use the residuals in the next analysis (i.e., the t-test), but in this case I receive only 16 residuals (16 residuals for each level of hour) instead of 240 residuals (1 for each observation). Is it possible to control for hour and gender in the t.test above?

    a2 <- aov_ez("id", "value",
                 obk.long,
                 between = c("gender"),
                 within = c("hour"),
                 observed = "gender")
    residuals(a2$lm)
    # X1         X2          X3          X4          X5
    # 1  -3.4583333 -3.4583333 -3.29166667 -2.95833333 -2.87500000
    # 2  -2.1250000 -2.1250000 -2.29166667 -1.95833333 -2.54166667
    # 3  -0.1250000 -0.1250000  0.04166667  0.37500000  0.45833333
    # 4  -0.8333333 -1.5416667 -1.37500000 -1.16666667 -0.70833333
    # 5  -0.1666667 -0.2083333  0.29166667 -0.83333333 -1.37500000
    # 6   2.8750000  3.2083333  2.37500000  2.70833333  2.79166667
    # 7   1.2083333  1.2083333  1.37500000  1.70833333  1.79166667
    # 8  -1.1666667 -0.5416667  0.29166667 -0.83333333 -0.04166667
    # 9  -0.5000000 -0.8750000 -0.70833333 -0.50000000 -1.37500000
    # 10  0.5416667  0.5416667 -0.29166667  0.04166667  1.45833333
    # 11 -1.1250000 -1.4583333 -0.62500000 -1.62500000 -1.20833333
    # 12  2.2083333  2.2083333  2.70833333  1.70833333  0.12500000
    # 13  0.8333333  1.1250000  0.29166667  3.16666667  2.29166667
    # 14 -0.1666667  0.1250000 -0.37500000 -1.50000000 -2.04166667
    # 15  0.5000000  0.1250000 -0.37500000  0.83333333  1.62500000
    # 16  1.5000000  1.7916667  1.95833333  0.83333333  1.62500000
$\endgroup$
  • 1
    $\begingroup$ What documentation do you refer to??? The contrast coefficients are always applied to the means as ordered in the grid. So your first impression is right. $\endgroup$ – rvl Aug 6 '18 at 18:56
  • 1
    $\begingroup$ Second, emmeans() computes predictions from the model, and those may not be the same as ordinary marginal means, as you are seeing happens. Read vignette(“basics”, package = “emmeans”) $\endgroup$ – rvl Aug 6 '18 at 19:00
  • $\begingroup$ @rvl: The package documentation in this link (on p. 16/17): cran.r-project.org/web/packages/afex/afex.pdf. Ah, that explains the slight difference. $\endgroup$ – Flo Aug 6 '18 at 19:04
  • $\begingroup$ Suggest reporting this to Henrik. Seems like a documentation blurp, or he’d meant to make phase an ordered factor. $\endgroup$ – rvl Aug 6 '18 at 19:29
  • 1
    $\begingroup$ @rvl Indeed a mistake in the documentation... $\endgroup$ – Henrik Aug 9 '18 at 8:37
2
$\begingroup$

Thanks a lot for reporting this. It is indeed a mistake in the documentation. This and all the following contrasts were wrong. I guess when I wrote the code the ordering of the factor was different or I just made a stupid error. Sorry for that.

The latest development version on github now does it correctly:

[...]
> r1 <- emmeans(a1, ~treatment +phase, model = "multivariate")
> r1
 treatment phase   emmean        SE df lower.CL upper.CL
 control   fup   4.333333 0.5511982 10 3.105187 5.561479
 A         fup   7.250000 0.6038074 10 5.904633 8.595367
 B         fup   7.291667 0.4611655 10 6.264126 8.319207
 control   post  4.083333 0.6277716 10 2.684571 5.482096
 A         post  6.500000 0.6876894 10 4.967733 8.032267
 B         post  6.625000 0.5252314 10 5.454711 7.795289
 control   pre   4.250000 0.7660323 10 2.543174 5.956826
 A         pre   5.000000 0.8391464 10 3.130265 6.869735
 B         pre   4.166667 0.6409086 10 2.738633 5.594700

Results are averaged over the levels of: gender, hour 
Confidence level used: 0.95 
> c1 <- list(
+   A_B_pre = c(rep(0, 6), 0, -1, 1),  # A versus B for pretest
+   A_B_comb = c(-0.5, 0.5, 0, -0.5, 0.5, 0, 0, 0, 0), # A vs. B for post and follow-up combined
+   effect_post = c(0, 0, 0, -1, 0.5, 0.5, 0, 0, 0), # control versus A&B post
+   effect_fup = c(-1, 0.5, 0.5, 0, 0, 0, 0, 0, 0), # control versus A&B follow-up
+   effect_comb = c(-0.5, 0.25, 0.25, -0.5, 0.25, 0.25, 0, 0, 0) # control versus A&B combined
+ )
> contrast(r1, c1)
 contrast      estimate        SE df t.ratio p.value
 A_B_pre     -0.8333333 1.0559027 10  -0.789  0.4483
 A_B_comb     2.6666667 0.7962586 10   3.349  0.0074
 effect_post  2.4791667 0.7624260 10   3.252  0.0087
 effect_fup   2.9375000 0.6694279 10   4.388  0.0014
 effect_comb  2.7083333 0.6519868 10   4.154  0.0020

Results are averaged over the levels of: gender, hour

You can install it via: devtools::install_github("singmann/afex@master")

Please let me know if you find other errors.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.