• Go to github.com/jensroes/bristol-ws-2025
• Click on: Code > Download ZIP > unzip directory on your machine.
• Open this project by double clicking on bristol-ws-2025.Rproj

• Unified framework of advanced statistical techniques that covers a wider range of data-analysis problems.
  • Linear regression and General linear models
  • Generalised linear models: transformation allows modelling of data-types including categorical, ordinal, count data.
  • Multi-level generalised linear models: data with hierarchical structure (repeated measures).
• Backbone of other advanced models: machine learning, time series models, path analysis, structural equation models, factor analysis, signal detection models.

We can use lm to implement such a model.

model <- lm(y ~ x, data = data)

We can simulate data that exactly meet model assumptions (normal distribution, equal variance, independence etc.) with known model parameters.

# Set parameter values
n <- 1000
beta_0 <- 500 
beta_1 <- 50
sigma <- 100

# Random data for group 1
group_1 <- rnorm(n = n/2, mean = beta_0, sd = sigma)
# Random data for group 2
group_2 <- rnorm(n = n/2, mean = beta_0 + beta_1, sd = sigma)

# Generate data
sim_data <- tibble(group_1 = group_1,
                   group_2 = group_2) 

# Preview data
glimpse(sim_data)

Rows: 500
Columns: 2
$ group_1 <dbl> 560, 495, 485, 546, 291, 402, 397, 450, 355, 522, 415, 347, 41…
$ group_2 <dbl> 615, 640, 514, 592, 564, 496, 632, 678, 510, 565, 532, 677, 63…

# Change data format
sim_data <- pivot_longer(sim_data, cols = c(group_1, group_2), names_to = "x", values_to = "y")

# Preview data
glimpse(sim_data)

Rows: 1,000
Columns: 2
$ x <chr> "group_1", "group_2", "group_1", "group_2", "group_1", "group_2", "g…
$ y <dbl> 560, 615, 495, 640, 485, 514, 546, 592, 291, 564, 402, 496, 397, 632…

# As a reminder, these are the parameter values
beta_0 <- 500 
beta_1 <- 50
sigma <- 100

lm uses maximum likelihood estimation to determine for which set of parameter values are the data most likely.

# Normal linear model of simulated data
model <- lm(y ~ x, data = sim_data)

coef(model) # Model coefficients

(Intercept)    xgroup_2 
      501.0        45.1 

sigma(model) # Standard deviation

[1] 96.9

Exercise: complete script normal_model_simulation.R

• By using linear models to estimate parameter values from data we make assumptions about the process that generates the data.
• For simulated data we know the process that generated the data.
  • normal distribution because rnorm samples normally distributed data
  • both groups have the same variance because sigma was the same.
• Even if we know that assumptions are met, parameter estimates obtained are not identical but similar to the true parameter values.
• In real life we don’t know the true parameter values or the underlying process that generates the data.

data <- read_csv("../data/martin-etal-2010-exp3a.csv")

# A tibble: 1,036 × 4
   item   ppt    rt nptype   
  <dbl> <dbl> <dbl> <chr>    
1    11    10  1055 conjoined
2    11     9  2010 conjoined
3    11    11   461 conjoined
4    11     6   977 conjoined
5    11     7  1152 conjoined
# ℹ 1,031 more rows

summarise(data, across(c(ppt, item), list(n = ~length(unique(.))), .names = "{col}"))

# A tibble: 1 × 2
    ppt  item
  <int> <int>
1    12    48

In standard repeated measures ANOVAs only one source of variance can be modelled at a time.

For example one would aggregate across items and by-participant neglecting the by-items variance.

data_ppt <- summarise(data, rt = mean(rt), .by = c(nptype, ppt))

# A tibble: 24 × 3
  nptype      ppt    rt
  <chr>     <dbl> <dbl>
1 conjoined     1 1382.
2 simple        1 1323.
3 conjoined     2 1077.
4 simple        2 1108.
5 conjoined     3 1162.
# ℹ 19 more rows

Each rt observation \(i \in 1 \dots N\) comes from a normal distribution with a mean \(\mu\) and a standard deviation \(\sigma^2\).

\[ \text{rt}_i \sim \mathcal{N}(\mu_i, \sigma^2)\\ \]

Using the linear regression equation, the mean \(\mu\) can be decomposed into

