7
$\begingroup$

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.

Loading data and libraries:
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)

Interaction of Contour and 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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
4
$\begingroup$

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$.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – skawano
    Feb 6, 2015 at 14:28
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Feb 7, 2015 at 2:58

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.