R自定义函数中emmeans报错：参考网格无p_t变量问题求助

错误原因 #

你在自定义函数mlm_sum中直接使用~ p_t | p_g作为emmeans的公式，但p_t是函数的参数名，并非模型mlm4中实际存在的变量名（模型里的变量是t和grp）。emmeans会把p_t当作模型变量去查找，自然找不到对应的变量，因此报错。

解决方案 #

需要将函数参数中传入的变量符号（比如t、grp）正确传递给emmeans，以下是两种常用修正方法：

方法1：使用rlang包的引用与反引用（推荐） #

通过rlang::enquo()捕获传入的变量符号，再用!!将其插入到emmeans的公式中，让emmeans识别到模型中的实际变量。

方法2：构造字符串公式 #

将变量名转为字符串，拼接成公式字符串后再转换为公式对象。

修正后的函数代码 #

方法1（rlang版本） #

library(rlang)
library(lme4)
library(emmeans)
library(ggplot2)
library(Matrix)

mlm_sum <- function(p_m, p_t, p_g) {
  # 捕获传入的变量符号
  t_var <- enquo(p_t)
  g_var <- enquo(p_g)
  
  sum <- summary(p_m)
  ci <- confint(p_m)
  
  # 方差-协方差矩阵
  mat <- Matrix::bdiag(VarCorr(p_m))
  
  # 诊断图
  pred_res <- data.frame(predicted=predict(p_m), residual=residuals(p_m))
  plot1 <- ggplot(pred_res, aes(x=predicted, y=residual)) + 
    geom_point() + 
    geom_hline(yintercept=0, lty=3)
  plot2 <- ggplot(pred_res, aes(x=residual)) + 
    geom_histogram(bins=20, color="black")
  plot3 <- ggplot(pred_res, aes(sample=residual)) + 
    stat_qq() + 
    stat_qq_line()

  # 估计边际均值：用!!插入捕获的变量符号
  emm1 <- emmeans(p_m, formula = !!t_var | !!g_var, 
                  at = list(!!t_var := c(-700, -365, 0, 365, 700)),
                  pbkrtest.limit = 43788)
  
  emm_plot1 <- emmip(p_m, formula = !!t_var | !!g_var,
                     at = list(!!t_var := c(-700, -365, 0, 365, 700)),
                     pbkrtest.limit = 43788)
  
  return(list(summary = sum, confint = ci, vcov_matrix = mat, 
              res_plot = plot1, hist_plot = plot2, qq_plot = plot3,
              emmeans = emm1, emm_plot = emm_plot1))
}

# 调用函数，传入模型中的实际变量名
mlm_sum(mlm4, t, grp)

方法2（字符串公式版本） #

library(lme4)
library(emmeans)
library(ggplot2)
library(Matrix)

mlm_sum <- function(p_m, p_t, p_g) {
  # 将变量转为字符串
  t_name <- deparse(substitute(p_t))
  g_name <- deparse(substitute(p_g))
  
  sum <- summary(p_m)
  ci <- confint(p_m)
  
  # 方差-协方差矩阵
  mat <- Matrix::bdiag(VarCorr(p_m))
  
  # 诊断图
  pred_res <- data.frame(predicted=predict(p_m), residual=residuals(p_m))
  plot1 <- ggplot(pred_res, aes(x=predicted, y=residual)) + 
    geom_point() + 
    geom_hline(yintercept=0, lty=3)
  plot2 <- ggplot(pred_res, aes(x=residual)) + 
    geom_histogram(bins=20, color="black")
  plot3 <- ggplot(pred_res, aes(sample=residual)) + 
    stat_qq() + 
    stat_qq_line()

  # 构造公式字符串并转换为公式对象
  emm_formula <- as.formula(paste0("~ ", t_name, " | ", g_name))
  # 构造at参数的列表
  at_list <- list()
  at_list[[t_name]] <- c(-700, -365, 0, 365, 700)
  
  emm1 <- emmeans(p_m, formula = emm_formula, 
                  at = at_list,
                  pbkrtest.limit = 43788)
  
  emm_plot1 <- emmip(p_m, formula = emm_formula,
                     at = at_list,
                     pbkrtest.limit = 43788)
  
  return(list(summary = sum, confint = ci, vcov_matrix = mat, 
              res_plot = plot1, hist_plot = plot2, qq_plot = plot3,
              emmeans = emm1, emm_plot = emm_plot1))
}

# 调用函数，传入模型中的实际变量名
mlm_sum(mlm4, t, grp)

额外提示 #

1. 你最初调用函数时写的mlm_sum(mlm4, t, g)，但模型mlm4中的分组变量是grp，需传入正确的变量名，否则会出现类似的变量找不到错误。
2. 返回列表时给每个元素命名（如summary = sum），后续提取结果更方便，比如result$summary就能直接拿到模型摘要。

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