2.1.1 Simple linear regression: Example

Suppose a researcher measured the mean length of utterance (MLU) in a group of young children (2-4 years) in a neighbourhood, and found that it depends on age (\(r=0.85\)). In other words, you can predict a particular child’s MLU from their age reasonably well, or at least better than if you don’t know it.

Now a new family moves into the neighbourhood, how would you go about making a guess for their child’s MLU? Without any information about the child’s age, your best guess can only bee the average MLU. But if you know their age, you can predict their MLU based on age (i.e. MLU as a function of age).

In a simple linear regression, the model function can be written as

\[ y = \alpha + {\beta}x \]

where \(\alpha\) is the intercept (where the line croses the y-axis), and \(\beta\) is the gradient or slope (how steep the line is). Data can be described by

\[ y_i = \alpha + {\beta}x_i + \epsilon_i \]

where \(\epsilon\) is the error or residual (the difference between predicted and observed value).

2.1.2 Assumptions of simple linear regression

• Linearity: the relationship between the prediction (\(x\)) and the outcome (\(y\)) is assumed to be linear.
• Normality of residuals: the residuals (errors) are assumed to be normally distributed.
• Homoscedasticity (Homogeneity of residual variance): The residuals are assumed to have a constant variance.
• Independence of residuals error terms: There is not a relationship between the residuals and the \(y\) variable.

2.1.3 Effect size of a linear regression: R-squared

The \(R^2\) is a statistical measure in a regression model that determines the proportion of variance in the dependent variable that can be explained by the independent variable (aka variance accounted for). In other words, the \(R^2\) shows how well the data fit the regression model (the goodness of fit).

In a simple linear regression, the \(R^2\) is simply the square of the Pearson correlation coefficient. The \(R^2\) can be interpreted in similar ways as the correlation coefficient: a \(R^2\) closer to 1 indicates that a larger proportion of the variability in the outcome canbe explained by the regression model, a number closer to 0 indicates less variance accounted for.

2.1.4 Do it in R: Simple linear regression

1. Prepare data33 You can download the data here.: data should be in short form: the two measures for the same observation should be in the same row.

# import data
simple_regression_exp <- read.csv('mlu.csv', header = TRUE)
head(simple_regression_exp)

##   participant      age       MLU
## 1           1 3.245548 11.704914
## 2           2 3.085281 12.035840
## 3           3 2.624829  7.582783
## 4           4 3.426735 11.920344
## 5           5 2.923928 11.801822
## 6           6 2.438216  7.292018

# visualise the data using a scatter plot
library(ggplot2)
ggplot(simple_regression_exp, aes(x=age, y=MLU)) + geom_point() + 
  geom_smooth(method=lm, se=FALSE, color = 'darkgrey') + 
  theme_light()

2. Compute the linear model

simple_regression_model <- lm(MLU ~ age, data = simple_regression_exp)
summary(simple_regression_model)

## 
## Call:
## lm(formula = MLU ~ age, data = simple_regression_exp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9116 -0.9599  0.2080  0.8839  3.0391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.6878     1.7405  -2.693   0.0118 *  
## age           4.9756     0.5784   8.603 2.39e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.502 on 28 degrees of freedom
## Multiple R-squared:  0.7255, Adjusted R-squared:  0.7157 
## F-statistic:    74 on 1 and 28 DF,  p-value: 2.39e-09

3. Interpret the summary output

The summary output shows 6 components, including:

• call: the function used to compute the regression model.
• residuals: a quick view of the distribution of the residuals, which by definition have a mean zero. Therefore, the median should not be far from zero, and the menimum and maximum should be roughly equal in absolute value.
• coefficients: the regression model’s \(\beta\) coefficients and their statistical significance.
• residual standard error (RSE), \(R^2\), and the F-statistics: metrics that are used to check how well the model fits to the data.

How to interpret the linear regression summary table:

The model function: prediction = intercept estimate + predictor estimate * predictor value. In this case, the model function is \(y = -4.6878 + 4.9756*\text{age}\).

