Note for me:

• Record the session if possible.

Plan for the next four sessions

• Revise linear model in R
  • Conceptual and statistical description of model with complete pooling
  • Linear model in R
• Intro to multilevel models
  • Identify hierarchies in data
  • Describe a fixed effects (no pooling) regression
  • Critically evaluate a fixed effects (no pooling) regression
  • Describe a multilevel model with random intercepts (partial pooling) conceptually and using statistical notation
• Multilevel models in R
  • Random intercept models in R
• Conceptual and statistical description of random slope models
  • Random slope models in R
• Model comparison
• Benefits and consequences of MLM
• Extensions to the MLM

Do i have to know how to use this approach?

• Yes

Item 3: Using a data-set of your own choice, perform and report a multilevel regression model to address a research question of relevance to psychology.

Plan for the next four sessions

• Session 1: Complete pooling model <- We are here

• Session 2: No pooling model

• Session 3: Partial pooling model (random intercepts)

• Session 4: Partial pooling model (random intercepts and slopes)

• in a nutshell: each session will progressively illustrate how we can handle data with ‘multiple levels’, starting from treating all data as identical, to treating each level of the data as being entirely separate, to sharing information across the levels

• sessions 3 and 4 are what you will be using in your analysis, sessions 1 and 2 are designed to help you understand why it matters.

Focusing on today: session outcomes

• After the session today students (you) are expected be able to:
  • Describe a linear model conceptually and using statistical notation
  • Run and interpret a linear model in R
  • Make predictions using a linear model in R
  • Critique the ‘one size fits all’ model in the presence of multiple ppts

Spot the difference

• Hierarchical linear models
• Mixed effects models
• Random effects models
• Random coefficient models
• Multilevel linear models

What are multilevel data?

• in disciplines like psychology and the wider social/behavioural sciences, data are often multilevel (aka hierarchical, nested, clustered, etc)
• essentially, this means that lower level units of analysis belong higher-level units of analysis
  • e.g., students may belong to specific school units (assuming both data points are captured)
• importantly, longitudinal data are also ‘multilevel’
  • e.g., repeated measures are nested within individuals

Packages

### packages to install / load

library(tidyverse)
library(lme4)
library(modelr)
library(lmerTest)

Linear mixed effects ‘4’

Reaction time and sleep deprivation

• Data from Belenky et al. (2003) in lme4
• 18 pps sleep deprived for 9 days
• RT taken on day 0 - 9 (10 measurements per pp)
• Does days sleep deprived affect RT?
• Data are in long format

glimpse(sleepstudy)

view(sleepstudy)

The data

• How would you analyse this data?

head(sleepstudy, 18)

##    Reaction Days Subject
## 1  249.5600    0     308
## 2  258.7047    1     308
## 3  250.8006    2     308
## 4  321.4398    3     308
## 5  356.8519    4     308
## 6  414.6901    5     308
## 7  382.2038    6     308
## 8  290.1486    7     308
## 9  430.5853    8     308
## 10 466.3535    9     308
## 11 222.7339    0     309
## 12 205.2658    1     309
## 13 202.9778    2     309
## 14 204.7070    3     309
## 15 207.7161    4     309
## 16 215.9618    5     309
## 17 213.6303    6     309
## 18 217.7272    7     309

The data

• How would you analyse this data?
  • Wright (1998) says there are four ways
    • ignore the fact that the data are clustered within a participant and treat them as statistically independent
    • calculate the average for each person and use those as the outcome (means as outcome approach)
    • create a dummy variable for each subject
    • multi level approach

The data

• visualise the data

ggplot(sleepstudy, aes(x = Days, y = Reaction)) +
  geom_point(size = 3) +
  facet_wrap(~ Subject) +
  theme_bw(base_size = 14)

The data

• How does RT change with sleep deprivation?

• Is this change the same for all participants?

• Their starting RT and their rate of change?

Reaction time and sleep deprivation

