28 апреля 2018

Интервальное оценивание параметров распределения



Автор: Владимир Шитиков

Два подхода к оценке доверительных интервалов 

Одна из главных задач статистической обработки заключается в оценке по выборочным данным таких параметров наблюдаемой случайной величины, как центр положения \(m\), интенсивность разброса \(s\) или форма кривой плотности распределения эмпирических частот. Для этого часто используется параметрический подход, который предполагает приближенную аппроксимацию эмпирического распределения наиболее подходящим теоретическим распределением (нормальным, логнормальным, биномиальным, гипергеометрическим или иным).




Если анализ природы исходных данных позволяет сделать допущение, что распределение случайной величины может быть описано нормальной моделью \(N(\mu, \sigma^2)\), то необходимо оценить только два параметра: математическое ожидание \(\theta_1 = \mu\) и дисперсию \(\theta_2 = \sigma^2\).  Однако при статистическом приближении не существует гарантированной точности и обычно задаются вопросом о несмещенности и надежности найденных оценок, для чего рассчитывают доверительные интервалы (ДИ) для параметров.
  • Параметрический доверительный интервал (подход А) определяется как границы диапазона, который с вероятностью \((1-\alpha)\times 100\text{%}\) накроет некоторый неизвестный параметр \(\theta\).
Чтобы найти границы ДИ,  подбирают такую  непрерывную центральную статистику \(G (\boldsymbol{X}, \theta)\) в рамках постулируемого теоретического распределения, которая строго монотонно связана с \(\theta\). Если известна функция плотности распределения этой статистики, то для любой доверительной вероятности \((1-\alpha)\) можно найти пару критических значений \(g_{H}\) и \(g_{B}\), которые далее позволят найти интервальную оценку  \(\theta\).

В случае предположения о нормальном распределении выборочных данных \(\boldsymbol{X}\) объемом \(n\) связь оцениваемых параметров с центральными статистиками основывается на теореме Р. Фишера (1925), которая утверждает, что случайные величины, определяющие эту зависимость, распределены как  \(F(\sqrt{n}(\overline{X} - \mu)/\sigma) = N(0, 1)\)  и \(F(nS^2/\sigma^2) = \chi^2 (n - 1)\)   соответственно,   где \(\overline{X}\) и S - выборочные среднее и стандартное отклонение (см. Ивченко, Медведев, 1984, с. 82-86). Тогда, используя в качестве критических значений тест-статистики квантили t-распределения Стьюдента \(t_{\alpha / 2, f}\), легко построить доверительный интервал для математического ожидания \(\mu\):
\[\overline{X} - t_{\alpha / 2, f} S/n^{0.5} < \mu < \overline{X} + t_{\alpha / 2, f} S/n^{0.5},\]
где \(\overline{X}\) и \(S\) - выборочные среднее и стандартное отклонение.

Поскольку для нормальных выборок несмещённая оценка среднеквадратичной вариации имеет распределение \(\chi^2\) с \((n-1)\) степенями свободы, то формула нахождения доверительных интервалов для дисперсии \(\sigma^2\) имеет вид:
\[\frac{(n - 1)S^2}{\chi^{2}_{1-\alpha / 2, n - 1}} \lt  \sigma^2 \lt  \frac{(n-1)S^2}{\chi^{2}_{\alpha / 2, n - 1}}\]
Прежде чем приступить к рассмотрению примеров, затронем проблему так называемых "малых выборок", требующую на наш взгляд особого внимания, поскольку такие выборки повсеместно встречаются на практике. Как известно, на заре человеческой цивилизации считали достаточно просто: "Один, два, три, много...". В дальнейшем континуум чисел "много" стал операционным полем для математической статистики. Однако в случае очень малых выборок (2-4 наблюдения) имеет место высокая статистическая неопределенность и фундаментальная теория не всегда может приводить нас к применимым для практики результатам (такое бывает со всеми естественно-научными законами, такими, как закон Гука, когда они сталкиваются с пограничными состояниями изучаемого объекта).

