boot.ci() возвращает странный доверительный интервал для примеров начальной загрузки

У меня есть две группы положительных вещественных чисел.

> dput(group1)
c(2.10753, 2.57251, 2.61687, 4.62551, 7.13166, 6347.63, 4.22139, 
10.7373, 2.11568, 2.71866, 4.09376, 10.9046, 109807, 5.87156, 
3.17082, 3.4703, 2.47262, 9.24319, 34.6945, 5.72567, 12.0134, 
108.33, 6.60707, 6.24304, 3.59048, 10.3174, 48.0265, 5.32097, 
3.77157, 6.67401, 22.633, 34.8186, 21.5315, 9.42882, 7.10627, ...)

> dput(group2)
c(4.88474, 65.4318, 128.101, 24.1271, 5.44262, 54.8987, 2.85175, 
14.1089, 172.23, 66.8563, 6.74067, 2.19603, 2.12985, 4.12735,
16.401, 3.22688, 15.6943, 4.32861, 36.4752, 7.33769, 75.855, 
62.7653, 35.1786, 3.71099, 29.0186, 34.4472, 19.1061, 2.75174, ...)

Группа 1 состоит из ~1000 значений, группа 2 из ~30000. Меня интересовало соотношение медиан между двумя группами, и я использовал следующую R-функцию для вычисления этого отношения для каждого из 2000 образцов начальной загрузки (для команды boot() см. Вывод функции ниже):

medianRatio <- function(x, i, noGroup1, noGroup2) {
    all <- x[i]
    currGroup1 <- all[1:noGroup1]
    currGroup2 <- all[c(noGroup1 + 1):length(all)]
    ratio <- median(currGroup1) / median(currGroup2)
    return(ratio)
}

Вызов из функции boot() выглядит следующим образом

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = c(group1, group2), statistic = medianRatio, R = noBs, 
    noGroup1 = length(group1), noGroup2 = length(group2))


Bootstrap Statistics :
    original      bias    std. error
t1*  1.08847 -0.08597889  0.05451763`

Среднее значение итогового распределения выборок начальной загрузки составляет 1,002, среднеквадратичное отклонение равно 0,054 (визуальный осмотр гистограммы подтверждает нормальное распределение около 1). Также:

диапазон (Group1_Group2.BS$t) [1] 0,823311 1,198469

Тем не менее, когда я запускаю boot.ci() на загрузочном объекте, сообщаемый доверительный интервал

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 2000 bootstrap replicates

CALL : 
boot.ci(boot.out = Group1_Group2.BS, type = "norm")

Intervals : 
Level      Normal        
95%   ( 1.068,  1.281 )  
Calculations and Intervals on Original Scale

Я не понимаю, что здесь происходит, так как сообщенный доверительный интервал даже не охватывает среднее (симметричное) распределение образцов начальной загрузки. Что мне не хватает?

2 ответа

Решение

Понимание того, как рассчитываются нормальные доверительные интервалы, может помочь прояснить ситуацию. Я нашел очень хорошее объяснение в ответе на этот вопрос.

Нормальные доверительные интервалы начальной загрузки построены вокруг наблюдаемого значения статистики с поправкой на смещение для устранения разницы между серединой статистики начальной загрузки и наблюдаемым значением. CI рассчитываются с использованием оценки статистических данных, представляющих интерес, из исходной выборки, вычитая смещение и прибавляя в 1,96 раза стандартную ошибку начальной загрузки. Если вы делаете это "вручную" со значениями из boot объект, вы увидите, что вы получите те же значения, что и от boot.ci,

1.08847 - -0.08597889 - 1.96*0.05451763
[1] 1.067594
1.08847 - -0.08597889 + 1.96*0.05451763
[1] 1.281303

Ваше смещение достаточно велико, чтобы сделать разницу между, скажем, процентильным КИ и нормальным КИ. Когда я думаю о доверительных интервалах начальной загрузки, я ожидал бы использовать среднее распределение начальной загрузки минус смещение, а не наблюдаемая статистика минус смещение. Я не потратил достаточно времени на размышления об этом, и при этом у меня нет моих начальных заметок, но вы можете проверить эту проблему немного более внимательно.

Вы делаете бутстрап неправильно. В вашей функции medianRatio вы предполагаете, что первая часть x[i] будет содержать значения из группы 1, а вторая часть из группы 2. Это неверно, потому что я просто образец с заменой от 1 до noGroup1+noGroup2. Вам нужно использовать параметр strata в команде загрузки. Я думаю, что следующее будет работать. Опция stype="f" указывает, что при загрузке будет сгенерирован вектор длины noGroup1 + noGroup2, который сообщает, сколько раз (0,1,2,...) каждое из наблюдений было выбрано для образца начальной загрузки. Этот код основан на примере в загрузочной документации.

medRatio <- function(x, f, noGroup1){
 gp1 <- 1:noGroup1
 m1 <- median(rep(x[gp1],f[gp1]))
 m2 <- median(rep(x[-gp1],f[-gp1]))
 return(m1/m2)
}
n1 <- length(group1)
n2 <- length(group2)
boot(c(group1,group2),medRatio,R=2000,stype = "f",strata = rep(1:2,c(n1,n2)),noGroup1=n1)
Другие вопросы по тегам