jamesbang

V1

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

🧐 lme4 | 多层线性模型小彩蛋解答

1写在前面

本期回答一下上期中小彩蛋部分的问题,如下:

Q: 不同的departmentbase不同,raise也不同,我们得出不同的αβ
可否等价为,先按照department分组,然后分别计算αβ
A: 不等价!

2用到的包

rm(list = ls())
library(tidyverse)
library(lme4)
library(modelr)
library(broom)
library(ggsci)
library(broom.mixed)

3示例数据

数据描述的是不同部门(department)的老师的收入(salary)情况。
具体可见上期

4多层线性模型

4.1 建模

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

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

4.2 加上predictions

df1 <- df %>% 
add_predictions(m1)

df1

4.3 可视化

p1 <- df1 %>%
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()

p1

5分组建立简单线性回归模型

5.1 批量建模

这里我们使用nest函数容纳一下建模的大量数据。 这里我就不做具体展示了,大家可以自己试一下。

m2 <- df %>%
group_by(department) %>%
nest() %>%
mutate(mdl = map(data, ~ lm(salary ~ 1 + experience, data=.))) %>%
mutate(fit = map(mdl, ~ .$fitted.values))

m2

5.2 加上predictions

df2 <- m2 %>% 
mutate(., data = map2(data, mdl, add_predictions)) %>%
select(., -mdl, -fit) %>%
unnest()

df2

5.3 可视化

p2 <- df2 %>%
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()

p2

6比较两种模型

6.1 比较pred值

差异不是很大,但还是不同的。🤣

df3 <- df1 %>%
dplyr::select(.,ids,pred) %>%
left_join(.,df2[,c(2,7)],by = "ids")

df3

6.2 可视化

library(patchwork)
p1 + p2

Note! 大家不要认为差异不大就随便选用建模方法,当我们纳入更多变量的时候,可能不同建模方式的差异就会显现啦!🥰


最后祝大家早日不卷!~

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

分类:

后端

标签:

后端

作者介绍

jamesbang
V1