Рассмотрим пример из книги Н. Цейтлина (2007, с. 150), где  лаборант выполнил определение концентрации магния в трех пробах раствора. Предположим, что наблюдения распределены по нормальному закону и оценим параметрическими методами доверительные интервалы математического ожидания и стандартного отклонения.

# Функция, возвращающая выборочные оценки параметров
#  нормального распределения:
fitdist <- function(x) {
 fit <- c(mean(x), sd(x),length(x)) 
 names(fit) <- c("Mean", "Std.Dev", "N")
 fit
}

#  Выборочные данные:
x <- c(5.433, 5.402, 5.367)
(fit <- fitdist(x))
##      Mean     Std.Dev           N 
## 5.4006667   0.0330202   3.0000000

#  Доверительные интервалы математического ожидания:
alfa <- 0.05
dp <- 1-alfa/2
df <-(fit[3]-1)
t <- qt(dp, df=df)
CI.mp <- c(fit[1]-t*fit[2]/sqrt(fit[3]), fit[1]+t*fit[2]/sqrt(fit[3]))
names(CI.mp) <- c("ДИнижн", "ДИверх")
##   ДИнижн       ДИверх 
## 5.318640     5.482693

#  Доверительные интервалы стандартного отклонения:
CI.sdp <- c(sqrt((df*fit[2]^2)/qchisq(dp, df=df)),
            sqrt((df*fit[2]^2)/qchisq(1-dp, df=df)))
names(CI.sdp) <- c("ДИнижн", "ДИверх")
##     ДИнижн      ДИверх 
## 0.01719224  0.20752317


Разумеется, точность оценки ДИ параметрическим методом напрямую зависит от того, насколько выполняется предположение о монотонном характере  связи выборочных оценок параметра с назначенной тест-статистикой. От этого абстрагирован второй подход, оценивающий доверительные интервалы с использованием имитационного моделирования.
  • Имитируемый доверительный интервал (подход Б) определяется как границы диапазона, в пределах которого варьирует \((1-\alpha)\times 100\text{%}\) от общего числа оценок параметра, найденных по бесконечной (очень большой) последовательности случайных выборок, извлекаемых из анализируемой совокупности.
Формально говоря, в этом определении речь уже не идет об "истинном" (но неизвестном!) значении параметра \(\theta\), поэтому имитируемые интервальные оценки принято иногда называть толерантными интервалами.

Выполним нахождение доверительных интервалов параметрическим бутстрепом. В случае нормального распределения он состоит из нескольких вполне очевидных шагов:
  1. Находят оценки параметров mean(x) и sd(x) для исходной выборки x;
  2. Из нормального распределения с параметрами, вычисленными в п. 1, извлекают случайную выборку xr объемом n;
  3. По бутстреп-выборке из п. 2 находят новые, несколько отличающиеся от исходных параметры mean(xr) и sd(xr)
  4. Пп. 2-3 повторяют достаточно большое число раз (например, 1000), что в итоге позволяет построить бутстреп-распределение выборочных оценок анализируемого параметра.
#  Функция, возвращающая оценки параметров для случайной
#  выборки из нормального распределения с параметрами  fit:
myboot0 <- function(fit){
  # генерация случайной выборки из заданного распределения:
  xr <- rnorm(fit[3], mean = fit[1], sd = fit[2])
  # подгонка параметров распределения под новые данные:
  fitr <- fitdist(xr)
  return(fitr[1:2])
}

Построим бутстреп-распределение 1000 выборочных значений оцениваемых параметров. Вычислим по бутстреп-выборке величину среднего и стандартного отклонения, а также их смещение относительно эмпирических значений. Найдем имитируемые доверительные интервалы  с использованием t-распределения и методом процентилей (Шитиков, Розенберг, 2014):