• How does RT change with sleep deprivation?

  • As days without sleep increases so does mean RT (and sd)

  • Apart from pp 335

• Is this change the same for all participants?

  • Starting RT differs (~ 200ms - 320ms)

  • Rate of change differs (some flatter/steeper than others)

What about using a RM ANOVA?

library(afex)
library(emmeans)

sleepstudy2 <- sleepstudy %>%
  mutate(
    Days_group = case_when(
      Days %in% 0:3 ~ "0-3",
      Days %in% 4:6 ~ "4-6",
      Days %in% 7:9 ~ "7-9"),
    Days_group = factor(Days_group, levels = c("0-3", "4-6", "7-9"))
  )

anova_model <- aov_car(Reaction ~ Days_group + Error(Subject / Days_group), data = sleepstudy2)

What about using a RM ANOVA?

library(afex)
library(emmeans)

anova_model

## Anova Table (Type 3 tests)
## 
## Response: Reaction
##       Effect          df    MSE         F  ges p.value
## 1 Days_group 1.72, 29.22 790.31 30.70 *** .284   <.001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

emm <- emmeans(anova_model, ~ Days_group)
pairs(emm)

##  contrast    estimate    SE df t.ratio p.value
##  X0.3 - X4.6    -35.7  8.56 17  -4.175  0.0017
##  X0.3 - X7.9    -68.0 10.10 17  -6.706  <.0001
##  X4.6 - X7.9    -32.3  7.08 17  -4.560  0.0008
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

What about using a RM ANOVA?

• How appropriate is this method?

• Error(Subject / Days_group) = a NESTED structure

  • “Days is a repeated measure within Subject”
  • Models the between-subject variability (Subject)
  • Models the within-subject variability over time (Subject:Days_group)

• ‘All models are regression models’

• Clear limitations of this approach

  • e.g., forcing groups from a continuous variable
  • assuming ppts act the same within each group

Normal linear model for one predictor

\[ y_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 x_i, \quad \text{for } i \in 1, \ldots, n. \]

• Two variables \(x\) and \(y\)
  • From our data and change or vary
  • \(x\) is the predictor (number of days without sleep)
  • \(y\) is the outcome (RT)
• \(y_i \sim N(\mu_i, \sigma^2)\)
  • Each observed value of \(y\) is a sample from a normal distribution
  • Whose mean is \(\mu_i\) and whose variance is \(\sigma^2\) (sd is \(\sigma\))
  • These are unknown, and we want to estimate them based on our sample
  • This may change depending on given levels of \(x\)

Normal linear model for one predictor

\[ y_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 x_i, \quad \text{for } i \in 1, \ldots, n. \]

• For \(x\) and \(y\) we have subscripts \(i\)
• Links back to our dataset
• Essentially \(i\) indexes a single row of sleepstudy
• \(x_i\) is the value of predictor variable for the \(i\)th observation
  • The number of days without sleep for the \(i\)th observation in the dataset
• \(y_i\) is the corresponding RT on that day
• We state that this formula is the same for all observations in the dataset
  • for \(i \in 1, \ldots, n.\)

Normal linear model for one predictor

\[ y_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 x_i, \quad \text{for } i \in 1, \ldots, n. \]

• The mean \(\mu_i\) is a sum of two coefficients \(b_0\) and \(b_1\) for a given level of \(x_i\)
  • Don’t change or are constants
    
  • \(b_0\) is the intercept
    • The average value of RT when \(x = 0\)
      
  • Average RT on Day 0 of the study
    • Note no subscript \(i\), so this is regarding the population

\[ \text{Reaction}_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 \,\text{Days}_i, \quad \text{for } i \in 1, \ldots, n. \]

Normal linear model for one predictor

\[ \text{Reaction}_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 \,\text{Days}_i, \quad \text{for } i \in 1, \ldots, n. \]

• \(b_1\) is the slope
  • How much the \(y\)-value increases (or decreases) for a 1-unit change in \(x\)
