# Running post hoc analysis after running Linear mixed model in R

I have a dataset that has 2 groups ("MICT" & "HIIT") and 3 time points ("pre" & "Trained + AcuteEx" & "Trained"). Below is the code that I ran to see if there are any significances in fixed effects.

lm3 <- function(x){
return(summary(lmer(x ~ Time+Group+Time:Group + (1|ID), data = htat122120)))
}


The output gives a summary like this:

I learned that a post-hoc analysis should be run if there are any significance in interaction (time x group) to find out the specific time point or a group that is showing a difference. However, as shown in the image file, I get 5 different fixed effects with different p-values that are: TimeTrained + AcuteEx, TimeTrained, GroupMICT, TimeTrained + AcuteEx:GroupMICT and TimeTrained:GroupMICT.

In this given example, there are significance in certain time, group, and interaction. What would be the correct interpretation? Should I run post-hoc analysis? If so, how can I run the post-hoc analysis?

Thank you so much.

You need to think through what each of your coefficients represents, given the treatment coding of your categorical predictors, the default in R. With that coding, the intercept represents the estimated baseline outcome when continuous predictors are at 0 and all categorical predictors are at their reference levels. Coefficients for individual predictors then represent the associated differences from that baseline. Coefficients for interactions represent the next level of differences, from what would be predicted simply from the corresponding individual-predictor coefficients. Each coefficient p-value is for a test of whether that coefficient equals 0. It's the probability that you would get at least that large a value due to random sampling if the true value was 0.

In your case, without continuous predictors, the Intercept is the estimated outcome for the HIIT group at Time=pre. That's clearly different from 0, based on the very small p-value. Be careful interpreting intercepts in general, however, as centering a continuous predictor or choosing a different reference for a categorical predictor might lead to an apparently "insignificant" intercept without affecting the model overall.

The coefficients for other values of Time are differences from that baseline value, within the HIIT reference group. Within that group, the value at Time = Trained + AcuteEx isn't different from that baseline, while that at Time = Trained is significantly lower.

The coefficient for the MICT group is the difference from the Intercept value (for the HIIT group) at Time=pre. By the usual standard of p < 0.05, the outcome for the MICT group is significantly lower than that for the HIIT group at that particular time.

The interaction terms are the differences from what you would have predicted otherwise, based on the Intercept and the individual coefficients. The first, for Time = Trained + AcuteEx and Group = MICT, isn't significantly different from 0 (p = 0.85 means you would have found at least that large a difference 85% of the time due to random sampling even if there really is no difference). Thus the response for that combination of Time and Group isn't significantly different from what you would have predicted based on the individual coefficients for that Time and Group.

In contrast, the interaction coefficient for Time = Trained and Group = MICT is significantly different from 0. The outcome for that combination is significantly higher than what you would have predicted based on the individual coefficients for that Time and Group.

Depending on the goal of your study, that might be all that you need to know. If you want to go farther, it can help to look at the actual estimated outcomes for each combination of predictor values, rather than the differences reported in your displays. Those estimates are:

Pre Trained + AcuteEx Trained
HIIT 2.55 2.48 1.49
MICT 1.15 1.16 1.44

Post-hoc tests of differences among these estimated outcomes are based on the standard errors of the linear combinations of coefficient estimates associated with those differences. As the coefficient estimates aren't independent of each other, you need to know the covariances among them to do that calculation, and apply the formula for the variance of a weighted sum of correlated variables. The variance-covariance matrix among the coefficient estimates is calculated as part of the model fit and is stored as part of the model object. In R, a vcov() function typically retrieves that matrix from the model.

Although it's important to understand these principles, you don't have to do these calculations by hand. Statistical software typically provides ways to summarize model outcome estimates and test differences of interest, based on those principles. As you are working in R, take the time to learn how to use the emmeans package, which greatly simplifies such calculations for a wide variety of model types.

Finally, you should recognize that how best to perform significance tests based on mixed models like yours is open to some dispute. This answer is a useful summary of the issues, with links for further study. Provided that you report exactly what you did you should be OK, but it's important to recognize the limitations as well as the strengths of the tools that you use.