How to Run Linear Mixed Effects Analysis for Pairwise Comparisons? A Tutorial and a Proposal for the Calculation of Standardized Effect Sizes

1.1.1. Analysis with t-test, ANOVA, and regression

The easiest-to-understand test for comparing two groups is the t-test. It basically compares the difference in mean score between the two groups and looks at how big this difference is relative to the variability in each group.

The t-test requires one observation per participant. This can pose a challenge for beginners as experimental data sets often contain several observations per participant. This is also true for the data in Table 1, where there are two observations per participant. A simple solution is to take the mean value per participant, as shown in Table 2.¹ Table 2 also includes the mean (M) and standard deviation (SD) of each group. This shows that the 65-year-olds generally had higher memory scores (M = 39.5) than the 75-year-olds (M = 34.5), a difference that can also be seen in Figure 1. We want to know whether the difference is statistically significant.

Figure 1

Because the data set is small, we can calculate the t-test by hand (Brysbaert, 2020), using the following equation (N = number of participants in each group):

To interpret the t-value, we need the degrees of freedom, which for a between-groups t-test is N_65yr – 1 + N_75yr – 1 = 10 – 1 + 10 – 1 = 18. This gives us a p-value of t(18) = 1.02, p = .323, two-tailed,² meaning that the difference is not significant at the .05 level.

An easier way to get the result is to use existing software. There are plenty of software packages doing the calculations for us. One of the richest is R. The t-test is part of many libraries in R and it is also part of the base packages, meaning that we do not need to install or load a package.

The standard way to set up the t-test for independent samples in R is to input the data in long format, with each observation on a separate row, as shown in Table 3. You can do this in R (see code below), but also in Excel. We will need the long format for most of the article.

Once you have imported the data set in the right format in R,³ you can run the following command:

The output gives you all the information you need, except that the better⁴ Welch t-test is used by default instead of the original Student t-test. Notice that the p-values is for a two-sided test as default. This will be case for all the packages we describe.

The p-value indicates the probability that we would observe this particular set of data if the null hypothesis were true. It can be seen as a quantification of the evidence against the null hypothesis: the smaller the p-value, the more evidence against the hypothesis of equal means. For the data set we analyzed, the evidence is not strong enough to assume that the null hypothesis is wrong, as the p-value is not smaller than .05. Because the p-value depends both on the difference between the groups, as well as on the sample size, it cannot be interpreted as a measure of effect size.

The most informative effect size consists of the mean and the standard deviation (SDs) in each condition (Baguley, 2009; Pek & Flora, 2018; Rights & Sterba, 2019; Vanhove, 2015). Therefore, it is crucial to always present this information in a clear and concise table within your statistical analysis section.

Additional information is provided by a standardized effect size. This metric expresses the size of a difference in relative terms. It allows you to measure the impact of your finding compared to findings from other studies, even if those studies worked with different materials and variables.

One of the most used standardized effect sizes in psychology is Cohen’s d. This is the difference between the means divided by the average standard deviation of the conditions.

You obtain Cohen’s d if you install and load the R package “effectsize” (Ben-Shachar et al., 2023). You do this as follows:

This gives the output:

The outcome not only tells you that the standardized effect size is d = .45, but also that its 95% confidence interval ranges from d = –.44 to d = 1.34. This basically means that with a sample size of 10 participants per group, you may obtain d-values going from a medium effect in favor of 75-year-olds to a very large effect in favor of 65-year-olds. This is a clear sign that your study is not very informative about the effect size at the population level, because of the small sample sizes (Brysbaert, 2019, 2021).

