From linear model to linear mixed model

• Linear model:
  • Subjects are independent.
  • Cannot handle missing data.

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

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

voice pitch data

d = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
str(d)

## 'data.frame':    84 obs. of  5 variables:
##  $ subject  : Factor w/ 6 levels "F1","F2","F3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gender   : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ scenario : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ attitude : Factor w/ 2 levels "inf","pol": 2 1 2 1 2 1 2 1 2 1 ...
##  $ 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))))

## [1] 39

visualize the data

library(ggplot2)
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.

Independence:

• 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
• Under the same scenario \(j\), people should also react similarly, not totally independent.
  • \(frequency_{1jk}, frequency_{2jk}, \ldots, frequency_{Ijk}\) should be correlated

data independence

library(tidyr)
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)) + 
  geom_smooth(method="lm")

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

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

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).

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

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

scenario independence

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

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

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

scenario independence

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

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

## 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
  • E.g., for repeated measures, the “group” 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 “group” effect is random if we can think of the levels we observe in that group 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

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
tryCatch({lmer(frequency ~ attitude, data = d)}, error=function(e) print(e))

## <simpleError: No random effects terms specified in formula>

• 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 + \beta_j + \alpha_i + \varepsilon_{ijk}\]

• \(i\): subject index
• \(j\): attitude index
• \(k\): scenario index
• \(\mu\): intercept
• \(\beta_j\): attitude coefficient of \(j^{th}\) group
• \(\alpha_i \sim N(0, \sigma^2_\alpha)\): random intercept for the \(i^{th}\) subject
• \(\varepsilon_{ijk} \sim N(0, \sigma^2)\):

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 intercept: \[\mu + \beta_j + \alpha_i\]
• Subject attitudepol (slope): \[\beta_j\]

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\)

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

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.

Getting 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.

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)
##            Df    AIC    BIC  logLik deviance  Chisq Chi 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 another covariate gender and its interaction with attitude

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: 778.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2334 -0.5105 -0.1551  0.4738  3.0875 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 605.4    24.61   
##  Residual             844.1    29.05   
## Number of obs: 83, groups:  subject, 6
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)          260.686     15.557    4.754  16.757 2.07e-05 ***
## attitudepol          -27.400      8.966   75.009  -3.056  0.00311 ** 
## genderM             -116.195     22.000    4.754  -5.282  0.00376 ** 
## attitudepol:genderM   16.190     12.763   75.018   1.269  0.20852    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp gendrM
## attitudepol -0.288              
## genderM     -0.707  0.204       
## atttdpl:gnM  0.202 -0.703 -0.286

Investigate scenario

ggplot(aes(1, frequency, col=scenario0), data=d) + 
  geom_boxplot() + 
  geom_jitter() + 
  facet_grid(attitude ~ scenario0)  + 
  theme_bw()

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

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

A larger mixed model

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

## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: frequency ~ attitude + gender + (1 | subject) + (1 | scenario)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##    807.1    821.6   -397.6    795.1       77 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2958 -0.6456 -0.0776  0.5448  3.5121 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  scenario (Intercept) 205.2    14.33   
##  subject  (Intercept) 417.0    20.42   
##  Residual             637.4    25.25   
## Number of obs: 83, groups:  scenario, 7; subject, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  256.847     13.827    8.500  18.576 3.53e-08 ***
## attitudepol  -19.722      5.547   70.920  -3.555 0.000677 ***
## genderM     -108.517     17.572    5.929  -6.176 0.000867 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp
## attitudepol -0.201       
## genderM     -0.635  0.004

A larger mixed model

coef(rs_intergen_subjscene_ml)

## $scenario
##   (Intercept) attitudepol   genderM
## 1    243.4859   -19.72207 -108.5173
## 2    263.3592   -19.72207 -108.5173
## 3    268.1322   -19.72207 -108.5173
## 4    277.2546   -19.72207 -108.5173
## 5    254.9319   -19.72207 -108.5173
## 6    244.8015   -19.72207 -108.5173
## 7    245.9618   -19.72207 -108.5173
## 
## $subject
##    (Intercept) attitudepol   genderM
## F1    243.3684   -19.72207 -108.5173
## F2    266.9443   -19.72207 -108.5173
## F3    260.2276   -19.72207 -108.5173
## M3    284.3536   -19.72207 -108.5173
## M4    262.0575   -19.72207 -108.5173
## M7    224.1292   -19.72207 -108.5173
## 
## attr(,"class")
## [1] "coef.mer"