\[ \mu_i = \beta_0 + \beta_1 \cdot \text{nptype}_i + \text{ppt}_i\\ \]

where \(\text{nptype}_i\) is dummy coded (by default as treatment contrast in alphabetical order)

\[ \text{nptype}_i = \left\{ \begin{array}{ll} 0, \text{ if nptype}_i = \texttt{conjoined}\\ 1, \text{ if nptype}_i = \texttt{simple}\\ \end{array} \right. \]

\(\beta_0\) is the average rt for conjoined NPs because

\(\beta_0 + \beta_1 \cdot 0 = \beta_0\),

the average rt for simple NPs is

\(\beta_0 + \beta_1 \cdot 1 = \beta_0 + \beta_1\),

and therefore \(\beta_1\) is the difference in rts between conjoined and simple NPs.

Finally \(\text{ppt}_i \sim \mathcal{N}(0, \sigma^2_\text{ppt})\).

Repeated measures ANOVA:

library(afex)
m <- aov_car(rt ~ Error(ppt/nptype), data = data_ppt)
summary(m)

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

              Sum Sq num Df Error SS den Df F value  Pr(>F)    
(Intercept) 28783750      1   587304     11  539.11 1.1e-10 ***
nptype          6217      1    13043     11    5.24   0.043 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Linear mixed-effects model:

library(lme4)
m <- lmer(rt ~ nptype + (1|ppt), data = data_ppt)
anova(m)

Type III Analysis of Variance Table with Satterthwaite's method
       Sum Sq Mean Sq NumDF DenDF F value Pr(>F)  
nptype   6217    6217     1    11    5.24  0.043 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exercises: complete script repeated_measures_model.R

and after that mixed_effects_model.R.

• Combination of fixed and random effects
• Allows to address dependency in the data: hierarchical (nested, multi-level) and crossed factor
• Fixed effects: systematic / deterministic effects (e.g. age, trial id, grouping variables)
• Random effects: non-systematic effects (e.g. ppts, items, class, country)
• Random effects can be nested: people within country; students within classroom within school within area etc.
• Sometimes the differences can be blurry: comparing differences between two test sites rather than treating test sites as repeated measures variable.

• Traditional repeated-measures ANOVA can only handle one source of variance (participant or item).
• Linear mixed-effects models can capture one and more (Baayen, Davidson, and Bates 2008).

m <- lmer(rt ~ nptype + ( 1 | ppt ) + ( 1 | item ), data = data)

summary(m)$varcor

 Groups   Name        Std.Dev.
 item     (Intercept)  40.6   
 ppt      (Intercept) 161.2   
 Residual             233.5   

sigma(m)

[1] 234

summary(m)$coef

             Estimate Std. Error t value
(Intercept)    1111.8       48.0    23.2
nptypesimple    -33.4       14.5    -2.3

(coefs <- m@beta) # extract estimates

[1] 1111.8  -33.4

# Predicted mean for np conjoined
coefs[1] + coefs[2] * 0

[1] 1112

# Predicted mean for np simple
coefs[1] + coefs[2] * 1

[1] 1078

# Check predictor contrasts (default: treatment)
contrasts(factor(data$nptype))

          simple
conjoined      0
simple         1

t-values \(> \mid 2 \mid\) correspond to \(\alpha = .05\) (Baayen 2008):

summary(m)$coef

             Estimate Std. Error t value
(Intercept)    1111.8       48.0    23.2
nptypesimple    -33.4       14.5    -2.3

Confidence intervals: 95% CI shows the range of hypotheses that can’t be ruled out for \(\alpha = .05\):

confint(m, "nptypesimple")

             2.5 % 97.5 %
nptypesimple -61.8  -4.93

Satterthwaite approximation for degrees of freedom (lmerTest):

m <- lmerTest::lmer(rt ~ nptype + ( 1 | ppt ) + ( 1 | item ), data = data)
summary(m)$coef

             Estimate Std. Error    df t value Pr(>|t|)
(Intercept)    1111.8       48.0  11.9    23.2 3.11e-11
nptypesimple    -33.4       14.5 977.9    -2.3 2.16e-02

Likelihood-ratio test:

m_null <- lmer(rt ~ 1 + ( 1 | ppt ) + ( 1 | item ), data = data)
anova(m_null, m)

Data: data
Models:
m_null: rt ~ 1 + (1 | ppt) + (1 | item)
m: rt ~ nptype + (1 | ppt) + (1 | item)
       npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)  