We can also calculate the t-test with a freely available, novice-friendly statistical package, such as jamovi (The jamovi project, 2023). We assume that entering data and clicking on buttons is clear for everyone with minimal computer knowledge (you also find all jamovi files with data entered and the right buttons clicked at https://osf.io/s3efz/; if you have jamovi installed on your computer, simply start the program and open the files). The output is shown in Figure 2.

Figure 2

The same dataset can be analyzed with an analysis of variance or ANOVA. For the example, we need an ANOVA with Age as the only between-subjects independent variable.

In R you get the outcome with the following commands:

It gives the output (see Brysbaert, 2020, for a discussion of how ANOVA tables are calculated and what the various entries mean):

Figure 3 gives the output of jamovi.

Figure 3

The p-values of the F-test in ANOVA and the t-test are the same, because the F and t-values are mathematically related to each other. When there are only two conditions, the F-value of the ANOVA equals the square of the t-value from the t-test [i.e., F(1,18) = 1.032 = t²(18) = 1.016²].

ANOVAs give us another measure of standardized effect size. It is the percentage of variance accounted for by the independent variable, called eta-squared (η²). In the jamovi output we see that the Group difference accounted for 5.4% of the variance between participants (η² = .054). To get the η² value in R, use the following commands:⁵

The notion “percentage of variance explained” is more familiar in linear regression analysis. There it is customary to see how much variance is explained by the predictors and call it R². As it happens, we can analyze our dataset with a linear regression analysis as well, using Group as a categorical (or nominal) variable meaning that there is no magnitude associated with it.

The R-commands are shown below. They are rather simple because on the basis of the participant names R automatically detects that the Group variable is a categorical variable.

The output of the linear regression gives the same information as the t-test or the ANOVA: The effect of Group is not significant (t = –1.02, p = .323) and accounts for 5.4% of the variance in the observations (R² = .05421 = η²).

To run the analysis in jamovi, we have to indicate that Group is a categorical variable by putting it under Factors and not under Covariates (which are continuous variables, like the age of the participants in years).⁶ Figure 4 shows the output.

Figure 4

Now that we know what our data look like in traditional analyses, we can compare them to a linear mixed effects (LME) analysis. The outcome will not be exactly the same, because in the traditional analyses we took the average performance of the participants, whereas in LME we take both measurements of each person into account. However, the pattern of findings should be very similar, because in our example we did not introduce big differences between the two measurements.

1.1.2. Analysis with LME

An advantage of LME over a t-test is that we do not have to average the data for each participant. We can just work with the raw values from Table 1. Thus, we do not have to calculate Table 2. This is not only convenient, but also better from a statistical point of view, because we are considering the full variability in the data set. When we average the scores of day 1 and day 2 for a t-test or ANOVA, we assume that there are no differences between these scores (or at least that they are not important for the question we are investigating).

A new element in LME analysis is that we must make a distinction between fixed effects and random effects. Fixed effects refer to the variables we manipulate and for which we are interested in the levels tested. In our example, that is the age of the participants (two groups). Random effects refer to variables that also influence the data, but for which we have no specific predictions or expectations. In our example, these are the participants. We accept that participants differ from each other in level of performance: Some participants have a higher working memory capacity than others. We have no further explanation for these differences. We simply assume that the participant scores are a random sample from a normally distributed population.

The difference in performance level between participants is called the participant intercept. To estimate the participant intercept, we need at least two observations per participant (if you have only one, the analysis reduces to a t-test, ANOVA, or regression). In an LME analysis of two independent groups we assume that the data are a combination of the fixed effect we manipulated (age group), a random intercept (performance level) for each participant, and residual measurement error.

The R package mostly used for LME is lme4 (Bates et al., 2022), often accessed through the lmerTest package (Kuznetsova et al., 2022), which provides p-values. For this analysis, we need a new input file in long format, which contains all the information from Table 1. It looks like Table 4.

Then you can run the following R code

The code (1 | Participant) indicates that we add random intercepts for the participants.

The R commands will provide you with the following output:⁷

From the output you can see that the critical test of the fixed effect Group is very similar to the one we found in the t-test analysis (t(18) = –1.02, p = .323). Under random effects, you also see the large variability due to differences in intercepts among participants (StdDev = 10.957). In fact, this estimate is close to the standard deviation observed in the raw scores (Table 2). This is because there are no other big sources of variance (differences between groups, days of testing).

The following command from the package effectsize can be used to obtain the standardized effect size η²:⁸

Jamovi gives the same output in an intuitive way (Figure 5). To obtain it, upload the package gamlj (Gallucci, 2022) from the jamovi library and run the mixed model, as shown below. You must enter Score as the dependent variable, Group as a categorical factor, and Participant as a cluster variable. In addition, under Random Effects, you must select the intercept of the participants.

Figure 5

If we look at the output, we see that the t- and p-values are reassuringly the same as before (t = 1.02, p = .323). We also see an entry called R-squared Marginal (=.0498). This values corresponds to η² and R² in the analyses we did before. So, all necessary information is present in jamovi.

However, there is an important difference between the default choices made in gamlj and in lme4. This can be seen in the value of the intercept in the fixed effects table. In gamlj, the intercept is 37, whereas in lme4 and lmerTest it is 39.5. In both analyses, the difference between the groups is –5. What is going on here? And why is it important?

The R packages lme4 and lmerTest use dummy coding as the default (like R in general). They assign 0 to the first condition (here 65yr) and 1 to the second condition (here 75yr). As a result, the mean score of the first condition corresponds to the intercept of the model (M_65yr = 39.5 + 0 * difference between conditions = 39.5) and the mean score of the second condition corresponds to the intercept plus the difference between the conditions (M_75yr = 39.5 + 1 * difference between conditions = 39.5 – 5 = 34.5).

In contrast, the jamovi package gamlj uses sum coding as default with contrast weights of –.5 and +.5. As a result, the intercept corresponds to the overall mean of the data set (i.e., (39.5 + 34.5)/2 = 37). The mean of condition 65yr corresponds to the intercept – 0.5 times the difference between conditions (= 37 + (–.5) * (–5) = 39.5) and the mean of condition 75yr corresponds to the intercept + 0.5 difference between conditions (= 37 + (+.5) * (–5) = 34.5).

You can get the same R output in gamlj by changing the default value of “Factors coding” from “simple” to “dummy”. You can also change the R code to use sum coding instead of dummy coding, as we will see below.

The distinction between dummy coding and sum coding may seem trivial now, but later in this article we will see that sum coding is needed to get useful interpretation for the main effects when your model includes an interaction between the independent variables!

LME offers a further advantage over traditional methods in dealing with missing data. LME can efficiently deal with missing values, provided they are not excessive and randomly distributed across conditions, so LME results obtained with missing data are very similar to those from full data sets. In contrast, t-tests and ANOVAs have more stringent requirements and require observations in all cells. This often leads to the exclusion of participants with missing data or the need for imputation techniques, as described by Austin et al. (2021).

Hopefully our example has shown you how LME relates to simpler tests. Three advantages of LME are: (1) you do not need to calculate means (you can simply use all values in long format like in Table 4), (2) the model deals efficiently with missing observations, and (3) the model can be extended to more than one random variable, as we will see in Part 2.

First, we will discuss two more designs to set the stage.

1.2.1. Analysis with t-test and ANOVA

For a t-test of paired observations, we need (1) a single score per participant per condition, and (2) the difference between the conditions for each participant. First, we take the mean of each participant at each age. Second, for each participant, we compute the difference between the two ages. Note that by computing the differences, performance difference between participants are removed. Table 6 and Figure 6 show what the data look like. They show that the average memory score of the participants was higher when they were 65 years old (M = 37.8) than when they were 75 years old (M = 34.5). Is the difference statistically significant?

Figure 6

The t-test for related or paired observations is calculated as follows (N = number of participants):

The test has N – 1 = 9 degrees of freedom and is significant: t(9) = 3.23, p = .010, two-tailed. The test is significant despite the small number of participants, because the effect size is large. We can calculate the effect size as follows:

The effect size expresses the expected average difference between the two conditions for a participant, relative to the variability in difference scores. It is traditionally called d_z.

The reason why the effect size and t-value are large is that most participants have greater working memory capacity at age 65 than at age 75. Therefore, the standard deviation of the difference scores is small (Table 6).

The effect size d_z builds on the assumption that differences in overall performance between participants are not important. Only the differences between conditions matter. Note that this assumption was the reason why a repeated-measures design was used in the first place. In some situations, however, it is interesting to also compare the condition effect to the overall variability in the sample. This is the case, for example, when we conduct a meta-analysis over studies that sometimes manipulated the variable within participants and sometimes between participants. So, we can define a second effect size, d_av, which takes all variability into account. Because this effect size includes all variability, the value of d_av is very similar to the d-value we obtained in our previous analysis of independent samples:

The difference between d_z and d_av is a headache for researchers running meta-analyses on topics in which both between-groups and repeated measures designs are used, because the statistics provided in the manuscripts often do not suffice to calculate both. You can help them by including this information in your text (see Fernández-Castilla et al., 2024, for guidelines about how to report data so that they can easily be included in meta-analysis).⁹

To run the t-test for related samples in R or jamovi, we use the format of Table 7.

The R code then is the following (if you have named your dataset Table_7).

This gives the following output:

The output tells us that the confidence interval of the effect size goes from d = .23 to d = 1.78; that is from a small effect size to a very large effect size. This informs us that the uncertainty about the effect size is high because of the small sample size used.

The jamovi output for the t-test is shown in Figure 7.

Figure 7

For the present discussion, it is important to understand what is done in a t-test for related conditions. First, the overall differences in performance between participants are left out of consideration by using the difference scores. Second, the t-test looks at whether the average difference between conditions differs from zero given the variability in difference scores. The latter is an important insight, because it means we not only consider the overall performance of the participants (their intercepts) but we also consider how participants vary in their differences or change between conditions. The change or difference between conditions is also called the slope. So, in LME we will have to enter both random intercepts and random slopes by participants.

Before discussing the LME analysis, we first analyze the data of Table 7 with an ANOVA for repeated measures to get a better feel for the data and the various statistics that can be calculated. In R we can do this with the code below (which incidentally shows that the basic R code is not always easy to remember; there are R packages that are more user-friendly):

The output of the same analysis in jamovi is shown in Figure 8.

Figure 8

The output of an ANOVA for repeated measures shows that a difference is made between within-subjects effects and between-subjects effects. The variance between subjects refers to the differences in overall performance between the participants (i.e., the participants’ intercepts). The variance within subjects refers to the age effect and the variability in the age effect (the participants’ slopes).

The jamovi output also gives two standardized effect sizes: η² (eta squared) and η²_p (eta squared partial).¹⁰ η² tells us how much of the total variance is explained by Age. This is only 2.5%. You obtain this value when you divide the Age sum of squares by the total of all three sums of squares in the ANOVA table, as can be seen in the equation below:

The value of η² is low because participants differ a lot in working memory capacity, as can be seen in Figure 6. This variability is captured by the between subjects sum of squares, which amounts to 2041, much more than the variability associated with the Age effect (54.5) and the variability in Age effect between participants (47.0).

Here, to some extent, η² is not very interesting, because we used a repeated measures design, exactly to deal with the problem that people differ in working memory capacity. We retested a group of participants 10 years later to see whether most of them would show a decline, independent of the capacity they had at the start of the study. For this, η²_p is more informative. You obtain it by dividing the Age sum of squares by all variance within-subjects:

The value η²_p gives the percentage of variance in the difference scores (or slopes) explained by Age. It tells us that more than half of the variability is accounted for by Age. This is a very big percentage, as psychological variables rarely explain more than 10% of variance (Bosco et al., 2015; Lovakov & Agadullina, 2021; Schäfer & Schwarz, 2019; Sladekova et al., 2023).

Further interesting is that η² gives the same information as d_av, whereas η²_p gives the same information as d_z. Indeed, you can translate one effect size into the other:¹¹

You may wonder why the equation of d_av includes a multiplication by 2, whereas that of d_z does not. This is because two independent distributions are compared in a between-groups comparison (the distribution of the 65yr and that of the 75yr), whereas only one distribution is considered in repeated measures (that of the difference scores between 65yr and 75yr). Still, the different relationship between d_av and η² and between d_z and η²_p is the origin of many misunderstandings in the literature. It is also the reason why a d effect size does not work in LME, because an effect can be between-conditions for one random variable but within-conditions for another random variable, as we will see shortly. Therefore, we argue it is better to stick to eta squared in LME rather than relative differences in means, just as we work with R² in linear regression.

You can get the values of η² and η²_p in R via the package effectsize:

1.2.2. Analysis with LME

Now that we have covered the basics of an analysis for repeated measures, let’s see what the LME analysis looks like. Remember that for this analysis we look at the entire dataset, not at the averages of two measurements (days tested). So, there will be some more variance in our dataset. Also remember that we must include both random intercepts and random slopes for the participants. We will start with the gamlj analysis in jamovi, as here we can use simple clicks to see what comes out. The input for LME must be in long format, as shown in Table 8.

We open mixed model in gamlj (under the heading Linear Models in jamovi). We enter Score as the dependent variable and Age as the independent variable. Because age is a categorical variable (consisting of two separate age groups), we put it under factor (and not under covariates). In addition, we must enter that the intercept and the Age slope are random variables by participants (each participant has their own performance level and shows variability in the size of the Age effect). If we do everything correctly, we get the output shown in Figure 9.

Figure 9

Reassuringly, the t-value for Age is similar to the t-test we ran: t(9) = –3.23, in line with the fact there were few differences between the two days of testing. Notice that the sign of the t-test is negative here whereas it was positive before. This is because the sign depends on how the conditions are compared to each other (65yr–75yr or 75yr–65yr). Therefore, it is necessary to always provide information about the means and SDs of the conditions when you report the outcome of a statistical analysis so that readers can check that you are drawing the right conclusions!

Further notice that the jamovi output gives information about the random effects. These tell you that the variability in intercepts is much larger than the variability in slopes and the residual noise of measurements.¹² Because jamovi coded Age as –.5 and .5, SD of the random slope can be interpreted as the variability in the score difference between the two age groups.

The jamovi output further tells us that R²_marginal = .0230, which is lower than the η² = .025 we obtained in the ANOVA. This is because the differences between the measurements are added to the total variability in the dataset. As a matter of fact, this is better because if we want to know the percentage of total variance explained, it is not a good idea to average values. Averaging reduces the variance observed in the original scores, so that the value of η² only applies to studies in which the same number of observations have been averaged (Brysbaert & Stevens, 2018).

To conduct the LME analysis in R, we use the code lines below. The instruction (Age | Participant) indicates that both random intercepts and random slopes must be included for participants.

To get the value of η² in LME, we can use the package r2glmm (Jaeger, 2022), which calculates marginal R², the equivalent of η². Because the package uses an algorithm pioneered by Nakagawa and Schielzeth (2013) and extended by Johnson (2014), the method is called nsj.

All in all, we have seen how we can analyze a repeated measures design with LME. It requires you to have more than 1 observation per condition per participant (often the case in research) and has the same advantages as we saw for the analysis of independent groups: (1) you do not need to calculate means (you can simply use all values in long format like in Table 8), (2) the model deals efficiently with missing observations, and (3) the model can be extended to more than one random variable (Part 2).

1.3.1. Analysis with ANOVA

Before starting an analysis, it is important to have a look at the data. Figure 10 shows mean performance as a function of Day and StimulusType. This shows that on Day 1 performance was higher for Stimulus type 2, whereas on Day 2 performance was better for Stimulus type 2. There was no big difference in performance between Day 1 and Day 2.

Figure 10

The most familiar analysis of an experimental design with several independent variables is ANOVA. For such an analysis, we have to average over measurement, because ANOVA requires a single observation per participant per condition.¹³ The outcome is shown in Table 10. This table also includes the means and the standard deviations of the conditions, showing the interaction between Day and Stimulus type.

To run a repeated measures analyses in R it is better to work with a long format, with 40 lines of data and extra columns for Day (Day1 vs. Day2) and StimulusType (Stim1 vs. Stim2). The code to do so can be found in the R file on osf. Once we have the right format, this is the code to run the repeated measures ANOVA in R (data from Table_10long that was made in R):

The output is:

The output confirms that there are no significant main effects, but that the interaction is significant: F(1,9) = 8.471, p = .0173.

The effectsize package (Ben-Shachar et al., 2023) allows us to easily calculate eta squared (similar to d_av) and partial eta squared (similar to d_z) for ANOVAs. First we calculate η². Remember that this is the percentage of variance explained by each effect relative to the total variance present in the data (including the variance in overall performance between participants).

The values of η² tend to be low, because they are calculated relative to the variability in the entire dataset. This is informative in some respects (e.g., to know how much variance would be explained in a between-subjects design when differences in overall performance between participants are important as well). However, it is uninformative in other respects. Researchers use a repeated measures design precisely to disregard overall differences in performance. So, an alternative measure is one that ignores the intercepts (i.e., average performance differences) of the participants.

One would hope that partial eta squared does this, like in the design with one independent variable we saw earlier, but this is not the case. Let us first see what values we obtain, using the package effectsize:

Jamovi gives the same information if you check the right boxes, as can be seen in Figure 11.

Figure 11

We can see that the values of η²_p are very high in multifactorial designs with repeated measures. Whereas the interaction between Day and Stimulus Type accounted for 7.91% of the total variance, it suddenly seems to account for a whopping 48.48% of the variance that matters. Even worse, the η²_p indices of the main effects and the interaction effect together seem to account for 12.43 + 4.00 + 48.48 = 65% of the variance. Indeed, in many multifactorial studies with repeated measures, the sum of η²_p values exceeds 1 (Norouzian & Plonsky, 2018), which is impossible if η²_p gives information about the percentage of variance accounted for by the factors we are investigating. What is going on?

The issue with η²_p is that the percentage of variance is not calculated against the total within-variance, but only against the error term associated with the effect itself (i.e., after taking away all other possible sources of variance). So, η²_p of Day is obtained as follows:

The same is true for the other two effects (Stimulus Type and the interaction). The calculation of η²_p is correct if we want η²_p to convey the same information as d_z (or the t-value for that matter) but the measure tells us little about the importance of a variable to explain variance in the phenomenon we are examining. Even worse, an effect that explains more variance can have a lower η²_p than an effect that explains less variance, because the denominators are different.

Because η²_p does not give information about the percentage of variance accounted for, no researcher really pays attention to η²_p values (Norouzian & Plonsky, 2018; Richardson, 2011). The value is duly reported, because editors and reviewers ask for it, but it actually puts everyone on the wrong foot. If we look at studies investigating correlations directly, we see that the vast majority of correlations between related variables in psychology are between r = .10 and r = .40 (Bartos et al., 2024; Bosco et al., 2015; Götz et al., 2022; Lovakov & Agadullina, 2021; Schäfer & Schwarz, 2019; Sladekova et al., 2023). This corresponds to R² or η² values between .10² = .01 and .40² = .16. Still, in multifactorial experiments with repeated measures, we regularly see η²_p values of more than .40 (suggesting a correlation of more than r = .63 = sqrt(.40)).

Partial eta squared has been made popular by SPSS (Richardson, 2011) and seems to be added to manuscripts as part of a self-imposed ritual rather than because it informs readers.

We think it is time to correct the situation. In experimental designs like the one in the example, what we want to know from an eta squared measure is how much meaningful variance it explains, similar to what happens when we calculate η²_p in an ANOVA with one within-participants factor (section 1.2.1). Looking at the ANOVA table, we need to compare the variance explained by an effect to the total variance observed within subjects, disregarding the variance between participants. The total variance within subjects is obtained by totaling all six sums of squares in the within subjects panel of the ANOVA table. We call the ratio of the sum of squares due to an effect relative to the total sums of squares within participants eta squared within (η²_w).

For the example we are discussing, this gives the following values:

Together, the three effects account for 30.3% of the variance within participants. The rest is noise (residual variance) in the within-participant data.

While ANOVA tables in R or jamovi make it easy to exclude between-participants variance, current implementations of LME packages lack such functionality. Fortunately, there is an alternative approach to neutralize differences between participants. It consists of centering the data per participant. We do so by subtracting the participant’s mean from their data points. The procedure is illustrated in Table 11.

To illustrate the procedure of calculating participant centered values, we take observation d1s1 from Participant 1. This observation has a value of 3. The mean value of Participant 1 over the four within conditions is 3.625. Subtracting the mean from the observed value gives the participant centered value of 3 – 3.625 = –.625. The same procedure is applied to all other values in Table 11.

Now we can look at what an ANOVA gives for participant centered values. Figure 12 shows the jamovi output.

Figure 12

Three things are noteworthy:

1. All the values in the Within Subjects Effects panel except η² are the same.

2. The sum of squares in the Between Subjects Effects panel is 0.¹⁴

3. The η² values in the Within Subjects Effects panel now match the η²_w values we calculated.

All information in Figure 12 is as it should be, because we reduced the sum of squares between subjects from 87.3 in the original analysis to 0 in the analysis with participant centered values, but everything else remained the same. In the next section we will see that we can use the same procedure to get the η²_w output from LME models.

When we have a significant interaction, we are likely to calculate the effect sizes of the most interesting pairs of conditions with the equations we have discussed under 1.2 (paired conditions). So, for the present study, we may want to test how large the difference between the Stimulus types is on Day 1 and on Day 2.

1.3.2. Analysis with LME

To run an LME analysis we need Table 9 in long format. This table will have 80 rows (10 participants * 2 days * 2 stimulus types * 2 measurements) and five columns (Participant, Day, StimulusType, Measurement, Score).

For the LME analysis, we have three fixed effects: the main effect of day, the main effect of stimulus type, and the interaction. In addition, we have four random effects: intercept per participant, slope of day per participant, slope of stimulus type per participant, and slope of the interaction effect per participant (remember that we need random slopes for all within effects).

We start with the jamovi output (Figure 13), because it is simpler.

Figure 13

For better visibility, we repeat the table showing the estimates of the fixed effect parameters.

We can easily verify that the statistics of the fixed effects are very similar to what we obtained in the ANOVA and t-tests on the averaged data.

In addition, we get information about the random effects.

The table with random effects tells us that in our toy example all random intercepts have variance associated with them, but that the variability of Stimulus type is rather small.¹⁵

Finally, we are informed that the fixed effects together explain 7.2% of the variance (R² marginal = .0725).

We discussed the jamovi output first, because lme4 and lmerTest have an unpleasant surprise in store for the uninformed user. The commands one would expect for the LME test are:

This gives the following output:

As you can see, the output of these commands does not at all correspond to what we have seen so far and what we would expect on the basis of the ANOVA! All effects seem to be significant, also the main effects, even though it is clear from Figure 10 that there are no differences between Day 1 and Day 2, or between Stimulus type 1 and Stimulus type 2 (as we also saw in the ANOVA analysis). Only the p-value for the interaction is the same. What is going on here?

The above analysis is an illustration of the observation that software packages sometimes include default choices that deviate from users’ intuitions (Brehm & Alday, 2022). As we saw before, the packages lme4 and lmerTest use the default dummy coded contrasts (0 and 1) for categorical variables. In the case of two levels, one level is chosen as the 0 value (in our example, these are Day 1 and Stimulus type 1) and the other as 1 (Day 2 and Stimulus type 2). What the programs then do is calculate whether there is a main effect in the 0 reference condition (i.e., is there a main effect of Day for Stimulus type 1, and is there a main effect for Stimulus type at Day 1?). As you can see in Figure 10, these specific main effects are indeed present, but they are absorbed by the wider interaction.

Dummy coded contrasts is relevant in regression-like situations, but for designs with orthogonal factors in experimental research, they must be avoided. In order for the analysis to correctly answer our research questions, the reference point must be in the middle of Figure 10: between Day 1 and Day 2, and between Stimulus type 1 and Stimulus type 2. This is achieved by sum contrast coding. So, we can try the following code:

This solves the issue of the t-values and p-values, but halves the regression weights of the main effects, because the default contrast weights in lme4 (and R in general) are –1 and +1. For the interaction, the effect estimate is even only one quarter of the difference in the data. Therefore, the jamovi gamlj option is more convenient. ¹⁶ You get it in R with following code:

Now the effects match what we expect.

Package r2glmm can be used to get the values of η².

To get the values of η²_w, we need to work with participant centered scores. You can calculate these values with Excel as we did in Table 11, but you can also use the following command in R to calculate participant centered scores in the long format table:

For the LME analysis with centered scores, one would think that it makes sense to take out the participant intercept as random variable, because there are no differences between participants anymore. However, current implementation of lme4 will then calculate separate random effects for each level of the centered variable, which may be even worse. So, better to leave the intercept in the model and ignore it. Having the intercept in the model does not make any difference in the calculation of η²_w, because that value is determined relative to the total variance, independent how it is allocated over components. So, the only thing that must be changed in the model is the dependent variable, which is the centered score instead of the original score:

This gives us the following values of η²_w:

The η²_w values of the LME analysis are lower than those of the ANOVA analysis (together explaining 17.5% of the within variance instead of 30.3%). This is because in the ANOVA analysis we eliminated within-subjects variance by taking the mean of two measurements. The η²_w values of the LME analysis are more informative because they remain the same regardless of the number of observations per participant per condition, while those of the ANOVA analysis typically increase as more observations are averaged. As a result, the standardized effect sizes based on ANOVA analysis may differ considerably between studies (Brysbaert & Stevens, 2018).¹⁷

While examining the various standardized effect sizes that have been proposed, we encountered a further problem with the η²_p values, as currently calculated by the R package effectsize (and probably by other packages as well). The issue is as follows. The effectsize package calculates η²_p based on the F-values in an ANOVA table. This is valid for ANOVA and simple regression, but not always for LME, because the degrees of freedom in LME are not well defined. This became clear when we calculated η²_p for the current example. The η²_p estimate of the interaction was very different if it was based on a model with random Day * StimulusType slopes (effect size is η²_p = .48) than if it was based on a model without random Day * StimulusType slopes (same effect size is η²_p = .20), even though the two models did not differ significantly in how well they fitted the data (see the accompanying R code for more information). Note also the unrealistically high estimate of the effect size for the interaction. A value of η²_w = .152 tells us that 15.2% of the variance within participants is explained by the interaction (equivalent to a correlation of r = .39), while a value of η²_p = .48 falsely suggests that the effect size would be similar to a correlation of r = .69. This is a further reason why researchers should not use η²_p values as estimates of effect sizes.

A second take-home message for LME analyses with interactions is to remember that you must make sure that the independent variables are centered on the mean (sum coding). This also applies to regression analysis if the analysis includes interactions between predictors. Otherwise, you are likely to misinterpret the outcome of your statistical analysis!

It is wise always to make a plot of your data, so that you can check whether the outcome of the analysis aligns with the data. Adding a plot in gamlj/jamovi is easy, as you only have to check a box to get a figure added. For a 2 × 2 design, you can use the package CGPfunctions. This is the code to use to get a picture like Figure 10.

2.2.1. Separate analysis by participants and by stimuli with LME

First, it may be interesting to remember that we can use LME analyses as an alternative to the t-tests by participants and stimuli described in 2.1. All we need to do is to run one analysis with participants as random variable and one analysis with stimuli as random variable.

For the analysis by participants (between-subjects variable, so no random slopes), this gives the following output (compare with the t-tests above):

Notice that η² is much lower than in the F1 ANOVA (where η² = .239). This is because in the ANOVA analysis we averaged over the 10 stimuli. The value of η² in LME gives the percentage of variance accounted for by Age in the ratings to the pictures. The value of η² in the ANOVA gives the percentage of variance accounted for by Age in average ratings over 10 pictures.

To some extent, η² is an underestimate of the effect size, because the same pictures were given to both groups of participants and we are not interested in differences between the pictures. So, it makes sense to calculate η²_w with stimulus centered scores. This is done as follows:

As expected, the value of η²_w is in-between the η² values obtained in the LME analysis and in the ANOVA. If we calculate η²_p with the package effect size, we get the same η²_p as for the ANOVA (η²_p = .239). It is an indication of the significant t-value obtained, but does not inform us about the size of the Age effect.

For the analysis by stimuli (repeated measure, so both random intercepts and random slopes), the LME analysis looks as follows:

Because the data set is the same in the analysis by stimuli as in the analysis by participants, the values of η² and η²_w are very much the same. For comparison, η² and η²_p in the F2 ANOVA were .067 and .264.

2.2.2. Combined LME analysis with both participants and stimuli as random variables

The analyses by participants and stimuli were interesting, because we could compare them with the F1 and F2 ANOVAs. However, a strength of LME is that it allows you to run combined analyses over both participants and stimuli, allowing you to generalize across both. If only participants are considered as random variable (and data are averaged across stimuli), we run the risk of drawing conclusions that are valid only for the specific stimuli used. If only stimuli are considered as random variable (and data are averaged across participants), we run the risk of drawing conclusions that are valid only for the specific participant group tested. In addition, the effect sizes will be exaggerated because they only apply to the average scores and not to the dependent variable itself.

As it happens, an LME analysis with both participants and stimuli as random variables is not more difficult to run than an LME analysis with one random variable. For the study of face ratings we have been analyzing, the following commands suffice:

(1 | Participant) indicates that we only include a random intercept for participants (between groups variable). (Age | Stimulus) indicates that we include both random intercepts and slopes for stimuli (repeated measure). This gives us the following output:

The output tells us that the effect of Age is not significant (t(14) = 1.549, p = .144). The evidence is not strong enough to assume that the effect generalizes to both new participants and new stimuli. We further see that there are three random effects (intercept of participant, intercept of stimulus, and variation of Age slope across stimuli), in addition to residual variance (the remaining variability in the data).

The same outcome is obtained with gamlj/jamovi (Figure 15), which has the advantage that the interface is more intuitive and that it makes convenient default choices (e.g., the use of –.5 and +.5 coding, which does not halve the difference between the means):

Figure 15

The package r2glmm informs us that η² = .032 and that η²_w = .053. These values agree quite well with the analysis by stimuli, because there were no large differences between participants.