pairwise ANOVA of subset of data

Question

I need to perform multiple pairwise ANOVA's in R, and correct the p-values using bonferroni. However I don't need to compare every CLASS to each other. Below is my data format and selcontrasts: the pairs of which I need to contrast the log10relquant. Does any of you know how I could execute this? I use the dplyr, lsmeans and broom packages.

SEX      EXPERIENCED    AGE  CLASS compound    relquant log10relquant

1 FEMALE          NO     1D     1F      C14 0.004012910     -2.396541
2 FEMALE          NO     1D     1F      C14 0.003759812     -2.424834
3 FEMALE          NO     1D     1F      C14 0.003838553     -2.415832
4 FEMALE          NO     1D     1F      C14 0.003582754     -2.445783
5   MALE          NO     1D     1M      C14 0.005099237     -2.292495
6   MALE          NO     1D     1M      C14 0.005379093     -2.269291

selcontrasts <- c("1F - 1M", "4F - 4M", "4EF - 4EM", 
                  "7F - 7M", "7EF - 7EM", # sex differences
              "1M - 4M", "4M - 7M", "1M - 7M", "1F - 4F", 
              "4F - 7F", "1F - 7F", # age differences
              "4M - 4EM", "7M - 7EM", "4F - 4EF", 
              "7F - 7EF" # social experience)

x=list(selcontrasts)

Currently I'm using this to pair the whole dataset (so to compair every class) instead of the selected contrasts:

pvalsage=data.frame(datagr %>% 
    do( data.frame(summary(contrast(lsmeans(
          aov(log10relquant ~ CLASS, data = .), ~ CLASS ),               
          method="pairwise",adjust="none"))) ))

To only do the selected contrasts in list x, I tried:

pvalsage=data.frame(datagr %>%  
    do(data.frame(summary(contrast(lsmeans(
        aov(log10relquant ~ CLASS, data = .),~ CLASS),
        method = x, adjust="none"))) ))

But I get the error:

 error in contrast.ref.grid(lsmeans(aov(log10relquant ~ CLASS, data = .),  : 
 Nonconforming number of contrast coefficients

Don’t use `”pairwise” in the `contrast()` function — give a named lust with the desired contrast coefficients. See the documentation. Also you probably should upgrade to the **emmeans** package which continues **lsmeans** — Russ Lenth, Dec 20 '17 at 02:37
I keep getting the error: 'nonconforming number of contrast coefficients' — Sam Vanbergen, Jan 14 '18 at 15:49
You have to give a named list of numerical vectors, each of whose length equals the number of means being contrasted. Look at the examples in the help page for `contrast`. For example, to compare the first two of six means, that list element might look like `con12 = c(-1,1,0,9,0,0)` — Russ Lenth, Jan 14 '18 at 16:49
Minus one times the first one, plus one times the second one, plus zero times each of the rest. — Russ Lenth, Jan 14 '18 at 17:19

Russ Lenth · Answer 1 · 2018-01-14T22:28:31.310

If I understand the question correctly (and I very well might not), there are really three factors involved: SEX (two levels), EXPERIENCED (two levels), and AGE (3 levels, namely 1, 4, and 7). And what is required is separate comparisons of the levels of each factor, for each combination of the other two.

If that is the case, then combining the three factors into one named CLASS just makes it a lot harder, because it makes it much harder to keep track of the levels of the factors separately. What's simpler is to fit a model that accounts for all three factors, estimate the means for each factor combination, and then do the required comparisons using by variables. Thus, for each dataset dat, you do:

require(emmeans)
mod = aov(log10relquant ~ SEX * EXPERIENCED * AGE, data = dat)
emm = emmeans(mod, ~ SEX * EXPERIENCED * AGE)
rbind(pairs(emm, by = c("EXPERIENCED", "AGE")),
      pairs(emm, by = c("SEX", "EXPERIENCED")),
      pairs(emm, by = c("SEX", "AGE")),
      adjust = "bonferroni")

I did not try to embed this in the functional-programming paradigm; I'll leave it to the OP to figure out those details.

Note: The emmeans package (estimated marginal means) is a continuation of lsmeans, which will be deprecated in the future. It works the same way.

PS -- Looking at the code in the question, I am concerned that the end results will not show the actual estimates being compared (the EMMs), only the comparisons; and the further implication in the naming that really, only P values are sought. This grates on me. I don't like to watch people go straight to statistical tests without even looking at the quantities being tested.

score 0 · Answer 2 · answered Jan 20 '18 at 20:49

You could do the pairwise contrast anyway and then filter the rows only containing your selcontrasts into a new dataframe followed by p.adjust= bonferroni with only the contrasts of your interest.

or you could write a mycontr.lmsc function and define selcontrasts and use that as method =

(Y)

pairwise ANOVA of subset of data

2 Answers2