Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

 
Добавить ответ в эту темуОткрыть тему
> Доверительные интервалы полиномиального распределения
nokh
сообщение 2.04.2020 - 09:24
Сообщение #1





Группа: Пользователи
Сообщений: 1202
Регистрация: 13.01.2008
Из: Челябинск
Пользователь №: 4704



Раньше считал (вероятно не совсем корректно) ДИ для долей всегда методами для биномиального распределения. Т.е., например, в ряду абсолютных частот 4-х категорий {13, 35, 18, 7} с суммой n=73 доля первой категории f=13/73*100%=17,8%. Для неё находил 95% ДИ методом Клоппера - Пирсона или в полседнее время методом Джеффриса (байесовский априорный интервал): [10,4; 27,7].
Сейчас решил посчитать ДИ для полиномиального распределения, думал, что раз информации больше, то они Уже будут. Ничего подобного. R-пакет DescTool считает одновременные ДИ для полиномиалного распределения функцией MultinomCI.
library(DescTools)
x<-c(13,35,18,7)
MultinomCI(x)
est lwr.ci upr.ci
[1,] 0.17808219 0.06849315 0.3006248
[2,] 0.47945205 0.36986301 0.6019947
[3,] 0.24657534 0.13698630 0.3691180
[4,] 0.09589041 0.00000000 0.2184330
По умолчанию считает ДИ методом Сайсона - Глаза по SAS-овскому алгоритму. Всё хуже, чем даже биномиальный Клоппер - Писон, который ругают за консервативность. Видно, что для 7 (9,6%) нижняя граница вообще ноль. Более адекватные результаты даёт только метод Уилсона:
> MultinomCI(x, method="wilson")
est lwr.ci upr.ci
[1,] 0.17808219 0.10713373 0.2812173
[2,] 0.47945205 0.36877454 0.5921840
[3,] 0.24657534 0.16204465 0.3564445
[4,] 0.09589041 0.04722895 0.1849564

Воросы:
1) Каким способом считаете вы?
2) Хочу попробовать сделать бутстреп. Думаю так: многократно пробублировать набор 4 типов в соотношении 13 : 35 : 18 : 7 и извлекать из него с возвратом случайные выборки размером n=73; для каждогго типа потом рассчитать ДИ методом процентилей. Корректно так будет организовать?

Сообщение отредактировал nokh - 2.04.2020 - 09:25
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 2.04.2020 - 09:49
Сообщение #2





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(nokh @ 2.04.2020 - 09:24) *
Раньше считал (вероятно не совсем корректно) ДИ для долей всегда методами для биномиального распределения. Т.е., например, в ряду абсолютных частот 4-х категорий {13, 35, 18, 7} с суммой n=73 доля первой категории f=13/73*100%=17,8%. Для неё находил 95% ДИ методом Клоппера - Пирсона или в полседнее время методом Джеффриса (байесовский априорный интервал): [10,4; 27,7].
Сейчас решил посчитать ДИ для полиномиального распределения, думал, что раз информации больше, то они Уже будут. Ничего подобного. R-пакет DescTool считает одновременные ДИ для полиномиалного распределения функцией MultinomCI.
library(DescTools)
x<-c(13,35,18,7)
MultinomCI(x)
est lwr.ci upr.ci
[1,] 0.17808219 0.06849315 0.3006248
[2,] 0.47945205 0.36986301 0.6019947
[3,] 0.24657534 0.13698630 0.3691180
[4,] 0.09589041 0.00000000 0.2184330
По умолчанию считает ДИ методом Сайсона - Глаза по SAS-овскому алгоритму. Всё хуже, чем даже биномиальный Клоппер - Писон, который ругают за консервативность. Видно, что для 7 (9,6%) нижняя граница вообще ноль. Более адекватные результаты даёт только метод Уилсона:
> MultinomCI(x, method="wilson")
est lwr.ci upr.ci
[1,] 0.17808219 0.10713373 0.2812173
[2,] 0.47945205 0.36877454 0.5921840
[3,] 0.24657534 0.16204465 0.3564445
[4,] 0.09589041 0.04722895 0.1849564

Воросы:
1) Каким способом считаете вы?
2) Хочу попробовать сделать бутстреп. Думаю так: многократно пробублировать набор 4 типов в соотношении 13 : 35 : 18 : 7 и извлекать из него с возвратом случайные выборки размером n=73; для каждогго типа потом рассчитать ДИ методом процентилей. Корректно так будет организовать?


Для долей надо "восстановить" выборку в виде объектов с признаками замеренными.

После этого или 1) перемешиваем признаки (рандомизация) и получаем "доверительный интервал для 0-гипотезы".

Или 2) извлекаем "с возвращением" перевыборки и считаем в каждой нужные нам доли, следим за "сходимостью процентилей долей в накопленной выборке", как только достигнутая точность нас устраивает прекращаем.