Significance of the predictor: read (the t-value and ) the p-value in the coefficients table. The null hypothesis in a regression is no relationship between \(x\) and \(y\). The alternative hypothesis is that there is a relationship between \(x\) and \(y\), or that the value of \(x\) can predict the value of \(y\). In this example, we find a significant predictor age (p<0.001).

In this case, we also have a significant intercept (p = 0.01), which basically means that when age = 0 (when the line crosses \(x=0\)), the MLU is significantly different from 0. This means very little to us in this case (what do you mean a new-born baby is predicted to have a non-zero utterance length?), but a significant intercepts may be interpretable in other cases.

Overall significance: read the F-statistics at the bottom of the output. The F-statistics asseess whether at least one predictor variable has a non-zero coefficient. In a simple linear regression, this is not really interesting because we just have one predictor in total. In fact, in this case, the F test is identical to the square of the t test $74 ^2 $.

The effect size: read the \(R^2\) (for simple regression) or the adjusted \(R^2\) (for multiple regression).

4. Check the assumptions (using diagnostic plots)

library(performance)
check_model(simple_regression_model)

5. Report the regression

Simple linear regression was used to predict children’s mean length of utterance (MLU) from their age. Age explained a significant amount of the variance in MLU, F(1,28)=74, p<0.001, R²=0.73. The regression coefficient (β=4.98) indicated that an increase in one year of age corresponded, on average, to an increase in MLU of 4.98 words.

2.2.1 Multiple linear regresion example

Let’s say the researcher in the last example also collected data on those children’s maternal education, measured in years of education (YoE); and they want to know whether MLU can be predicted as a function of age and maternal education.

2.2.2 Do it in R: Multiple linear regression

1. Prepare data44 You can download the data here.
2. Run the model (and check assumptions)

multiple_regression_exp <- read.csv('mlu-multi.csv', header = TRUE)
head(multiple_regression_exp)

##   participant      age       MLU  YoE
## 1           1 3.245548 11.704914 13.1
## 2           2 3.085281 12.035840  7.5
## 3           3 2.624829  7.582783 12.3
## 4           4 3.426735 11.920344  9.5
## 5           5 2.923928 11.801822 10.3
## 6           6 2.438216  7.292018 13.8

# compute the linear regression model
## Note that I put down age + YoE because in this example we are not interested in the interaction between age and YoE.
## If you want the interaction to be included in the model, put down age * YoE instead.
multiple_regression_model <- lm(MLU ~ age + YoE, data = multiple_regression_exp)
summary(multiple_regression_model)

## 
## Call:
## lm(formula = MLU ~ age + YoE, data = multiple_regression_exp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9530 -0.9383  0.2373  0.8541  3.0798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.87977    1.89934  -2.569    0.016 *  
## age          4.95104    0.59472   8.325  6.2e-09 ***
## YoE          0.02455    0.08812   0.279    0.783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.527 on 27 degrees of freedom
## Multiple R-squared:  0.7263, Adjusted R-squared:  0.706 
## F-statistic: 35.82 on 2 and 27 DF,  p-value: 2.532e-08

3. Model selection

From the summary output we can see that YoE is not a significant predictor in the model.

A good model not only needs to fit the data well, it also needs to be parsimonious (i.e. not more complex than necessary). If you are choosing beteween a simple model with one predictor and a complex model with, say, 10 predictors, the complex model needs to provide a much better fit to the data in order to be a better model. If it can’t, then the simpler model should be preferred.

In this case, we can compare whether the model with just one predictor age (the one we got from the simple regression) or the model with both age and YoE as predictors is better. For this we can use the anova()55 Note that the function we used to actually conduct an ANOVA is called aov(), and the anova() function is used to compare regression models. function.

# compare difference models using anova()
# model one (one predictor: age)
model_1 <- lm(MLU ~ age, data = multiple_regression_exp)
# model two (two predictors: age and YoE)
model_2 <- lm(MLU ~ age + YoE, data = multiple_regression_exp)

# compare model one with model two
anova(model_1, model_2)

## Analysis of Variance Table
## 
## Model 1: MLU ~ age
## Model 2: MLU ~ age + YoE
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     28 63.168                           
## 2     27 62.987  1   0.18105 0.0776 0.7827