m_null    4 14319 14339  -7155    14311                      
m         5 14316 14340  -7153    14306  5.28  1      0.022 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Random intercepts: size of np-type effect is constant.
• Random slopes: size of np-type effect varies by participants and items.
• Correlation between random intercepts and random slopes: e.g. np-type effect might be weaker for participants / items with longer rts.
• Maximal random effects structure (Barr et al. 2013): for repeated-measures designs with crossed participants and items, we need random intercepts for participants and items with by-participants and by-items random slopes.

# Random intercepts model
m <- lmer(rt ~ nptype + ( 1 | ppt ) + ( 1 | item ), data = data)

# Maximal random effects structure 
m_max <- lmer(rt ~ nptype + ( nptype | ppt ) + ( nptype | item ), data = data)

# Maximal random effects model 
summary(m_max)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: rt ~ nptype + (nptype | ppt) + (nptype | item)
   Data: data

REML criterion at convergence: 14288

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.505 -0.596 -0.147  0.449  5.002 

Random effects:
 Groups   Name         Variance Std.Dev. Corr 
 item     (Intercept)   2749.8   52.44        
          nptypesimple   522.1   22.85   -1.00
 ppt      (Intercept)  27159.7  164.80        
          nptypesimple    52.4    7.24   -1.00
 Residual              54359.2  233.15        
Number of obs: 1036, groups:  item, 48; ppt, 12

Fixed effects:
             Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)    1111.9       49.3   11.5   22.58    7e-11 ***
nptypesimple    -33.6       15.0  125.7   -2.24    0.027 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
nptypesimpl -0.310
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

• Convergence failures (and other problems like floor and ceiling effects) can be (easily) addressed in Bayesian models.
• Bayesian model converge by definition (as the number of iterations approaches \(\infty\)).

Frequentist mixed-effects models

library(lme4)
fit_lmer <- lmer(rt ~ nptype + ( nptype | ppt ) + ( nptype | item ), data = data)

Bayesian mixed-effects models

library(brms)
fit_brm <- brm(rt ~ nptype + ( nptype | ppt ) + ( nptype | item ), data = data)

Two different philosophical frameworks of generalising from sample to the population: What do the data tell us about “the truth”?

• Null-hypothesis significance testing
• Bayesian inference

• Classical / frequentist statistics; p-value based statistics
• p-value: How plausible are the data (or something more extreme) if we assume that there’s no effect?
• Evaluating the probability of data (e.g. t-value, F-statistic) assuming that the null hypothesis \(H_0\) is true.
• If the data are extreme enough, \(H_0\) is concluded to be implausible (rejected).
• If \(H_0\) is implausible, we assume \(H_1\) (alternative hypothesis).
• Statistically significant means data are unexpected under \(H_0\).

\[ Pr(\text{data} \mid H_0 = \text{TRUE}) \]

• Expressing uncertainty about a hypothesis expressed in probability distributions
• Updating one’s belief (prior, example later) about a hypothesis (e.g. difference between conditions, interactions etc) as the new information becomes available.

\[ Pr(H \mid \text{data}) \]

\[ \underbrace{Pr(\theta \mid \text{data})}_{\text{posterior}} \propto \underbrace{Pr(\theta)}_{\text{prior}} \cdot \underbrace{Pr(\text{data} \mid \theta)}_{\text{likelihood}} \]

• Likelihood: probability of data given the model parameter value(s)
• Prior: what do we already know about model parameter(s)
• Posterior: our updated belief in the parameter value(s) after seeing the data

Instead of concluding

\[ Pr(\text{data} \mid H_0 = \text{TRUE}) \]

we can answer

\[ Pr(H \mid \text{data}) \]

It tells us what we want to know!

• What’s the probability of an effect given the data (and what we already knew before)?
• Summary stats of a Bayesian model are often very similar to what people think frequentist quantities (p-values, CIs) mean (Nicenboim and Vasishth 2016):
• What’s the relative evidence for one model (e.g., \(H_0\)) vs. another (e.g., \(H_1\))?
• What interval contains an unknown parameter value with .95 probability?