Random Slopes

ggplot(aes(attitude, frequency, color = subject), data=d) + 
  geom_point() + 
  geom_smooth(method="lm", aes(group=1)) + 
  facet_wrap(~subject ) + 
  theme_bw()

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

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

Random Slopes in R

politeness.model.rs = lmer(frequency ~  attitude + gender + (1 + attitude | subject) + 
    (1 | scenario), REML = FALSE, data = d)
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) + (1 |  
##     scenario)
##    Data: d
## 
##      AIC      BIC   logLik deviance df.resid 
##    811.1    830.4   -397.5    795.1       75 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2799 -0.6554 -0.0803  0.5368  3.5015 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  scenario (Intercept) 205.099  14.321       
##  subject  (Intercept) 395.367  19.884       
##           attitudepol   1.257   1.121   1.00
##  Residual             637.064  25.240       
## Number of obs: 83, groups:  scenario, 7; subject, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  257.804     13.682    7.977  18.843 6.73e-08 ***
## attitudepol  -19.715      5.565   58.458  -3.543 0.000785 ***
## genderM     -110.432     17.530    5.929  -6.300 0.000781 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp
## attitudepol -0.152       
## genderM     -0.641  0.003

Random Slopes in R

coef(politeness.model.rs)

## $scenario
##   (Intercept) attitudepol   genderM
## 1    244.4408   -19.71467 -110.4322
## 2    264.3135   -19.71467 -110.4322
## 3    269.0864   -19.71467 -110.4322
## 4    278.2085   -19.71467 -110.4322
## 5    255.8865   -19.71467 -110.4322
## 6    245.7769   -19.71467 -110.4322
## 7    246.9166   -19.71467 -110.4322
## 
## $subject
##    (Intercept) attitudepol   genderM
## F1    243.7405   -20.50762 -110.4322
## F2    266.7312   -19.21134 -110.4322
## F3    260.1365   -19.58316 -110.4322
## M3    285.4191   -18.15766 -110.4322
## M4    263.9136   -19.37020 -110.4322
## M7    226.8841   -21.45803 -110.4322
## 
## attr(,"class")
## [1] "coef.mer"

Is random slope necessary

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

rs_gen_subjscene_reml = lmer(frequency ~ attitude + gender + (1 | subject) + (1 | 
    scenario), 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) + (1 | scenario)
## rs_gen_subjscene_con_reml: frequency ~ attitude + gender + (1 + attitude | subject) + (1 | 
## rs_gen_subjscene_con_reml:     scenario)
##                           Df    AIC    BIC  logLik deviance  Chisq Chi Df
## rs_gen_subjscene_reml      6 807.10 821.61 -397.55   795.10              
## rs_gen_subjscene_con_reml  8 811.07 830.42 -397.54   795.07 0.0275      2
##                           Pr(>Chisq)
## rs_gen_subjscene_reml               
## rs_gen_subjscene_con_reml     0.9863

• 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

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.

Another example sleep study

sleepStudy <- read.table("https://www.stat.tamu.edu/~sheather/book/docs/datasets/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.

Simple longitudinal data

• Repeated measures data consist of measurements of a response (and, perhaps, some covariates) on several experimental (or observational) units.
• 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.

Data exploratory

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

Fit linear mixed model in R

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

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Reaction ~ Days + (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.09   24.740       
##           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.6662579
## 309    211.0065  1.8475828
## 310    212.4449  5.0184061
## 330    275.0956  5.6529547
## 331    273.6653  7.3973914
## 332    260.4446 10.1951153
## 333    268.2455 10.2436615
## 334    244.1725 11.5418620
## 335    251.0714 -0.2848731
## 337    286.2955 19.0955699
## 349    226.1950 11.6407002
## 350    238.3351 17.0814910
## 351    255.9829  7.4520288
## 352    272.2687 14.0032993
## 369    254.6806 11.3395026
## 370    225.7922 15.2897506
## 371    252.2121  9.4791309
## 372    263.7196 11.7513157
## 
## 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