R语言merlin包mlrcs模型预测/对比实现及merlin语法问题

一、问题背景 #

使用mlrcs拟合含随机效应的限制性立方样条模型后，无法调用predict()函数获取预测值；直接替换为merlin()函数时，因语法错误报错。需实现模型预测值获取，以及不同time水平下location的差异对比分析。

二、mlrcs模型的手动预测方法 #

mlrcs没有内置predict方法，可通过手动计算实现预测，步骤如下：

1. 提取模型系数
2. 构造包含所有time和location组合的新数据集
3. 生成限制性立方样条设计矩阵，结合系数计算预测值

代码示例：

# 提取模型系数
coefs <- coef(mod)

# 构造预测用数据集（覆盖所有time、location组合，随机效应取0时用任意fish）
new_dat <- expand.grid(time = unique(dat$time),
                       location = unique(dat$location),
                       fish = unique(dat$fish)[1])

# 生成限制性立方样条设计矩阵
rcs_mat <- rms::rcspline.eval(new_dat$time, knots = attr(mod$terms, "knots"), inclx = FALSE)
# 合并设计矩阵与其他变量（对应模型的固定效应项）
new_dat_design <- cbind(1, rcs_mat, new_dat$location, new_dat$time * new_dat$location)

# 计算固定效应预测值（不含随机效应，对应组水平趋势）
new_dat$pred_fixed <- new_dat_design %*% coefs[c("rcs():1", "rcs():2", "rcs():3", "location", "time:location", "_cons")]

三、merlin函数的正确语法 #

merlin不支持random参数，需通过f()函数在模型公式中定义随机效应，正确语法如下：

library(merlin)
library(rms)

# 重新拟合含随机截距的限制性立方样条模型
mod2 <- merlin(
  model = value ~ 1 + rcs(time, 3) + location + time:location + f(fish, model = "iid"),
  data = dat,
  family = "gaussian" # 因变量为连续型，指定高斯分布族
)

summary(mod2)

说明：f(fish, model = "iid")对应原模型的~1|fish，表示为fish定义独立同分布的随机截距；需显式指定family确保与mlrcs模型拟合逻辑一致。

四、merlin模型的预测与差异对比 #

1. 获取预测值 #

merlin支持predict()方法，可指定是否包含随机效应：

# 构造覆盖所有time、location、fish组合的预测数据集
pred_dat <- expand.grid(time = 1:6,
                        location = c(0,1),
                        fish = unique(dat$fish))

# 固定效应预测（组水平趋势，不含随机效应）
pred_dat$pred_fixed <- predict(mod2, newdata = pred_dat, type = "response", include.random = FALSE)

# 包含随机效应的预测（个体水平趋势）
pred_dat$pred_random <- predict(mod2, newdata = pred_dat, type = "response", include.random = TRUE)

2. 不同time下location的差异对比 #

使用merlin的lincom()函数，检验每个time点的location差异：

# 遍历所有time点，输出location差异的统计检验结果
for(t in unique(dat$time)){
  # 构造对比式：location=1与location=0的预测值差异 = location + t*location
  contrast <- paste0("location + ", t, "*location")
  cat("Time =", t, ":\n")
  print(lincom(mod2, contrast))
  cat("\n")
}

也可通过整理预测结果直接计算差异：

library(dplyr)
library(tidyr)

# 计算每个time下location的预测值差异
pred_summary <- pred_dat %>%
  group_by(time, location) %>%
  summarise(pred_mean = mean(pred_fixed), .groups = "drop") %>%
  pivot_wider(names_from = location, values_from = pred_mean) %>%
  mutate(diff = `1` - `0`)

print(pred_summary)

五、预测趋势绘图示例 #

ggplot(pred_dat, aes(x = time, y = pred_fixed, color = factor(location))) +
  geom_line(size = 1) +
  geom_point(data = dat, aes(y = value), alpha = 0.5) +
  labs(x = "Time", y = "Value", color = "Location") +
  theme_bw()

内容的提问来源于stack exchange，提问作者Dr.FishGirl