From the anova() output, we can see that adding YoE as a predictor did not significantly improve the model fit. Thus we should reject model two and stick to the simpler model one.

2.3.1 Random intercepts, random slopes, and the maximal random effect structure

In plain words, a model with random effects allow the model predictions to systematically vary based on a random factor. Now there are different ways that those systematic variation can happen, and they are random intercepts and random slopes.

A model with random intercepts allow each level of the random factor to have a different intercept, while keeping the slopes constant. In contrast, random slopes allow each level of the random factor to have a different slope, while all predictions share the same intercept. A random factor can have random intercepts, random slopes for any number of the fix factors, or both. Here are some example data and what the model predictions look like with by-participant random intercepts, by-participant random slopes, or both.

Who decides what random factors should be in a model? You do. Just like it’s your responsibility to design the experiment (thus decide what the fixed factors are), it’s your responsibility to decide what random factors to include, and whether those should be random intercepts and/or slopes. What random factors to include usually isn’t too difficult: you’ll have a fairly well idea as long as you understand your experiment design. For example, in an experiment where 50 participants each read 100 English sentences and their reading times are measured, it’s a good idea to include participant and sentence (experimental item) as two random factors; as there are likely to be by-participant and by-sentence/item differences in the data that are systematic (e.g. some people read faster than others; some sentences takes longer to read than others.).

But how do I decide whether participant and item should have random intercepts or/and random slopes? You can always plot your data and decide. It sounds stressful, but luckily there are some guidelines we can follow. Barr et al. (2012) suggested that we should always begin with what’s called the maximal random effect structure: including both random intercepts and random slopes for all random factors. And unless the maximal model has problems88 Usually non-convergence or singular-fit problems. Don’t worry R will warn you if those happens.
And of course another possibility is that your computer isn’t powerful enough to handle the maximal random effect structure (and you don’t have access to a more powerful one), which surprisingly isn’t too uncommon in my experience. (in which case plot your data and see what terms you can take out), we should stick to it.

2.3.2 Coding categorical factors

If you have a categorical factor in a regression (fixed or random), you’ll need to code it (convert it to numbers) before the computer can make sense of it. Sometimes this is straightforward, for example if you have two groups of patients, one underwent treatment and the other didn’t, it’s easy to say that the group that did not undergo treatment is some sort of default/reference (and should receive the value 0), and the group that had treatment is where something actually happened (and should receive the value 1). What I just described is called treatment coding, or dummy coding.

But a lot of the time there is no reference-change relationship between the levels of your factor: male vs. female, north vs. south, for example. There is no default and all levels in those factors are equal, they’re just categorically different. In this case we usually use what’s called sum coding or deviation coding, which makes all levels of the factor equally far away from the reference (0). If you have two levels in a factor, they might be coded as -0.5, +0.5, or -1, +1.

There are a ton of other coding schemes but they’re outside the scope of this class.

2.3.3 Do it in R: Mixed-effects linear regression

Suppose the data from the last example is collected from two different neighbourhoods: one in the north and one in the south. Suppose that we are not interested in the effect of region (north vs. south) on MLU, it just happens that we had to collect data from two regions because there weren’t enough children in one region. We can include region (north vs. south) as a random factor.

1. Prepare data99 You can download the data here.

mixed_regression_exp <- read.csv('mlu-mix.csv', header=TRUE)
head(mixed_regression_exp)

##   participant      age       MLU region
## 1           1 3.245548 11.704914  south
## 2           2 3.085281 12.035840  south
## 3           3 2.624829  7.582783  north
## 4           4 3.426735 11.920344  south
## 5           5 2.923928 11.801822  south
## 6           6 2.438216  7.292018  north

2. Code categorical factors.

mixed_regression_exp$region <- as.factor(mixed_regression_exp$region)
contrasts(mixed_regression_exp$region) <- contr.sum(2)
contrasts(mixed_regression_exp$region)

##       [,1]
## north    1
## south   -1

3. Run the regression model (lme4 package).

The formula for a mixed-effects linear model with a maximal random effect structure is:

lmer(DV ~ IV1 * IV2 + (1 + IV1 * IV2 | RV1) + (1 + IV1 * IV2 | RV2), data = data)

where 1 indicate random intercepts (0 for no random intercepts), and the fixed factor terms in the randon effect structure indicate random slopes.

# the summary table of mixed effect models does not tell you the p value by default
# This is because p-values can be sometimes misleading in a mixed model
# But for this class that's out of our scope
# we can force the p-values to show using the lmerTest package
library(lme4)
library(lmerTest)
mixed_model <- lmer(MLU ~ age + (1 + age | region), data = mixed_regression_exp)
summary(mixed_model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MLU ~ age + (1 + age | region)
##    Data: mixed_regression_exp
## 
## REML criterion at convergence: 100
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17690 -0.49386 -0.00459  0.52190  1.92606 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  region   (Intercept) 5.850962 2.41888       
##           age         0.008988 0.09481  -1.00
##  Residual             1.568340 1.25233       
## Number of obs: 30, groups:  region, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   3.0603     3.1498  1.5322   0.972    0.459  
## age           2.3770     0.8872 18.5487   2.679    0.015 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## age -0.875

4. Check the assumptions (same as linear regression).

library(performance)
check_model(mixed_model)

5. Report the results: fixed and random effect structures, and effects.

Mixed-effect linear regression with one fixed factor age as well as by-region random intercepts and random slopes revealed a sigficiant main effect of age on the mean length of utterance in children (\(\beta\)=2.35, t=2.68, p=0.007).

2.4.1 GLMM example: Mixed-effects logistic regression

In this class I’ll just talk about one case of GLMMs that might be handy for us: the binomial distribution and the logistic regression. Your data follows a binomial distribution if it’s binary (e.g. yes vs. no, success vs. failure); this can be common in surveys and some experimental measures.

There is no time to delve into how logistic regression works1010 Basically the linear predictor is linked to the predictions (here probability of success) using a linking function, in this case the inverse logit function., let’s just take an example and learn how to do it in R.

Suppose an experimental phonologist want to know whether a phonological rule in a language A is gradually disappearing (if applying such a rule depends on the speaker’s age). 50 participants varying in age participated in an experiment where they were each asked to read out 20 words where the rule could apply. Data was recorded such that if the speaker applied the rule, the response was coded as 1, otherwise the response was coded as 0.

1. Prepare data1111 You can download the data here.

logit_exp <- read.csv('phonology.csv', header = TRUE)
head(logit_exp)

##   subject age word response
## 1       1  40    1        1
## 2       1  40    2        1
## 3       1  40    3        1
## 4       1  40    4        1
## 5       1  40    5        0
## 6       1  40    6        1

2. Run the logistic regression model

# here we are only including by-participant random intercepts
# the optimizer is there to optimize the way the computer finds the model parameters (coefficients), the details are outside of this class's scope
logistic_model <- glmer(response ~ age + (1 | subject), 
  data = logit_exp, 
  family = 'binomial',
  glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
summary(logistic_model)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: response ~ age + (1 | subject)
##    Data: logit_exp
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1008.9   1023.6   -501.5   1002.9      997 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7810  0.2092  0.4178  0.5561  1.2839 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  subject (Intercept) 0.01252  0.1119  
## Number of obs: 1000, groups:  subject, 50
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.41072    0.39618  -6.085 1.17e-09 ***
## age          0.09897    0.01109   8.925  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## age -0.979

3. Check the assumptions:

library(performance)
check_model(logistic_model)

4. Report the results

Mixed-effect logistic regression with one fixed factor age as well as by-participant random intercepts revealed a significant effect of age on whether a speaker applies the phonological rule in question (\(\beta\)=0.1, z=8.93, p<0.001).

In this case the researcher might also be interested in exactly how fast this phonological rule is disappearing, so reporting the following may be helpful (although you need to do some calculations based on the results):

It was found that with one year of age decreased, the odds of applying the phonological rule decreased by 10.4%1212 \(e^{\beta}\) (95% CI [.08, .13])1313 \(e^{{\beta}-1.96*\text{sd}}\), \(e^{{\beta}+1.96*\text{sd}}\).