0

有些软件包允许您使用偏差校正和加速引导程序确定 95% 置信区间。我试图完全理解这个过程并编写 R 代码来重现 DescTools::MeanDiffCI 函数产生的相同结果。我已经创建了合成数据，并使用该函数通过 BCa 调整返回 95% CI。

但是，当我复制它并生成重新采样时，我陷入了如何应用校正的困境。我会很感激任何指示。

library(dplyr)

set.seed(54321)

N = 40
c1 <- rnorm(N, mean = 100, sd = 25)
c2 <- rnorm(N, mean = 100, sd = 50)
g1 <- rnorm(N, mean = 120, sd = 25)
g2 <- rnorm(N, mean = 80, sd = 50)
g3 <- rnorm(N, mean = 100, sd = 12)
g4 <- rnorm(N, mean = 100, sd = 50)
gender <- c(rep('Male', N/2), rep('Female', N/2))
dummy <- rep("Dummy", N)
id <- 1: N


wideData <- 
  tibble::tibble(
    Control1 = c1, Control2 = c2,
    Group1 = g1, Group2 = g2, Group3 = g3, Group4 = g4,
    Dummy = dummy,
    Gender = gender, ID = id)


myDat  <- 
  wideData %>%
  tidyr::gather(key = Group, value = Measurement, -ID, -Gender, -Dummy)

# Determine 95% CI using bca correction
library(DescTools)

Control1Data= 
  my.data %>%
  filter(Group== 'Control1') %>%
  select(Measurement)

ctlSample= as.vector(Control1Data$Measurement)

Group1Data= 
  my.data %>%
  filter(Group== 'Group1') %>%
  select(Measurement)
grpSample= as.vector(Group1Data$Measurement)

# Bootstrap
MeanDiffCIResults= MeanDiffCI(x= grpSample, y = ctlSample, method= 'bca', R= 5000)

### MANUAL RESAMPLE ###

ctlNo= length(ctlSample)
grpNo= length(grpSample)

# Create one large matrix since we are sampling WITH replacement
resampNo= 5000

ctlResampMatrix= matrix(sample(x= ctlSample, size= ctlNo * resampNo, replace = T), ncol= resampNo)
grpResampMatrix= matrix(sample(x= grpSample, size= grpNo * resampNo, replace = T), ncol= resampNo)

mean(grpMeansResamp - ctlMeansResamp)
quantile(grpMeansResamp - ctlMeansResamp, probs= c(0.05, 0.95))

Thee 函数产生的 CI 为 (7.724227, 31.010417)，而我的代码在应用 BCa 之前返回 (9.697589, 28.802789)。

4

0 回答 0