Предположим, у вас есть популяция размером N = 10000, в которой за один день произошло 6500 событий.
pop <- rep(c(0,1), times=c(3500, 6500))
table(pop)
#pop
# 0 1
#3500 6500
Теперь предположим, что вы можете выбрать 100 из этих (0,1) событий без замены .
set.seed(123)
N <- 10000
n <- 100
sam <- data.frame(id=1:n, event=sample(pop, size=n), prob=n/N)
table(sam$event)
# 0 1
#30 70
Таким образом, мы получили 70 из 100. Максимальная оценка вероятности общего числа событий в популяции составляет просто 70/100 x 10 000 = 7 000. Стандартная ошибка этой оценки определяется как
sqrt((N-n)/N * N^2 * var(sam$event)/n)
#[1] 473.71
. 95-процентный доверительный интервал составляет [6101 - 7898], который охватывает истинное население в общей сложности 6500 человек. Но 1 из 20 дней вы, скорее всего, получите плохой образец.
R пакетов? Не очень нужно для этого эксперимента. Для более сложных схем выборки я могу думать только о пакете survey , но могут быть и другие.
Теперь, если вы делали это неоднократно, скажем, в течение 10 дней, вы ' получаю оценку за каждый день. Оценка общего количества, согласно статистике, занимающейся частыми исследованиями, будет представлять собой сумму x N / n и оценку SE, рассчитанную аналогичным образом. Например, предположим, что вы обнаружили 3, 4, 5, 11, 6, 8, 14, 8, 17 и 19 «положительных» событий из выборок 100 в течение десяти дней подряд:
events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)
Что означает «Отрицательное» или не происходящее событие:
events0 <- 100 - events1
Вектор (0,1) событий можно построить следующим образом, используя rep.
events <- rep(rep(c(0,1), each=10), times=c(events0, events1))
Давайте определим n и N как число событий в вашей десятидневной выборке и в десятидневной популяции соответственно.
n <- 100 * 10
N <- 10000 * 10
Количество "положительных" событий в вашей десятидневной выборке равно:
sum(events1)
#[1] 95
И MLE в десятидневной популяции:
(T <- sum(events1) * N / n)
[1] 9500
Стандартная ошибка этой десятидневной оценки:
SE <- sqrt((N-n)/N * N^2 * var(events)/n); SE
[1] 923.0409
При 95% ДИ :
T + c(-1,1) * 1.96*SE
[1] 7690.84 11309.16
Байесовец может сказать, что каждый день должен быть переоценен или обновлен на основе оценки предыдущего дня, но я думаю, что результаты будут довольно схожими.
Байесовец использовал бы правило Байеса и использовал бы Униформу (0,1) в качестве разумного до распределения для пропорции "положительных" событий для t он десятидневный период. Unif (0,1) совпадает с Beta (1,1). Опытный статистик (частый или байесовский) признал бы, что бета-распределение сопряжено с биномиальным распределением. Таким образом, байесовский алгоритм будет использовать распределение бета (1 + X, 1 + NX) для доли «положительных» событий за десятидневный период, где X - общее количество «положительных» событий (95), а N - Размер выборки (1000). Обратите внимание, что среднее значение бета (альфа, бета) = альфа / (альфа + бета).
В R:
n <- rep(100, 10)
events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)
X <- sum(events1)
N <- sum(n)
pmean = (1+X)/(2+N); pmean
#[1] 0.09580838
CI = qbeta(c(.025,.975), 1+X, 1+N-X); CI # 95% credible interval
#[1] 0.07837295 0.11477134
Таким образом, в течение десятидневного периода доля положительные события будут составлять 9,58% от всех событий, и существует 95% вероятность того, что истинная пропорция находится между 7,84% и 11,48%. Экстраполируя на все население, мы можем сказать, что 9,58% из 100 000 событий или 9 581 событие были бы положительными. Это, как я уже сказал, очень похоже на частый подход.
Мета-анализ
Теперь эти два метода эффективно объединяют все десять дней в одну большую выборку и оценки доли положительных событий или общего количества положительных событий в общей популяции. Может оказаться более интуитивно понятным комбинировать результаты каждого дня более подходящим образом, основываясь на весах, как это делается в метаанализе.
Если p [k] является оценочной пропорцией в день k, а se [k] является ее стандартной ошибкой, то объединенная оценка дается как p_hat = sum (p [k] * w [k]) / sum (w [k]), где w [k] = (1 / se [k]) ^ 2, а стандартная ошибка равна 1 / sqrt (sum (w [k]).
In R:
N <- rep(100000, 10) # Population and 10 day period
n <- rep(100, 10)
events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)
events0 <- n - events1
p <- NULL; SE <- NULL; w <- NULL
for(k in seq_along(events1)){
events <- c(rep(0, events0[k]), rep(1, events1[k]))
p[k] <- sum(events1[k]) / n[k]
SE[k] <- sqrt((N[k]-n[k]) / N[k] * var(events)/n[k])
w[k] <- 1 / (SE[k])^2
}
p_hat <- sum(p*w)/sum(w); p_hat
#[1] 0.06997464
SE_p <- 1 / sqrt(sum(w)); SE_p
#[1] 0.007943816
(p_hat + c(-1,1) * 1.96 * SE_p)
#[1] 0.05440476 0.08554452
Таким образом, около 7% всех событий будут положительными с 95% доверительным интервалом (5,44% - 8,55%), который не слишком отличается от двух грубых методов, использованных выше. меньшая (и, возможно, более точная) оценка из-за искаженного характера десятидневной выборки.