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

1写在前面

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

Q: 不同的department的base不同，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! 大家不要认为差异不大就随便选用建模方法，当我们纳入更多变量的时候，可能不同建模方式的差异就会显现啦！🥰

---

最后祝大家早日不卷!~

---

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