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:
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.
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
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.↩
Ignoring any differences or artifacts that may arise from the differences in the design itself, such as order effects, etc.↩