Доверительные интервалы в R: t.test(), confint() и prop.test()

Что на самом деле означает «95% уверенности»

Доверительный интервал превращает точечную оценку («выборочное среднее равно 5.61») в честный диапазон («истинное среднее правдоподобно лежит между 5.29 и 5.93»). Но фраза 95% уверенности — одна из самых неверно прочитываемых в статистике, так что разберёмся сразу.

95% описывают процедуру, а не интервал. Представьте, что вы прогоняете своё исследование снова и снова: каждый раз новая выборка, каждый раз новый интервал. Интервалы прыгали бы, и около 95% из них накрывали бы истинное значение; 5% промахивались бы. Ваш единственный реальный интервал — одна реализация этого процесса. Чего 95% не означают: «вероятность 95%, что истинное среднее находится между этими числами». Истинное среднее является фиксированным (неизвестным) числом — оно не бродит внутрь и наружу интервалов; варьируются как раз интервалы.

На практике дружелюбное прочтение вполне годится: интервал — это диапазон значений, совместимых с вашими данными. Просто знайте, какое утверждение вы вправе сделать, когда на вас надавят.

Доверительный интервал для среднего, простой путь

Каждый t.test() несёт доверительный интервал; критерий можно запустить хотя бы ради того, чтобы его собрать:

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

result <- t.test(times)
print(result$conf.int)
print(result$estimate)

Среднее 5.61, 95-процентный интервал примерно от 5.29 до 5.93. Интервал делает то, чего не может голое среднее: он показывает, сколько точности реально покупают десять наблюдений.

Тот же интервал вручную

Формулу за этим интервалом стоит один раз собрать самому, потому что она демистифицирует любой доверительный интервал, который вам встретится: оценка ± критическое значение × стандартная ошибка.

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

n <- length(times)
m <- mean(times)
se <- sd(times) / sqrt(n)          # standard error of the mean
t_crit <- qt(0.975, df = n - 1)    # t critical value, 2.5% in each tail

margin <- t_crit * se
print(t_crit)
print(margin)
print(c(m - margin, m + margin))

Сверьте с выводом t.test() — тот же интервал до последней цифры. Три подвижные части:

• Стандартная ошибка sd(x)/sqrt(n) измеряет, насколько колеблется выборочное среднее (различие sd и SE разобрано в описательных статистиках).
• qt(0.975, df) — критическое значение t: для 95-процентного покрытия вы оставляете по 2.5% в каждом хвосте, отсюда 0.975. При df = 9 оно составляет около 2.26 — толще нормального значения 1.96, налог за малую выборку, взимаемый за оценку sd по тем же данным.
• Предельная ошибка — их произведение, то самое «±», которое выносят в заголовки.

Смена уровня: 90%, 99% и компромисс

conf.level управляет покрытием, а вместе с ним и шириной:

times <- c(5.1, 4.9, 6.2, 5.8, 5.5, 5.9, 6.1, 5.2, 5.4, 6.0)

print(t.test(times, conf.level = 0.90)$conf.int)
print(t.test(times, conf.level = 0.95)$conf.int)
print(t.test(times, conf.level = 0.99)$conf.int)

Больше уверенности стоит большей ширины: 99-процентный интервал должен быть достаточно широким, чтобы оказываться прав 99 раз из 100, поэтому он растягивается; 90-процентный уже, но промахивается вдвое чаще, чем 95-процентный. Бесплатного обеда нет, есть только регулятор. 95% — чистое соглашение: умолчание, которое стоит сохранять, если нет причины его менять, а не закон природы.

Доверительный интервал для доли: prop.test()

Допустим, 47 из 120 опрошенных пользователей осваивают новую функцию. Каков правдоподобный диапазон истинной доли освоения?

result <- prop.test(47, 120)
print(result$estimate)
print(result$conf.int)

