为何自定义R置信区间函数与t.test()计算结果不一致？

问题背景 #

我编写了一个计算95%置信区间的R函数，但运行结果与base包中的t.test()输出不一致：

测试数据 #

testData <- structure(list(group = c("Group1", "Group1", "Group1", "Group1", 
                                    "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", 
                                    "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", 
                                    "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", 
                                    "Group1", "Group1", "Group1", "Group1", "Group1", "Group1", "Group1"
), year = c(2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
            2015, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 
            2016, 2016, 2016, 2016, 2016, 2017, 2017, 2017, 2017, 2017, 2017, 
            2017), category = c("cat1", "cat1", "cat1", "cat1", "cat1", "cat1", 
                                "cat1", "cat1", "cat1", "cat1", "cat2", "cat2", "cat2", "cat2", 
                                "cat2", "cat2", "cat2", "cat2", "cat2", "cat2", "cat2", "cat2", 
                                "cat2", "cat2", "cat2", "cat3", "cat3", "cat3", "cat3", "cat3", 
                                "cat3", "cat3"), value = c(15.1382663715558, 38.7804544564934, 
                                                           46.8153764828161, 167.414619767484, 147.819242182614, 163.289605383038, 
                                                           97.4909154781249, 76.4990140823147, 10.2998099118541, 106.829837472452, 
                                                           47.9470225625797, 117.481510505374, 103.353651531038, 82.8258992025231, 
                                                           75.8617413682001, 0.895652854035013, 158.506322595117, 153.09256856583, 
                                                           223.536384788365, 75.748851191101, 46.9191391269587, 1.05445490408603, 
                                                           34.2440937279552, 12.5493519758163, 81.9894639436096, 102.38603104988, 
                                                           11.8608226647822, 16.0662436435422, 0.883484884196097, 58.1467542647205, 
                                                           145.495946136843, 106.259860732627)), row.names = c(NA, -32L), class = c("data.table", 
                                                                                                                                    "data.frame"))

自定义函数 #

calculateCI <- function(value){
  
  avg <- mean(value)
  s <- sqrt(var(value))
  n <- length(value)
  
  error <- qnorm(0.975)*s/sqrt(n)
  
  lower <- avg - error
  upper <- avg + error 
  
  return(list(lowerCI = lower, 
              upperCI = upper))
  
}

运行结果 #

• 自定义函数输出：

  $lowerCI
  [1] 58.49955
  
  $upperCI
  [1] 99.4681
  

• t.test(testData$value)输出：

  One Sample t-test
  
  data:  testData$value
  t = 7.5573, df = 31, p-value = 1.615e-08
  alternative hypothesis: true mean is not equal to 0
  95 percent confidence interval:
    57.66815 100.29950
  sample estimates:
  mean of x 
    78.98382 
  

差异原因 #

核心差异在于分位数的选择：

• 自定义函数使用了正态分布的分位数qnorm(0.975)（对应值约为1.96），这仅适用于总体标准差已知、或样本量极大（n>30）的场景。
• t.test()默认使用t分布的分位数qt(0.975, df = n-1)，这里样本量n=32，自由度df=31，对应的分位数约为2.0395。

t分布的尾部比正态分布更厚，小样本场景下t分位数更大，计算出的误差区间更宽，因此t.test()的置信区间比自定义函数的结果更宽，这就是二者结果不一致的原因。

修正后的函数 #

将qnorm替换为qt，并指定自由度为n-1，即可得到与t.test()一致的结果：

calculateCI <- function(value){
  
  avg <- mean(value)
  s <- sqrt(var(value))
  n <- length(value)
  df <- n - 1
  
  error <- qt(0.975, df = df)*s/sqrt(n)
  
  lower <- avg - error
  upper <- avg + error 
  
  return(list(lowerCI = lower, 
              upperCI = upper))
  
}

运行修正后的函数，输出结果将与t.test()的置信区间完全一致：

$lowerCI
[1] 57.66815

$upperCI
[1] 100.2995

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