This example follows a blog post by Mattan S. Ben-Shachar: https://shouldbewriting.netlify.app/posts/2019-10-21-accounting-for-within-and-between-subject-effect The key idea is that group-level data does not always reflect individual-level processes.

Let’s work with the now-classic typing speed example. We take a group of 5 typists, and measure the speed of their typing (words per minute), and the rate of typing errors (errors per 100-words). Looking at the data we might get something like this:

Let’s estimate a simple linear model:

## 
## Call:
## lm(formula = errors ~ speed, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6704 -1.2742 -0.1443  1.4222  4.6534 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.490e-16  1.497e-01   0.000        1    
## speed       -9.040e-01  1.466e-01  -6.164 6.39e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.833 on 148 degrees of freedom
## Multiple R-squared:  0.2043, Adjusted R-squared:  0.1989 
## F-statistic:    38 on 1 and 148 DF,  p-value: 6.391e-09

Note how the model suggests a negative relationship between speed and errors: the faster people type, the fewer errors they make.

Now let’s break down the data by typist:

As we can see, we have two sources of variation that can be used to explain or predict the rate of errors:

  1. Overall, faster typists make less mistakes (group-level pattern).
  2. When typing faster, typists make more mistakes (individual-level pattern).

Let’s model the typist as a fixed effect.

## 
## Call:
## lm(formula = errors ~ speed + as.factor(ID), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.36716 -0.53360 -0.00212  0.42014  1.98377 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.2384     0.1327  -1.796 0.074583 .  
## speed            1.4000     0.1190  11.762  < 2e-16 ***
## as.factor(ID)2  -4.4978     0.2543 -17.688  < 2e-16 ***
## as.factor(ID)3  -0.7338     0.2163  -3.393 0.000892 ***
## as.factor(ID)4   2.6441     0.2150  12.296  < 2e-16 ***
## as.factor(ID)5   3.7795     0.1961  19.276  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7166 on 144 degrees of freedom
## Multiple R-squared:  0.8817, Adjusted R-squared:  0.8776 
## F-statistic: 214.6 on 5 and 144 DF,  p-value: < 2.2e-16

Now we see that typing speed is positively correlated with errors, controlling for typist.

Aside: We can also model these using liner mixed-effects models.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: errors ~ speed + (1 | ID)
##    Data: data
## 
## REML criterion at convergence: 355.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3055 -0.7421  0.0012  0.5969  2.7401 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 10.3393  3.2155  
##  Residual              0.5136  0.7167  
## Number of obs: 150, groups:  ID, 5
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept) -1.948e-13  1.439e+00  3.963e+00    0.00        1    
## speed        1.384e+00  1.187e-01  1.454e+02   11.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## speed 0.000

What if we want to capture both the within- and between-group patterns?

We can model these using liner mixed models, but first we need to split our predictor (speed) into two variables, each representing a different source of variance - each typist’s average typing speed, and the deviation of each measurement from the typist’s overall mean:1

library(dplyr)
data <- data %>% 
  group_by(ID) %>% 
  mutate(speed_M = mean(speed),
         speed_E = speed - speed_M) %>% 
  ungroup()

head(data)
## # A tibble: 6 x 5
##      ID  speed errors speed_M speed_E
##   <int>  <dbl>  <dbl>   <dbl>   <dbl>
## 1     1 -0.773 -1.74   -0.188 -0.585 
## 2     1 -0.144 -0.703  -0.188  0.0438
## 3     1 -0.686 -1.73   -0.188 -0.498 
## 4     1  0.560  1.17   -0.188  0.748 
## 5     1  0.214  0.316  -0.188  0.402 
## 6     1  0.179  0.392  -0.188  0.367

Let’s fit a liner mixed model and see how we can detect both patterns correctly.

fit <- lmer(errors ~ speed_M + speed_E + (1 | ID), data = data)
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: errors ~ speed_M + speed_E + (1 | ID)
##    Data: data
## 
## REML criterion at convergence: 346.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3132 -0.7550  0.0006  0.6020  2.7570 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1.9029   1.3794  
##  Residual             0.5135   0.7166  
## Number of obs: 150, groups:  ID, 5
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  3.515e-15  6.197e-01  3.000e+00   0.000    1.000    
## speed_M     -1.600e+00  6.928e-01  3.000e+00  -2.309    0.104    
## speed_E      1.400e+00  1.190e-01  1.440e+02  11.762   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) sped_M
## speed_M 0.000        
## speed_E 0.000  0.000

As we can see, the slope for speed_M is negative (-1.6), reflecting the group-level pattern where typists who are overall faster have fewer errors; whereas the slope for speed_E is positive (1.4), reflecting the individual-level pattern where faster typing leads to more errors.

We can access the estimated deviation between each subject average typing speed and the overall average:

## $ID
##   (Intercept)
## 1  -0.7955464
## 2  -0.8960398
## 3   1.2740378
## 4  -0.9118556
## 5   1.3294040
## 
## with conditional variances for "ID"

Question: Do we need a random slope?

## boundary (singular) fit: see ?isSingular
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: errors ~ speed + (1 + speed | ID)
##    Data: data
## 
## REML criterion at convergence: 355.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3056 -0.7419  0.0023  0.5954  2.7289 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr
##  ID       (Intercept) 1.032e+01 3.213124     
##           speed       5.225e-05 0.007229 1.00
##  Residual             5.136e-01 0.716659     
## Number of obs: 150, groups:  ID, 5
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 5.357e-03  1.438e+00 3.957e+00   0.004    0.997    
## speed       1.384e+00  1.187e-01 1.378e+02  11.652   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## speed 0.027 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Answer: probably not.

Question: Could we capture the between-effect with a fixed-effects model?

## 
## Call:
## lm(formula = errors ~ speed_M, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2864 -1.0748  0.0028  0.9154  3.7784 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.702e-16  1.196e-01    0.00        1    
## speed_M     -1.600e+00  1.338e-01  -11.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.465 on 148 degrees of freedom
## Multiple R-squared:  0.4915, Adjusted R-squared:  0.4881 
## F-statistic: 143.1 on 1 and 148 DF,  p-value: < 2.2e-16

Answer: It is different from the -0.9 estimated in the model m1. The above model 1 ( m1 ) returns a biased estimate, which is a “weighted average” of the within- and between-effects.

Note that standard errors differ, because the variance in the grouping structure is more accurately taken into account by the mixed-effects model.

When are individual-level the same as group-level patterns?

Experiments!
Or to be more precise, when we control the values of the independent variable. Why is this so? Because we control the values of the independent variable, the independent variable cannot be split into different sources of variance: there is either variance between subjects (the variable is manipulated in a between-subjects design) or there is variance within subjects (the variable is manipulated in a within-subjects design), but never both. Thus, although there can be huge heterogeneity in the way subjects present an effect, the average individual-level effect will be the same as the group-level effect (depending on the design).2


  1. Read more in: Hoffman, L. (2015). Time-varying predictors in models of within-person fluctuation. In Longitudinal analysis: Modeling within-person fluctuation and change (pp. 327-392). Routledge.

  2. Ignoring any differences or artifacts that may arise from the differences in the design itself, such as order effects, etc.