在函数中使用svycontrast时含反引号与I()的对比问题

问题描述 #

编写实现survey包峰度估算的svykurt函数时，在svycontrast中使用bquote传递I(x^2)这类带反引号的表达式，触发错误：

Error in eval(e[[2L]], where) : object 'I(x^2)' not found

问题根源 #

bquote的插值语法.无法直接识别带反引号的变量名。在.(I(x^2))中，R会尝试寻找名为I(x^2)的独立对象，而非引用momnts结果中对应的列名，导致变量未找到的错误。

解决思路与修改方案 #

1. 用字符串构建对比表达式（推荐） #

放弃bquote，直接通过字符串拼接生成对比公式，再用parse(text=...)转换为R可识别的表达式，这样能正确保留带反引号的列名。

修改核心代码部分：

# 用字符串拼接对比表达式，替换为实际变量名
mu4_str <- paste0("-3 * (", var_name, ")^4 + 6 * (", var_name, ")^2 * `I(", var_name, "^2)` - 4 * (", var_name, ") * `I(", var_name, "^3)` + `I(", var_name, "^4)`")
sigma2_str <- paste0("`I(", var_name, "^2)` - (", var_name, ")^2")

centrl_momnts <- svycontrast(
  momnts,
  list(
    mu4 = parse(text = mu4_str)[[1]],
    sigma2 = parse(text = sigma2_str)[[1]]
  )
)

如果使用glue包，拼接会更直观：

library(glue)
mu4_str <- glue("-3 * ({var_name})^4 + 6 * ({var_name})^2 * `I({var_name}^2)` - 4 * ({var_name}) * `I({var_name}^3)` + `I({var_name}^4)`")
sigma2_str <- glue("`I({var_name}^2)` - ({var_name})^2")

2. 修正原函数的缺失值处理逻辑 #

原函数中手动过滤x的缺失值（x <- x[!is.na(x)]）会破坏抽样设计的权重和结构，正确的做法是依赖svymean的na.rm参数来处理缺失值，无需手动修改变量。

完整修改后的函数 #

svykurt <- function(
    x,
    design,
    na.rm = FALSE,
    excess = TRUE
) {
  
  if (!inherits(design, "survey.design"))
    stop("design is not a survey design")
  
  var_name <- as.character(x)[2]
  
  # 直接构建原始矩的公式
  momnts_fmla <- as.formula(paste0(
    "~",
    var_name,
    " + I(",
    var_name,
    "^2) + I(",
    var_name,
    "^3) + I(",
    var_name,
    "^4)"))
  
  momnts <- svymean(
    momnts_fmla,
    design,
    na.rm = na.rm
  )
  
  # 字符串拼接中心矩的对比表达式
  mu4_str <- paste0("-3 * (", var_name, ")^4 + 6 * (", var_name, ")^2 * `I(", var_name, "^2)` - 4 * (", var_name, ") * `I(", var_name, "^3)` + `I(", var_name, "^4)`")
  sigma2_str <- paste0("`I(", var_name, "^2)` - (", var_name, ")^2")
  
  centrl_momnts <- svycontrast(
    momnts,
    list(
      mu4 = parse(text = mu4_str)[[1]],
      sigma2 = parse(text = sigma2_str)[[1]]
    )
  )
  
  # 计算峰度（支持超额峰度选项）
  kurtosis <- svycontrast(centrl_momnts, list(kurt = ~ mu4 / sigma2^2))
  if (excess) {
    kurtosis <- svycontrast(kurtosis, list(excess_kurt = ~ kurt - 3))
  }
  
  return(kurtosis)
}

测试运行 #

data(api)
dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
svykurt(x = ~api00, dclus1, na.rm = TRUE)

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