# Addressing "NOTE: Results may be misleading due to involvement in interactions" warning with Tukey post-hoc comparisons in lsmeans R package

Background:
I am using linear mixed-effects models (LMMs) in order to determine how the interaction between two fixed effects influences measures of a response variable. Since I am working with a dataset in which there are multiple samples from multiple individuals that could violate the assumption of independence of data points, I am treating "individual" as a random effect. Thus, the generic model I am working with is:

lmer(y ~ Factor1*ContinuousVariable1 + (1|Ind), dataset, REML=T)


Note: for my actual dataset, I used a likelihood ratio test to determine whether I needed to also nest the multiple trials within individual [i.e., lmer(y ~ Factor1*ContinuousVariable1 + (1|Ind/Trial), dataset)], and failed to reject the null hypothesis that this "fuller" model contributed significantly to accounting for additional variation in the data.

Problem to solve:
Determine whether the results from my Tukey's post-hoc comparisons are reliable, given the interactions included in my LMM model.

library(car) # for Soils dataset
data(Soils)
library(lme4) # for lmer()
library(lsmeans) # for remaining functions

Example code:
## Create the LMM
## "Na" is a numeric continuous response variable
## "Contour" is a factor, with character categories, and is treated as a fixed effect
## "P" is an integer variable, is treated as a fixed effect, and differs across the Contour groups
## "Group" is a a numerical factor and is treated as a random effect

Na.LMER <- lmer(Na ~ Contour*P + (1|Group), Soils, REML=T)
Na.LMER

Linear mixed model fit by REML ['lmerMod']
Formula: Na ~ Contour * P + (1 | Group)
Data: Soils
REML criterion at convergence: 190.4919
Random effects:
Groups   Name        Std.Dev.
Group    (Intercept) 2.514
Residual             1.063
Number of obs: 48, groups:  Group, 12
Fixed Effects:
(Intercept)    ContourSlope      ContourTop               P  ContourSlope:P    ContourTop:P
7.104951        4.381251       -0.260527       -0.006811       -0.026952       -0.006258

### Conduct Tukey's post-hoc comparisons
Na.Tukey <- lsmeans(Na.LMER, pairwise~Contour, adjust="tukey")


NOTE: Results may be misleading due to involvement in interactions

Na.Tukey

$lsmeans Contour lsmean SE df lower.CL upper.CL Depression 5.973118 1.289466 8.15 3.008857 8.937380 Slope 5.875929 1.286895 8.08 2.913697 8.838160 Top 4.672639 1.294933 8.24 1.701416 7.643863 Confidence level used: 0.95$contrasts
contrast             estimate       SE   df t.ratio p.value
Depression - Slope 0.09718976 1.821763 8.11   0.053  0.9984
Depression - Top   1.30047917 1.827450 8.19   0.712  0.7635
Slope - Top        1.20328941 1.825636 8.16   0.659  0.7925

P value adjustment: tukey method for a family of 3 means


So this is where the question comes in.
Since I received the warning message ("NOTE: Results may be misleading due to involvement in interactions"), I want to verify whether I can reliably use the p-values output from lsmeans() to determine which contrasts were different from each other. So how can I tell whether the interactions from this particular dataset could be problematic for interpreting the results from the Tukey's post-hoc comparisons.

Here is what I have tried to investigate this issue.
Based on the recommendations by Professor Russell Lenth (developer of the lsmeans R package), I used additional functions from the lsmeans R package to investigate what's going on with the data.

### First, here are the F-tests of the fixed effects of the LMM.
anova(Na.LMER)

Analysis of Variance Table
Df  Sum Sq Mean Sq F value
Contour    2  0.5696  0.2848  0.2520
P          1 10.4083 10.4083  9.2093
Contour:P  2  6.7070  3.3535  2.9672


Does the Contour:P interaction seem relatively strong?

Next, I'm going to evaluate whether this interaction is important by determining to what extent the values of P varies across the Contour groups, using lsmip().

Na.lsm <- lsmeans(Na.LMER, ~Contour|P, at=list(P = c(75, 100, 200, 300, 400)))
Na.lsm

P =  75:
Contour       lsmean       SE    df   lower.CL  upper.CL
Depression  6.594094 1.399580 10.70  3.5029413  9.685246
Slope       8.953983 1.562754 13.53  5.5913341 12.316631
Top         5.864180 1.511863 12.76  2.5917619  9.136598

P = 100:
Contour       lsmean       SE    df   lower.CL  upper.CL
Depression  6.423808 1.355688  9.64  3.3876590  9.459957
Slope       8.109909 1.429365 10.79  4.9562943 11.263524
Top         5.537432 1.391548 10.16  2.4433848  8.631479

P = 200:
Contour       lsmean       SE    df   lower.CL  upper.CL
Depression  5.742665 1.286244  8.08  2.7814923  8.703838
Slope       4.733616 1.354120  9.32  1.6863856  7.780847
Top         4.230440 1.384598 10.01  1.1459415  7.314939

P = 300:
Contour       lsmean       SE    df   lower.CL  upper.CL
Depression  5.061522 1.396923 10.63  1.9738402  8.149204
Slope       1.357323 1.960472 21.77 -2.7109112  5.425557
Top         2.923449 2.025495 24.22 -1.2549312  7.101829

