# 第 56 章 多层线性模型

## 56.2 案例

create_data <- function() {
df <- tibble(
ids = 1:100,
department = rep(c("sociology", "biology", "english", "informatics", "statistics"), 20),
bases = rep(c(40000, 50000, 60000, 70000, 80000), 20) * runif(100, .9, 1.1),
experience = floor(runif(100, 0, 10)),
raises = rep(c(2000, 500, 500, 1700, 500), 20) * runif(100, .9, 1.1)
)

df <- df %>% mutate(
salary = bases + experience * raises
)
df
}
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.2
## Warning: package 'readr' was built under R version 4.2.2
## Warning: package 'purrr' was built under R version 4.2.2
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.2
## Warning: package 'lubridate' was built under R version 4.2.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lme4)
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.2.3
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
##     expand, pack, unpack
library(modelr)
## Warning: package 'modelr' was built under R version 4.2.3
library(broom)
## Warning: package 'broom' was built under R version 4.2.3
##
## Attaching package: 'broom'
##
## The following object is masked from 'package:modelr':
##
##     bootstrap
library(broom.mixed)

df <- create_data()
df
## # A tibble: 100 × 6
##      ids department   bases experience raises salary
##    <int> <chr>        <dbl>      <dbl>  <dbl>  <dbl>
##  1     1 sociology   43472.          6  1861. 54636.
##  2     2 biology     49262.          5   545. 51989.
##  3     3 english     56677.          4   482. 58605.
##  4     4 informatics 67269.          1  1599. 68868.
##  5     5 statistics  79641.          8   544. 83994.
##  6     6 sociology   41590.          7  2184. 56876.
##  7     7 biology     48098.          9   541. 52972.
##  8     8 english     57433.          7   536. 61181.
##  9     9 informatics 64145.          4  1794. 71323.
## 10    10 statistics  83165.          3   458. 84541.
## # ℹ 90 more rows

## 56.3 线性模型

$\hat{y} = \alpha + \beta_1x_1 + ... + \beta_nx_n$

# Model without respect to grouping
m1 <- lm(salary ~ experience, data = df)
m1
##
## Call:
## lm(formula = salary ~ experience, data = df)
##
## Coefficients:
## (Intercept)   experience
##       59509         1033
broom::tidy(m1)
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   59509.     2382.     25.0  2.75e-44
## 2 experience     1033.      467.      2.21 2.93e- 2
df %>% modelr::add_predictions(m1)
## # A tibble: 100 × 7
##      ids department   bases experience raises salary   pred
##    <int> <chr>        <dbl>      <dbl>  <dbl>  <dbl>  <dbl>
##  1     1 sociology   43472.          6  1861. 54636. 65710.
##  2     2 biology     49262.          5   545. 51989. 64676.
##  3     3 english     56677.          4   482. 58605. 63643.
##  4     4 informatics 67269.          1  1599. 68868. 60542.
##  5     5 statistics  79641.          8   544. 83994. 67777.
##  6     6 sociology   41590.          7  2184. 56876. 66743.
##  7     7 biology     48098.          9   541. 52972. 68810.
##  8     8 english     57433.          7   536. 61181. 66743.
##  9     9 informatics 64145.          4  1794. 71323. 63643.
## 10    10 statistics  83165.          3   458. 84541. 62609.
## # ℹ 90 more rows
# Model without respect to grouping
df %>%
ggplot(aes(x = experience, y = salary)) +
geom_point() +
geom_line(aes(x = experience, y = pred)) +
labs(x = "Experience", y = "Predicted Salary") +
ggtitle("linear model Salary Prediction") +
scale_colour_discrete("Department")

## 56.4 变化的截距

$\hat{y_i} = \alpha_{j[i]} + \beta x_i$

