From linear model to linear mixed model

• Linear model:
  • Subjects are independent.
  • Cannot handle missing date (missing in Y).

\[Y = X^\top\beta + \varepsilon\]

• Linear mixed model:
  • Account for correlated structure between subjects.
  • Can accommodate missing data (missing in Y).

voice pitch data

d = read.csv("https://caleb-huo.github.io/teaching/data/politeness/politeness_data.csv")
str(d)

## 'data.frame':    84 obs. of  5 variables:
##  $ subject  : chr  "F1" "F1" "F1" "F1" ...
##  $ gender   : chr  "F" "F" "F" "F" ...
##  $ scenario : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ attitude : chr  "pol" "inf" "pol" "inf" ...
##  $ frequency: num  213 204 285 260 204 ...

d$attitude <- as.factor(d$attitude)
d$scenario0 <- as.factor(d$scenario)

• subject: (6 subjects)
• gender: (F/M)
• scenario: different scenario (7 categories)
• attitude:
  • inf: informal situations (e.g., explaining to a friend why you’re late)
  • pol: politeness (e.g., giving an excuse for being late to a professor)
• frequency: outcome variable, frequency of voice pitch

voice pitch data

head(d)

##   subject gender scenario attitude frequency scenario0
## 1      F1      F        1      pol     213.3         1
## 2      F1      F        1      inf     204.5         1
## 3      F1      F        2      pol     285.1         2
## 4      F1      F        2      inf     259.7         2
## 5      F1      F        3      pol     203.9         3
## 6      F1      F        3      inf     286.9         3

with(d, table(subject, gender))

##        gender
## subject  F  M
##      F1 14  0
##      F2 14  0
##      F3 14  0
##      M3  0 14
##      M4  0 14
##      M7  0 14

which(apply(d,1,function(x) any(is.na(x)))) ## one subject contails missing data (missing outcome variable)

## [1] 39

visualize the data

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

ggplot(aes(attitude, frequency, color = subject), data=d) + facet_grid(~subject) + geom_boxplot() + geom_jitter() + theme_bw()

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

## Warning: Removed 1 rows containing missing values (geom_point).

Fit linear model

• Outcome variable: pitch (frequency)
• primary predictor: attitude (politeness)
• \(\varepsilon\): deviations from our predictions due to random factors that we cannot control experimentally

\[frequency \sim \mu + politeness + \varepsilon\]

Linear model in R

lmFit <- lm(frequency ~ attitude, data = d)
summary(lmFit)

## 
## Call:
## lm(formula = frequency ~ attitude, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -103.488  -62.122    9.044   51.178  105.044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   202.59      10.08  20.107   <2e-16 ***
## attitudepol   -18.23      14.34  -1.272    0.207    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.3 on 81 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.01958,    Adjusted R-squared:  0.007475 
## F-statistic: 1.618 on 1 and 81 DF,  p-value: 0.2071

Note:

• \(1\) observation deleted due to missingness
• attitude is not a significant predictor

Problem with linear model

\[frequency_{ijk} \sim \mu + politeness_{k} + \varepsilon_{ijk},\] where \(i\) is a subject index, \(j\) is a scenario index, and \(k\) is the politeness index.

Dependent structure:

• Multiple responses from the same subject \(i\) cannot be regarded as independent from each other.
  • \(frequency_{i1k}, frequency_{i2k}, \ldots, frequency_{iJk}\) should be correlated

data independence

subj_spread <- d %>% spread(attitude, frequency)
subj_spread$scenario0 <- as.factor(subj_spread$scenario)
head(subj_spread)

##   subject gender scenario scenario0   inf   pol
## 1      F1      F        1         1 204.5 213.3
## 2      F1      F        2         2 259.7 285.1
## 3      F1      F        3         3 286.9 203.9
## 4      F1      F        4         4 276.8 250.8
## 5      F1      F        5         5 252.4 231.9
## 6      F1      F        6         6 230.7 181.2

Subject dependence

ggplot(data=subj_spread, aes(pol, inf)) + theme_bw() +
  geom_point(aes(col=subject)) + 
  stat_ellipse(aes(col=subject)) +
  geom_smooth(method="lm")

## Warning: Removed 1 rows containing non-finite values (stat_ellipse).

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

Random effects:

• Fixed effects:
  • Categorical variables with well-defined categories (e.g., gender, race)
• Random effects:
  • We only sample a subset of the entire population
    • The 6 participants in the data are a subset of the entire population
    • Each participant has repeated measurement.
    • The “subject” effect \(\alpha_i\) for subject \(i\) is best thought of as random: \(\alpha_i \sim N(0, \sigma_0^2)\)

When to use random effects?

