Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Доверительный интервал долей
Форум врачей-аспирантов > Разделы форума > Медицинская статистика
плав
Поскольку в другой теме было много споров о разных ДИ для долей и огрномный список цитат, решил - для иллюстрации - провести вычислительный эксперимент.
Итак. Были смоделированы популяции в которой содержится х объектов одного класса и 1-х объектов другого класса (х менялась от 1 до 10%). Из этой популяции брались случайные выборки размером 40 объектов. Оценивалось количество объектов одного и другого класса в выборке и рассчитывались доверительные интервалы по Клопперу-Пирсону, Агрести-Коулу и по распределению Пуассона (значения менее 0 заменялись на нулевые).
Таких выборок бралось 10 000 и затем рассчитывался вероятность покрытия популяционного значения доверительным интервалом и средняя ширина доверительного интервала. Что в результате (это небольшой размер выборки и малая вероятнсть:
Ширина 95%ДИ Покрытие
pi__ КП__ АК__ Пу__ КП__ АК__ Пу__
1% 0,104 0,118 0,110 0,993 0,993 0,993
2% 0,119 0,130 0,126 0,992 0,951 0,992
3% 0,132 0,141 0,141 0,994 0,970 0,994
4% 0,146 0,152 0,156 0,979 0,979 0,979
5% 0,157 0,162 0,169 0,986 0,952 0,986
6% 0,168 0,170 0,181 0,991 0,970 0,991
7% 0,177 0,178 0,192 0,981 0,981 0,981
8% 0,186 0,185 0,203 0,988 0,965 0,988
9% 0,195 0,192 0,214 0,953 0,974 0,970
10% 0,203 0,198 0,223 0,972 0,962 0,972
При малых значениях популяционной вероятности (менее 8%) интервал Клоппера-Пирсона является более узким, при больших - боле узкий интервал Агрести-Коула.
Покрытие прыгает, почти всюду больше номинального уровня 95%, однако среднее покрытие для КП - 98,3%, для Агрести - 97,0% и для Пуассона - 98,5%. Агрести-Коула немного ближе к номинальному уровню.
Однако если смотреть на данные реально, принципиальных различий между этими тремя методами нет. В большинстве случаев они дают одинаковые результаты и, как и показано в других работах, КП немного более консервативен, а АК немного более широкий при малых значениях популяционной вероятности.
плав
Теперь сделал при большем размере выборки - 150 объектов
Результат
1% 0,040 0,044 0,041 0,982 0,982 0,982
2% 0,052 0,054 0,053 0,989 0,968 0,989
3% 0,061 0,062 0,063 0,978 0,964 0,978
4% 0,069 0,069 0,071 0,968 0,963 0,968
5% 0,076 0,075 0,079 0,965 0,960 0,965
6% 0,082 0,080 0,086 0,965 0,948 0,965
7% 0,087 0,085 0,092 0,966 0,951 0,966
8% 0,093 0,090 0,098 0,969 0,954 0,969
9% 0,097 0,094 0,103 0,972 0,961 0,972
10% 0,102 0,098 0,108 0,971 0,960 0,971
Опять же при малых значениях популяционной вероятности интервал Клоппера-Пирсона уже, при больших (более 5%) - уже интервал Агрести-Коула. Среднее покрытие для Клоппера-Пирсона 97,2%, для Агрести - 96,1% и для пуассоновской аппроксимации - 97,2%. Опять интервал Агрести-Коула ближе к номианльному уровню, хотя его и не достигает.
В целом это означает, что при небольших значениях популяционной вероятности интервал покрытия соответствующий 95% для метода Клоппера-Пирсона получается при использовании 92% интервала.
Опять же серьезных различий между методами нет, но не один из них не является идеальным и даже нельзя выбрать однозначно предпочительный (ибо по выборке мы не знаем популяционную вероятность).
Для сравнения я сделал 10 000 выборок из нормальной популяции и оценил доверительный интервал по известной формуле M+/-t*m. Для 95% интервала вероятность накрытия при размере выборки 40 объектов 0,9536, для 150 объектов - 0,9476 - т.е. теоретические 95%. Более того, если выборка из прямоугольного распределения (явно не нормального), то интервал накрытия 0,9501 для размера выборки 40 и 0,9473 для размера выборки 150. Иными словами, опять интервал равен номинальному.
Если популяция с "тяжелыми хвостами" (смесь двух популяций с разными стандартными отклонениями), то для 40 объектов интервал покрытия 0,9511, для 150 - 0,9525. Все равно близко к номинальному уровню.
Собственно это и подтверждает идею о том, что разные методы близки для количественных переменных, а для качественных всегда остаются проблемы.
nokh
Исследуются частоты клеток с ядерными и цитоплазматическими нарушениями. Количество проанализированных клеток составляет 500, 1000 или чаще 2000. По некоторым показателям выборки крайне неоднородны (размах от промилле до процентов), что подтверждается G-критерием. Нужно найти доверительный интервал. Если объединить все данные внутри выборки и найти ДИ, скажем, по Клопперу-Пирсону, то это будет некорректно: выборка гетерогенная. Если найти ДИ бутстрепом индивидуальных частот - не будет учтена информация, что эти частоты оценены с разной точностью (по 500-2000 клеток). Как лучше всего поступить в этом случае?
плав
Цитата(nokh @ 24.06.2009 - 06:30) *
Исследуются частоты клеток с ядерными и цитоплазматическими нарушениями. Количество проанализированных клеток составляет 500, 1000 или чаще 2000. По некоторым показателям выборки крайне неоднородны (размах от промилле до процентов), что подтверждается G-критерием. Нужно найти доверительный интервал. Если объединить все данные внутри выборки и найти ДИ, скажем, по Клопперу-Пирсону, то это будет некорректно: выборка гетерогенная. Если найти ДИ бутстрепом индивидуальных частот - не будет учтена информация, что эти частоты оценены с разной точностью (по 500-2000 клеток). Как лучше всего поступить в этом случае?

Я чего-то не понял. Если частоты в выборках разные, то объединять их нельзя, поскольку они пришли из разных популяций. Соответственно, оценивается диапазон возможных "родительских" популяций для каждой выборки, например ДИ Клоппера-Пирсона (или Пуассона, если промилле). А то, что ширина интервалов будет разной, ну так это естественно...
nokh
Цитата(плав @ 3.07.2009 - 20:26) *
Я чего-то не понял. Если частоты в выборках разные, то объединять их нельзя, поскольку они пришли из разных популяций. Соответственно, оценивается диапазон возможных "родительских" популяций для каждой выборки, например ДИ Клоппера-Пирсона (или Пуассона, если промилле). А то, что ширина интервалов будет разной, ну так это естественно...

Имел в виду другое. Разные выборки я не объединяю, просто их несколько и для каждой хочу найти ДИ. Но проверяя каждую из этих выборок на внутреннюю однородность G-критерием (G2) часто обнаруживаю значимую гетерогенность. Получается, что раз выборка гетерогенная, то объединять клетки от всех особей в пределах выборки некорректно и правильнее рассчитывать ДИ по индивидуальным частотам входящих в выборку особей. Но значения этих частот получены с разной точностью.
Например, в выборке из 15 особей:
10 особей имели по 2 клетки с ядерными аномалиями на 2000 проанализированных клеток (т.е. по 0,1%),
3 особи - по 5 клеток из 1000 проанализированных (по 0,5%),
2 особи - по 10 клеток из 500 проанализированных (по 2%).
Я могу рассчитать ДИ для частоты каждой входящей в выборку особи, но не знаю как правильно рассчитать ДИ средней частоты клеток с аномалиями для всей этой выборки. (Если забыть о количестве проанализированных от каждой особи клеток и искать ДИ бутстрепом индивидуальных частот, то имеющаяся информация будет недоиспользована). Может здесь подойдет какая-то техника из мета-анализа?
плав
Цитата(nokh @ 3.07.2009 - 22:07) *
Имел в виду другое. Разные выборки я не объединяю, просто их несколько и для каждой хочу найти ДИ. Но проверяя каждую из этих выборок на внутреннюю однородность G-критерием (G2) часто обнаруживаю значимую гетерогенность. Получается, что раз выборка гетерогенная, то объединять клетки от всех особей в пределах выборки некорректно и правильнее рассчитывать ДИ по индивидуальным частотам входящих в выборку особей. Но значения этих частот получены с разной точностью.
Например, в выборке из 15 особей:
10 особей имели по 2 клетки с ядерными аномалиями на 2000 проанализированных клеток (т.е. по 0,1%),
3 особи - по 5 клеток из 1000 проанализированных (по 0,5%),
2 особи - по 10 клеток из 500 проанализированных (по 2%).
Я могу рассчитать ДИ для частоты каждой входящей в выборку особи, но не знаю как правильно рассчитать ДИ средней частоты клеток с аномалиями для всей этой выборки. (Если забыть о количестве проанализированных от каждой особи клеток и искать ДИ бутстрепом индивидуальных частот, то имеющаяся информация будет недоиспользована). Может здесь подойдет какая-то техника из мета-анализа?

Оптимальный вариант - оценка по пуассоновой регрессии со случайным фактором - особь. Зависимая переменная - число клеток, offset - количество проанализированных клеток.
nokh
Благодарю, попробую.
la.vi.na.
Здравствуйте хочу вернуться к этой теме. У меня подобное цитогенетическое исследование на растениях. До сих пор я считала различия исходя количества клеток как N., в силу того, что признак нарушение может проявиться "1" или не проявиться "0" в клетке, то есть признак является дихотомическим именно по отношению к клетке. Ведь разные частоты могут быть и в пределах одной особи (в разных корешках одного растения) в равной степени подверженных исследуемому фактору.
В связи с этим имеются вопросы:
1. Как будет все-таки оценивать значимость различий между выборками. если использовать предложенный мета-анализ, то пожалуйста, напишите как он реализуется в Statistica 6.0. или Эксель
2. Я оцениваю разность долей между опытом и контролем, использую метод Фишера через фи преобразование или Хи квадрат, с поправкой Йейтса по таблицам сопряженности 2х2. Итак, получаю в некоторых вариантах значимые различия на уровне p<0,05. Для графического отображения подсчитываю ДИ для каждой частоты методом Клоппера-Пирсона для этого же уровня значимости и обнаруживаю, что в некоторых случаях, где по критерию были значимые различия, ДИ перекрываются. Чему верить? Рецензент попросил отметить ошибки или ДИ на графике, а разве ошибки изображают в виде усов? Как поступить?
Буду признательна за ответы!
nokh
Цитата(la.vi.na. @ 23.11.2011 - 17:12) *
1. Как будет все-таки оценивать значимость различий между выборками. если использовать предложенный мета-анализ, то пожалуйста, напишите как он реализуется в Statistica 6.0. или Эксель
2. Я оцениваю разность долей между опытом и контролем, использую метод Фишера через фи преобразование или Хи квадрат, с поправкой Йейтса по таблицам сопряженности 2х2. Итак, получаю в некоторых вариантах значимые различия на уровне p<0,05. Для графического отображения подсчитываю ДИ для каждой частоты методом Клоппера-Пирсона для этого же уровня значимости и обнаруживаю, что в некоторых случаях, где по критерию были значимые различия, ДИ перекрываются. Чему верить? Рецензент попросил отметить ошибки или ДИ на графике, а разве ошибки изображают в виде усов? Как поступить?
Буду признательна за ответы!

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

2. С частотами сложность одна - что использовать в качестве единицы анализа: (1) Клетку или (2) Индивида (у вас, как я понял, это даже не растение а корешок).

(1) Клетка. Для того чтобы работать с клеткой нужно быть уверенным в полной однородности материала. Для этого составляется большая таблица сопряжённости 2хn, где 2 - это входы для абсолютных частот клеток с нарушениями и без, а n - разные индивиды. Анализируется хи-квадратом Пирсона или лучше отношением правдоподобия (G-критерий=критерий G-квадрат). Если индивиды не отличаются статистически значимо друг от друга (выборка однородна) - можно объединять данные по клеткам от всех индивидов в одну кучу. Т.е. для всей выборки будут только 2 числа: клетки с нарушениями и без нарушений. По ним и рассчитать долю с клеток с нарушениями. В этом случае рассчитывать ДИ для выборки следует, например, по Клопперу-Пирсону (есть и другие формулы, но я тоже считаю по К-П). Их и приводить на графике. Сравнивать такую выборку с другой такой же однородной выборкой следует также через анализ таблиц сопряжённости, как вы и делаете.

(2) Индивид. Если тест на однородность выборки не проходят, значит существует дополнительный источник изменчивости, который необходимо учитывать - индивидуальная вариабельность внутри выборки. Объединять всё нельзя. В этом случае можно поступить так: рассчитать доли клеток с нарушениями для каждого индивида, фи-преобразовать (преобразование арксинуса), рассчитать среднее арифметическое и 95%-ный ДИ для него (нормальная аппроксимация), а затем ретрансформировать полученные значения обратно в доли (или в %). Их и приводить на графике. ДИ будет тем сильнее асимметричен, чем сильнее средняя частота отклоняется от 0,5. Сравнивать такие неоднородные выборки следует с использованием в качестве отдельных значений фи-преобразованных частот, например, t-критерием Стьюдента или его модификацией Вэлча для неравных дисперсий (в зависимости от результатов проверки на равенство дисперсий F-критерием Снедекора).
Точно не помню, но такая схема выбора критериев для проверки (без упоминаний ДИ) давалась кажется в рекомендациях ВОЗ по тестам на мутагенность. Выделил "можно", т.к. можно и иначе: через логлинейную модель или даже через обобщённую линейную модель - в этом случае вся информация будет использована по максимуму и это будет самый правильный и мощный анализ. Но иначе пока не пробовал - обходился и так.

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

И ещё по графикам. По какой-то существующей традиции данные, представленные частотами принято представлять на графиках в виде столбчатых диаграмм. В этом случае столбики закрашивают разными цветами, и приводят только верхнюю границу ДИ, тогда как нижнюю или не приводят вовсе (думая что она равна верхней или ничего не думая) или она закрашивается штриховкой. Не повторяйте ошибок этих, зачастую именитых, неучей. В подавляющем большинстве случаев для средних частот (как бы их не считали) ДИ асимметричны, а значит нужно искать софт, в котором это можно задать. В Statistica через кнопочный интерфейс этого сделать нельзя (может можно программировать на Statistica Basic); я пользуюсь бесплатной версией KyPlot и получается очень красиво:
p2004r
Цитата(nokh @ 3.07.2009 - 21:07) *
Например, в выборке из 15 особей:
10 особей имели по 2 клетки с ядерными аномалиями на 2000 проанализированных клеток (т.е. по 0,1%),
3 особи - по 5 клеток из 1000 проанализированных (по 0,5%),
2 особи - по 10 клеток из 500 проанализированных (по 2%).


а что не отражает такой бутстреп?
Код
vyb1<-c(c(1,1),rep(0,2000-2))
vyb2<-c(rep(1,5),rep(0,1000-5))
vyb3<-c(rep(1,10),rep(0,500-10))

exp <- function () {
    c(replicate(10, sum(sample(vyb1, 2000, replace = TRUE))),
      2*replicate(3, sum(sample(vyb2, 1000, replace = TRUE))),
      4*replicate(2, sum(sample(vyb3, 500, replace = TRUE))))
}

hist(colMeans(replicate(100000, exp())))


ну или что там считать надо с выборкой?

Код
> replicate(10, exp())
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    2    0    2    3    3    4    3    1    2     2
[2,]    1    4    2    4    1    1    1    1    3     2
[3,]    1    1    2    3    1    3    3    4    0     2
[4,]    0    0    3    5    1    0    0    3    2     3
[5,]    5    4    0    4    1    4    1    2    2     2
[6,]    1    3    0    2    2    4    3    3    1     0
[7,]    2    0    3    2    2    1    0    0    1     2
[8,]    3    1    4    3    0    2    1    2    1     1
[9,]    2    2    2    3    0    0    4    3    3     1
[10,]    1    4    3    3    3    0    0    0    0     3
[11,]    8   16   10    4    8   12   18   14    4    12
[12,]   12    6   10   14   20   14   18   20   18    18
[13,]    8   16   18    6   10   16   10   10    8     4
[14,]   56   52   56   68   44   16   44   52   72    36
[15,]   48   40   36   36   24   52   32   48   36    36
nokh
Цитата(p2004r @ 23.11.2011 - 21:24) *
а что не отражает такой бутстреп?

Да наверное нормально отражает, просто я с ресэмплингом пока никак.
Цитата(p2004r @ 23.11.2011 - 21:24) *
ну или что там считать надо с выборкой?

95% ДИ. Кстати можно по этим выборкам рассчитать ДИ с поправкой на смещение (bias corrected accelerated)?
И что означает "exp()"?
p2004r
Цитата(nokh @ 23.11.2011 - 19:54) *
Да наверное нормально отражает, просто я с ресэмплингом пока никак.

95% ДИ. Кстати можно по этим выборкам рассчитать ДИ с поправкой на смещение (bias corrected accelerated)?
И что означает "exp()"?


это я функцию "эксперимент" объявил, что бы код не загромождать (что бы он понятнее стал, и перегрузил экспоненту smile.gif

Код
> expr <- function () {c(replicate(10, sum(sample(vyb1, 2000, replace = TRUE))), 2*replicate(3, sum(sample(vyb2, 1000, replace = TRUE))), 4*replicate(2, sum(sample(vyb3, 500, replace = TRUE))))}

> result<-colMeans(replicate(100000, expr()))

> summary(result)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  3.667   7.733   8.600   8.669   9.533  15.070

> sort(result)[2500]
[1] 6.2
> sort(result)[100000-2500]
[1] 11.33333


по 2.5% отрезал с каждой стороны
p2004r
Цитата(nokh @ 23.11.2011 - 19:54) *
Кстати можно по этим выборкам рассчитать ДИ с поправкой на смещение (bias corrected accelerated)?


надо заводить boot объект и скармливать его например bmem.ci.bca причем вместе с результатами складного ножа

мне почему то boot не нравится использовать, но наверное придется smile.gif
nokh
Благодарю за готовый код! Теперь когда понадобится - буду использовать. Именно BCa - не критично, просто из модного бутстрепа он более моден laugh.gif . Появится время - сравню ДИ по бутстрепу с ДИ по регресии.
p2004r
Цитата(nokh @ 23.11.2011 - 20:47) *
Благодарю за готовый код! Теперь когда понадобится - буду использовать. Именно BCa - не критично, просто из модного бутстрепа он более моден laugh.gif . Появится время - сравню ДИ по бутстрепу с ДИ по регресии.


оно оказывается вот как выглядит

Код
> bmem.ci.bc
function (par.boot, par0, cl = 0.95)
{
    se.boot <- apply(par.boot, 2, sd, na.rm = TRUE)
    estimate <- par0
    p <- ncol(par.boot)
    ci <- NULL
    for (i in 1:p) {
        ci <- rbind(ci, bmem.ci.bc1(par.boot[, i], par0[i], cl))
    }
    cbind(estimate, se.boot, ci)
}
<environment: namespace:bmem>

> bmem.ci.bc1
function (x, b, cl = 0.95)
{
    n <- length(x)
    z0 <- qnorm(sum(x < b, na.rm = TRUE)/n)
    alpha <- (1 - cl)/2
    alpha <- c(alpha, 1 - alpha)
    alpha <- sort(alpha)
    alpha1 <- alpha
    alpha <- pnorm(2 * z0 + qnorm(alpha))
    dig <- max(2L, getOption("digits"))
    np <- length(alpha)
    qs <- quantile(x, alpha, na.rm = TRUE)
    names(qs) <- paste(if (np < 100)
        formatC(100 * alpha1, format = "fg", width = 1, digits = dig)
    else format(100 * alpha1, trim = TRUE, digits = dig), "%",
        sep = "")
    qs
}
<environment: namespace:bmem>


они имели в виду под результатом бутстрепа что то не из boot smile.gif

вот такое выходит после загрузки library(bmem)

Код
>  bmem.ci.norm(as.data.frame(result),mean(c(rep(2,10),rep(5,3)*2,rep(10,2)*4)),cl=.95)
       estimate  se.boot     2.5%    97.5%
result 8.666667 1.315171 6.088979 11.24435

>  bmem.ci.bc(as.data.frame(result),mean(c(rep(2,10),rep(5,3)*2,rep(10,2)*4)),cl=.95)
       estimate  se.boot 2.5%    97.5%
result 8.666667 1.315171  6.2 11.33333
la.vi.na.
nokh спасибище за развернутый ответ!

[quote name='nokh' date='23.11.2011 - 20:59' post='12357']
"можно и иначе: через логлинейную модель или даже через обобщённую линейную модель"

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

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


nokh
Цитата(la.vi.na. @ 24.11.2011 - 08:02) *
... Где-то слышала,что можно логит-регрессию к таким данным применить, но сама не разобралась как правильно составить матрицу и интерпретировать результаты...

Тоже не разбирался, пока вполне обходился описанными подходами.
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Форум IP.Board © 2001-2025 IPS, Inc.