为何封装survey包循环逻辑的函数仅返回sexual_orientation最后层级结果？

目标 #

我尝试使用R语言的survey包，针对加利福尼亚州各郡，计算变量sexual_orientation（异性恋、同性恋/女同性恋、双性恋、其他）四个层级的比例及置信区间，以处理复杂抽样设计。预期输出如下（感谢@jpsmith提供的svyby循环代码）：

county        mean         ci_l       ci_u sex_orient
0.A      A 0.911193691 0.8520671385 0.94812694          0
1.A      A 0.036149402 0.0149796639 0.08466537          1
2.A      A 0.042363554 0.0196451241 0.08897020          2
3.A      A 0.010293353 0.0018801543 0.05430515          3
0.B      B 0.923758035 0.8636429959 0.95863968          0
1.B      B 0.036042649 0.0141337843 0.08885182          1
2.B      B 0.038706656 0.0164219903 0.08851038          2
3.B      B 0.001492660 0.0002056946 0.01074521          3

问题 #

但当我将循环代码封装为函数时，仅得到sexual_orientation最后一个层级（“其他”）的比例结果：

0.A      A 0.010293353 0.0018801543 0.05430515                  0
1.A      A 0.010293353 0.0018801543 0.05430515                  1
2.A      A 0.010293353 0.0018801543 0.05430515                  2
3.A      A 0.010293353 0.0018801543 0.05430515                  3
0.B      B 0.001492660 0.0002056946 0.01074521                  0
1.B      B 0.001492660 0.0002056946 0.01074521                  1
2.B      B 0.001492660 0.0002056946 0.01074521                  2
3.B      B 0.001492660 0.0002056946 0.01074521                  3

复现代码 #

library(survey)
library(dplyr)

# 创建数据集
set.seed(130)
probs <- c(.94, .03, .02, .01)
n <- 2000

sogi <- data.frame(t(rmultinom(n, 1, prob = probs)))
names(sogi) <- c("straight","gay","bisexual","other")
df_sogi <- sogi %>% mutate(sexual_orientation =
                             as.factor(case_when(
                               straight %in% 1 ~ 0,
                               gay %in% 1 ~ 1,
                               bisexual %in% 1 ~ 2,
                               other %in% 1 ~ 3))) %>%
  mutate(pw = runif(n, 0.2, 3.5)) %>%
  mutate(county = as.factor(rep(LETTERS[1:10], each = n/length(LETTERS[1:10])))) %>%
  mutate(ind = row_number()) %>%
  select(ind, county, sexual_orientation, pw)

table(df_sogi$county, df_sogi$sexual_orientation)

# 构建复杂抽样设计
df_svy_z1 <- svydesign(id=~ind, data = df_sogi, weights = ~pw, strata = ~county)

尝试的解决方案 #

方案1：传入数据集的函数

get_table <- function(df) {
  svCounty <- list()
  for(i in sort(unique(df[ , "sexual_orientation"]))) {
    svCounty[[i]]  <- svyby(~ I(sexual_orientation == i),
                           design = df_svy_z1,
                           by = ~ county,
                           FUN = svyciprop,
                           vartype = "ci",
                           method = "xlogit") 
    
    svCounty[[i]]$sexual_orientation <- i
    names(svCounty[[i]])[2] <- "mean"
  }
  return(svCounty)
}

do.call(rbind, get_table(df_sogi))[order(do.call(rbind, get_table(df_sogi))$county), ]

方案2：循环遍历性取向值的函数

svCounty <- list()
get_table <- function(i) {
  svCounty[[i]] <- svyby(~ I(sexual_orientation == i),
                         design = df_svy_z1,
                         by = ~ county,
                         FUN = svyciprop,
                         vartype = "ci",
                         method = "xlogit")
  svCounty[[i]]$sex_orient <- i
  names(svCounty[[i]])[2] <- "mean"
  
  svCounty[[i]]
}

rbindlist(lapply(sort(unique(df_sogi$sexual_orientation)), get_table))

修正方法 #

问题根源是公式的延迟求值特性：循环中的I(sexual_orientation == i)不会立即计算i的值，而是等到最后才读取i的最终值（即性取向的最后一个层级）。以下是两种修正方式：

方式1：动态生成公式

通过reformulate动态创建包含当前i值的公式，避免延迟求值问题：

get_table <- function(df) {
  svCounty <- list()
  so_vals <- sort(unique(df$sexual_orientation))
  
  for(i in so_vals) {
    # 动态生成对应i的公式
    formula <- reformulate(paste0("I(sexual_orientation == ", i, ")"))
    svCounty[[as.character(i)]] <- svyby(formula,
                                         design = df_svy_z1,
                                         by = ~ county,
                                         FUN = svyciprop,
                                         vartype = "ci",
                                         method = "xlogit")
    
    svCounty[[as.character(i)]]$sexual_orientation <- i
    names(svCounty[[as.character(i)]])[2] <- "mean"
  }
  
  # 合并结果并按郡和性取向排序
  result <- do.call(rbind, svCounty)
  result[order(result$county, result$sexual_orientation), ]
}

# 调用函数
get_table(df_sogi)

方式2：用local()捕获当前循环变量

在循环内部用local()创建独立环境，捕获当前的i值：

get_table <- function(df) {
  svCounty <- list()
  so_vals <- sort(unique(df$sexual_orientation))
  
  for(i in so_vals) {
    svCounty[[i]] <- local({
      current_i <- i  # 捕获当前循环的i值
      res <- svyby(~ I(sexual_orientation == current_i),
                   design = df_svy_z1,
                   by = ~ county,
                   FUN = svyciprop,
                   vartype = "ci",
                   method = "xlogit")
      res$sexual_orientation <- current_i
      names(res)[2] <- "mean"
      res
    })
  }
  
  result <- do.call(rbind, svCounty)
  result[order(result$county, result$sexual_orientation), ]
}

# 调用函数
get_table(df_sogi)

修正方案2（基于lapply的版本）

去掉全局的svCounty列表，直接在函数内返回当前i对应的结果，同时确保i被正确捕获：

get_table <- function(i) {
  # 动态生成公式避免延迟求值
  formula <- reformulate(paste0("I(sexual_orientation == ", i, ")"))
  res <- svyby(formula,
               design = df_svy_z1,
               by = ~ county,
               FUN = svyciprop,
               vartype = "ci",
               method = "xlogit")
  res$sex_orient <- i
  names(res)[2] <- "mean"
  res
}

# 调用并排序结果
result <- rbindlist(lapply(sort(unique(df_sogi$sexual_orientation)), get_table))
result[order(result$county, result$sex_orient), ]

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