Код
> d <- c(rep(1, 13), rep(2, 35), rep(3, 18), rep(4, 7))
> table(d)
d
1  2  3  4
13 35 18  7
> length(d)
[1] 73
> dd <- replicate(100000, table(sample(factor(d), replace=T)))
> str(dd)
int [1:4, 1:100000] 12 36 15 10 12 37 16 8 17 40 ...
- attr(*, "dimnames")=List of 2
  ..$ :8322456 [1:4] "1" "2" "3" "4"
  ..$ : NULL
> sapply(1:nrow(dd), function(i) quantile(dd[i,], probs=c(0.025,0.5,0.975)))
      [,1] [,2] [,3] [,4]
2.5%     7   27   11    3
50%     13   35   18    7
97.5%   20   43   25   12



PS

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

Сообщение отредактировал p2004r - 2.04.2020 - 10:06


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
100$
сообщение 2.04.2020 - 18:35
Сообщение #3





Группа: Пользователи
Сообщений: 902
Регистрация: 23.08.2010
Пользователь №: 22694



Цитата
Каким способом считаете вы?


1. Wilson's CI, которые для данного случая якобы дают более адекватный результат, представляют собой CI для каждого из параметров полиномиального распределения как для биномиального, невзирая на то, что элементарных исходов у эксперимента уже не 2, а 4. Поэтому Уилсоном мы не пользуемся.

2. В оригинальной статье Сайсон* - Глаза (1995) первый из обсуждаемых методов построения доверительного множества для многомерного параметра основан на существовании и оценке такого неотрицательного целочисленного гиперпараметра C, для которого одновременно справедлива система уравнений (6). Поэтому MultinomCI() вовсе не обязаны напоминать о таковых для биномиальных параметров.

Резюме: считайте MultinomCI() с установками по умолчанию.

2.
Цитата
Хочу попробовать ... Корректно так будет организовать?


Поскольку мультиномиальное распределение асимметрично, то Эфронов ДИ - это худшее, что можно измыслить в данной ситуации. Патамушта:
а) смещение исходной выборки все равно не устраняет;
б) непонятно, какие процентили вообще искать: то ли 2.5 - 97.5, то ли 3 - 98, то ли 2 - 97.
_____________________________

* Кристина Сайсон - женщина.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
nokh
сообщение 4.04.2020 - 05:52
Сообщение #4





Группа: Пользователи
Сообщений: 1202
Регистрация: 13.01.2008
Из: Челябинск
Пользователь №: 4704



Благодарю за мнения и код! Попробую всё-таки ещё свой вариант, интересно будет сравнить с результатом р2004r.
По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию. Последнее для меня очень затратно по времени написания кодов, но может за самоизоляцию и получится (как ни странно, сейчас времени вообще нет: в НИИ дана команда сидеть дома и писать статьи на год вперёд))), а в универе народ у кого занятий много вообще вешается с этой дистанционкой...)

Сообщение отредактировал nokh - 4.04.2020 - 05:54
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
nokh
сообщение 4.04.2020 - 05:53
Сообщение #5





Группа: Пользователи
Сообщений: 1202
Регистрация: 13.01.2008
Из: Челябинск
Пользователь №: 4704



Не нуда запастил, сорри... Стёр.

Сообщение отредактировал nokh - 4.04.2020 - 05:55
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 4.04.2020 - 08:15
Сообщение #6





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(nokh @ 4.04.2020 - 05:52) *
Благодарю за мнения и код! Попробую всё-таки ещё свой вариант, интересно будет сравнить с результатом р2004r.
По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию. Последнее для меня очень затратно по времени написания кодов, но может за самоизоляцию и получится (как ни странно, сейчас времени вообще нет: в НИИ дана команда сидеть дома и писать статьи на год вперёд))), а в универе народ у кого занятий много вообще вешается с этой дистанционкой...)


Возьмите готовые доверительные интервалы отсюда https://cran.r-project.org/web/packages/Ternary/index.html


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
100$
сообщение 4.04.2020 - 19:25
Сообщение #7





Группа: Пользователи
Сообщений: 902
Регистрация: 23.08.2010
Пользователь №: 22694



Цитата(nokh @ 4.04.2020 - 05:52) *
...По поводу Сайсон - Глаза ничего не читал, но мне решительно не понравился ноль в качестве нижней границы. Получается так: по набору в 73 объекта частота почти 9,5%, а нижняя граница ноль. Причём не 0.0001, что и так нереалистично мало, а вообще 0.00000000. Т.е. по-сути, метод говорит, что несмотря на то, что в выборке у меня оказалось почти 10%, если я продолжу процесс извлечения выборок, то в 95% выборок не обнуружу ни одного объёта такой категории. Не верю. Поэтому более склонен довериться моделированию.


Все это - следствие гримас и ухмылок смещения малой выборки.

Удвоим вектор частот:

> x<-2*x

> x
[1] 26 70 36 14

> MultinomCI(x)
est lwr.ci upr.ci
[1,] 0.17808219 0.09589041 0.2622103
[2,] 0.47945205 0.39726027 0.5635801
[3,] 0.24657534 0.16438356 0.3307034
[4,] 0.09589041 0.01369863 0.1800185
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

Добавить ответ в эту темуОткрыть тему