逻辑多层模型多概率分组可视化实现技术求助

问题概述 #

需要拟合并可视化一个逻辑多层模型：

• 因变量：二元就业状态（empl，0=失业/1=就业）
• 自变量：连续型互联网连接水平（isoc）
• 核心结构：以教育水平（educ，1=低技能/2=中等技能/3=高技能）为分组的随机截距+随机斜率，含其他控制变量
• 工具：lme4::glmer()

样本数据与初始代码 #

加载依赖并生成样本数据：

library(lme4)
library(lmerTest)
library(tidyverse)
library(sjPlot)
library(sjmisc)

set.seed(1212)
d <- data.frame(empl=c(1,1,1,0,1,0,1,1,0,1,1,1,0,1,0,1,1,1,1,0), 
      isoc=runif(20, min=-0.2, max=0.2), 
      educ=sample(1:3, 20, replace=TRUE))

样本数据预览：

empl         isoc educ
1     1  0.078604108    1
2     1  0.093667591    3
3     1 -0.061523272    2
4     0  0.009468908    3
5     1 -0.169220134    2
6     0 -0.038594789    3
7     1  0.170506490    1
8     1 -0.098487991    1
9     0  0.073339737    1
10    1  0.144211813    3
11    1 -0.133510687    1
12    1 -0.045306606    3
13    0  0.124211903    3
14    1 -0.003908486    3
15    0 -0.080673475    3
16    1  0.061406993    3
17    1  0.015401951    2
18    1  0.073351501    2
19    1  0.075648137    2
20    0  0.041450192    1

初始模型与可视化问题 #

拟合正确的随机效应模型：

m <- glmer(empl ~ isoc + (1 + isoc | educ),
            data=d,
            family=binomial("logit"),
            nAGQ = 0)

使用sjPlot::plot_model()生成分面子图时，结果无法解读：

plot_model(m, type="pred",
           terms=c("isoc [all]", "educ"),
           show.data=TRUE)

尝试修改模型（存在语法错误与共线性风险）后，仍无法得到正确子图：

# 错误模型：语法错误，educ未正确纳入固定效应
m <- glmer(empl ~ isoc + educ (1 + isoc | educ),
            data=d,
            family=binomial("logit"),
            nAGQ = 0)

解决方案 #

1. 修正模型结构 #

若需将educ作为固定效应纳入，必须将其转换为分类因子，避免连续变量处理引发的问题：

# 将educ转换为因子，明确分组含义
d$educ <- factor(d$educ, levels=1:3, labels=c("低技能", "中等技能", "高技能"))

# 拟合正确的模型（固定效应+随机截距/斜率）
m <- glmer(empl ~ isoc + educ + (1 + isoc | educ),
            data=d,
            family=binomial("logit"),
            nAGQ = 0)

注：若仅需educ作为随机效应分组，初始模型结构正确，无需额外加入固定效应。

2. 用sjPlot生成正确可视化 #

调整plot_model()参数，确保按教育水平分面展示0-1区间的预测概率：

plot_model(m, type="pred",
           terms=c("isoc", "educ"),  # 去掉[all]，让函数自动处理自变量范围
           show.data=TRUE,
           ylim=c(0,1),  # 强制y轴在0-1区间
           facet.grid=TRUE,  # 按educ分面
           axis.title=c("互联网连接水平", "就业概率"),
           title="不同教育水平下互联网连接与就业概率的关系")

3. 手动计算预测值（ggplot2方案） #

若sjPlot仍有问题，可手动生成预测数据，用ggplot2绘制更灵活的可视化：

# 生成覆盖isoc全范围的预测网格
pred_grid <- expand_grid(
  isoc=seq(min(d$isoc), max(d$isoc), length.out=100),
  educ=unique(d$educ)
)

# 计算包含随机效应的预测概率
pred_grid$prob <- predict(m, newdata=pred_grid, type="response", re.form=NULL)

# 绘制分面子图
ggplot(pred_grid, aes(x=isoc, y=prob)) +
  geom_line(color="#2c3e50", linewidth=1) +
  geom_point(data=d, aes(x=isoc, y=empl), alpha=0.6, color="#e74c3c") +
  facet_wrap(~educ) +
  ylim(0,1) +
  labs(x="互联网连接水平", y="就业概率") +
  theme_bw()

内容的提问来源于stack exchange，提问作者ghost_trevor