Критерий для анализа сверхмалых выборок, выборки из 3 значений. Применение стат. анализа Опции

[Vitek_22]

Салют! Столкнулся с проблемой сравнения двух выборок, в каждой из которых по 3 значения. Это результаты иммуноблоттинга (определение концентрации целевого белка в пробе) очень ценных образцов, полученных от трансгенных животных. Но их - образцов, было всего 3 для каждой из групп (2 группы: интактная и подвергнутая воздействию исследуемого фактора). Покопавшись в литературе, нашёл статьи, где есть такие выборки и авторы как ни в чём не бывало используют t-критерий Стьюдента для сравнения средних. Нашёл статьи, где используют U-критерий Манна-Уитни... Скажем так, в биологии 3 образца - это нормально для публикации, если речь идёт об особо ценном и сложнополучаемом биоматериале (к примеру как у меня, когда животные практически не дают потомства). Т.е. представить эти данные можно и не стыдно. Но вот как сравнить, как показать, что эти выборки отличаются статистически значимо, иными словами, что наш исследуемый фактор значимо повлиял на концентрацию целевого белка?
Почитал ещё о таком методе, как ресамплинг или бутстреп, когда объём выборки искусственно увеличивают. Ну, не знаю насколько это правильно... также не нашёл софт и чёткого понимания как это сделать у меня нет.
Работаю в проге Statistica 12

Посоветуйте, как всё же обработать эти данные. Вот пример исходных цифр:
Выборка 1:
221,60112
305,217725
295,251684

Выборка 2:
371,3313
397,452722
437,212724

[ИНО]

Хороший вопрос! Документация пакета boot (сверяя с которым, я и пришел к вышеозначенному заключению относительно типа ДИ в PAST) лаконично отвечает на него следующее: "The formulae on which the calculations are based can be found in Chapter 5 of Davison and Hinkley (1997)."
Ну и библиографическая ссылка в конце: Davison, A.C. and Hinkley, D.V. (1997) Bootstrap Methods and Their Application, Chapter 5.Cambridge University Press. У меня такой книжки нету.

Но попадались статьи, в которых некий "basic bootstrap CI" нещадно поносили и хвалили BCa. Особо не вникал (учитывая, что, повторяю, саму глубинную идею бутстрепа так и не смог постичь), просто пользовал по возможности BCa. Скажем так, впал в грех слепого следования авторитетам

Немного помоделировал (осбо в это дело не умею, так что, если что не так, тапками не кидайте)

Загрузка "бигдаты" ТС (пожалуй, стоит уже в рабочем пространстве сохранить - как ни крути, классический датасет получился, ирисы Фишера отдыхают) :

Код

x<-c(221.60112, 305.217725, 295.251684)
y<-c(371.3313, 397.452722, 437.212724)

Генерация 10000 пар выборок из нормальных распределений с похожими параметрами, построение 95% ДИ методом Уэлча и проверка, попадает ли в них истинная разница матожиданий:

Код

diff<-274-402
res<-rep(0, 10000)
for (i in 1:10000)
{
x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(y), sd(y))
ci<- t.test(x_sim, y_sim)$conf.int
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
#Доля попаданий:
sum(res)/10000

Вывод: истинный доверительный уровень на смоделированных данных получился чуть выше номинального, но разница не существенна

Теперь возьмем выборки, каждую из смеси двух нормальных распределений с одинаковыми мю но разными сигмами. Страшненькая формула генератора получена методом тыка, так чтобы стандартные отклонения все еще оставались похожими на оригинал. Если умеете генерировать смеси с заданной дисперсией, замените на свою элегантную:)

Код

diff<-274-402
res<-rep(0, 10000)
for (i in 1:10000)
{
x_sim<-c(rnorm(2, mean(x), sd(x)*1/3), rnorm(1, mean(x), sd(x)*2/3*2.5))
y_sim<-c(rnorm(2, mean(y), sd(y)*1/3), rnorm(1, mean(y), sd(y)*2/3*2.5))
ci<- t.test(x_sim, y_sim)$conf.int
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
sum(res)/10000

А здесь уже не все так хорошо: метод Уэлча дает примерно 99% интервал вместо желаемого 95%. И туда уже обычно ноль таки попадает!

Справится ли BCa лучше?

Код

library(confintr)
diff<-274-402
res<-rep(0, 1000)
for (i in 1:1000)
{
x_sim<-c(rnorm(2, mean(x), sd(x)*1/3), rnorm(1, mean(x), sd(x)*2/3*2.5))
y_sim<-c(rnorm(2, mean(y), sd(y)*1/3), rnorm(1, mean(y), sd(y)*2/3*2.5))
ci<- ci_mean_diff(x_sim, y_sim, type="boot", boot_type="bca")$interval
if(diff>ci[1]&diff<ci[2]) res[i]<-1
}
sum(res)/1000

Если дождетесь выдачи, сообщите мне. Я не смог, несмотря на то, что урезал количество симуляций до 1000. Можно, конечно, еще и количество псевдовыборок на каждом шаге урезать (параметр R в ci_mean_diff()), но тогда уж совсем как-то неубедительно получится. Философский вопрос: существует ли в мире более медленный язык программирования чем R?

P.S. Таки нашел в Шитиков,Розенберг "Рандомизация и бутстреп:статистический анализ в биологии и экологиис использованием R", что что за бэйзик-интевал такой: это тот же процентильный, но перевернутый относительно среднего. Теоретические основы данной процедуры не понял. Но ,в третий раз повторяю: я и основы бутсрепа вообще плохо понимаю. Но судя по тому, как извращаются с этими интервалами, начиная от переворота в случае бэйзика этого, заканчивая расчетом неких параметров смещения и акселерации в случае BCa, с "восстановлением распределения путем замены одних вариант другими" на самом деле все совсем не так просто, как это иногда преподносят. И интуиция тут работает неважно (моя, по крайней мере).

Сообщение отредактировал ИНО - 17.06.2022 - 18:01

[100$]

Вывод: истинный доверительный уровень на смоделированных данных получился чуть выше номинального, но разница не существенна

При параметре ф-ции t.test(..., var.equal=TRUE)
оценка доверительной вероятности составляет 94,88%, при var.equal=FALSE, да, зашкаливает за 96%.

Справится ли BCa лучше?Если дождетесь выдачи, сообщите мне

Эт можно.

С параметром R=999 в ci_mean_diff() при общем числе повторностей=1000 мой комп ковырялся 3,5 мин. Разумеется, при параметре по умолчанию R=9999 ковырялся бы на порядок дольше. Я бы тоже не дождался.

Мой результат
LOW______ UP _____ L________ PROB
-147.78218 -114.44736 33.33482 0.77600

где L- длина полученного интервала (она загуливает совершенно диким образом), PROB=sum(res)/1000 - оценка доверительной вероятности. Она находится в диапазоне ,765 - ,778.

Не хочется превращать мероприятие в форум для программистов, однако для исходных выборок x и y эмпирическая корреляция равна cor(x,y)=.73041, что, конечно, надо учитывать при моделировании системы случайных величин. Конструкция

Код

x_sim<-rnorm(3, mean(x), sd(x))
y_sim<-rnorm(3, mean(y), sd(y))

смотрится просто неприлично.



Сообщение отредактировал 100$ - 18.06.2022 - 01:30