How to analyze interdependent interaction terms of lmer model

Question

Assume I test a number of patients repeatedly over time to see how a certain treatment changes their skin conductance in response to a certain colour (cond) after 2 months, 4 months, ... etc. I test the skin conductance on the palm, wrist or arm simultaneously to see whether the place matters.

I have a longitudinal dataset of the form:

SC = skin conductance: dependent variable (metric)
time = fixed factor with 6 levels (time1, time2, ..., time6) for each measurement in time
place = fixed factor with 3 levels (place A, place B, place C)
cond = fixed factor with 2 levels (cond1=blue, cond2=red)
ID = random factor with subject ID (because not each subject could be tested at each time point)

I put the data into the following form:

data$time  <- relevel(data$time, "time1")
data$cond  <- relevel(data$time, "blue")
data$place <- relevel(data$time, "placeA")

and do some model comparison:

model_1  <- lmer(sc ~ time * cond * place + (1 | ID), data)
model_1x <- update(model_1, REML = F)

model_2  <- update(model_1, .~. - time:cond:place)
model_2x <- update(model_2, REML = F)
anova(model_2x, model_1x)
# three-way interaction is not significant (p=0.93) --> leave three-way 
#  interaction out and continue with model_2

# remove two-way interactions:
model_3  <- update(model_2, .~. - time:cond)
model_3x <- update(model_3, REML = F)
model_4  <- update(model_2, .~. - time:place)
model_4x <- update(model_4, REML = F)
model_5  <- update(model_2, .~. - cond:place)
model_5x <- update(model_5, REML = F)
anova(model_3x, model_2x) # interaction is significant (p=0.003)**
anova(model_4x, model_2x) # interaction is not significant (p=0.46)
anova(model_5x, model_2x) # interaction is significant (p=0.039) *

In the end I end up with the final model being sc ~ time + cond + place + time:cond + cond:place + (1|ID), data)

Here are my questions:

• How can I look into the interaction terms? If e.g. place would have been dropped from the model altogether, looking at the following post hoc test would take the average over the factor place, right?

  posthoc_test <- glht(model_final, c("condred == 0", "condred + time2:condred == 0", 
                    ..., "condred + time6:condred == 0"))
  

Because, however, the factor place is still in the model, the above post hoc test is applied only to the baseline level of place, i.e. placeA

What if I want to see how time and cond interact, regardless of place? Is this even possible given that place is itself "captured" in the interaction cond:place?

• How do I report the significant interactions? Is the p-value I get from the model selection procedure (i.e., p=0.003 for time:cond) also the one I can later report for time:cond?

P.S.
By now I'm almost convinced that the answer to question 2 is, that the values I get from the model comparison are the ones I can report for time:cond in general. In short: it doesn't matter that I got the values from a model comparison instead of getting them from a preformulated ANOVA.

Please correct me if I'm wrong!

Answer 1

First, to look into the time:cond term I would check least-square (LS) means of sc at the 6 unique levels of time within each level of cond, and vice versa, using the final model. I would keep significant place and cond:place in calculating the LS means because removing these might give misleading values. This can be done using the lsmeans package:

lsmeans(model_final, ~ time|cond, cov.reduce = FALSE) ## or  ~ cond|time

Then I would repeat the same for the cond:place interaction.

Second, if you want to obtain "exact" p-values, I suggest you use lmerTest package. (The Kenward-Roger approximation it uses gives you exact degrees of freedom and hence p-values for the F-test if the data are balanced and have no no missing values.)