Assumptions

The assumptions of Repeated Measures ANOVA are:

1. Absence of outliers
   
2. Normality of the response variable
3. Constant variance in time or sphericity

Absence of outliers

The boxplots and density plots in the Data page show no substantial number of outliers in the data, satisfying this assumption.

Normality of the response variable

To assess the normality of the response variable Weight, we show QQ-plots for each of the eight factor combinations. No deviation from normality is noted in any of the plots.

Sphericity

Sphericity refers to the variance of the differences in time in the response variable. To test for it visually, we take the difference of the response variable across time periods. Since there are four time periods, Week 0, 2, 4 and 8, we have to calculate this difference for 4 choose 2 or six periods. We draw a boxplot for each period corresponding to differences in the response variable for the period, and compare the height or IQR of the boxplots to see whether variance is constant across periods. In this case, we see that variance is not constant, as a couple of boxplots have short, squat IQR while others have tall, extended IQR, pointing to unequal variance in time or non-sphericity.

We can also test for sphericity using Mauchly’s test. To do so, we use the anova_test() function from the rstatix library. A p-value below 0.05 means rejection of the null hypothesis that variance is constant across differences of the response variabel in time. Here, Mauchly’s test confirms the earlier finding from the boxplots that the data is not spheric, and therefore that this assumption of Repeated Measures ANOVA is not satisfied.

anova_test(data = WeightLossData,
           dv = Weight, # Response variable
           wid = Participant, # Random effect
           within = Week)$`Mauchly's Test for Sphericity`

  Effect     W        p p<.05
1   Week 0.046 8.56e-63     *

So, is absence of sphericity a show-stopper in Repeated Measures ANOVA? No. Non-sphericity can be corrected in one of two ways. The first is via corrections in the degrees of freedom in the ANOVA calculation, resulting in adjustments to F-values, after fitting the model. These post hoc corrections are beyond the scope of this tutorial, and have been increasingly superseded by a second way to account for non-sphericity. This second way is to correct for non-sphericity by specifying the model covariate structure prior to fitting the model. This is the method described in this tutorial.

Before we continue, we must understand that, regardless of the method employed, non-sphericity needs to be corrected lest the resulting parameter estimates increase Type I errors (‘false positives’), leading us to reject the null hypothesis when we ought to admit it.

Correlated residuals

Recall that for each participant in the study four measures of weight in time were taken, and therefore they are correlated. Let’s check whether this is true. To do so we model the data with an Ordinary Least Squares (OLS) regression as if the observations were all independent and check the residuals for correlation or a time-series structure. The ACF and PACF plots for the residuals show that there are significant lags or correlations. The tall spike at lag1 on the PACF and the tapering sinusoidal on the ACF hint strongly at an AR(1) time-series structure in the residuals, pointing to the unsuitability of an OLS model that does not account for time or repeated measures.

Residual ACF / PACF plot

Before moving on, we need to understand why correlated residuals are problematic. One reason is that they skew our parameter estimates, returning models with less reliable coefficients. A second reason is that patterns in the residuals point to poor model fit. For if the model explained the variation in the data, it would not leave behind significant patterns or discernible variation in the residuals. It would leave instead white noise.

Code

attach(WeightLossData)

### Fit OLS regression model
model.OLS = lm(Weight ~ Regimen + Week + Regimen*Week)

### Extract and plot ACF/PACF of residuals
acf(residuals(model.OLS)) ; pacf(residuals(model.OLS)))

Covariate structure

To fit a Repeated Measures ANOVA model in R, we use the gls() function from the nmle library. A first, essential step in a fitting the model is to determine which time correlation structure or covariate structure best fits the data. To do so, we fit the model with a few covariate structures and chose the best one based on model fit measured by AIC, BIC and LogLikelihood. A list of the correlation structure classes available with the gls() function is found here.

Code

In specifying the model, note two things. One, the ‘subject’ or experimental unit for which repeated measures are taken is specified in the corr parameter. Two, the experimental units, the participants, are modelled as a Random Effect, denoted by the numeral one followed by bar, in the 1 | Participant specification. Why? One reason is that we are dealing with a sample of participants rather than the population. Another reason is that we are not interested in the individual change in weight of each participant, but rather we are interested in accounting for the variability of the participants as a collective or group. To do so in the model, they are specified as a Random Effect.

### Covariate structure with no correlation
mod.nostruct = gls(Weight ~ Regimen + Week + Regimen*Week,
                   corr=corSymm(, form= ~ 1 | Participant),
                   weights = varIdent(form = ~ 1 | Week))

### Covariate structure with 'compound symmetry' or constant, non-zero correlation
mod.compsym = gls(Weight ~ Regimen + Week + Regimen*Week,
                  corr=corCompSymm(, form= ~ 1 | Participant))

### Covariate structure with time-series AR(1) or correlation of order 1
mod.ar1 = gls(Weight ~ Regimen + Week + Regimen*Week,
              corr=corAR1(, form= ~ 1 | Participant))

anova(mod.nostruct, mod.compsym, mod.ar1) ### Compare the three models

corrSymm() specifies a covariate structure with no structure or no correlation along measures. This is identical to specifying an OLS model that does not take time repeated measures into account. corrCompSymm specifies a covariate structure with compound symmetry or constant, non-zero correlation along measures.

corAR1 specifies a covariate structure with correlation of order one (1) along measures, meaning that correlations decay to the power of time difference. The greater the time difference along measures, the smaller the correlation. For example, if the correlation between measures at times 1 and 2 is, \(0.80^{2-1=1}\), then the correlation between measures at times 1 and 3 is \(0.80^{3-1=2} = 0.64\), and the correlation between measures at times 1 and 4 is \(0.80^{4-1=3} = 0.51\). This decay in correlation along measures is usually a more realistic assumption that either no correlation or constant correlation. The anova() function compares the three models.

Model comparison

In the model comparison, the AR(1) covariance structure model returns both lower AIC and BIC scores than the two competitors, suggesting that it offers the best fit among the models tested. The differences among the models are not stark, but confirm the earlier finding from the OLS model ACF/PACF plots that the residuals have an AR(1) structure. Hence, we select the AR(1) covariate structure as the best fit for the data.

             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
mod.nostruct     1 18 2220.687 2292.170 -1092.344                        
mod.compsym      2 10 2283.248 2322.960 -1131.624 1 vs 2 78.56038  <.0001
mod.ar1          3 10 2209.081 2248.794 -1094.541                        

Model diagnostics

The ANOVA table for the AR(1) model shows that all terms, including the interaction term, are significant at 95% confidence, pointing to the statistically significant impact of the factors of time and regimen, as well as the interaction of the two, on weight.

Model diagnostics are satisfactory, satisfying the assumptions of constant variance and normality of the residuals. We are good to go and interpret the results.

AR(1) model ANOVA

Model fit diagnostics