set.seed(123)
nr = 1000
bootpar <- replicate(nr, myboot0(fit))
print(bootpar[, 1:5])
##               [,1]       [,2]       [,3]       [,4]       [,5]
## Mean    5.40912044 5.42174303 5.38425565 5.41319490 5.40017811
## Std.Dev 0.03764918 0.03080718 0.02900863 0.02757333 0.01619625

# Среднее, ошибка, смещение среднего
c(mean(bootpar[1, ]), sqrt(var(bootpar[1, ])), mean(bootpar[1, ]) - fit[1])
## 5.4010902122 0.0188399199 0.0004235456

#  Доверительные интервалы с использованием t-распределения
mean(bootpar[1, ]) + qt(dp, nr-1)*sqrt(var(bootpar[1, ]))*c(-1, 1)
## [1] 5.364120 5.438061

#  Доверительные интервалы методом процентилей
CI.mb <- quantile(bootpar[1, ], prob = c(0.025,0.5,0.975))
##     2.5%      50%    97.5% 
## 5.362222 5.401146 5.440371 

Покажем на графике эмпирическую функцию плотности распределения выборочных средних.  Здесь красными пунктирными линиями представлены доверительные границы, вычисленные параметрическим методом А, а синими – методом процентилей Б.

hist(bootpar[1,], freq = FALSE, breaks=8,col="grey88", 
       xlim = c(5.25, 5.55),
       main="Гистограмма и ядерная плотность оценки среднего")
lines(density(bootpar[1, ]), lwd = 2, col = "blue")
abline(v = CI.mb, col = "blue", lty = 2, lwd = 2)
abline(v = CI.mp, col = "red", lty = 2, lwd = 2)
abline(v = fit[1], col = "red", lty = 2, lwd = 2)

Аналогичным способом выполним оценку доверительных интервалов и построение графика для стандартного отклонения:

(CI.sdb <- quantile(bootpar[2,], prob = c(0.025,0.5,0.975)))
##        2.5%         50%       97.5% 
## 0.005573054 0.026825969 0.064515259 

hist(bootpar[2,], freq = FALSE, breaks = 8, col = "grey88", xlim = c(0, 0.21),
     main="Гистограмма и ядерная плотность оценки станд. отклонения")
lines(density(bootpar[2, ]), lwd = 2, col = "blue")
abline(v = CI.sdb, col = "blue", lty = 2, lwd = 2)
abline(v = CI.sdp, col = "red", lty = 2, lwd = 2)
abline(v = fit[2], col = "red", lty = 2, lwd = 2)



Очевидно, что в условиях столь малой выборки доверительные интервалы, вычисленные параметрическим и имитационным методами, сильно различаются, поскольку основаны на двух разных исходных концепциях.  Если для параметрического интервала ставится условие, что максимально возможное значение неизвестного параметра с заданной вероятностью не выйдет за пределы доверительной области, то имитационный интервал основывается на подсчете частот экстремальных оценок параметров.

Однако здесь мы не склонны столь безоговорочно становиться на сторону адептов бутстреп-метода, которые считают: "Если на одних и тех же данных параметрические и непараметрические методы приводят к различным выводам, то необходимо принять на веру результаты рандомизации и бутстрепа". В репрезентативной выборке малого объёма наблюдения могут случайно отстоять друг от друга как очень близко, так и очень далеко. Учитывая это, например, при объёме n = 2 правый квантиль распределения Стьюдента "перестраховывается" до значения \(t_{0.975, 1} = 12.7\), а величина \(1/\chi^2_{0.025, 1}\) достигает \(1018\)! Но можно лишь бездоказательно рассуждать, насколько практически целесообразной в случае малых выборок является такая перестраховка, обусловленная теорией, и являются ли параметрические ДИ оптимальными (т.е. обладающими свойствами минимальной ширины).

Другая проблема заключается в том, что классический бутстреп обычно предполагает нормальное распределение оценок имитируемого параметра, в то время, как теоретические предпосылки Фишера постулируют t-распределение их оценок для математического ожидания и \(\chi^2\) для дисперсии. В этом заключается еще одно принципиальное различие между подходами А и Б к нахождению доверительных интервалов (некоторые интересные подробности см. в сообщении).