• A “subject” effect is random if we can think of the levels we observe for that subject to be samples from a larger population.
• Example: we are investigating the GPA of UF undergraduate students by surveying a small group of students.
  • Gender and race are considered fixed effect since their categories are well-defined.
  • Department is considered as random effect, since our survey may not covering all departments in UF.

Mixed model

• Linear model is “fixed-effects only” model
  • one or more fixed effects
  • a general error term \(\varepsilon\)
  • error term is considered no structure.
• Mixed model is a mixture of fixed effect and random effect

Benefit of random effect models

• Allows us to resolve this non-independence data structure for each subject
• Allows missing data (missing outcome variable)

Fit linear mixed model in R

• In R, lmer() function in lme4 package will do the job

alm <- lm(frequency ~ attitude, data = d)
suppressMessages(library(lme4))
suppressMessages(library(lmerTest)) ## this package helps provide a p-value

lmer(frequency ~ attitude, data = d)

• We have to specify random effect when using lmer funciton.

Fit linear mixed model in R

rs_subj_reml = lmer(frequency ~ attitude + (1 | subject), data = d)
summary(rs_subj_reml)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: frequency ~ attitude + (1 | subject)
##    Data: d
## 
## REML criterion at convergence: 804.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2953 -0.6018 -0.2005  0.4774  3.1772 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 3982     63.10   
##  Residual              851     29.17   
## Number of obs: 83, groups:  subject, 6
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  202.588     26.151   5.150   7.747 0.000501 ***
## attitudepol  -19.376      6.407  76.003  -3.024 0.003399 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## attitudepol -0.121

How to write the linear mixed model mathmatically

\[Y_{ijk} = \mu + \gamma \mathbb{I}(k=pol) + \alpha_i + \varepsilon_{ijk}\]

• \(i\): subject index
• \(j\): scenario index
• \(k\): attitude index
• \(\mu\): intercept
• \(\gamma\): attitude coefficient of the polite group
• \(\alpha_i \sim N(0, \sigma^2_\alpha)\): random effect for the \(i^{th}\) subject
• \(\varepsilon_{ijk} \sim N(0, \sigma^2)\):
• Total error term: \(\alpha_i + \varepsilon_{ijk}\)
  • so the error for \(i1k\) and \(i2k\) are dependent
• If there is no \(\alpha_i\), the above model will reduce to linear regression model.

Further explore rs_subj_reml

coef(alm)

## (Intercept) attitudepol 
##    202.5881    -18.2320

coef(rs_subj_reml)

## $subject
##    (Intercept) attitudepol
## F1    241.1352   -19.37584
## F2    266.8920   -19.37584
## F3    259.5540   -19.37584
## M3    179.0262   -19.37584
## M4    155.6906   -19.37584
## M7    113.2306   -19.37584
## 
## attr(,"class")
## [1] "coef.mer"

AIC(rs_subj_reml)

## [1] 812.7023

Further explore rs_subj_reml

• Subject total intercept: \[\mu + \alpha_i\]
  • vary among different subjects
• Attitude effect \[\beta_j\]
  • same among different subjects

Random Intercepts

• linear model \[frequency \sim 1 + politeness + \varepsilon\]

• linear model with random effect (random intercept) \[frequency \sim 1 + politeness +(1|subject) + \varepsilon\] \((1|subject)\) is the R syntax for a random intercept.

• intercepts are different for each subject

• “\(1\)” stands for the intercept in R

• this formula indicates multiple responses per subject

• resolves the non-independence multiple responses by the same subject

• the formula still contains a general error term \(\varepsilon\)

Visualization for the random intercept

## Warning: Removed 1 rows containing missing values (geom_point).

Fit linear mixed model in R (REML = FALSE)

rs_subj_ml = lmer(frequency ~ attitude + (1 | subject), REML = FALSE, data = d)
summary(rs_subj_ml)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: frequency ~ attitude + (1 | subject)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##    826.5    836.2   -409.3    818.5       79 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3067 -0.6051 -0.2042  0.4793  3.1931 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 3308.7   57.52   
##  Residual              839.9   28.98   
## Number of obs: 83, groups:  subject, 6
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  202.588     23.905   6.214   8.475 0.000122 ***
## attitudepol  -19.372      6.365  77.003  -3.043 0.003199 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## attitudepol -0.131

REML option

REML represents restricted maximum likelihood

• if REML=FALSE, maximum likelihood estimation will be used (biased estimation)
• if REML=TRUE, its estimation does not base on a maximum likelihood, but instead uses a likelihood function calculated from a transformed set of data (unbiased estimation)

Compare to linear model

summary(lm(frequency ~ attitude, data = d))

