如何在R Markdown中用table1包生成带95%置信区间的统计表格

问题描述 #

我正在使用R Markdown，基于melanoma数据集，通过table1包开展场景模拟。请问能否在自定义渲染中，为连续变量在均值和标准差下方以[下限CI, 上限CI]的格式添加95%置信区间？相关代码如下：

library(boot) 

melanoma2 <- melanoma
 
# Factor the basic variables that
# we're interested in

    melanoma2$status <- 
      factor(melanoma2$status, 
             levels=c(2,1,3),
             labels=c("Alive", # Reference
                      "Melanoma death", 
                      "Non-melanoma death"))

    labels <- list(
        variables=list(sex="Sex",
                       age="Age (years)",
                       ulcer="Ulceration",
                       thickness="Thicknessᵃ (mm)"),
        groups=list("", "", "Death"))

# Remove the word "death" from the labels, since it now appears above

    levels(melanoma2$status) <- c("Alive", "Melanoma", "Non-melanoma")

    strata <- c(list(Total=melanoma2), split(melanoma2, melanoma2$status))

    my.render.cont <- function(x) {
        with(stats.apply.rounding(stats.default(x), digits=2), c("",
            "Mean (SD)"=sprintf("%s (&plusmn; %s)", MEAN, SD)))
    }
    my.render.cat <- function(x) {
        c("", sapply(stats.default(x), function(y) with(y,
            sprintf("%d (%0.0f %%)", FREQ, PCT))))

}

解决方案 #

完全可以实现这个需求，只需要修改自定义的连续变量渲染函数my.render.cont，在函数内计算均值的95%置信区间，并将其格式化为[下限CI, 上限CI]的样式添加到输出行中。

以下是修改后的完整代码：

library(boot)
library(table1) # 确保加载table1包

melanoma2 <- melanoma
 
# 转换感兴趣变量的因子类型
melanoma2$status <- 
  factor(melanoma2$status, 
         levels=c(2,1,3),
         labels=c("Alive", # 参考组
                  "Melanoma death", 
                  "Non-melanoma death"))

labels <- list(
    variables=list(sex="Sex",
                   age="Age (years)",
                   ulcer="Ulceration",
                   thickness="Thicknessᵃ (mm)"),
    groups=list("", "", "Death"))

# 调整分组标签，移除重复的"death"
levels(melanoma2$status) <- c("Alive", "Melanoma", "Non-melanoma")

strata <- c(list(Total=melanoma2), split(melanoma2, melanoma2$status))

# 修改连续变量渲染函数，添加95%置信区间
my.render.cont <- function(x) {
    # 计算基础统计量（均值、标准差）
    base_stats <- stats.apply.rounding(stats.default(x), digits=2)
    # 计算均值的95%置信区间
    ci <- t.test(x)$conf.int
    ci_rounded <- round(ci, 2)
    ci_str <- sprintf("[%s, %s]", ci_rounded[1], ci_rounded[2])
    
    # 组合输出内容
    c("",
      "Mean (SD)"=sprintf("%s (&plusmn; %s)", base_stats$MEAN, base_stats$SD),
      "95% CI"=ci_str)
}

my.render.cat <- function(x) {
    c("", sapply(stats.default(x), function(y) with(y,
        sprintf("%d (%0.0f %%)", FREQ, PCT))))
}

# 生成表格
table1(strata, labels=labels, render.continuous=my.render.cont, render.categorical=my.render.cat)

关键说明 #

1. 置信区间计算：这里用t.test()函数直接计算均值的95%置信区间，适用于正态分布的连续变量；如果数据偏态，也可以用boot包做bootstrap置信区间，示例代码如下：

   # 用bootstrap计算置信区间的替代方法
   boot_mean <- function(data, indices) {
       mean(data[indices], na.rm=TRUE)
   }
   boot_result <- boot(data=x, statistic=boot_mean, R=1000)
   ci_boot <- boot.ci(boot_result, type="bca")$bca[4:5]
   ci_str <- sprintf("[%s, %s]", round(ci_boot[1],2), round(ci_boot[2],2))
   

2. 格式调整：通过sprintf()将置信区间格式化为[下限, 上限]的样式，和均值标准差行并列展示；
3. 表格生成：最后调用table1()函数时，指定自定义的渲染函数即可生效。

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