Рассмотрим, как изменяется соотношение между оценками доверительных интервалов при увеличении объема выборки. Увеличим размер эмпирической совокупности до 50 элементов, включив в нее набор случайных значений, извлеченных из нормального распределения c параметрами \(N(5.4, 0.033)\). Найдем зависимость разности ДИ, вычисленных параметрическим и имитационным методом, от объема выборки для \(n\) в интервале от 3 до 50.  Пусть в анализе нас будет интересовать только ширина полуинтервала от среднего до верхнего граничного значения.

# Функция вычисляет ширину верхней части доверительного интервала
# среднего и СО параметрическим и имитационным методом для выборки xcur:
getEst <- function(xcur) {
  fit.cur <- fitdist(xcur) 
  df <-(fit.cur[3]-1)
  t <- qt(dp, df=df)
  CI.mp <- t*fit.cur[2]/sqrt(fit.cur[3])
  CI.sdp <- sqrt((df*fit.cur[2]^2)/qchisq(1-dp, df=df)) - fit.cur[2]
  bootpar <- replicate(1000, myboot0(fit.cur)) 
  CI.mb <- quantile(bootpar[1,], prob = dp) - fit.cur[1]
  CI.sdb <- quantile(bootpar[2,], prob = dp) - fit.cur[2]
  res <- c(n = fit.cur[3], CI.mp, CI.mb, CI.sdp, CI.sdb)
  names(res) <- c("n", "CI.mp", "CI.mb", "CI.sdp", "CI.sdb")
  res
}
set.seed(123) 
fit <- fitdist(x)

# Формируем выборку из 50 случайных значений
# из того же распределения, что и эмпирические данные:
x50 <- c(x, rnorm(47, mean = fit[1], sd = fit[2]))
np <- c(3, 4, 5, 7, 9, 12, 15, 20, 25, 30, 35, 40, 45, 50)
set.seed(123)
MasRes <- sapply(np, function (N) getEst(x50[1:N]))
MasResT <- as.data.frame(t(MasRes))

plot(MasResT$n, MasResT$CI.mp, type = "l", xlab = "Объем выборки", 
     ylab = "Верхная часть ДИ", col = 2, lwd = 2, ylim = c(0,0.1))
lines(MasResT$n, MasResT$CI.sdp, col = 2, lwd = 2, lty = 2) 
lines(MasResT$n, MasResT$CI.mb, col = 4, lwd = 2) 
lines(MasResT$n, MasResT$CI.sdb, col = 4, lwd = 2, lty = 2) 
legend("topright", col = c(2, 2, 4, 4), lty= c(1, 2, 1, 2), lwd = 2,
legend = c("Mean parametric", "SD parametric",  "Mean bootstrap", "SD bootstrap"))



Очевидно, что при увеличении объема эмпирической выборки уже до n = 15 доверительные интервалы математического ожидания становятся практически неразличимыми. Доверительные границы для стандартного отклонения также, как правило, сближаются по мере того, как (при \(n \gt 30\))  \(\chi^2\)-распределение асимптотически приближается к нормальному.

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


Доверительные полосы и доверительные эллипсы

Естественным является стремление рассматривать оценки доверительных интервалов не по отдельности для каждого параметра, а как результат совместного распределения двух случайных величин. Например, легко увидеть прямую аналогию между функцией плотности нормального распределения и обычной нелинейной функцией регрессии с двумя параметрами \(\theta_1\) и \(\theta_2\), которые с некоторой ошибкой \(\epsilon\) подгоняются по эмпирической выборке:

\[p = f(\theta_{1}, \theta_{2}) = \frac{1}{\sigma \sqrt{2\pi}} exp(\frac{-(x-\mu)^2}{2 \sigma^2}) +  \epsilon\]
Тогда, используя понятие доверительной полосы (confidence band) для всей линии регрессии (Себер, 1980), мы можем оценить доверительную область, которая с вероятностью \((1-\alpha)\times 100\text{%}\) накроет теоретическую аппроксимирующую функцию распределения.

