Сравнение частично зависимых выборок, 95% ДИ и ресэмплинг-техниками Опции

[nokh]

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

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

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

Пример 2. Работа по ветеринарии. Из разных популяций взяли выборки животных и определили в них 1) животных без явных нарушений по комплексу показателей, 2) со слабыми отклонениями от нормы (типа предпатология) и 3) сильными отклонениями, указывающими на разную патологию. Посчитали %, сравнили популяции между собой. Теперь стало нужно определить какие показатели и насколько отклоняются от нормы в ту или другую сторону в популяциях и о чём это говорит. Поскольку норма из книжек/статей не очень хороша и к тому же "плывёт" во времени, было решено скомпоновать свою норму: из всех популяций всех здоровых животных объединили в одну группу, охарактеризовали, в планах - расчёт референтных ("референсных") интервалов (кстати для их расчётов по EPC28A3C / C28-A3c всё нужное есть в MedCalc и R ("referenceIntervals" и др.)). Теперь нужно сравнить каждую популяцию с этой нормой, чтобы оценить по каким показателям и в какую сторону отклонения, рассчитать дельту в %. Таким образом имеем две оценки, как некие самостоятельные характеристики:
1) популяционные оценки показателей (медианы не для всех годились, беру средние с ДИ бутстрэпом) и
2) значения нормы, но рассчитанные частично по животным, входящим в эти популяции.
Сначала тоже хотел сравнить по 95% ДИ, а потом решил получить р-значения и вышел на литературу по сравнению partially overlapping samples, ссылки на пару статей - внизу. В ней есть много подходов и формул, но не увидел сопоставления 95% ДИ и вполне логичного сравнения средних в таких группах техникой Монте-Карло. В связи с чем начал сомневаться, читать про гипотезы, которые проверяют ресэмплинг-техники и т.п. В этом втором примере у меня есть полная информация по степени и структуре перекрытия выборок, и сравнить по формулам из статей в принципе можно, но муторно и не универсально, хочу универсального подхода для кода в R)))

А вы что по этому поводу думаете? Буду признателен за советы/критику/парустатей.

https://www.tqmp.org/RegularArticles/vol18-1/p055/p055.pdf
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9042157/

Сообщение отредактировал nokh - 25.01.2024 - 23:13

[comisora]

Доброго дня.

Первая задача мне напоминает single case analysis или time series, при котором мы имеет только две точки для анализа. Я пишу про две точки, так как в данном случае оценивается организация в целом и происходящие в ней процессы. Ничего плохого в допущении об независимости наблюдений при опросе не вижу, особенно в случае, когда информация полностью анонимная. Если есть возможность, сведите Ваши данные к тестированию тренда. В рамках "продуктивных методов анализа" время опроса в организации определить как зависимую переменную, а ответы - как независимые. Далее посмотреть характер паттерна ответов.

По поводу второй задачи. Если Вы объединили данные и получили параметры распределения, которое Вас устраивает, то нагенерируйте вагон референтных данных, сравните с ними Ваши выборки, усредните тестовые статистики. В синтетических данных реальных животных не будет.

Вообще, если известны id испытуемых, то реально использовать mixed model. В коде изложил результаты беглого поиска по проблеме частично парным выборкам.

CODE

library(lmerTest)
library(robustrank)
library(IncomPair)
library(contingencytables)

# Зависимая переменная
y <- c(
1, 2, 3, 3, 2, 1, 2, 4, 1, 2,
3, 4, 3, 4, 3, 2, 1, 4, 3, 2
)

# Независимая переменная
x <- rep(c('a', 'b'), each = 10)

# Номер
id <- paste0(
'id',
c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
1, 2, 3, 4, 5, 6, 11, 12, 13, 14)
)

# Объединяем данные
df <- data.frame(id = factor(id), x = factor(x), y = y)

# Строим частично смешанную модель
fit <- lmer(y ~ x + (1|id), data = df)

# Плохо задокументированный пакет - пользоваться с осторожностью
robustrank::pm.wilcox.test(
Xpaired = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) > 1))],
Ypaired = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) > 1))],
Xextra = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) == 1))],
Yextra = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) == 1))]
)
robustrank::mw.mw.2.perm(
X = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) > 1))],
Y = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) > 1))],
Xprime = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) == 1))],
Yprime = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) == 1))],
.corr = -.5
)

# Более понятная библиотека
IncomPair::permb(
xp = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) > 1))],
yp = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) > 1))],
xu = df$y[df$x == 'a' & df$id %in% names(which(table(df$id) == 1))],
yu = df$y[df$x == 'b' & df$id %in% names(which(table(df$id) == 1))],
r = .5
)

# Бонусом
contingencytables::JonckheereTerpstra_test_rxc(contingencytables::table_7.7)