## 
## Call:
## lm(formula = frequency ~ attitude, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -103.488  -62.122    9.044   51.178  105.044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   202.59      10.08  20.107   <2e-16 ***
## attitudepol   -18.23      14.34  -1.272    0.207    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.3 on 81 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.01958,    Adjusted R-squared:  0.007475 
## F-statistic: 1.618 on 1 and 81 DF,  p-value: 0.2071

Without adjusting for random effect, attitude won’t be a significant predictor.

Alternative ways to get p-values (likelihood ratio test)

• likelihood ratio test: comparing full model and reduced model
• get the p-value for reduced variable
• need to use REML = FALSE, otherwise, anova will re-fit the reml model such that REML = FALSE.

rs_subj_ml = lmer(frequency ~ attitude + (1 | subject), REML = FALSE, data = d)
rs_null_ml = lmer(frequency ~ 1 + (1 | subject), REML = FALSE, data = d)

anova(rs_null_ml, rs_subj_ml)

## Data: d
## Models:
## rs_null_ml: frequency ~ 1 + (1 | subject)
## rs_subj_ml: frequency ~ attitude + (1 | subject)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## rs_null_ml    3 833.25 840.51 -413.62   827.25                        
## rs_subj_ml    4 826.51 836.18 -409.25   818.51 8.7444  1   0.003106 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Add covariates

summary(lmer(frequency ~ attitude + gender +  (1 | subject), data = d))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: frequency ~ attitude + gender + (1 | subject)
##    Data: d
## 
## REML criterion at convergence: 786.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3619 -0.5305 -0.1724  0.4647  3.2260 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 603.9    24.57   
##  Residual             850.9    29.17   
## Number of obs: 83, groups:  subject, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  256.691     15.226    4.378  16.859 3.78e-05 ***
## attitudepol  -19.410      6.407   76.018  -3.030  0.00334 ** 
## genderM     -108.205     21.063    4.009  -5.137  0.00677 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp
## attitudepol -0.210       
## genderM     -0.692  0.004

Randon intercept and random slope

• Random intercept only:

\[Y_{ijk} = (\mu + \alpha_i) + \gamma \mathbb{I}(k=pol) + \varepsilon_{ijk}\]

• Random intercept and random slope:

\[Y_{ijk} = (\mu + \alpha_i) + (\gamma + \eta_i) \mathbb{I}(k=pol) + \varepsilon_{ijk}\]

• \((\alpha_i, \eta_i)^\top \sim N(0, \Sigma)\): random intercept and random slope for the \(i^{th}\) subject
• \(\varepsilon_{ijk} \sim N(0, \sigma^2)\):

Visualization for the random slope

## boundary (singular) fit: see ?isSingular

## Warning: Removed 1 rows containing missing values (geom_point).

Random Slopes in R

politeness.model.rs = lmer(frequency ~  attitude + gender + (1 + attitude | subject), REML = FALSE, data = d)

## boundary (singular) fit: see ?isSingular

summary(politeness.model.rs)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: frequency ~ attitude + gender + (1 + attitude | subject)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##    820.3    837.2   -403.2    806.3       76 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3848 -0.5270 -0.2221  0.5179  3.2462 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  subject  (Intercept) 359.95   18.972       
##           attitudepol   1.57    1.253   1.00
##  Residual             839.41   28.973       
## Number of obs: 83, groups:  subject, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  257.503     12.396    6.274  20.773 5.06e-07 ***
## attitudepol  -19.408      6.384   63.951  -3.040 0.003423 ** 
## genderM     -109.829     17.167    6.021  -6.398 0.000677 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp
## attitudepol -0.205       
## genderM     -0.692  0.004
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Random Slopes in R

coef(politeness.model.rs)

## $subject
##    (Intercept) attitudepol   genderM
## F1    244.1807   -20.28814 -109.8295
## F2    266.1358   -18.83824 -109.8295
## F3    259.8311   -19.25459 -109.8295
## M3    283.4821   -17.69271 -109.8295
## M4    263.7729   -18.99429 -109.8295
## M7    227.6144   -21.38216 -109.8295
## 
## attr(,"class")
## [1] "coef.mer"

Is random slope necessary

rs_gen_subjscene_con_reml = lmer(frequency ~ attitude + gender + (1 + attitude | 
    subject), REML = TRUE, data = d)

## boundary (singular) fit: see ?isSingular

rs_gen_subjscene_reml = lmer(frequency ~ attitude + gender + (1 | subject), REML = TRUE, data = d)

anova(rs_gen_subjscene_reml, rs_gen_subjscene_con_reml)

## refitting model(s) with ML (instead of REML)