# Model with varying intercept
m2 <- lmer(salary ~ experience + (1 | department), data = df)
m2
## Linear mixed model fit by REML ['lmerMod']
## Formula: salary ~ experience + (1 | department)
##    Data: df
## REML criterion at convergence: 1911.831
## Random effects:
##  Groups     Name        Std.Dev.
##  department (Intercept) 14055
##  Residual                3499
## Number of obs: 100, groups:  department, 5
## Fixed Effects:
## (Intercept)   experience
##     59921.8        936.3
broom.mixed::tidy(m2, effects = "fixed")
## # A tibble: 2 × 5
##   effect term        estimate std.error statistic
##   <chr>  <chr>          <dbl>     <dbl>     <dbl>
## 1 fixed  (Intercept)   59922.     6320.      9.48
## 2 fixed  experience      936.      130.      7.20
broom.mixed::tidy(m2, effects = "ran_vals")
## # A tibble: 5 × 6
##   effect   group      level       term        estimate std.error
##   <chr>    <chr>      <chr>       <chr>          <dbl>     <dbl>
## 1 ran_vals department biology     (Intercept)  -10420.      781.
## 2 ran_vals department english     (Intercept)   -2442.      781.
## 3 ran_vals department informatics (Intercept)   11860.      781.
## 4 ran_vals department sociology   (Intercept)  -15808.      781.
## 5 ran_vals department statistics  (Intercept)   16810.      781.
df %>%
ggplot(aes(
x = experience, y = salary, group = department,
colour = department
)) +
geom_point() +
geom_line(aes(x = experience, y = pred)) +
labs(x = "Experience", y = "Predicted Salary") +
ggtitle("Varying Intercept Salary Prediction") +
scale_colour_discrete("Department")

## 56.5 变化的斜率

$\hat{y_i} = \alpha + \beta_{j[i]}x_i$

# Model with varying slope
m3 <- lmer(salary ~ experience + (0 + experience | department), data = df)
m3
## Linear mixed model fit by REML ['lmerMod']
## Formula: salary ~ experience + (0 + experience | department)
##    Data: df
## REML criterion at convergence: 2071.328
## Random effects:
##  Groups     Name       Std.Dev.
##  department experience 2376
##  Residual              8226
## Number of obs: 100, groups:  department, 5
## Fixed Effects:
## (Intercept)   experience
##       60121          852
broom.mixed::tidy(m3, effects = "fixed")
## # A tibble: 2 × 5
##   effect term        estimate std.error statistic
##   <chr>  <chr>          <dbl>     <dbl>     <dbl>
## 1 fixed  (Intercept)   60121.     1502.    40.0
## 2 fixed  experience      852.     1104.     0.772
broom.mixed::tidy(m3, effects = "ran_vals")
## # A tibble: 5 × 6
##   effect   group      level       term       estimate std.error
##   <chr>    <chr>      <chr>       <chr>         <dbl>     <dbl>
## 1 ran_vals department biology     experience   -1876.      402.
## 2 ran_vals department english     experience    -335.      289.
## 3 ran_vals department informatics experience    2261.      389.
## 4 ran_vals department sociology   experience   -2604.      393.
## 5 ran_vals department statistics  experience    2553.      352.
df %>%
ggplot(aes(
x = experience, y = salary, group = department,
colour = department
)) +
geom_point() +
geom_line(aes(x = experience, y = pred)) +
labs(x = "Experience", y = "Predicted Salary") +
ggtitle("Varying slope Salary Prediction") +
scale_colour_discrete("Department")

## 56.6 变化的斜率 + 变化的截距

$\hat{y_i} = \alpha_{j[i]} + \beta_{j[i]}x_i$ 具体来说，教师$$i$$，所在的学院$$j$$, 他的入职的起始收入表示为 ($$\alpha_{j[i]}$$)，年度增长率表示为($$\beta_{j[i]}$$).

# Model with varying slope and intercept
m4 <- lmer(salary ~ experience + (1 + experience | department), data = df)
m4
## Linear mixed model fit by REML ['lmerMod']
## Formula: salary ~ experience + (1 + experience | department)
##    Data: df
## REML criterion at convergence: 1906.88
## Random effects:
##  Groups     Name        Std.Dev. Corr
##  department (Intercept) 14398.1
##             experience    464.3  -0.23
##  Residual                3313.2
## Number of obs: 100, groups:  department, 5
## Fixed Effects:
## (Intercept)   experience
##     59758.9        983.8
broom.mixed::tidy(m4, effects = "fixed")
## # A tibble: 2 × 5
##   effect term        estimate std.error statistic
##   <chr>  <chr>          <dbl>     <dbl>     <dbl>
## 1 fixed  (Intercept)   59759.     6471.      9.24
## 2 fixed  experience      984.      242.      4.06
broom.mixed::tidy(m4, effects = "ran_vals")
## # A tibble: 10 × 6
##    effect   group      level       term        estimate std.error
##    <chr>    <chr>      <chr>       <chr>          <dbl>     <dbl>
##  1 ran_vals department biology     (Intercept)   -8913.     1054.
##  2 ran_vals department english     (Intercept)   -1369.     1687.
##  3 ran_vals department informatics (Intercept)   10660.     1209.
##  4 ran_vals department sociology   (Intercept)  -17902.     1257.
##  5 ran_vals department statistics  (Intercept)   17524.     1154.
##  6 ran_vals department biology     experience     -442.      220.
##  7 ran_vals department english     experience     -203.      259.
##  8 ran_vals department informatics experience      302.      245.
##  9 ran_vals department sociology   experience      522.      257.
## 10 ran_vals department statistics  experience     -179.      213.
df %>%
ggplot(aes(
x = experience, y = salary, group = department,
colour = department
)) +
geom_point() +
geom_line(aes(x = experience, y = pred)) +
labs(x = "Experience", y = "Predicted Salary") +
ggtitle("Varying Intercept and Slopes Salary Prediction") +
scale_colour_discrete("Department")