• \(b_1 x_i\) (weighted)
  • Remember \(x_i\) is the predictor variable (Days without sleep) for the \(i\)th observation
  • How much the \(y\)-value increases (or decreases) multiplied by the value of the predictor variable for the \(i\)th observation
  • In this study: How much RT increases or decreases multiplied by the number of days without sleep for the \(i\)th observation

Normal linear model for one predictor

\[ \text{Reaction}_i \sim N(\mu_i, \sigma^2), \quad \mu_i = b_0 + b_1 \,\text{Days}_i, \quad \text{for } i \in 1, \ldots, n. \]

• In this model we aim to estimate:

• The mean and variance (\(\sigma^2\), or sd \(\sigma\)) of the normal distribution from which our data are a sample

• We assume this same model for all participants

  • One model fits all

• Is this a good assumption to make?

One model fits all?

ggplot(sleepstudy, aes(x = Days, y = Reaction)) +
  geom_point(size = 3) +
  facet_wrap(~ Subject) +
  theme_bw(base_size = 14)

• Nevertheless, let’s analyse the data

  • Linear model in R with one continuous predictor

  • Estimate \(b_0\) and \(b_1\) to estimate the population mean (and standard deviation) for a given value of \(x\)

Linear Model in R

lm(outcome ~ predictor, data = dataset)

### simplest of linear regressions

m1 <- lm(Reaction ~ Days, data = sleepstudy)
summary(m1)

## 
## Call:
## lm(formula = Reaction ~ Days, data = sleepstudy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -110.848  -27.483    1.546   26.142  139.953 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  251.405      6.610  38.033  < 2e-16 ***
## Days          10.467      1.238   8.454 9.89e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.71 on 178 degrees of freedom
## Multiple R-squared:  0.2865, Adjusted R-squared:  0.2825 
## F-statistic: 71.46 on 1 and 178 DF,  p-value: 9.894e-15

Linear Model in R

• How well did the model fit the data?
  • Plot it

### model 1 plot

sleepstudy %>% 
  ggplot(aes(x = Days, y = Reaction)) + 
  stat_smooth(method = "lm", se = F) +
  geom_point(alpha = 0.5) +
  scale_x_continuous(breaks = 0:9)+
  labs(y = "Reaction Time") + 
  theme_classic()

Fit

• Residuals: difference between actual and predicted values
  • Look for symmetric distribution centered on 0
  • Model assumptions checked with plot(m1)

Fit

plot(m1, which = 2)

Fit

plot(m1, which = 1)

Fit

• spread of residuals is not equal (heteroscedastic)
• indicators of violations of normality of residuals
• some outliers -> one size may not fit all

Linear model in R: Intercept

• Coefficients: Intercept
  • Estimate (Intercept) is the intercept
  • The population mean (\(\mu_i\)) of RT when \(x = 0\)
  • Example: 251 ms
• Std. Error
  • Average distance that the observed values fall from the intercept / slope
  • How “wrong” the estimate is in units of the outcome variable

Linear model in R: Intercept

• Confidence Interval

confint(m1)

##                  2.5 %    97.5 %
## (Intercept) 238.360753 264.44946
## Days          8.023855  12.91072

• a measure of the uncertainty of our estimate

• range of values of the coefficients that would not be rejected in a null hypothesis test

• all possible suspects

• Residual standard error (actually the population standard deviation, \(\sigma\))

sigma(m1)

## [1] 47.71472

• Example: When Days without sleep = 0, population mean is 251.41 ms with sd of 47.71

\[\text{Reaction}_i \sim N(251.41_{} , 47.71^2) \quad\]

Only when Days = 0

Linear model in R: Slope

• Coefficients: Days
  • Estimate for Days is the slope (\(b_1\))
  • How much the \(y\)-value increases (or decreases) for a 1-unit change in \(x\)
  • a linear function
  • Example: RT increases by 10.5 ms when Days without sleep increases by 1
• Std. Error
  • Average distance that the predicted values fall from the regression line
