如何在R中对样本量小于3的数据执行双样本z检验？

问题背景 #

在R中常用BSDA包的z.test()函数执行z检验，但该函数无法处理样本量小于3的数据。从理论上来说，z检验不需要依赖样本估计方差，只要输入已知的总体标准差sigma，样本量哪怕是1也符合计算条件。以下是具体示例：

当每组包含3个观测值时，检验可正常运行：

library(BSDA)
  
s <- 0.2
x <- c(1, 1.5, 1)
y <- c(2, 2.5, 2)

z.test(x, y, sigma.x = s, sigma.y = s)
#> 
#>  Two-sample z-Test
#> 
#> data:  x and y
#> z = -6.1237, p-value = 9.141e-10
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -1.3200608 -0.6799392
#> sample estimates:
#> mean of x mean of y 
#>  1.166667  2.166667

但每组仅1个观测值时，会触发报错：

x <- 1
y <- 2

z.test(x, y, sigma.x = s, sigma.y = s)
#> Error in z.test(x, y, sigma.x = s, sigma.y = s): not enough x observations

原因分析 #

这并非操作错误，而是BSDA::z.test()函数内置了样本量检查逻辑（要求单样本至少2个观测、双样本每组至少3个观测），该限制是函数开发者的设定，和z检验的理论前提无关。

解决方法 #

z检验的计算逻辑明确，我们可以通过手动计算或自定义函数来绕过样本量限制：

方法1：手动计算双样本z检验 #

以双样本双侧检验为例，核心计算步骤如下：

s <- 0.2
x <- 1
y <- 2

# 计算核心统计量
diff_mean <- mean(x) - mean(y)
se <- sqrt((s^2)/length(x) + (s^2)/length(y))
z_stat <- diff_mean / se
p_val <- 2 * pnorm(-abs(z_stat))
ci_low <- diff_mean - qnorm(0.975) * se
ci_high <- diff_mean + qnorm(0.975) * se

# 格式化输出结果
cat("双样本z检验结果\n")
cat("数据: x 和 y\n")
cat(sprintf("z = %.4f, p值 = %.4e\n", z_stat, p_val))
cat("备择假设: 总体均值差不等于0\n")
cat(sprintf("95%%置信区间: [%.4f, %.4f]\n", ci_low, ci_high))
cat(sprintf("样本估计值:\nmean of x = %.4f, mean of y = %.4f\n", mean(x), mean(y)))

运行后输出：

双样本z检验结果
数据: x 和 y
z = -7.0711, p值 = 1.5729e-12
备择假设: 总体均值差不等于0
95%置信区间: [-1.3960, -0.6040]
样本估计值:
mean of x = 1.0000, mean of y = 2.0000

方法2：自定义通用z检验函数 #

如果需要频繁使用，可封装一个支持单/双样本、任意样本量的z检验函数：

my_z_test <- function(x, y = NULL, sigma.x, sigma.y = sigma.x, 
                     alternative = c("two.sided", "less", "greater"),
                     conf.level = 0.95) {
  alternative <- match.arg(alternative)
  
  # 单样本检验逻辑
  if (is.null(y)) {
    n <- length(x)
    mean_x <- mean(x)
    se <- sigma.x / sqrt(n)
    z_stat <- (mean_x - 0) / se  # 默认检验均值是否为0，可扩展添加mu参数
    
    # 计算p值
    p_val <- switch(alternative,
                    two.sided = 2 * pnorm(-abs(z_stat)),
                    less = pnorm(z_stat),
                    greater = 1 - pnorm(z_stat))
    
    # 计算置信区间
    ci <- mean_x + c(-1, 1) * qnorm((1 + conf.level)/2) * se
    
    # 输出结果
    cat("单样本z检验结果\n")
    cat(sprintf("数据: x\nz = %.4f, p值 = %.4e\n", z_stat, p_val))
    cat(sprintf("备择假设: 总体均值不等于0\n"))
    cat(sprintf("%.0f%%置信区间: [%.4f, %.4f]\n", conf.level*100, ci[1], ci[2]))
    cat(sprintf("样本均值: %.4f\n", mean_x))
  } else {
    # 双样本检验逻辑
    n_x <- length(x)
    n_y <- length(y)
    mean_x <- mean(x)
    mean_y <- mean(y)
    diff_mean <- mean_x - mean_y
    se <- sqrt((sigma.x^2)/n_x + (sigma.y^2)/n_y)
    z_stat <- diff_mean / se
    
    # 计算p值
    p_val <- switch(alternative,
                    two.sided = 2 * pnorm(-abs(z_stat)),
                    less = pnorm(z_stat),
                    greater = 1 - pnorm(z_stat))
    
    # 计算置信区间
    ci <- diff_mean + c(-1, 1) * qnorm((1 + conf.level)/2) * se
    
    # 输出结果
    cat("双样本z检验结果\n")
    cat(sprintf("数据: x 和 y\nz = %.4f, p值 = %.4e\n", z_stat, p_val))
    cat(sprintf("备择假设: 总体均值差不等于0\n"))
    cat(sprintf("%.0f%%置信区间: [%.4f, %.4f]\n", conf.level*100, ci[1], ci[2]))
    cat(sprintf("样本估计值:\nmean of x = %.4f, mean of y = %.4f\n", mean_x, mean_y))
  }
}

使用该函数处理小样本场景：

s <- 0.2
x <- 1
y <- 2

my_z_test(x, y, sigma.x = s, sigma.y = s)

输出结果与手动计算一致，且支持单样本、不同备择假设等场景。

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