Выборочная доля 0.39, 95-процентный интервал примерно от 0.30 до 0.49 — так что «около 40% освоения» честно только с этим ореолом в ±9 пунктов. prop.test() использует лучшее приближение, чем учебная формула p ± 1.96 × sqrt(p(1−p)/n) (это интервал типа Уилсона с поправкой на непрерывность), и это важнее всего вблизи 0 или 1: учебный интервал может вылезти за пределы [0, 1], а этот не может. Для очень малых количеств binom.test(47, 120)$conf.int даёт точную версию.

Интервалы для коэффициентов модели: confint()

Подогнанные модели получают интервалы через одну обобщённую функцию:

fit <- lm(mpg ~ wt, data = mtcars)
print(confint(fit))
print(confint(fit, level = 0.99))

Каждая строка ограничивает один коэффициент линейной регрессии: 95-процентный интервал наклона при wt идёт примерно от −6.5 до −4.2 мили на галлон на тысячу фунтов. Это информативнее p-значения из сводки: он говорит, что эффект не просто ненулевой, а составляет минимум около 4 миль на галлон и, возможно, около 6.5.

Для логистической glm() функция confint(fit) работает так же, но возвращает границы в логарифмах шансов; проэкспоненцируйте, чтобы получить интервалы отношений шансов — exp(confint(fit)), — и помните, что опорным значением «нет эффекта» становится 1 вместо 0.

Размер выборки: закон квадратного корня

Ширина уменьшается как sqrt(n), и у этого есть запоминающееся следствие: вчетверо больше данных — вдвое уже интервал. Посмотрите на это на смоделированных данных: то же распределение, одна выборка вчетверо больше (зерно делает результат воспроизводимым; см. случайные числа):

set.seed(42)
small <- rnorm(50, mean = 100, sd = 15)
large <- rnorm(200, mean = 100, sd = 15)   # 4 times the n

ci_small <- t.test(small)$conf.int
ci_large <- t.test(large)$conf.int

print(diff(ci_small))   # width with n = 50
print(diff(ci_large))   # width with n = 200 - about half

Ширины укладываются в отношение, близкое к 2:1 (выборочный шум не даёт ему быть точным). Из-за квадратного корня точность дорожает: первые 100 наблюдений сужают интервал сильнее, чем следующие 300 вместе взятые, а сокращение вдвое уже узкого интервала всегда стоит вчетверо больше того, что вы заплатили до сих пор.

Интервал и p-значение: два взгляда на один критерий

Доверительный интервал и проверка гипотезы — это одна и та же информация в разных одеждах. 95-процентный интервал содержит ровно те значения параметра, которые двусторонний критерий на уровне 0.05 не смог бы отвергнуть. Отсюда:

• Интервал для разности средних не содержит 0 ⇔ t-критерий даёт p < 0.05.
• Интервал для отношения шансов не содержит 1 ⇔ p-значение коэффициента ниже 0.05.

Когда у вас есть интервал, обычно у вас есть и лучшая сводка: он выносит тот же вердикт о значимости плюс даёт величину эффекта в реальных единицах. «p = 0.03» говорит, что различие существует; «95% ДИ: от 0.2 до 7.6» говорит, что оно существует и может быть как ничтожным, так и огромным, — а это часто и есть тот вывод, который имеет значение.

Что вы уносите с собой

• 95% описывают долгосрочную долю попаданий процедуры, а не вероятность того, что именно этот интервал поймал истину.
• Среднее: t.test(x)$conf.int или вручную как среднее ± qt(0.975, n−1) × SE.
• Доля: prop.test(x, n)$conf.int; коэффициенты модели: confint(fit) (для отношений шансов в glm экспоненцируйте).
• Выше уверенность = шире интервал; 95% это соглашение, а не закон.
• Ширина уменьшается как sqrt(n): вчетверо больше данных сужают интервал вдвое.
• Интервал содержит каждое значение, которое критерий уровня 0.05 не отверг бы, — и, в отличие от p-значения, показывает величину эффекта.

Дальше: инструментарий симуляции, стоящий за всем этим, — случайные числа через rnorm, runif, sample и set.seed.