How to check for homoscedasticity in a mixed effect model with longitudinal data in r

I am new to using mixed effect models and all the information online has me quite confused. I hope someone can help me.

I have data of patients who did the same test at three different timepoints. I want to use a mixed effect model (lme4 package) to look at effect of time on change in score. I added an example data set below. I have 3 variables. (1) The person who did the test (2) the change in score (3) if the change in score was between timepoints 1 and 2 (t12) or 1 and 3 (t13). I fitted the mixed effect model below with 'Timepoint' as a fixed effect and 'Patient' as a random effect (hoping this is correct).

model <- lmer(Change_score ~ Timepoint + (1|Patient), data=Example_data)

If im correct I now have the model but im struggling to plot the standardized residuals and do the Levene's test to see if the data has homoscedasticity.

Example dataset:

structure(list(Patient = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5), Timepoint = structure(c(1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), levels = c("t12", "t13"), class = "factor"), 
    Change_score = c(10, 0, 7, 15, -3, 13, 9, -4, 2, NA)), row.names = c(NA, 
-10L), class = c("tbl_df", "tbl", "data.frame"))

To plot the standardized residuals I found some information online suggesting I should do use the code below. This does give me a plot but I don't understand at all what im seeing on the plot and I feel like this is not what im looking for. I want to visualise the residuals from the model to asses homoscedasticity.

plot(model, resid(., type="pearson") ~ fitted(.), abline=0)

I tried doing the Levene's test (see code below) but due to missing values I get the error "Error in model.frame.default(form) : variable lengths differ (found for 'Example_data$Timepoint')" and I can't find a way to do the analysis with missing data.

Example_data$Timepoint <- as.factor(Example_data$Timepoint)
leveneTest(residuals(model) ~ Example_data$Timepoint)

I hope someone can help me with some sugestion on how to visualise the residuals and test for homoscedasticity in my data.

Solution

Your proximal problem is that you have an NA in your data set (thank you for giving an appropriate example!). The default na.action is to omit these values from the residuals etc. etc.. If you re-run your model setting the argument na.action = na.exclude, your test should work:

car::leveneTest(residuals(model) ~ factor(Example_data$Timepoint))

Levene's Test for Homogeneity of Variance (center = median)
      Df F value  Pr(>F)  
group  1  3.6071 0.09932 .
       7

You might also find performance::check_model(model) and library(DHARMa); plot(simulateResiduals(model)) useful (although note that the section option does not play nicely with na.action=na.exclude; you'd have to switch back to the default na.action=na.omit setting).