I have collected data looking at the effect of three dietary treatment groups on protein synthesis across three time points. Because I have missing data, I have elected to conduct a linear mixed model analysis on the data as opposed to a two-way repeated measures ANOVA. Specifically, I am using the 'nlme' package in R (version 4.0.2) on a mac computer [x86_64-apple-darwin17.0 (64-bit)].
I am particularly interested in looking at the interaction between diet (group) and time point (phase)on the dependent variable (in my case, fractional synthesis rate, fsr). Thus, I ran the following line of code:
model5 <- lme(fsr~group + phase + group*phase,data=data, random=~1|subID, na.action=na.omit, method="ML")
and get the following output (truncated for brevity):
Linear mixed-effects model fit by maximum likelihood
Data: data
AIC BIC logLik
-26.46228 -3.424489 24.23114
Random effects:
Formula: ~1 | subID
(Intercept) Residual
StdDev: 0.1827453 0.1092466
Fixed effects: fsr ~ group + phase + group * phase
Value Std.Error DF t-value p-value
(Intercept) 1.3264977 0.08164724 32 16.246694 0.0000
groupSM 0.0109315 0.12426358 19 0.087970 0.9308
groupWM 0.0763769 0.11951932 19 0.639034 0.5304
phase2 0.1161971 0.06281472 32 1.849838 0.0736
phase3 0.3330094 0.06281472 32 5.301454 0.0000
groupSM:phase2 0.0350255 0.09878970 32 0.354546 0.7253
groupWM:phase2 0.0151018 0.09183417 32 0.164447 0.8704
groupSM:phase3 -0.0940546 0.09546707 32 -0.985205 0.3319
groupWM:phase3 -0.0371642 0.08920422 32 -0.416620 0.6797
If I then look at Type III tests of Fixed effects using:
anova.lme(model5, type="marginal", adjustSigma = F)
I get the following output:
numDF denDF F-value p-value
(Intercept) 1 32 310.53538 <.0001
group 2 19 0.26934 0.7667
phase 2 32 16.97659 <.0001
group:phase 4 32 0.62842 0.6457
So, there is no significant interaction between group and phase, but there is a main effect of phase.
My Question: How can I look specifically at what phases are different? There were three time points in my study, so I want to look at the pairwise comparisons (1 v 2, 1 v 3, 2 v 3) to see where there was an effect. Here is what I have tried:
tukey_mod5<- summary(glht(model5,linfct=mcp(phase="Tukey")))
...which gives me the following output:
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = fsr ~ group + phase + group * phase, data = data,
random = ~1 | subID, method = "ML", na.action = na.omit)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 0.11620 0.05791 2.006 0.111
3 - 1 == 0 0.33301 0.05791 5.750 <0.001 ***
3 - 2 == 0 0.21681 0.05839 3.713 <0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Does this seem like a logical (i.e. correct) approach to take? I will admit that I am new to mixed models and may be going about it incorrectly. If so, I would appreciate being steered in the right direction!
