Критерий для анализа сверхмалых выборок, выборки из 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

[comisora]

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

CODE

## Создадим данные
a <- c(221.60112, 305.217725, 295.251684)
b <- c(371.3313, 397.452722, 437.212724)
ab <- c(a, b)
n <- length(ab)
g <- rep(c('a', 'b'), time = c(length(a), length(b)))

## Трансформация к диапазону (0; 1)
vec <- scales::rescale_max(ab,
from = c(
1/max(ab)/n,
max(ab)+max(ab)/n
)
)

## Бета-регрессия
fit <- betareg::betareg(vec~g)
emm <- emmeans::emmeans(fit, ~g)

## Средние и ДИ для групп
res <- data.frame(emm, row.names = 'g')[c('emmean', 'asymp.LCL', 'asymp.UCL')]
cev <- apply(res, 2, function(x) scales::rescale_max(x,
to = c(
1/max(ab)/n,
max(ab)+max(ab)/n
),
from = c(0, 1)))

## Разница средних и ДИ
ctr <- data.frame(confint(pairs(emm)))[c('estimate', 'asymp.LCL', 'asymp.UCL')]
rtc <- apply(ctr, 2, function(x) scales::rescale_max(x,
to = c(
1/max(ab)/n,
max(ab)+max(ab)/n
),
from = c(0, 1)))

## Расчёт средних в группах
set.seed(32167)
mean.a <- Hmisc::smean.cl.boot(a, B = 9999)
mean.b <- Hmisc::smean.cl.boot(b, B = 9999)
means <- rbind(cev, mean.a, mean.b)
colnames(means) <- c('mu', 'lwr', 'upr')

## Расчёт разницы средних и ДИ
set.seed(32167)
ci.bca <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'bca')
ci.per <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'perc')
ci.nor <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'norm')
ci.bas <- confintr::ci_mean_diff(a, b, type = 'bootstrap', boot_type = 'basic')
fun_extr <- function(x) { #Для единообразного извлечения
estimate <- x[['estimate']]
interval <- x[['interval']]
result <- c(estimate, interval)
return(result)
}
diffs <- rbind(
beta = rtc,
bca = fun_extr(ci.bca),
perc = fun_extr(ci.per),
norm = fun_extr(ci.nor),
basic = fun_extr(ci.bas)
)
colnames(diffs) <- c('mu', 'lwr', 'upr')

## Выводим результат
summary(fit)
print(means, digits = 5)
print(diffs, digits = 5)

> Call:
betareg::betareg(formula = vec ~ g)

Standardized weighted residuals 2:
Min 1Q Median 3Q Max
-1.7892 -1.0236 0.2118 0.9816 1.5627

Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.1480 0.1607 0.921 0.357
gb 1.1632 0.2525 4.606 4.1e-06 ***

Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 50.91 29.16 1.746 0.0808 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type of estimator: ML (maximum likelihood)
Log-likelihood: 8.148 on 3 Df
Pseudo R-squared: 0.7817
Number of iterations: 68 (BFGS) + 2 (Fisher scoring)
> mu lwr upr
a 273.88 233.94 313.82
b 401.80 369.21 434.40
mean.a 274.02 221.60 305.22
mean.b 402.00 371.33 437.21
> mu lwr upr
beta -127.92 -179.45 -76.395
bca -127.98 -191.06 -86.850
perc -127.98 -180.40 -78.143
norm -127.98 -180.30 -74.747
basic -127.98 -177.81 -75.553

Судя по результатам, средние и доверительные интервалы, разница средних, которые получены через бета-модель, аналогичны таковым при вычислении через бутстреп. Теперь попробуем повторить анализ для смеси распределений бета-регрессией, разными вариантами бутстрепа. Форма вывода - как у 100$.

CODE

##Дополнительные функции #fix
fun_int <- function(x) {
int <- x[['interval']]
return(int)
}
fun_res <- function(data, diff) {
len <- length(data)
res <- vapply(X = 1:n, FUN.VALUE = logical(1),
function(i) {
cond <- diff>data[[i]][1] & diff<data[[i]][2]
return(cond)
}
)
pr <- sum(res)/len
sums <- rowMeans(simplify2array(data))
res <- list(low = sums[1], up = sums[2], l = abs(sums[1]-sums[2]), prob = pr)
return(res)
}

x <- c(221.60112, 305.217725, 295.251684)
y <- c(371.3313, 397.452722, 437.212724)
mu_x <- mean(x)
mu_y <- mean(y)
len_x <- length(x)
sd_x <- sd(x)
sd_y <- sd(y)
len_y <- length(y)
diff <- mu_x-mu_y
len <- len_x+len_y
g <- rep(c('x', 'y'), time = c(len_x, len_y))
idx <- c('asymp.LCL', 'asymp.UCL')

n <- 1000
res <- vector(mode = 'logical', length = n)

set.seed(32167)
x_sim <- replicate(n,
c(
rnorm(round(len_x*2/3), mu_x, sd_x*1/3),
rnorm(round(len_x*1/3), mu_x, sd_x*2/3*2.5)
)
)
y_sim <- replicate(n,
c(
rnorm(round(len_x*2/3), mu_y, sd_y*1/3),
rnorm(round(len_x*1/3), mu_y, sd_y*2/3*2.5)
)
)
xy_sim <- rbind(x_sim, y_sim)
vec_sim <- apply(xy_sim, 2, function(i) scales::rescale_max(i,
from = c(
1/max(i)/len,
max(i)+max(i)/len
)
))

ci_beta <- lapply(
X = 1:n, FUN = function(i) {
fit <- betareg::betareg(vec_sim[,i]~g)
ci <- as.numeric(confint(pairs(emmeans::emmeans(fit, ~g)))[idx])
ic <- scales::rescale_max(ci,
to = c(
1/max(xy_sim[,i])/len,
max(xy_sim[,i])+max(xy_sim[,i])/len
),
from = c(0, 1)
)
return(ic)
})
ci_bca <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'bca', R = n)
ci <- fun_int(fit)
return(ci)
})
ci_perc <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'perc', R = n)
ci <- fun_int(fit)
return(ci)
})
ci_norm <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'norm', R = n) #fix
ci <- fun_int(fit)
return(ci)
})
ci_basic <- lapply(
X = 1:n, FUN = function(i) {
fit <- confintr::ci_mean_diff(x_sim[,i], y_sim[,i],
type = 'bootstrap', boot_type = 'basic', R = n) #fix
ci <- fun_int(fit)
return(ci)
})
vct <- ls(pattern = 'ci_')

res_ci <- t(simplify2array(lapply(1:length(vct), function(i) {
#do.call(fun_res, list(data = get(vct[i]), diff = get('diff')))
fun_res(data = get(vct[i]), get('diff')) #fix
})))
row.names(res_ci) <- vct
print(res_ci, digits = 5) # fix

low up l prob
ci_basic -173.1 -81.769 91.328 0.885
ci_bca -176.03 -81.595 94.435 0.773
ci_beta -172.44 -79.885 92.552 0.866
ci_norm -174 -80.332 93.666 0.876
ci_perc -172.38 -81.104 91.271 0.814

Сообщение отредактировал comisora - 18.06.2022 - 21:58