• Estimate the posterior uncertainty over parameter values based on the observed data (Kruschke 2014).
• For this we need

1. a probability model (e.g. a normal distribution and a set of predictors)
2. priors (on e.g. means, difference, variability)

• brms uses probabilistic sampling to determine the posterior uncertainty about the parameter value.
• This involve calculating the weighted probability of the likelihood of the data given a proposed parameter value in the (potentially high dimensional) parameter space.
• This can take time.

• Prior knowledge about plausible parameter values.
• This knowledge is expressed as probability distributions (e.g. normal distributions).
• Help probabilistic sampler by limiting the parameter space.
• Small data samples are sensitive to prior information which makes intuitively sense.
• Otherwise data typically overcome the prior (automatic Ockham’s razor).
• Less common: test data against a (prior) effect suggested by the literature.

• brms uses default priors that can be changed as appropriate.
• Some defaults are good, (flat) priors are worth specifying.

Check defaults used earlier:

prior_summary(fit_brm)

Exercise: complete script mixed_effects_model_brms_with_priors.R (this time with priors)

• To test a hypothesis about parameter values (e.g. difference between groups, interactions, change in outcome):
  • Summary of posterior parameter value
  • Region of practical equivalence (ROPE)
  • Savage-Dickey Bayes Factor: for nested models
• Cross-validation (leave-one-out): comparison of different probability models

Get posterior of slope \(\beta\) (difference between conditions)

beta <- as_draws_df(fit_brm) %>% pull(b_nptypesimple)

length(beta)

[1] 4000

beta[1:5]

[1] -34.0 -45.2 -33.0 -37.6 -23.7

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

sigma <- as_draws_df(fit_brm) %>% pull(sigma)
delta <- beta / sigma

abs(mean(delta)) # very small effect (see Cohen's d)

[1] 0.142

• Determine ROPE to assess the uncertainty of effect size (Makowski, Ben-Shachar, and Lüdecke 2019)
• Calculate the standardised effect size \(\delta\) from posterior: \(\delta = \frac{\beta}{\sigma}\)
• ROPE of [-0.1, 0.1] is the range of negligible effect sizes (Cohn 1988; Kruschke 2018).

rope(delta, range = c(-0.1, 0.1)) 

# Proportion of samples inside the ROPE [-0.10, 0.10]:

inside ROPE
-----------
26.34 %    

• The ROPE value indicates the extent to which the posterior cannot rule out a negligible effect.
• A meaningful effect size should have a small proportion of posterior samples within the ROPE.

rope(delta, range = c(-0.1, 0.1)) 

# Proportion of samples inside the ROPE [-0.10, 0.10]:

inside ROPE
-----------
26.34 %    

Exercise: complete script hypothesis_testing_rope.R

• Regularising (skeptical) priors prevent the model to get overexcited by the sample (his words: McElreath 2020).
• If in doubt use weakly informative priors like student-t distributions that allow for larger tails, e.g. student_t(3, 0, 100)

• Explore hypothesis() in the exercise: hypothesis_testing_bf_1.R
• Bonus: explore the Savage-Dickey method in hypothesis_testing_bf_2.R

H: there is no difference between conjoined and simple NPs

hypothesis(fit_brm, "nptypesimple = 0")

Hypothesis Tests for class b:
          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (nptypesimple) = 0    -33.2      17.2    -66.3     0.75         NA        NA
  Star
1     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

H: simple NPs are faster than conjoined NPs

hypothesis(fit_brm, "nptypesimple < 0")

Hypothesis Tests for class b:
          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (nptypesimple) < 0    -33.2      17.2    -61.4    -5.88         36      0.97
  Star
1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

H: simple NPs are more than 10 msecs faster than conjoined NPs

hypothesis(fit_brm, "nptypesimple < -10")

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (nptypesimple)-(-10) < 0    -23.2      17.2    -51.4     4.12       11.1
  Post.Prob Star
1      0.92     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

“Models are devices that connect theories to data. A model is an instantiation of a theory […]” (Rouder, Morey, and Wagenmakers 2016)

• Our models describe how we understand reality.
• Model allows us to describe reality qualitatively (with their parameters) and quantitatively (parameter value estimates).
• Interpretation of parameters depends on data-modeling context.
• At minimum, distribution families (probability models) depend on data type.

