Mixed effects model for repeated measures approach in R

Question

I'm using mixed effects models for repeated measures (MMRM) in R with the nlme package for the first time as part of a research project and have read lots of posts here to learn about the intricacies of these models. I am hoping to verify that I'm using the right approach — and ask a couple of questions that I haven't seen on here. Your feedback and insights would be much appreciated!

The data is from a trial with repeated measures on a survey that was measured from Baseline to Week 12 at uneven time points. The intervention was given at Day 0. The timepoints of interest for the study outcome are Week 1, Week 3, and Week 12. It's a small dataset (n < 20) and it has missing values for individuals that dropped out of the study at different timepoints (Week 1, Week 6, Week 9) or who randomly missed a week.

I've created some mock data that reflects the shape and roughly the values of the data:

ID	time	category	gender	score	change
P1	bl00	0	1	33	0
P1	wk01	0	1	29	-4
P1	wk03	0	1	23	-10
P1	wk12	0	1	24	-9
P2	bl00	1	0	40	0
P2	wk01	1	0	20	-20
P2	wk03	1	0	25	-15
P2	wk12	1	0	25	-15
P3	bl00	0	0	60	0
P3	wk01	0	0	21	-39
P3	wk03	0	0	NA	NA
P3	wk12	0	0	NA	NA

ID is the participant ID, time is the time point of the measurement (bl00 = baseline, wk01 = Week 1, etc), category is a binary categorical variable, gender is binary of M/F, score is the score for that time point, and change is the change in score from baseline to that time point. I'm treating ID, time, category, and gender as factors.

A couple of questions:

Question 1.

I plotted the mean change over the four timepoints and found that it is not totally linear. It drops down dramatically to week 1 and then increases a bit towards week 12. So I transformed the visit to days (bl00 = 0, wk01 = 7, etc) and added it as a quadratic effect for time using splines. Then I tested model performance based on AIC and BIC:

data$days_ns <- ns(data$days, 3, intercept=FALSE)

null <- nlme::lme(
  score ~ 1, random = ~1|ID, 
  data, 
  na.action=na.exclude, method="ML"
)

linear <- nlme::lme(
  score ~ time, random = ~1|ID, 
  data, na.action=na.exclude, method="ML"
)

quad <- lme(
  score ~ days_ns, 
  random = list(ID = pdDiag(~ days_ns)), data, 
  na.action=na.exclude, method="ML"
)

Each model was a significant (p < 0.001) improvement from the previous one:

Model	df	AIC	BIC	logLik
null	1	6	336.8193	348.0465
linear	2	6	279.7095	290.9367
quad	3	9	264.3566	281.1974

But with the quad model the interpretability goes out the window since the mean inferred change in the summary table is no longer related to the different timepoints. Does it make sense to go with the slightly worse performing linear model in favor of interpretability of the model output? Also, is there something wrong in my approach to the models themselves?

Question 2.

Given that nothing is glaringly wrong here (fingers crossed)... When it comes to describing and interpreting the output of the model using the linear MMRM model:

Ultimately, I used the model:

linear <- lme(score ~ time + category + gender, random = ~1|ID, data, method="REML", na.action=na.exclude)

• Does this look correct? And is it right to say: "participant ID was added as a random effect, which allows for varying intercepts and slopes for each participant"?
• I ran an anova on the model output using aov <- anova(linear.model). Is it appropriate to say, "There was a large, statistically significant effect of time on scores using the MRMM model (F(3, 32) = 32.35 , p <.0001)" and in describing the approach: "The P value associated with the time effect in the MMRM model obtained through type II analysis was used to assess the overall change in scores"?
• Lastly, based on this post and the method developed by Nakagawa & Schielzeth (2013), I used the r.squaredGLMM(x) function to calculate Pseudo-R-squared values for the fixed and fixed + random effects of the MMRM model. According the result, 83% of the variance can be explained by fixed + random effects, and 60% is explained by fixed effects alone. Does this reflect a lot of variance for each participant? If the amount of variance explained by fixed vs. fixed + random were the same, does it mean that the random effects should not be included and a different model should be used?