1
$\begingroup$

I have performed lmer on my data and, to be generic, my model has 3 fixed effects (A,B,C) and 3 nested random effects (d,e,f) following the response variable y: y~A * B * C+(1|d/e/f).

I also ran lsmeans(mymodel, specs=~ A|B:C) and I'm interested in seeing the radio contrasts instead of the linear / quadratic contrasts. For example, if A1 and A2 are levels of fixed effect A and B1 and B2 are levels of fixed effect B (C1 & C2 are levels of fixed effect C), I'm interested in calculating the estimate and pvalue of (B2-B1)/B1 at A1 and C1.

How can I accomplish this? I usually define my own contrast by making a matrix but that is a "linear" approach.


So as a follow up question, the pvalue is naturally very large because I'm limiting the difference to the absolute value. Is there a way to get the actual factor difference (A-B/A) instead of (|A-B|/A).

I'm using mvt for multiple comparison but is this the least stringent one? I would rather have false positives than false negatives.

$\endgroup$
4
$\begingroup$

It is doable! You have to do it in stages, though. Using the warpbreaks data to illustrate, I'll do such comparisons of wool at each tension:

> warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)

> library(emmeans)
> (warp.emm = emmeans(warp.lm, ~ wool | tension))

tension = L:
 wool   emmean       SE df lower.CL upper.CL
 A    44.55556 3.646761 48 37.22325 51.88786
 B    28.22222 3.646761 48 20.88992 35.55453

tension = M:
 wool   emmean       SE df lower.CL upper.CL
 A    24.00000 3.646761 48 16.66769 31.33231
 B    28.77778 3.646761 48 21.44547 36.11008

tension = H:
 wool   emmean       SE df lower.CL upper.CL
 A    24.55556 3.646761 48 17.22325 31.88786
 B    18.77778 3.646761 48 11.44547 26.11008

Confidence level used: 0.95 

The second stage is to set up the numerators and denominators:

> (custom.emm = contrast(warp.emm, list(`A-B` = c(1, -1), A = c(1, 0))))

tension = L:
 contrast  estimate       SE df t.ratio p.value
 A-B      16.333333 5.157299 48   3.167  0.0027
 A        44.555556 3.646761 48  12.218  <.0001

tension = M:
 contrast  estimate       SE df t.ratio p.value
 A-B      -4.777778 5.157299 48  -0.926  0.3589
 A        24.000000 3.646761 48   6.581  <.0001

tension = H:
 contrast  estimate       SE df t.ratio p.value
 A-B       5.777778 5.157299 48   1.120  0.2682
 A        24.555556 3.646761 48   6.734  <.0001

Now, use regrid() to fake a log transformation, needed because differences of logs are logs of ratios. Thus, we get the differences and back-transform them:

> pairs(regrid(custom.emm, transf = "log"), type = "response")

Non-positive response predictions are flagged as non-estimable
tension = L:
 contrast     ratio         SE df t.ratio p.value
 A-B / A  0.3665835 0.09688534 48  -3.797  0.0004

tension = M:
 contrast     ratio         SE df t.ratio p.value
 A-B / A     nonEst         NA NA      NA      NA

tension = H:
 contrast     ratio         SE df t.ratio p.value
 A-B / A  0.2352941 0.18695686 48  -1.821  0.0748

Tests are performed on the log scale 

Note that we need to re-run custom.emm with the first contrast reversed to get the result for medium tension, because the log transformation choked on the negative difference.

A couple of misc. notes: I used the new emmeans package, which is just out and succeeds lsmeans. If you still like the term "least-squares means", you may use lsmeans() instead of emmeans() and it will do the same thiong with different labels attached. Also, this does also work with the existing lsmeans package, except the final results are labeled a bit differently.

Addendum

Here's a way to get them all at once:

> custom.emm1 = contrast(warp.emm, list(`|A-B|` = c(1,-1), A = c(1,0)))
> custom.emm2 = contrast(warp.emm, list(`|A-B|` = c(-1,1), A = c(1,0)))
> custom.emm = custom.emm1[1:2] + custom.emm2[3:4] + custom.emm1[5:6]
> custom.emm
 contrast tension  estimate       SE df t.ratio p.value
 |A-B|    L       16.333333 5.157299 48   3.167  0.0161
 A        L       44.555556 3.646761 48  12.218  <.0001
 |A-B|    M        4.777778 5.157299 48   0.926  1.0000
 A        M       24.000000 3.646761 48   6.581  <.0001
 |A-B|    H        5.777778 5.157299 48   1.120  1.0000
 A        H       24.555556 3.646761 48   6.734  <.0001

P value adjustment: bonferroni method for 6 tests 

> test(pairs(regrid(custom.emm, tran = "log"), by = "tension", type = "response"), 
       by = NULL, adjust = "mvt")
 contrast  tension     ratio         SE df t.ratio p.value
 |A-B| / A L       0.3665835 0.09688534 48  -3.797  0.0012
 |A-B| / A M       0.1990741 0.23724291 48  -1.354  0.4477
 |A-B| / A H       0.2352941 0.18695686 48  -1.821  0.2057

P value adjustment: mvt method for 3 tests 
Tests are performed on the log scale

The Tukey method won't be allowed here because these aren't recognized as pairwise comparisons; but the mvt method is the moral equivalent of the Tukey adjustment. Using the + operator to do rbind() is a new feature in emmeans.

By the way, don't put whitespace around the - in |A-B| -- it will turn the labels into |A / B| / A

$\endgroup$
  • $\begingroup$ Thank you for your fast reply! So if I wanted to gt the medium tension, I should reverse the contrast matrix so that it would be (B-A) and then rbind this with the results from pairs(custom.emm) to get the pvalue adjusted for multiple comparison? Or is should I run the test() again on the combined matrix? $\endgroup$ – NoaHo Oct 24 '17 at 20:05
  • $\begingroup$ Yes, something like that. You need to do something like rbind(pairs1[c(1,3)], pairs2[2]) with the two results from pairs() $\endgroup$ – rvl Oct 24 '17 at 23:29

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.