Доверительную полосу для функции нормального распределения можно получить, если многократно извлекать из генеральной совокупности различные выборки одинакового размерома, оценивать для них параметры \(\mu\) и \(\sigma\) и строить по ним интегральные кривые распределения. Тогда \(100\times\alpha\text{%}\) таких линий окажутся за пределами области, где находится большинство (\((1-\alpha)\times 100\text{%}\)) линий, а кривая, ограничивающая эту область, будет называться огибающей (envelope).

# 1. Функция для нахождения p-квантилей случайной выборки из 
#  нормального распределения с параметрами  fit:
myboot <- function(fit, p) {
  # генерация случайной выборки из заданного распределения:
  xr <- rnorm(fit[3], mean = fit[1], sd = fit[2])
  # подгонка параметров распределения под новые данные:
    fitr <- fitdist(xr)
    hc5r <- qnorm(p, mean = fitr[1], sd = fitr[2])
  return(hc5r)
}

# 2. Функция, возвращающая значения вероятностей для случайной
#  выборки из нормального распределения с параметрами  fit:
myboot2 <- function(fit, newxs){
  # генерация случайной выборки из заданного распределения:
  xr <- rnorm(fit[3], mean = fit[1], sd = fit[2])
  # подгонка параметров распределения под новые данные:
  fitr <- fitdist(xr)
  # прогноз вероятностей под новые данные
  pyr <- pnorm(newxs, mean = fitr[1], sd = fitr[2])
  return(pyr)
}

x <- c(5.433, 5.402, 5.367)
fit <- fitdist(x)
set.seed(1234) # установка генератора случайных чисел
nrep = 200
D <- data.frame(val=sort(x), frac = ppoints(x, 0.5))

# Новые данные для построения плавной кривой:
newxs <- (seq(min(D$val)-1, max(D$val)+1, length.out = nrep))

# Получение матрицы для построения nrep кривых:
boots <- replicate(nrep, myboot2(fit, newxs))
require(reshape2)
bootdat <- data.frame(boots)
bootdat$newxs <- newxs
bootdat <- melt(bootdat, id = 'newxs')

# Извлечение доверительных интервалов:
cis <- apply(boots, 1, quantile, c(0.025, 0.975))
rownames(cis) <- c('lwr', 'upr')

# Добавление в итоговую таблицу подогнанных значений:
pdat <- data.frame(newxs, py = pnorm(newxs, mean = fit[1], sd = fit[2]))

# Добавление доверительных интервалов:
pdat <- cbind(pdat, t(cis))

# Вывод графика с использованием пакета  ggplot2:
library(ggplot2)
ggplot()+
  labs(x = 'Концентрация Mg', y = 'Накопленная вероятность')  + 
  xlim(c(5.1, 5.7))+ 
  geom_line(data = bootdat, aes(x = newxs, y = value, 
    group = variable), col = 'steelblue', alpha = 0.5) +
  geom_line(data = pdat, aes(x = newxs, y = py), 
    col = 'red', size=1.2) +
  geom_line(data = pdat, aes(x = newxs, y = lwr), 
    linetype = 'dashed', size = 1.2) +
  geom_line(data = pdat, aes(x = newxs, y = upr), 
    linetype = 'dashed', size = 1.2) +
  geom_point(data = D, aes(x = val, y = frac), col = 'green') +
  theme_bw()


Кривой красного цвета представлена теоретическая функция нормального распределения с эмпирическими значениями параметров \(N(5.4, 0.033)\), а голубым цветом - 200 кривых с параметрами, полученными бутстрепом. Пунктирными линиями представлены 95%-ные огибающие, которые позволяют оценить доверительные интервалы концентрации магния не только при \(p = 0.5\) (т.е. для математического ожидания), но и для любой другой вероятности.

