jamesbang

V1

2022/09/09阅读:13主题:雁栖湖

🧐 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,则每个departmentbase是一样的,并有相同的年度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 建模

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

截距α代表着bases,随department变化,每个学院对应一个值,称之为随机效应项
模型中既有固定效应项又有随机效应项,因此称之为混合效应模型(mixed)。
这里,老师i所在departmentj, 描述为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 建模

我们先假设不同的departmentbase相同,但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 建模

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

这里可以解释为具体来说,教师i,所在的departmentj, 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: 不同的departmentbase不同,raise也不同,我们得出不同的αβ
可否等价为,先按照department分组,然后分别计算αβ
A: 不等价!


最后祝大家早日不卷!~

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

分类:

后端

标签:

后端

作者介绍

jamesbang
V1