在R中按sn_id分组执行Cochran-Armitage检验（catt函数）报错求助

问题描述 #

我手头有个带is_severe、encoding、sn_id字段的数据集，需要按每个唯一的sn_id分组，跑Cochran-Armitage趋势检验（用catt函数实现）。单独调用CATT(data$is_severe, data$encoding)没问题，但试了by(df, df$sn_id, CATT(df$is_severe, df$encoding))分组运行时，直接报错could not find function 'FUN'。

后来我改用了归档包HapEstXXR里的catt函数（代码如下），手动测函数是正常的，但套到分组数据里还是不行，求解决办法。

使用的catt函数代码 #

catt <-
  function(y, x, score = c(0, 1, 2)) {
    miss <- unique(c(which(is.na(y)), which(is.na(x))))
    n.miss <- length(miss)
    if(n.miss > 0) {
      y <- y[-miss]
      x <- x[-miss]
    }
    if(!all((y == 0) | (y == 1))) 
      stop("y should be only 0 or 1.")
    if(!all((x == 0) | (x == 1) |(x == 2))) 
      stop("x should be only 0, 1 or 2.")
    ca <- x [y == 1]
    co <- x [y == 0]
    htca <- table(ca)
    htco <- table(co)
    A <- matrix(0, 2, 3)
    colnames(A) <- c(0, 1, 2)
    rownames(A) <- c(0, 1)
    A[1, names(htca)] <- htca
    A[2, names(htco)] <- htco
    ptt <- prop.trend.test(A[1, ], colSums(A), score = score)
    res <- list("2x3-table" = A, 
                chisq = as.numeric(ptt$statistic), 
                df = as.numeric(ptt$parameter), 
                p.value = as.numeric(ptt$p.value), 
                n.miss = n.miss)
    return(res)
}

数据集示例与dput信息 #

数据示例 #

is_severe     encoding    sn_id
1           1        1   chr1 14907
2           1        1   chr1 14930
...
12          0        1   chr1 69511

dput结构 #

structure(list(X = 0:4, CHROM = c("chr1", "chr1", "chr1", "chr1", 
"chr1"), POS = c("14907", "14930", "15211", "15274", "16378"), 
    REF = c("A", "A", "T", "A", "T"), ALT = c("G", "G", "G", 
    "T", "C"), AF_VALUE = c(0.5, 0.5, 0.5, 1, 0.5), hetro.homo = c("hetro", 
    "hetro", "hetro", "homo", "hetro"), id_ = c("i_peb107_270", 
    "i_peb107_270", "i_peb107_270", "i_peb107_270", "i_peb107_270"
    ), Phenotype = c("Severe", "Severe", "Severe", "Severe", 
    "Severe"), is_severe = c(1L, 1L, 1L, 1L, 1L), encoding = c(1L, 
    1L, 1L, 2L, 1L), sn_id = c("chr1 14907", "chr1 14930", "chr1 15211", 
    "chr1 15274", "chr1 16378")), row.names = c(NA, 5L), class = "data.frame")

问题原因 #

你用by()的时候犯了两个关键错误：

1. 第三个参数传的是CATT(df$is_severe, df$encoding)——这是直接调用函数的结果，不是函数本身，by()需要的是一个能处理每个分组的函数对象。
2. 没有使用分组后的子数据集的列，而是直接调用原数据集的列，完全失去了分组的意义。

解决方法 #

方法1：修正by()的用法 #

把catt调用包装成匿名函数，让每个分组的子数据作为参数传入，取对应列计算：

# 注意用你定义的小写catt函数
result_by <- by(df, df$sn_id, function(sub_data) {
  catt(sub_data$is_severe, sub_data$encoding)
})
# 查看分组检验结果
result_by

方法2：用dplyr分组（更易读易分析） #

如果习惯tidyverse风格，用dplyr可以把检验结果整理成结构化的数据框，方便后续筛选、可视化：

library(dplyr)

result_df <- df %>%
  group_by(sn_id) %>%
  summarise(
    # 保存完整检验结果
    catt_res = list(catt(is_severe, encoding)),
    # 提取关键统计量到单独列
    chisq = catt_res[[1]]$chisq,
    p_value = catt_res[[1]]$p.value,
    df = catt_res[[1]]$df,
    missing_count = catt_res[[1]]$n.miss,
    .groups = "drop"
  )

# 查看整理后的结构化结果
print(result_df)

额外注意事项 #

• 提前清理数据：确保is_severe只有0/1，encoding只有0/1/2，避免函数报错：

df_clean <- df %>%
  filter(is_severe %in% c(0,1), encoding %in% c(0,1,2))

• 处理极端分组：如果某个sn_id对应的子数据里，is_severe全是0或全是1，prop.trend.test会无法计算，建议提前过滤这类分组：

df_clean <- df_clean %>%
  group_by(sn_id) %>%
  filter(n_distinct(is_severe) == 2) %>%
  ungroup()

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