Отметим также, что все три выборочные наблюдения хорошо ложатся на функцию нормального распределения (что всегда характерно для малых выборок), следовательно, с этой точки зрения в корректности параметрических ДИ сомневаться не приходится.


Другим вариантом представления совместного распределения оценок параметров нормального распределения является построение доверительных эллипсов, что часто применяется для  анализа связи между двумя коэффициентами модели множественной регрессии. В системе R имеется ряд методов построения доверительных эллипсов, один из которых реализован в пакете car:


set.seed(123)
bootpar <- replicate(1000, myboot0(fit))
library(car)
dataEllipse(bootpar[1,], bootpar[2,], xlab="Среднее",
    ylab = "Станд. отклонение", cex = .5, levels = c(.9, .95, .99), robust = TRUE,
    ellipse.label = c(.9, .95, .99), lty = 2, fill = TRUE, 
    fill.alpha = 0.1)

Сравнение групповых средних

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

Рассмотрим второй пример (Цейтлин, 2007, с. 150) основанный на наблюдениях за перепадом давлений (мм. рт. ст.), завихрителя вихревого пылеулавливающего аппарата при углах установки лопаток 15°, 25° и 35°. Делаем предположение, что все три выборки подчиняются нормальному закону распределения.

Рассчитаем вначале параметрические доверительные интервалы. Поскольку мы собираемся проверять гипотезу о равенстве групповых средних для всех возможных пар, то выполним коррекцию доверительной вероятности по формуле Бонферрони.

D15 <- c(192, 201, 207, 212)
D25 <- c(167, 171, 183, 194)
D35 <- c(142, 150, 158, 164)
dp <- 1 - (p.adjust(rep(0.05,3), "bonferroni")/2)[1]
## [1] 0.925

Для оценки параметров будем использовать функцию fitdist() из пакета fitdistrplus, в которой реализован метод максимизации функции правдоподобия:

library(fitdistrplus)
fit <- fitdist(D15, "norm")
Warning message:
In dnorm(c(192, 201, 207, 212), mean = 205.69609375,
sd = -2.30368341212433,  :   NaNs produced

summary((fit))  # Оптимальные параметры распределения
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##        estimate Std. Error
## mean 203.000000   3.724916
## sd     7.449832   2.633913
## Loglikelihood:  -13.70852   AIC:  31.41704   BIC:  30.18963 
## Correlation matrix:
##               mean            sd
## mean  1.000000e+00 -2.178503e-09
## sd   -2.178503e-09  1.000000e+00

Мы получили диагностическое сообщение и неверно вычисленное значение sd! Видимо, метод максимального правдоподобия некорректно работает на малых выборках (\(n \lt 5\)). Придется вернуться к предыдущему варианту функции fitdist(), рассчитывающей выборочные параметры с использованием функций mean() и sd():
fit.D15 <- fitdist(D15)
fit.D25 <- fitdist(D25)
fit.D35 <- fitdist(D35)
df.par <- t(data.frame(fit.D15, fit.D25, fit.D35))
df.par <- as.data.frame(df.par)
df.par$exp <- as.factor(c("D15", "D25", "D35"))
t.sd <- qt(dp, df = df.par$N - 1)*df.par$Std.Dev/sqrt(df.par$N) 
df.par$CI.mpl <- df.par$Mean - t.sd
df.par$CI.mph <- df.par$Mean + t.sd
df.par$CI.sdpl <- sqrt((df.par$N - 1)*df.par$Std.Dev^2)/
      qchisq(dp, df=(df.par$N - 1))
df.par$CI.sdph <- sqrt((df.par$N - 1)*df.par$Std.Dev^2)/
      qchisq(1 - dp, df=(df.par$N - 1))
print(df.par[, c(5,1,6,7,2,8)],2)
##           CI.mpl   Mean   CI.mph  CI.sdpl   Std.Dev  CI.sdph
## fit.D15 194.7232 203.00 211.2768 2.157919  8.602325 31.56758
## fit.D25 166.9824 178.75 190.5176 3.068039 12.230427 44.88147
## fit.D35 144.2880 153.50 162.7120 2.401735  9.574271 35.13429

Изобразим полученные интервалы на графике:

library(Hmisc)
errbar(x = c(1:3), y = df.par$Mean, yplus = df.par$CI.mph, 
       yminus = df.par$CI.mpl, ylab = "Перепад давления", 
       pch = 15, xaxt = "n", xlab = "", xlim = c(0.5, 3.5),
       main = "95% доверительные интервалы среднего")
axis(1, c(1:3), labels=substr(row.names(df.par),5,7))
errbar(x = c(1:3), y = df.par$Std.Dev, yplus = df.par$CI.sdph, 
       yminus = df.par$CI.sdpl, ylab = "Перепад давления", 
       pch = 15, xaxt = "n", xlab = "", xlim = c(0.5, 3.5),
       main = "95% доверительные интервалы станд. отклонения")
axis(1, c(1:3), labels = substr(row.names(df.par), 5, 7))


Рассчитаем теперь имитационные доверительные интервалы с использованием того же параметрического бутстрепа с 333 повторностями для каждой выборки:

set.seed(123)
boot.D15 <- replicate(333, myboot0(fit.D15))
boot.D25 <- replicate(333, myboot0(fit.D25))
boot.D35 <- replicate(333, myboot0(fit.D35))
df.boot <- data.frame(t(boot.D15), exp = rep("D15", 333))
df.boot <- rbind(df.boot, data.frame(t(boot.D25), exp = rep("D25", 333)))
df.boot <- rbind(df.boot, data.frame(t(boot.D35), exp = rep("D35", 333)))

Выведем полученные результаты на график и построим для каждой группы доверительные эллипсы:

library(ggplot2)
p2 <-qplot(data = df.boot, x = Mean, y = Std.Dev, colour = exp) +
           stat_ellipse(level = dp, type = "t", size = 1) +
           stat_ellipse(level = dp, type = "norm", linetype = 2, size = 1) +
           theme_bw()
 p2 + geom_point(data = df.par, aes(x = Mean, y = Std.Dev, shape = exp), 
                 colour = "black", size = 4.5)


Для каждой группы наблюдений на графике жирными символами представлено положение оценок параметров (среднего и стандартного отклонения). Пунктиром показаны доверительные эллипсы при аппроксимации двухмерным нормальным распределением, а сплошной линией - при аппроксимации многомерным распределением Стьюдента. Поскольку доверительные интервалы нигде не пересекаются, то очевидно, что следует отвергнуть нулевую гипотезу о статистической гомогенности групп.

Отметим, что проекции доверительных интервалов на ось абсцисс хорошо соответствуют значениям параметрических доверительных интервалов для групповых средних. В то же время, верхняя граница доверительного интервала стандартного отклонения, достигающая 44 мм рт. ст., может выглядеть на практике недостаточно реалистичной. По крайней мере, из 999 сгенерированных случайных выборок ни для одной из них оценка стандартного отклонения не превысила 25 мм.


Благодарности

Выражаю искреннюю признательность своему многолетнему знакомому, математику-статистику Н. А. Цейтлину (г. Гамбург) за плодотворное обсуждение и ценные замечания.


Список литературы
  • Ивченко Г. И., Медведев Ю. И. (1984) Математическая статистика. М.: Высш. школа, 840 с.
  • Себер Дж. (1980) Линейный регрессионный анализ. М.: Мир, 456 с.
  • Цейтлин Н. А. (2007) Из опыта аналитического статистика.  M.: Солар, 906 с.  PDF
  • Шитиков В. К., Розенберг Г. С. (2014) Рандомизация и бутстреп: статистический анализ в биологии и экологии с использованием R.  Тольятти: Кассандра, 314 с. (PDF, данные и скрипты доступны на сайте авторов)

2 комментария :

Анонимный комментирует...

Большое спасибо за статью и за ссылку на книгу Цейтлина!

Unknown комментирует...

Качественный контент"

Отправить комментарий