使用sjPlot绘制对角随机协方差结构混合模型时遇错误求助

问题根源分析 #

你遇到的Error in gregexpr(pattern = "\|", re)[[1]] : subscript out of bounds错误，主要有两个可能原因：

1. 变量名拼写错误：plot_model的terms参数中使用了"V1","V2"，但你的数据和模型中对应的变量是IV1和IV2，这会导致变量匹配异常；
2. sjPlot对复杂随机结构的解析问题：当nlme模型使用pdDiag定义包含交互项的随机效应时，sjPlot内部的正则表达式解析逻辑可能出现冲突。

分步解决方案 #

方案1：修正变量名并更新sjPlot版本 #

首先修正变量名拼写错误，同时确保sjPlot和依赖包是最新版本，这能解决大部分兼容性问题：

# 安装最新版依赖包
install.packages(c("sjPlot", "nlme", "sjlabelled"))

# 重新加载包并运行完整流程
library(sjPlot)
library(nlme)

set.seed(42)  
n <- 90
data <- data.frame ( id = rep(1:6, each = 15), 
                  DV=rnorm(n),
                  IV1=rnorm(n),
                  IV2=rnorm(n),
                  IV3=rnorm(n),
                  IV4=rnorm(n),
                  IV5=rnorm(n)
                  )

model<- lme(fixed= DV~ IV1*IV2+ IV3 + IV4 + IV5,
            control=list(opt="nlminb"), data=data,
            random = list(id=pdDiag(~1+IV1*IV2)),
            na.action=na.omit) 

# 修正terms参数的变量名
plot_model(model, type="pred", terms=c("IV1","IV2"), mdrt.values="meansd")

方案2：手动用emmeans计算预测值后绘图 #

如果更新包后仍报错，可绕过sjPlot的自动解析，用emmeans手动计算预测值，再用ggplot绘图，这是更稳定的替代方案：

library(emmeans)
library(ggplot2)

# 生成均值±标准差的自变量取值（对应mdrt.values="meansd"）
iv1_vals <- with(data, c(mean(IV1)-sd(IV1), mean(IV1), mean(IV1)+sd(IV1)))
iv2_vals <- with(data, c(mean(IV2)-sd(IV2), mean(IV2), mean(IV2)+sd(IV2)))

# 计算预测边际均值及置信区间
emm_results <- emmeans(model, ~ IV1*IV2, at = list(IV1 = iv1_vals, IV2 = iv2_vals))
emm_df <- as.data.frame(emm_results)

# 绘制预测曲线
ggplot(emm_df, aes(x = IV1, y = emmean, color = factor(round(IV2,2)), group = factor(round(IV2,2)))) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = factor(round(IV2,2)), color = NULL), alpha = 0.2) +
  labs(x = "IV1", y = "Predicted DV", color = "IV2", fill = "IV2") +
  theme_minimal()

方案3：调整随机结构（可选） #

如果研究设计允许，可尝试将随机结构改为pdIdent（同方差对角结构），这能避免sjPlot的解析冲突，但注意这会改变模型的随机效应假设：

model <- lme(fixed= DV~ IV1*IV2+ IV3 + IV4 + IV5,
             control=list(opt="nlminb"), data=data,
             random = list(id=pdIdent(~1+IV1*IV2)),  # 改用pdIdent
             na.action=na.omit) 

plot_model(model, type="pred", terms=c("IV1","IV2"), mdrt.values="meansd")

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