使用sae包direct函数时纳入权重致估计量级异常的问题

我在用R语言的sae包进行小区域估计时，发现当在sae::direct()函数中指定sweight权重参数后，得到的估计值量级出现明显异常。以下是完整可复现代码，同时附上了emdi包的计算结果作为对比，求问问题原因？

packages <- c("emdi", "sae", "sf", "sp", "ggplot2", "dplyr", "SUMMER")

for (pkg in packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
    message(sprintf("✅ Installed package: %s", pkg))
  } else {
    message(sprintf("✔ Package already installed: %s", pkg))
  }
  
  # Load the package
  library(pkg, character.only = TRUE)
  message(sprintf("📦 Loaded package: %s", pkg))
}

# Load data sets
data("eusilcA_pop")
data("eusilcA_smp")
data("eusilcA_popAgg")
data("eusilcA_prox")


p <- 8500000 # total population
eusilcA_popAgg$total_pop <- eusilcA_popAgg$ratio_n * p # calculate pop per district
pop <- eusilcA_popAgg %>%select ("Domain","total_pop")%>%filter(Domain %in% eusilcA_smp$district)

# We make sure that district  level are ordered in the same way in the two dataframe
names(pop)=c('district', 'poptot')
pop$district=droplevels(pop$district)
# alphabetical ordering
eusilcA_smp$district= factor( eusilcA_smp$district, levels = sort(levels( eusilcA_smp$district)))
# simple check 
levels(eusilcA_smp$district)
levels(pop$district)


sae.DIR <- sae::direct(y = eusilcA_smp$eqIncome, 
                       dom =eusilcA_smp$district ,
                       #sweight = eusilcA_smp$weight, 
                       domsize = pop[,c(1,2)]) |>
 select(Domain, Direct, SD)


# Direct estimator using SAE
direct_sae <- sae::direct(y = eqIncome, dom = district, data = eusilcA_smp, domsize=pop)

# Direct means by district
ans1=eusilcA_smp%>%
 group_by(district) %>% 
   summarise(
    directmanual = mean(eqIncome),
    .groups = "drop"
  ) |> 
  as.data.frame()


result1 <- sae::direct(y=eqIncome, dom=district, sweight=weight,
                   domsize=pop, data= eusilcA_smp)

head(result1)

result3 <-sae::direct(y=eqIncome, dom=district, domsize=pop[,c(1,2)],
                  data=eusilcA_smp)

head(result3)

# Alternative with emdi:::direct 
library(emdi)
library(laeken)
data("eusilcA_smp")

# Example 1: With weights and naive bootstrap
emdi_direct <-emdi::direct(
  y = "eqIncome", smp_data = eusilcA_smp,
  smp_domains = "district", weights = "weight", threshold = 11064.82,
  var = TRUE, boot_type = "naive", B = 50, seed = 123, X_calib = NULL,
  totals = NULL, na.rm = TRUE
)

head(emdi_direct$ind$Mean)

运行结果 #

head(result1)
                Domain SampSize    Direct        SD       CV
3            Amstetten       33  44.70926  8.009790 17.91528
4                Baden       40  67.41675 10.483811 15.55075
67             Bludenz       17  35.49880  8.695522 24.49526
39      Braunau am Inn       29  35.85324  6.747762 18.82051
68             Bregenz       34 105.09176 17.895180 17.02815
5  Bruck an der Leitha       27  69.86268 13.163270 18.84163

> result3 <-sae::direct(y=eqIncome, dom=district, domsize=pop[,c(1,2)],
+                   data=eusilcA_smp)
>
> head(result3)
                Domain SampSize   Direct        SD       CV
3            Amstetten       33 15201.15 1141.5632 7.509718
4                Baden       40 22921.69 1007.5916 4.395799
67             Bludenz       17 12069.59 1127.2980 9.339984
39      Braunau am Inn       29 12190.10  876.9662 7.194085
68             Bregenz       34 35731.20 1928.0847 5.396082
5  Bruck an der Leitha       27 23753.31 1207.9449 5.085375

问题原因 #

sae::direct()函数的sweight参数默认是抽样权重的倒数，而非常规的抽样权重。你传入的eusilcA_smp$weight是标准抽样权重（代表每个样本对应的总体数量），直接传入会被当作倒数处理，导致估计值量级大幅缩小。

解决方法 #

有两种方式修正这个问题：

• 方法1：将权重取倒数后传入sweight
• 方法2：使用weight_type = "sampling"参数明确指定传入的是抽样权重

方法1代码示例

result1_fixed <- sae::direct(y=eqIncome, dom=district, 
                             sweight=1/weight,  # 对权重取倒数
                             domsize=pop, data= eusilcA_smp)

方法2代码示例

result1_fixed <- sae::direct(y=eqIncome, dom=district, 
                             sweight=weight, 
                             weight_type = "sampling",  # 明确权重类型
                             domsize=pop, data= eusilcA_smp)

修正后，估计值的量级会和不带权重的result3以及emdi包的结果一致。

内容的提问来源于Stack Exchange，提问作者Silchara