## 56.7 信息池

### 56.7.1 提问

df %>%
group_by(department) %>%
group_modify(
~ broom::tidy(lm(salary ~ 1 + experience, data = .))
)
## # A tibble: 10 × 6
## # Groups:   department [5]
##    department  term        estimate std.error statistic  p.value
##    <chr>       <chr>          <dbl>     <dbl>     <dbl>    <dbl>
##  1 biology     (Intercept)   51297.      937.     54.7  1.81e-21
##  2 biology     experience      399.      208.      1.92 7.13e- 2
##  3 english     (Intercept)   58929.     2046.     28.8  1.65e-16
##  4 english     experience      687.      324.      2.12 4.81e- 2
##  5 informatics (Intercept)   69934.     1327.     52.7  3.55e-21
##  6 informatics experience     1420.      284.      4.99 9.40e- 5
##  7 sociology   (Intercept)   40988.     1047.     39.2  7.13e-19
##  8 sociology   experience     1715.      227.      7.57 5.33e- 7
##  9 statistics  (Intercept)   77470.     1630.     47.5  2.24e-20
## 10 statistics  experience      771.      315.      2.45 2.49e- 2

### 56.7.2 信息共享

• 完全共享
broom::tidy(m1)
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   59509.     2382.     25.0  2.75e-44
## 2 experience     1033.      467.      2.21 2.93e- 2
complete_pooling <-
broom::tidy(m1) %>%
dplyr::select(term, estimate) %>%
tidyr::pivot_wider(
names_from = term,
values_from = estimate
) %>%
dplyr::rename(Intercept = (Intercept), slope = experience) %>%
dplyr::mutate(pooled = "complete_pool") %>%
dplyr::select(pooled, Intercept, slope)

complete_pooling
## # A tibble: 1 × 3
##   pooled        Intercept slope
##   <chr>             <dbl> <dbl>
## 1 complete_pool    59509. 1033.
• 部分共享
fix_effect <- broom.mixed::tidy(m4, effects = "fixed")
fix_effect
fix_effect$estimate[1] fix_effect$estimate[2]
var_effect <- broom.mixed::tidy(m4, effects = "ran_vals")
var_effect
# random effects plus fixed effect parameters
partial_pooling <- var_effect %>%
dplyr::select(level, term, estimate) %>%
tidyr::pivot_wider(
names_from = term,
values_from = estimate
) %>%
dplyr::rename(Intercept = (Intercept), estimate = experience) %>%
dplyr::mutate(
Intercept = Intercept + fix_effect$estimate[1], estimate = estimate + fix_effect$estimate[2]
) %>%
dplyr::mutate(pool = "partial_pool") %>%
dplyr::select(pool, level, Intercept, estimate)

partial_pooling
partial_pooling <-
coef(m4)\$department %>%
tibble::rownames_to_column() %>%
dplyr::rename(level = rowname, Intercept = (Intercept), slope = experience) %>%
dplyr::mutate(pooled = "partial_pool") %>%
dplyr::select(pooled, level, Intercept, slope)