• Gaussian: data come from normal distribution

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = gaussian())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = bernoulli())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = poisson())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = zero_inflated_poisson())

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)
• Cumulative: ordinal (ordered categorical) data (Bürkner and Vuorre 2019) (Likert scales)

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = cumulative()) # or acat, sratio

• Gaussian: data come from normal distribution
• Bernoulli: binomial outcomes (yes / no; correct / incorrect)
• Poisson: count data (number of words / mistakes) (Winter and Bürkner 2021)
• Zero-inflated Poisson: count data with lots of zeros (number of typos per word)
• Cumulative: ordinal (ordered categorical) data (Bürkner and Vuorre 2019) (Likert scales)
• Lognormal: zero bound data with positive skew; also skewed / shifted (log)normal, Wiener Diffusion models etc. [for RT data; Matzke and Wagenmakers (2009)]

fit <- brm(outcome ~ predictor + (1|participant), data = data, family = lognormal())

So far we used a Gaussian model which showed a poor fit to the data.

Alternative models: skewed normal, shifted (log)normal, Wiener Diffusion models, ex-Gaussian, mixture of lognormal etc. (see Matzke and Wagenmakers 2009; Roeser et al. 2024).

\[y \sim \mathcal{N}(\mu, \sigma^2)\]

• Data generating process can be described as a normal distribution with mean \(\mu\) and standard deviation \(\sigma\).

\[ \log(y) \sim \mathcal{N}(\mu, \sigma^2) \]

• Often used to address positive skew; e.g. response times (Baayen 2008)
• Log-values are only defined if \(y \in [0, \infty]\)
• De-emphasize large values over small values.
• Models the percentage change and not absolute differences:
  • Rather than saying “keystroke intervals are 40 msecs slower”, you say “keystrokes are 7% slower”.
  • Difference of 40 msecs can be huge for fast typing but negligible for pauses around 10 secs.

\[ \log(y - \exp(\text{ndt})) \sim \mathcal{N}(\mu, \sigma^2) \]

• Variation of the lognormal distribution: a zero shift would be equivalent to lognormal
• Includes a horizontal shift.
• Useful for modeling rt data because they are typically not just
  • right-skewed
  • non-negative
• But also
  • positively shifted: there is a minimum time required for any response

Accounts for a minimum rt: The shift parameter represents the minimum possible rt, accounting for the fact that no response can be instantaneous.

Fit two other probability models, namely a lognormal and a shifted lognormal model:

• Complete script mixed_effects_models_brms_lognormal.R
• Complete script mixed_effects_models_brms_shifted_lognormal.R

• How accurately does the model predict new data.
• \(R^2\) et al.: more complex models are generally better ( \(R^2\) can’t go down).
• Overfitting: models with more parameters can make unreasonable predictions.
• Generalisations to data that were used to fit the model are over optimistic (Gelman, Hill, and Vehtari 2020).
• We will look at both fit and predictive performance of the model.

• Cross validation: performance of a model on new data to remove the overfitting problem.
• Leave-One-Out (LOO) Information Criterion (LOO-IC):
  • Train model on \(N-1\) observations
  • Predict remaining data point from training model.
  • Repeat process \(N\) times to predict every observation from a model of the remaining data.
• Adding up prediction results gives an estimate of expected log-predictive density (\(elpd\)); i.e. approximation of results that would be expected for new data.
• loo() uses the probability calculations to approximate LOO-IC (i.e. Pareto smoothed importance sampling) (Vehtari, Gelman, and Gabry 2015, 2017).

• brms isn’t much more complicated than lme4 and comes with more flexibility.
• Priors require some thinking.
• Bayesian models allow straight forward interpretation of parameter values without an arbitrary threshold for significance.

• Bürkner tutorials: Bürkner (2018); Bürkner and Vuorre (2019); Bürkner (2019); Winter and Bürkner (2021)
• Vasishth tutorials: Sorensen, Hohenstein, and Vasishth (2016); Nicenboim and Vasishth (2016); Vasishth and Nicenboim (2016)
• Excellent book-length intros: Gelman, Hill, and Vehtari (2020), McElreath (2020), Kruschke (2014), Lambert (2018), Lee and Wagenmakers (2014)
• Convince yourself that Bayesian inference is the right answer: Clayton (2021)