![alt](https://i-blog.csdnimg.cn/blog_migrate/431e2d46210f405d95e791b69ba56060.png)
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
)情况。
具体可见上期。
![alt](https://i-blog.csdnimg.cn/blog_migrate/33199946e4cdb30d1e253c8e00c9d2aa.png)
4多层线性模型
4.1 建模
m1 <- lmer(salary ~ experience + (1 + experience | department), data = df)
m1
broom.mixed::tidy(m1, effects = "ran_vals")
![alt](https://i-blog.csdnimg.cn/blog_migrate/bab55efebac742723f4863b13a9946ab.png)
4.2 加上predictions
df1 <- df %>%
add_predictions(m1)
df1
![alt](https://i-blog.csdnimg.cn/blog_migrate/400dfa8522c949e6ff30e808023f18a5.png)
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
![alt](https://i-blog.csdnimg.cn/blog_migrate/6b6cfaf1235d38a4ba8ea1fed95e6c69.png)
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
![alt](https://i-blog.csdnimg.cn/blog_migrate/7d124a60ed2dd74344670b834b600c4e.png)
5.2 加上predictions
df2 <- m2 %>%
mutate(., data = map2(data, mdl, add_predictions)) %>%
select(., -mdl, -fit) %>%
unnest()
df2
![alt](https://i-blog.csdnimg.cn/blog_migrate/4a358d2e08725150d0147dba502421f4.png)
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
![alt](https://i-blog.csdnimg.cn/blog_migrate/7d63a3ca1fc9bc128bfe90921c5fb7b8.png)
6比较两种模型
6.1 比较pred值
差异不是很大,但还是不同的。🤣
df3 <- df1 %>%
dplyr::select(.,ids,pred) %>%
left_join(.,df2[,c(2,7)],by = "ids")
df3
![alt](https://i-blog.csdnimg.cn/blog_migrate/1a3e887d1dec44662686f3c6654a7a52.png)
6.2 可视化
library(patchwork)
p1 + p2
![alt](https://i-blog.csdnimg.cn/blog_migrate/9439c187e5565947177ce45bcc81bf35.png)
Note! 大家不要认为差异不大就随便选用建模方法,当我们纳入更多变量的时候,可能不同建模方式的差异就会显现啦!🥰
![](https://i-blog.csdnimg.cn/blog_migrate/91d02c5eaa8ef681a50e636966384ea4.png)
点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
![alt](https://i-blog.csdnimg.cn/blog_migrate/0d97df03d585ae22202511c03905aeaa.png)
本文由 mdnice 多平台发布