I apologise for the long question, I have made it as general as possible and I hope by keeping it as one question (with sample data) rather than splitting it up it may be useful to others fitting similar models.
Although I have a reasonable background in statistics this is my first experience with mixed effects models and unfortunately there seems to be a lot less clarity in the resources I have read compared to other statistical approaches.
The Problem
I have an observational study in which a number of sports people are monitored for a number of training sessions. The players mood state is measured each day as well as their training work rate (for the sake of simplicity). I am looking to examine the relationship between these to variable examine the hypothesis that mood state influences training work rate.
The Data
A sample of the data using R looks as follows:
sample <- read.table(text = " Session ID Response Predictor
1: 1 1 60.47 0.012
2: 2 1 42.53 -0.50
3: 3 1 64.45 0.01
4: 1 2 67.64 0.01
5: 2 2 51.77 0.04
6: 3 2 68.84 0.09
7: 1 3 79.80 -0.05
8: 2 3 46.95 0.43
9: 3 3 83.3 -0.05 ", h = T)
sample$Session <- factor(sample$Session)
sample$ID <- factor(sample$ID)
Where session, ID, Response and Predictor refer to training session, individual players, training work rate and mood state (normalised) respectively.
Question 1: The wellness score has been reported as a Z-score as seems to be the norm in similar studies (see here and here) but my intuition says that this might be unnecessary (as the mixed effect model should capture these between player differences) and perhaps a bad idea (as this papers have dropped all players that were not deemed to have normally distributed predictor observations). What is the best approach here?
The Methods
I'm setting the model up so that Session and ID are the random effects with predictor being the fixed effect. I am intereperating these random effects as being crossed random effects according following this blog. The model I fit is then:
model.1 <- lmer(Response ~ Predictor + (1 | ID) + (1 | Session), data=sample, REML = FALSE)
Question 2: Does this coding do what I think it should? And if so is this the best/only approach that would make sense within a mixed effects framework?
The Analysis
Assuming the model has been set up correctly, my main issue is determine whether there is a relationship between predictor and response. "The R Book" (p708) seemed to take a backwards stepwise regression approach using likelihood ratio test in a similar situation. In my case that would mean comparing my model to the null model as follows:
model.null <- lmer(Response ~ 1 + (1 | ID) + (1 | Session), data=sample, REML = FALSE)
anova(model.1, model.null)
Another approach I stumbled across was to use the stargazer package with:
stargazer(model.1, type = "text",
digits = 3,
star.cutoffs = c(0.05, 0.01, 0.001),
digit.separator = "")
Both of these approaches return the publication required p-values but I am unsure of the legitimacy of these values.
Question 3: Paired with appropriate diagnostics would either of these approaches be reasonable? What better measure could I find to determine if the relationship is significant?
Extension
Bit of a bonus question but now I'm considering the situation where I wish to extend my analysis to consider multiple measures of work rate and examine which if any are related to mood state. Would it make sense to set the previous predictor mood state as the response and set work_rate_1*work_rate_2*...*work_rate_n as the predictors applying a backwards stepwise process? Or should each of the work rate's be fit as their own model with some sort of correction factor e.g. Bonferroni? In the first case this post discusses the importance of considering random slopes in addition to random intercepts. Would random slopes apply here (only in the case of >1 predictor)?