## Data: d
## Models:
## rs_gen_subjscene_reml: frequency ~ attitude + gender + (1 | subject)
## rs_gen_subjscene_con_reml: frequency ~ attitude + gender + (1 + attitude | subject)
##                           npar    AIC    BIC  logLik deviance  Chisq Df
## rs_gen_subjscene_reml        5 816.34 828.43 -403.17   806.34          
## rs_gen_subjscene_con_reml    7 820.31 837.24 -403.15   806.31 0.0282  2
##                           Pr(>Chisq)
## rs_gen_subjscene_reml               
## rs_gen_subjscene_con_reml      0.986

• it appears that we don’t need to include random slope for condition in the model, because the likelihood ratio test is not significant
• however, others would argue that we should keep our models maximal

Beyond linear mixed model (optional)

If the outcome variable is

• Continuous: use a linear regression model with mixed effects
• Binary: use a logistic regression model with mixed effects
• Count: use a Poisson/negative binomial regression model with mixed effects

Function lmer is used to fit linear mixed models, function glmer is used to fit generalized (non-Gaussian) linear mixed models.

Longitudinal data

• Longitudinal data are repeated measures data in which the observations are taken over time.
• We wish to characterize the response over time within subjects and the variation in the time trends between subjects.

Longitudinal data – sleep study

sleepStudy <- read.table("https://caleb-huo.github.io/teaching/data/sleep/sleepstudy.txt", header=T, row.names = 1)
sleepStudy$Subject <- as.factor(sleepStudy$Subject)

head(sleepStudy, n=10)

##    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

str(sleepStudy)

## 'data.frame':    180 obs. of  3 variables:
##  $ Reaction: num  250 259 251 321 357 ...
##  $ Days    : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...

Sleep deprivation data

• This laboratory experiment measured the effect of sleep deprivation on cognitive performance.
• There were 18 subjects, chosen from the population of interest (long-distance truck drivers), in the 10 day trial. These subjects were restricted to 3 hours sleep per night during the trial.
• On each day of the trial each subject’s reaction time was measured. The reaction time shown here is the average of several measurements.
• These data are balanced in that each subject is measured the same number of times and on the same occasions.

Data exploratory

ggplot(aes(Days, Reaction, col=Subject), data = sleepStudy) + geom_point() + geom_smooth(method="lm") + facet_wrap(~Subject)

## `geom_smooth()` using formula 'y ~ x'

Fit linear mixed model in R

• Examine the relationship of Reaction time with respect to Days (slope)
• Each subject has repeated measurement of Reaction time
• Consider subject specific random intercept only

fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
summary(fm1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 251.4051     9.7467  22.8102   25.79   <2e-16 ***
## Days         10.4673     0.8042 161.0000   13.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Fit linear mixed model in R

• Examine the relationship of Reaction time with respect to Days (slope)
• Each subject has repeated measurement of Reaction time
• Consider both subject specific random intercept and subject specific random slope

fm1 <- lmer(Reaction ~ Days + (1 + Days|Subject), sleepstudy)
#fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) ## same as this one
summary(fm1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (1 + Days | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1743.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9536 -0.4634  0.0231  0.4634  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 612.10   24.741       
##           Days         35.07    5.922   0.07
##  Residual             654.94   25.592       
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  251.405      6.825  17.000  36.838  < 2e-16 ***
## Days          10.467      1.546  17.000   6.771 3.26e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.138

coef

coef(fm1)

## $Subject
##     (Intercept)       Days
## 308    253.6637 19.6662617
## 309    211.0064  1.8476053
## 310    212.4447  5.0184295
## 330    275.0957  5.6529356
## 331    273.6654  7.3973743
## 332    260.4447 10.1951090
## 333    268.2456 10.2436499
## 334    244.1725 11.5418676
## 335    251.0714 -0.2848792
## 337    286.2956 19.0955511
## 349    226.1949 11.6407181
## 350    238.3351 17.0815038
## 351    255.9830  7.4520239
## 352    272.2688 14.0032871
## 369    254.6806 11.3395008
## 370    225.7921 15.2897709
## 371    252.2122  9.4791297
## 372    263.7197 11.7513080
## 
## attr(,"class")
## [1] "coef.mer"

References:

Credits give to:

• https://web.stanford.edu/class/psych252/section_2013/Mixed_models_tutorial.html
• Winter, B. (2013). Linear models and linear mixed effects models in R with linguistic applications. arXiv:1308.5499. https://arxiv.org/pdf/1308.5499.pdf
• http://lme4.r-forge.r-project.org/slides/2011-03-16-Amsterdam/2Longitudinal.pdf
• http://statweb.stanford.edu/~jtaylo/courses/stats203/notes/fixed+random.pdf