partial_pooling
##         pooled       level Intercept     slope
## 1 partial_pool     biology  50845.60  541.3142
## 2 partial_pool     english  58389.97  780.7522
## 3 partial_pool informatics  70418.78 1286.2142
## 4 partial_pool   sociology  41857.11 1505.3033
## 5 partial_pool  statistics  77282.90  805.1879
• 不共享
no_pool <- df %>%
dplyr::group_by(department) %>%
dplyr::group_modify(
~ broom::tidy(lm(salary ~ 1 + experience, data = .))
)
no_pool
## # A tibble: 10 × 6
## # Groups:   department [5]
##    department  term        estimate std.error statistic  p.value
##    <chr>       <chr>          <dbl>     <dbl>     <dbl>    <dbl>
##  1 biology     (Intercept)   51297.      937.     54.7  1.81e-21
##  2 biology     experience      399.      208.      1.92 7.13e- 2
##  3 english     (Intercept)   58929.     2046.     28.8  1.65e-16
##  4 english     experience      687.      324.      2.12 4.81e- 2
##  5 informatics (Intercept)   69934.     1327.     52.7  3.55e-21
##  6 informatics experience     1420.      284.      4.99 9.40e- 5
##  7 sociology   (Intercept)   40988.     1047.     39.2  7.13e-19
##  8 sociology   experience     1715.      227.      7.57 5.33e- 7
##  9 statistics  (Intercept)   77470.     1630.     47.5  2.24e-20
## 10 statistics  experience      771.      315.      2.45 2.49e- 2
un_pooling <- no_pool %>%
dplyr::select(department, term, estimate) %>%
tidyr::pivot_wider(
names_from = term,
values_from = estimate
) %>%
dplyr::rename(Intercept = (Intercept), slope = experience) %>%
dplyr::mutate(pooled = "no_pool") %>%
dplyr::select(pooled, level = department, Intercept, slope)

un_pooling
## # A tibble: 5 × 4
## # Groups:   level [5]
##   pooled  level       Intercept slope
##   <chr>   <chr>           <dbl> <dbl>
## 1 no_pool biology        51297.  399.
## 2 no_pool english        58929.  687.
## 3 no_pool informatics    69934. 1420.
## 4 no_pool sociology      40988. 1715.
## 5 no_pool statistics     77470.  771.

### 56.7.3 可视化

library(ggrepel)

un_pooling %>%
dplyr::bind_rows(partial_pooling) %>%
ggplot(aes(x = Intercept, y = slope)) +
purrr::map(
c(seq(from = 0.1, to = 0.9, by = 0.1)),
.f = function(level) {
stat_ellipse(
geom = "polygon", type = "norm",
size = 0, alpha = 1 / 10, fill = "gray10",
level = level
)
}
) +
geom_point(aes(group = pooled, color = pooled)) +
geom_line(aes(group = level), size = 1 / 4) +
# geom_point(data = complete_pooling, size = 4, color = "red") +
geom_text_repel(
data = . %>% filter(pooled == "no_pool"),
aes(label = level)
) +
scale_color_manual(
name = "information pool",
values = c(
"no_pool" = "black",
"partial_pool" = "red" # ,
# "complete_pool" = "#A65141"
),
labels = c(
"no_pool" = "no share",
"partial_pool" = "partial share" # ,
# "complete_pool" = "complete share"
)
) #+
## Warning: Using size aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use linewidth instead.
## This warning is displayed once every 8 hours.
## Call lifecycle::last_lifecycle_warnings() to see where this warning was
## generated.
# theme_classic()

## 56.8 更多

• 解释模型的含义
lmer(salary ~ 1 + (0 + experience | department), data = df)
# vs
lmer(salary ~ 1 + experience + (0 + experience | department), data = df)
lmer(salary ~ 1 + (1 + experience | department), data = df)
## Linear mixed model fit by REML ['lmerMod']
## Formula: salary ~ 1 + (1 + experience | department)
##    Data: df
## REML criterion at convergence: 1927.814
## Random effects:
##  Groups     Name        Std.Dev. Corr
##  department (Intercept) 15656
##             experience   1070    -0.44
##  Residual                3313
## Number of obs: 100, groups:  department, 5
## Fixed Effects:
## (Intercept)
##       66047
# vs
lmer(salary ~ 1 + (1 | department) + (0 + experience | department), data = df)
## Linear mixed model fit by REML ['lmerMod']
## Formula: salary ~ 1 + (1 | department) + (0 + experience | department)
##    Data: df
## REML criterion at convergence: 1928.025
## Random effects:
##  Groups       Name        Std.Dev.
##  department   (Intercept) 14391
##  department.1 experience   1064
##  Residual                  3314
## Number of obs: 100, groups:  department, 5
## Fixed Effects:
## (Intercept)
##       60018