P = 400:
Contour       lsmean       SE    df   lower.CL  upper.CL
Depression  4.380379 1.651907 17.57  0.9037052  7.857053
Slope      -2.018970 2.841216 33.67 -7.7950921  3.757152
Top         1.616457 2.914268 36.01 -4.2938885  7.526803

Confidence level used: 0.95


### Plotting the interactions

Na.lsmip <- lsmip(Na.lsm, Contour~P)

### It seems like the levels of Contour vary at different values of P (especially for Slope), but I'm going to use pairs() to verify this using pairwise comparison at each value of P.

pairs(Na.lsm)
P =  75:
contrast             estimate       SE    df t.ratio p.value
Depression - Slope -2.3598888 2.097862 12.15  -1.125  0.5175
Depression - Top    0.7299139 2.060232 11.74   0.354  0.9335
Slope - Top         3.0898026 2.174381 13.15   1.421  0.3589

P = 100:
contrast             estimate       SE    df t.ratio p.value
Depression - Slope -1.6861012 1.970019 10.22  -0.856  0.6784
Depression - Top    0.8863760 1.942755  9.90   0.456  0.8928
Slope - Top         2.5724773 1.994865 10.47   1.290  0.4308

P = 200:
contrast             estimate       SE    df t.ratio p.value
Depression - Slope  1.0090489 1.867637  8.70   0.540  0.8539
Depression - Top    1.5122246 1.889851  9.04   0.800  0.7122
Slope - Top         0.5031757 1.936686  9.67   0.260  0.9636

P = 300:
contrast             estimate       SE    df t.ratio p.value
Depression - Slope  3.7041990 2.407248 16.78   1.539  0.2988
Depression - Top    2.1380732 2.460493 18.21   0.869  0.6660
Slope - Top        -1.5661258 2.818879 23.01  -0.556  0.8447

P = 400:
contrast             estimate       SE    df t.ratio p.value
Depression - Slope  6.3993492 3.286534 28.79   1.947  0.1439
Depression - Top    2.7639218 3.349889 30.81   0.825  0.6905
Slope - Top        -3.6354273 4.070070 34.89  -0.893  0.6481

P value adjustment: tukey method for a family of 3 means


### Based on the pairs() output, it doesn't seem like Contour groups vary at these incremental values of P.

Since the Contour groups do not seem to vary at different levels of P, does that mean that the interaction strength is not that strong? and thus, I am okay to ignore the warning message that "NOTE: Results may be misleading due to involvement in interactions"?

I would appreciate any feedback about interpreting these results, and whether there are additional analyses that I should be conducting in order to address my concern. If there is any additional information that would be helpful in tackling this problem, please let me know.

Thank you for your time!

UPDATE (2/6/15): I had a minor typo at the beginning, in which the first line of code read "Dens.LMER <- lmer(...)". The lmer product should have been named "Na.LMER", which was used in the remaining code. Thus, the Dens.LMER product that rvl mentions is equivalent to Na.LMER. I apologize for the inconvenience.

My view is that the $F$ test of statistical significance of the interaction effect is less important than the subjective nature of the interaction, as evidenced by the plot. The plot tells me that it is reasonably sensible to compare the overall averages of Depression and Top, but it'd be silly to compare those averages with the overall average of Slope -- whether or not these comparisons are statistically significant. Basically, I'd say to avoid doing comparisons that don't make sense -- so my advice is do not ignore the warning note in this case. If the curve for Top were fairly parallel with the other two, that's when you could ignore it.

In general, I suggest looking at enough plots that you can tell what's going on, and then restrict your post-hoc testing to things that are sensible.

Since P is continuous, you're really fitting straight lines (they look curved because you chose unequally spaced points). You can compare the slopes of these lines:

R> lstrends(Dens.LMER, pairwise ~ Contour, var = "P")

$lstrends Contour P.trend SE df lower.CL upper.CL Depression -0.00681143 0.004901195 39.68 -0.01671957 0.003096714 Slope -0.03376293 0.010533875 41.88 -0.05502295 -0.012502911 Top -0.01306992 0.010499548 41.97 -0.03425936 0.008119525 Confidence level used: 0.95$contrasts
contrast               estimate         SE    df t.ratio p.value
Depression - Slope  0.026951501 0.01161827 42.00   2.320  0.0639
Depression - Top    0.006258486 0.01158716 41.81   0.540  0.8520
Slope - Top        -0.020693015 0.01487290 41.99  -1.391  0.3545

P value adjustment: tukey method for a family of 3 tests


The comparison between the shallowest and largest slopes has an adjusted $P$ value of about $.06$.

• Thank you so much @rvl for the thorough answer! Am I correct that although the output from pairs() and lstrends() suggest that the Contour factor levels do not vary at different values of P and do not have significantly different levels (alpha = 0.05), respectively, the plot from lsmip() illustrates quite starkly that Slope has a much different interaction with P than the other Contour factor levels. Is this correct? If so, is there a better way to evaluate whether Na differs amongst the Contour levels, while accounting for variation from Group, but ALSO correct for the variation of P? Commented Feb 6, 2015 at 14:28
• The output I show essentially tests the significance of what you see in the graph-- that the slopes are unequal. You could do another contrast on the lstrends results comparing the Slope trend with the average of the other two. I'd suggest using the Scheffe correction as it's based on exploration. If I understand your last question correctly, no. The whole point is there is an interaction, so you can't ignore P while comparing the Contour means. Commented Feb 7, 2015 at 2:58