lme4 | 这个线性模型对你来说可能更合理

1写在前面

在进行数据分析时，我们可能经常会遇到分层的数据结构，指每一次观察属于某个特定的组。
比如考察老师的教学成果，而这些老师属于某个班，班又属于某个学校。

2用到的包

rm(list = ls())
library(tidyverse)
library(lme4)
library(modelr)
library(broom)

3示例数据

这里我们模拟一个数据，数据描述的是不同部门(department)的老师的收入(salary)情况.

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(broom.mixed)
df <- create_data()

---

Note! 参数依次为：
✅ "ids" → 教师代号
✅ "department → 部门
✅ "experience" → 经验
✅ "raises" → 增长率
✅ "salary" → 收入

4线性模型

4.1 建模

假设老师的salary主要取决于他的experience，则每个department的base是一样的，并有相同的年度raise。
那么，对salary的评估就是一个简单线性模型：

m1 <- lm(salary ~ experience, data = df)
m1
broom::tidy(m1)

---

接着我们再加入预测值，对比一下。

df %>% modelr::add_predictions(m1)

---

4.2 可视化

df %>%
  add_predictions(m1) %>%
  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_color_npg()

这种线性模型的建模方法过于粗暴，对于每个老师，不管来自哪个department，系数α和β都是是一样的，是固定的，因此这种简单线性模型也称之为固定效应模型。
显然，这样的建模方法并不合理。

5变化的截距

5.1 建模

我们先假设不同的department的base不同，但raise是相同的。
这时，截距会随所在department不同而变化，统计模型写为：

截距α代表着bases，随department变化，每个学院对应一个值，称之为随机效应项。
模型中既有固定效应项又有随机效应项，因此称之为混合效应模型(mixed)。
这里，老师i所在department为j, 描述为j[i], 其所在department的α就表述为 。

---

m2 <- lmer(salary ~ experience + (1 | department), data = df)
m2

broom.mixed::tidy(m2, effects = "ran_vals") # "ran_vals" = random effect values

---

5.2 可视化

df %>%
  add_predictions(m2) %>%
  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_color_npg()

这个时候，我们就能看到department不同带来的salary的差别啦！

6变化的斜率

6.1 建模

我们先假设不同的department的base相同，但raise是不同的，与之前正好相反。 这时，斜率会随所在department不同而变化，统计模型写为:

这里，α对所有老师而言是固定不变的，而β会随department不同而变化，不同department对应不同β。

m3 <- lmer(salary ~ experience + (0 + experience | department), data = df)
m3
broom.mixed::tidy(m3, effects = "ran_vals")

---

6.2 可视化

df %>%
  add_predictions(m3) %>%
  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_color_npg()

7变化的斜率和截距

7.1 建模

我们先假设不同的department的base不同，而且raise也是不同的。 这时，斜率和截距都会随所在department不同而变化，统计模型写为:

这里可以解释为具体来说，教师i，所在的department为j, base表示为 ，raise表示为 .

m4 <- lmer(salary ~ experience + (1 + experience | department), data = df)
m4

broom.mixed::tidy(m4, effects = "ran_vals")

---

7.2 可视化

df %>%
  add_predictions(m4) %>%
  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_color_npg()

8. 小彩蛋

Q: 不同的department的base不同，raise也不同，我们得出不同的α和β。
可否等价为，先按照department分组，然后分别计算α和β。
A: 不等价！

---

最后祝大家早日不卷!~

---

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

本文由 mdnice 多平台发布