🧐 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: 不等价！
最后祝大家早日不卷!~
腾讯云开发者
