如何筛选同CBSERIAL且同US2021A_SOCP的双人样本以分析夫妻收入差

我希望分析职业相同与不同的夫妻之间的收入差异，目前拥有包含职业（US2021A_SOCP）和家庭编码（CBSERIAL）的个体级数据，数据结构如下：

structure(list(YEAR = structure(c(2021L, 2021L, 2021L, 2021L, 
2021L, 2021L), 
    SAMPLE = structure(c(202101L, 202101L, 202101L, 202101L, 
    202101L, 202101L), 
    CBSERIAL = structure(c(2021010000026, 2021010000031, 2021010000063, 
    2021010000067, 2021010000100, 2021010000227), label = "Original Census Bureau household serial number", 
    HHWT = structure(c(13, 51, 17, 61, 15, 46), label = "Household weight", 
    CLUSTER = structure(c(2021000000011, 2021000000021, 2021000000031, 
    2021000000041, 2021000000051, 2021000000061), label = "Household cluster for variance estimation", 
    STRATA = structure(c(80001, 80001, 120001, 170001, 50001, 
    160001), label = "Household strata for variance estimation", 
    GQ = structure(c(3L, 3L, 3L, 3L, 3L, 4L), labels = structure(0:6, names = c("Vacant unit", 
    "Households under 1970 definition", "Additional households under 1990 definition", 
    "Group quarters--Institutions", "Other group quarters", "Additional households under 2000 definition", 
    "Fragment")), label = "Group quarters status",   
 US2021A_DIVISION = structure(c("6", 
    "6", "6", "6", "6", "6"), labels = c(`New England (Northeast region)` = "1", 
    `Middle Atlantic (Northeast region)` = "2", `East North Central (Midwest region)` = "3", 
    `West North Central (Midwest region)` = "4", `South Atlantic (South region)` = "5", 
    `East South Central (South region)` = "6", `West South Central (South Region)` = "7", 
    `Mountain (West region)` = "8", `Pacific (West region)` = "9"
    ), label = "Division code based on 2010 Census definitions", var_desc = "", 
 US2021A_PUMA = structure(c("00800", 
    "00800", "01200", "01700", "00500", "01600"), label = "Public use microdata area code (PUMA) based on 2010 Census definition (areas with population of 100,000 or more, use with ST for unique code)", 

PERNUM = structure(c(1, 1, 1, 
    1, 1, 1), label = "Person number in sample unit", 
    PERWT = structure(c(13, 51, 17, 61, 15, 46), label = "Person weight", 
    US2021A_AGEP = structure(c("85", "67", "74", "16", "83", 
    "19"), 
 US2021A_MAR = structure(c("1", 
    "2", "2", "5", "3", "5"), 
    
    US2021A_PINCP = structure(c("0015000", "0004800", "0036000", 
    "0000000", "0007200", "0008000"), label = "Total person's income (signed, use ADJINC to adjust to constant dollars)", var_desc = ""), 
    US2021A_SOCP = structure(c("BBBBBB", "BBBBBB", "BBBBBB", 
    "BBBBBB", "BBBBBB", "412031"), 
 PERWT_SP = structure(c(NA_real_, 
    NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), label = "Person weight [of Spouse's location in household]", 
    US2021A_AGEP_SP = structure(c("", "", "", "", "", ""),  label = "Age [of Spouse's location in household]", 
US2021A_SEX_SP = structure(c("", 
    "", "", "", "", ""), labels = c(Male = "1", Female = "2"), label = "Sex [of Spouse's location in household]", 

  US2021A_ESR_SP = structure(c("", 
    "", "", "", "", ""), labels = c(`Civilian employed, at work` = "1", 
    `Civilian employed, with a job but not at work` = "2", Unemployed = "3", 
    `Armed forces, at work` = "4", `Armed forces, with a job but not at work` = "5", 
    `Not in labor force` = "6", `N/A (less than 16 years old)` = "B"
    ), 
US2021A_PERNP_SP = structure(c("", 
    "", "", "", "", ""), label = "Total person's earnings (use ADJINC to adjust to constant dollars) [of Spouse's location in household]", var_desc = ""), 
    US2021A_PINCP_SP = structure(c("", "", "", "", "", ""), label = "Total person's income (signed, use ADJINC to adjust to constant dollars) [of Spouse's location in household]", var_desc = ""), 
    US2021A_SOCP_SP = structure(c("", "", "", "", "", ""), label = "Standard Occupational Classification (SOC) codes for 2018 and later based on 2018 SOC codes [of Spouse's location in household]", 
 class = c("haven_labelled", 
    "vctrs_vctr", "character"))), row.names = c(NA, -6L), class = c("tbl_df", 
"tbl", "data.frame"))

我需要筛选出同一CBSERIAL（家庭）下的夫妻两人组，且两人的US2021A_SOCP相同。尝试过用for循环遍历，但setequal函数不符合需求。

使用tidyverse的dplyr包可以高效处理这类分组匹配问题，无需手动循环，以下提供两种适配不同数据情况的方法：

方法一：利用现有配偶字段（优先选择） #

如果数据中包含US2021A_SOCP_SP这类直接记录配偶职业的字段，直接基于该字段处理：

library(dplyr)

# 加载数据（假设数据对象名为df）
df <- structure(...) # 替换为你的数据结构代码

# 筛选夫妻二人组并判断职业是否相同
couple_matched <- df %>%
  # 统计每个家庭的人数，确保是二人家庭
  add_count(CBSERIAL, name = "household_size") %>%
  # 筛选条件：已婚（根据数据标签，US2021A_MAR=1为已婚）、家庭人数为2、职业编码有效、配偶职业非空
  filter(household_size == 2, 
         US2021A_MAR == "1", 
         US2021A_SOCP != "BBBBBB", 
         US2021A_SOCP_SP != "") %>%
  # 标记夫妻职业是否相同
  mutate(same_occupation = US2021A_SOCP == US2021A_SOCP_SP) %>%
  # 保留核心分析字段
  select(CBSERIAL, US2021A_PINCP, US2021A_PINCP_SP, US2021A_SOCP, US2021A_SOCP_SP, same_occupation)

# 查看结果
print(couple_matched)

方法二：自连接匹配夫妻（适用于配偶字段缺失的情况） #

如果配偶字段不全，通过自连接同一家庭的两个个体来配对：

library(dplyr)

df_clean <- df %>%
  # 清理无效数据：移除缺失职业编码、未成年人、非已婚个体
  filter(US2021A_SOCP != "BBBBBB", 
         as.numeric(US2021A_AGEP) >= 18, 
         US2021A_MAR == "1") %>%
  # 给每个家庭内的个体分配序号，用于后续配对
  group_by(CBSERIAL) %>%
  mutate(person_idx = row_number()) %>%
  ungroup()

# 自连接同一家庭的两个个体，生成夫妻对
couple_pairs <- df_clean %>%
  inner_join(df_clean, by = "CBSERIAL", suffix = c("_self", "_spouse")) %>%
  # 确保配对的是不同个体，且只保留唯一配对（避免重复）
  filter(person_idx_self < person_idx_spouse) %>%
  # 判断夫妻职业是否相同
  mutate(same_occupation = US2021A_SOCP_self == US2021A_SOCP_spouse) %>%
  # 保留核心分析字段
  select(CBSERIAL, 
         US2021A_PINCP_self, US2021A_PINCP_spouse, 
         US2021A_SOCP_self, US2021A_SOCP_spouse, 
         same_occupation)

# 查看结果
print(couple_pairs)

注意事项 #

• 请根据数据实际的婚姻状况编码调整US2021A_MAR的筛选值（比如部分数据中2可能代表已婚，需核对字段标签）；
• BBBBBB属于无效职业编码，需提前过滤，避免干扰分析；
• 如果需要分析收入差异，可以基于same_occupation分组，计算两组的收入均值、差值或进行统计检验。

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