• Confidence Interval confint(m1)
  • If we repeated the study many times, 95% of the time the true value of the slope would fall between these numbers
• Interpretation
  • When Days without sleep = 1:
    \(\mu = 251.41 + 10.47 \times 1 = 261.88\) ms
    
  • The standard deviation remains constant at 47.71 ms across all values of the predictor

Linear model in R

• t value
  • \(\dfrac{\text{Estimate}}{\text{Std. Error}}\)
  • Ratio of the estimate to the error in the model
• p-value
  • Probability that we would have obtained our results (or results more extreme) if the null hypothesis is true
    • i.e., \(t\)-distribution centered over 0 with 178 (\(n - 2\)) degrees of freedom
  • Probability of getting our results if in fact there is no relationship between RT and Days without sleep
  • Probability ranges between 0 and 1
    • 0 = “won’t happen”, 1 = “will happen”
  • Very low p-value here says it’s very unlikely we’d get these results if in the population there were no relationship between RT and Days without sleep

Linear model in R

• Let’s make predictions using our model
  • I’m just showing predictions for the first 3 Days for 4 participants
• Note how RTs are the same for each day for each person

sleepstudy %>%
  add_predictions(m1, var = "Predictions") %>%
  filter(Days %in% c(1, 2, 3) & Subject %in% c(308, 309, 310, 332)) %>%
  select(Days, Subject, Predictions)

##    Days Subject Predictions
## 1     1     308    261.8724
## 2     2     308    272.3397
## 3     3     308    282.8070
## 4     1     309    261.8724
## 5     2     309    272.3397
## 6     3     309    282.8070
## 7     1     310    261.8724
## 8     2     310    272.3397
## 9     3     310    282.8070
## 10    1     332    261.8724
## 11    2     332    272.3397
## 12    3     332    282.8070

Linear model in R

• Ignoring the fact that there are groups (different ppts) in the data, the model looks ok

Linear model in R

• What assumption is our model making about RT and days without sleep, particularly regarding each of the pps?
  • That everyone’s RT is the same at day 0 (intercept)
  • That subsequent days without sleep affects everyone’s RT by the same amount (slope)
  • All data points are independent
• Are these realistic assumptions to make?
• Does or should one model fit all?
• This one size fits all modeling is sometimes called complete pooling
• Simple linear regression uses complete pooling

Linear model in R

sleepstudy %>% 
  ggplot(aes(x = Days, y = Reaction)) + 
  geom_abline(intercept = coef(m1)[1],
              slope = coef(m1)[2],
              color = 'blue') +
  geom_point(alpha = 0.5) +
  scale_x_continuous(breaks = 0:9) +
  facet_wrap( ~ Subject) +
  labs(y = "Reaction Time") + 
  theme_classic()

Two primary points

• Estimates of the normal linear model provide the same RT for all participants at Day 0 (intercept)
  • Is this a sensible assumption?
• The model assumes the same (linear) effect of Days without sleep for all participants (slope)
  • Is this a sensible assumption?

Another point for good measure

• this dataset was balanced (equal number of L1 in L2)
• if it was unbalanced, groups with more observations would weigh more heavily because the model assumes all observations are independent (but they are not)

• ppt B contributes more because they have more trials

Session outcomes

• Describe a linear model conceptually and using mathematical notation
  
• Run and interpret a linear model in R
  
• Make predictions using a linear model in R
  
• Critique the “one size fits all” model in the presence of multiple participants

Activity

• Using the (mlmRev::Oxboys) data
• How many subjects are there?
• How many Occasions occurred?
• is this a balanced dataset?
• Plot the data so that we can see how height varies as a function of occasion for each subject
• Run a ‘complete pooling’ model testing whether height varies as a linear function of occasion, first make Occasion a numeric variable
• plot a model of the model derived slope over the spread of raw data
• does the model meet the assumptions of normality and spread of residuals
• Plot how the trends would look if superimposed on raw data scatter for each ppt, using the coef from m1