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

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

4 страниц V   1 2 3 > » 

comisora
Отправлено: 18.02.2024 - 12:49





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


Цитата(ИНО @ 18.02.2024 - 11:13) *
"в коллективе стало больше знатоков Достоевского благодаря проведенной руководством в интервале между опросами просветительской работе (а не потому что половина сотрудников, не пожелавшая учить Достоевского, уволилась smile.gif)" недопустимо.

Если в результате "просветительской работы" половина уволилась и доля незнающих Достоевского упала, то задача руководства по "ликвидации безграмотности" решена crazy.gif . В конце концов если прийдут новые уже знающие Достоевского или желающие его изучать, то предприятие только выиграет.
  Форум: Медицинская статистика · Просмотр сообщения: #28865 · Ответов: 24 · Просмотров: 3114

comisora
Отправлено: 18.02.2024 - 12:14





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


Цитата(ИНО @ 18.02.2024 - 11:34) *
Зачем GAM если есть только две точки для каждого животного. Чем обычные линейный модели не угодили? Кривую все равно по двум точкам не построить.

Идея с подстановкой нуля умершим меня тоже посещала, но была отброшена как неправомерная. Почему 0, а не 0,001 или -1?

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


Из коробки один синтаксис позволяет строить lm/glm, betareg и заодно прикрутить случайный эффект. На двух точка кривую, конечно, не построить, но их соединить прямой линией библиотека может.
Исходя из условий задачи, 0 - это какое-то время. Я допустил, что 0 - это железный показатель того, что "пациент мёртв" и никакой активности не может быть по определению. Отрицательное значение времени (исходя из задачи) ситуация невозможная. Создавать колонку статуса "жив/мертв" смысла нет, если мы знаем, что 0 - это смерть. Считаю, что подставлять мы имеем право. Другое дело, что придётся подумать, как этот 0 учесть, чтобы он проблем не создавал. Первое, что приходит в голову - Zero inflated beta (ZOIB) regression. Если у всех результат измерения укалдывается в 60 с, то мы легко перейдём к данной модели.

Еще лучше знать дату смерти особи, чтобы задачу свести к анализу выживаемости. Вторая модель (на живих) уже сделана.
  Форум: Медицинская статистика · Просмотр сообщения: #28864 · Ответов: 28 · Просмотров: 10551

comisora
Отправлено: 18.02.2024 - 00:07





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


Цитата(nokh @ 4.02.2024 - 21:43) *
Благодарю участников за мнения. Как-то безрадостно пока...
Ну а что всё-таки думаете по поводу 95% ДИ? (Пока второй пример отложим).

Вот есть больница, пусть это генеральная совокупность. Оценили показатели в интересующих выборках персонала, рассчитали 95% ДИ. После мероприятий через год смотрим снова и снова частоты с 95% ДИ. Если не перекрываются, значит отличия значимы. Анонимность, по-моему здесь роли не играет, в том плане, что делает сравнение невозможным. Я когда аспирантом полёвок в лесах ловил никакого их реестра не было, да и потом не метили никак. Кстати хороший пример: экологи изучают вид на какой-то территории какое-то время. Ведётся мониторинг чего-то и всё анонимно. Играет роль время, которое прошло между двумя исследованиями. Ну, например, если 5 лет прошло, то различия скорее можно объяснить просто изменением контингента. Ну а если эти выборки процентов на 90-95% перекрываются, т.е. это почти те же люди и прошёл всего только год - логично приписать различия проведённой работе с персоналом. Мне видится, что сопоставление ДИ одинаково применимо как к независимым зависимым, так и к зависимым (в том числе - частично) выборкам. Или нет?

Как и ранее считаю, что пользоваться такой логикой допустимо. Впрочем, есть предложение выполнить на полном наборе MCA и отталкиваться от трансформированных данных.
Статья - для сведения.
Прикрепленные файлы
Прикрепленный файл  Greenland_et_al_2016_Statistical_tests__P_values__confidence_intervals__and_power.pdf ( 452,95 килобайт ) Кол-во скачиваний: 11
 
  Форум: Медицинская статистика · Просмотр сообщения: #28861 · Ответов: 24 · Просмотров: 3114

comisora
Отправлено: 17.02.2024 - 23:31





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


Цитата(Vitek_22 @ 17.02.2024 - 15:25) *
Что значат эти "остатки модели"?

Посмотрите учебные материалы по ссылке.

Цитата(Vitek_22 @ 17.02.2024 - 15:25) *
Разницу вычислять - не очень, т.к. разный размер выборок (у части животных симптоматика прогрессировала быстрее и они погибли раньше повторного анализа). как в этом случае произвести вычитание?


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

В коде ниже указана попытка сделать анализ повторных измерений при допущении различных трансформаций. Лучшей оказалась модель с предварительным логарифмированием зависимой переменной. На графике видно, что группы изначально различались, к концу исследования эти различия усилились. На мой взгляд, если различия объясняются эффектом лечения, то результат (сила эффекта) достаточно скромный. Чтобы в этом убедиться, нужно или усилить эффект от лечения или увеличить количество наблюдений.

CODE
library(mgcv)
txt <- "id value group time
1 31 control first
2 40 control first
3 19 control first
4 55 control first
5 27 control first
6 3 control first
7 40 control first
8 29 control first
9 17 control first
10 46 control first
11 55 treatment first
12 55 treatment first
13 38 treatment first
14 55 treatment first
15 55 treatment first
16 55 treatment first
17 27 treatment first
18 40 treatment first
1 0,7 control second
2 1,4 control second
3 0,5 control second
4 3 control second
5 0,5 control second
6 1,5 control second
7 1,6 control second
11 3 treatment second
12 6,6 treatment second
13 2,4 treatment second
14 2,1 treatment second
15 3 treatment second
16 1,4 treatment second
17 4 treatment second
18 9,7 treatment second"
df <- read.table(text = txt, header = TRUE, dec = ',')
#df <- df[-which(df$id %in% c(8, 9, 10)),]
df$group <- factor(df$group, levels = c('control', 'treatment'))
df$time <- factor(df$time, levels = c('first', 'second'))
df$val1 <- log(df$value)
df$val2 <- df$value/60
fitA <- gam(value~time*group + s(id, bs = 're'),
family = gaussian(link = 'identity'),
data = df)
fitB <- gam(value~time*group + s(id, bs = 're'),
family = gaussian(link = 'log'),
data = df)
fitC <- gam(value~time*group + s(id, bs = 're'),
family = gaussian(link = 'inverse'),
data = df)
fit1 <- gam(val1~time*group + s(id, bs = 're'),
family = gaussian(link = 'identity'),
data = df)
fit2 <- gam(val2~time*group + s(id, bs = 're'),
family = betar(link = 'logit'),
data = df)
mods <- ls(pattern = 'fit')
modSum <- do.call(cbind,
lapply(seq_along(mods),
function(i) {
mod <- summary(get(mods[i]))
result <- rbind(mod$dev.expl, mod$r.sq)
row.names(result) <- c('DevExpl', 'RsqAdj')
colnames(result) <- mods[i]
return(result)
}
)
)

pred <- predict(fit1,
newdata = expand.grid(group = levels(df$group),
time = levels(df$time)),
newdata.guaranteed = TRUE,
type = 'response',
exclude = 's(id)',
se.fit = TRUE)
pred$cil <- pred$fit - pred$se.fit*qnorm(0.975)
pred$ciu <- pred$fit + pred$se.fit*qnorm(0.975)
pred$mu <- exp(pred$fit)
pred$lwr <- exp(pred$cil)
pred$upr <- exp(pred$ciu)
pred <- data.frame(pred)
pred <- cbind(pred,
expand.grid(group = levels(df$group), time = levels(df$time))
)
pl <- ggplot(pred, aes(x = time, y = mu, ymax = upr, ymin = lwr,
shape = group)) +
geom_pointrange(position = position_dodge(0.1)) +
theme_bw() + theme(strip.text = element_blank()) +
scale_y_log10()


Цитата(ИНО @ 17.02.2024 - 21:07) *
Единственный недостаток - тормоза, но сомневаюсь, что Ваши выборки столь велики, чтобы Вы из заметили.

Это не баг, это - фича (с).
Эскизы прикрепленных изображений
Прикрепленное изображение
 

Прикрепленные файлы
Прикрепленный файл  txt.txt ( 3,94 килобайт ) Кол-во скачиваний: 11
 
  Форум: Медицинская статистика · Просмотр сообщения: #28860 · Ответов: 28 · Просмотров: 10551

comisora
Отправлено: 17.02.2024 - 13:25





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


Цитата(Vitek_22 @ 16.02.2024 - 19:14) *
Есть две группы животных, подвергнуты лечению и контроль. И для каждой группы было проведено 2 анализа в разных временных точках. Причём, во второй временной точке выборки подчиняются нормальному распределению (Шапиро-Уилка), но всё равно отклонение от нормальности субкритическое так сказать, а во второй - нет. Снижение времени - показатель прогрессии заболевания.

Выборки проверять на соответствие нормальному распределению не надо. В линейной модели предполагается нормальное распределение остатков. В примере показано, как можно создать выборки, по которым тест Шапиро-Уилка отклонит нулевую гипотезу в двух из трех подгрупп, но не по остаткам модели. В приложении - статьи по процедуре проверке нормальности.

CODE
set.seed(324)
n <- 1000
nu <- 1.5
sd <- 1
x <- gl(n = 3, k = n, labels = c('a', 'b', 'c'))
a <- gamlss.dist::rSN2(n = n, mu = 160, sigma = 1, nu = nu)
b <- rnorm(n = n, mean = 140, sd = 1)
c <- gamlss.dist::rSN2(n = n, mu = 120, sigma = 1, nu = 1/nu)
y <- c(a, b, c)
fit <- lm(y~x)
res <- resid(fit)
lst <- list(a, b, c, y, res)
round(
sapply(seq_along(lst), function(i) {shapiro.test(lst[[i]])[['p.value']]}),
digits = 3
) #[1] 0.000 0.847 0.000 0.000 0.634
dev.new(noRStudioGD = TRUE)
layout(matrix(c(1,2,3,4,5,5), nrow = 3, ncol = 2, byrow = TRUE))
plot(density(a))
plot(density(b))
plot(density©)
plot(density(y))
plot(density(res))


По Вашим данным можно сделать следующее:
0. Пояснить следующее: измерение через месяц - это реальное время или величина изменения относительно первоначального замера?
1. Взять разницу между состояниями "до-после" и сравнивать её между группами любым тестом, который Вам нравится (Стьюдента/Уэлча, Манна-Уитни, точный тест Фишера-Питмана и т.п.).
2. Ознакомиться с выполнением смешанного дисперсионного анализа по ссылке за авторством nokh.
3. Попытаться выполнить трансформацию (логарифмирование) и посмотреть, что будет.
Эскизы прикрепленных изображений
Прикрепленное изображение
 

Прикрепленные файлы
Прикрепленный файл  1471_2288_12_81.pdf ( 790,19 килобайт ) Кол-во скачиваний: 8
Прикрепленный файл  p302.pdf ( 1,96 мегабайт ) Кол-во скачиваний: 14
 
  Форум: Медицинская статистика · Просмотр сообщения: #28855 · Ответов: 28 · Просмотров: 10551

comisora
Отправлено: 28.01.2024 - 11:33





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


Доброго дня.

Первая задача мне напоминает 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)
  Форум: Медицинская статистика · Просмотр сообщения: #28819 · Ответов: 24 · Просмотров: 3114

comisora
Отправлено: 13.12.2023 - 20:53





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


2 ИНО

https://packagemanager.posit.co/client/#/repos/cran/setup?snapshot=2023-12-13 - посмотрите этот проект. Я лично не пользовался, но указание на "snapshot" позволяют предполагать функционал, аналогичный MRAN.
  Форум: Медицинская статистика · Просмотр сообщения: #28805 · Ответов: 3 · Просмотров: 1555

comisora
Отправлено: 2.07.2023 - 19:13





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


2 ИНО

Как я понял, нужно выбрать тип распределения для categorical depended variable, раскидать по нужным полям фиксированные и случайные эффекты и попытаться запустить модель. Однако это может не сработать так как распределения предназначены(?) для моделирования дихотомических переменных. Можно выкрутиться, если знать условные нормы для интерлейкинов или разбить данные на интересующие области, хотя дихотомизация не всегда удачный выбор для анализа.

На мой взгляд, порядковая регрессия является способом решения данной задачи, так как все интервальные переменные - это порядковые данные (обратное, разумеется, не только лишь верно). В данной публикации рассмотрен пример моделирования клеток CD4 при ВИЧ-инфекции (настоящие дискретные количественные данные). Чем непараметрический аналог параметрической ANOVA с повторными измерениями?

Раз 100$ затронул MANOVA, то приведу ссылку на реализацию, которая не требует специальных допущений, что в случае с интерлейкинами выглядит неплохо.
  Форум: Медицинская статистика · Просмотр сообщения: #28658 · Ответов: 28 · Просмотров: 10551

comisora
Отправлено: 1.07.2023 - 13:25





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


2 Varta

В Вашем случае можно использовать порядковую регрессию со смешанными эффектами, так как данная модель обобщает критерии Манна-Уитни и Краскела-Уоллиса. Попробуйте программу Jamovi и модуль GAMLj к ней.
  Форум: Медицинская статистика · Просмотр сообщения: #28654 · Ответов: 28 · Просмотров: 10551

comisora
Отправлено: 20.02.2023 - 13:57





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


2 passant

Добрый день. Мой ответ не касается Вашего исходного вопроса, но при прочтении Ваших пояснений мне вспомнились тесты the Bartels rank test of randomness и библиотека для Change point analysis. Может на Вашу задачу следует посмотреть с этой стороны?

Цитата(passant @ 20.02.2023 - 13:11) *
В какой момент - при какой доле в выборке - эти данные становятся значимыми?


Есть ситуация, где нужна 100% уверенность, поэтому нужно предпринимать все возможные меры по профилактике и решению образовавшихся проблем. В таком случае вопрос о статистической значимости смысла не имеет. Если есть возможность построить зависимость условных расходов от доли условного сигнала (хотя бы на уровне допущения), то нужно найти ту величину условных расходов и соответствующую долю условного сигнала, которая будет Вашим ориентиром. Расходы и сигнал - любые метрики успешности процесса, хоть количество расходов на каждый неправильно заполненный бланк в отчёте маркетолога формата А4.

На мой взгляд, готовность платить или принимать решение (Cost-effectiveness thresholds) зависит исключительно от Вашей задачи.
  Форум: Медицинская статистика · Просмотр сообщения: #28093 · Ответов: 25 · Просмотров: 24018

comisora
Отправлено: 20.11.2022 - 23:53





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


Цитата(ИНО @ 20.11.2022 - 22:29) *
comisora, эта гибрид из логичстической модели и ZPT, применяется. наоборот, когда в выборке нулей овер дофига. Наверное, можно расковырять готовую функцию, и оторвать половину гибридной модели, но получим все тот же pospoisson.


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

CODE

n = 10
lambda = 0.5

set.seed(32167)
x <- qpois(p = runif(n, min = dpois(0, lambda)), lambda) #rpospois()

set.seed(32167)
y <- qpois(p = runif(n, min = 0), lambda)

set.seed(32167)
z <- rpois(n, lambda)

print(x);print(y);print(z)

[x] 1 1 1 1 1 2 2 1 1 2
[y] 0 0 0 0 0 1 2 0 0 2
[z] 0 0 0 0 0 1 2 0 0 2



Я пока не придумал как адаптировать Ваш код для доказательства/опровержения утверждения, что Ваш генератор и rpospois() даёт одно и тоже. Пока только демонстрирую разницу между rpospois() и rpois() с невозможностью арифметическими операциями свести одно распределение к другому. Даже если мы возьмём и создадим последовательность Пуассона без нулей, то она не всё равно не совпадёт при большом количестве наблюдений и lambda менее 10. Дальнейшее увеличение lambda приводит к расхождению рядов генерации.

CODE


n = 10000
lambda = 9
set.seed(32167)
x <- qpois(p = runif(n, min = dpois(0, lambda)), lambda)
set.seed(32167)
y <- qpois(p = runif(n, min = 0), lambda)
set.seed(32167)
z <- rpois(n, lambda)
table(x-y);table(x-z);table(y-z)


  Форум: Медицинская статистика · Просмотр сообщения: #27871 · Ответов: 81 · Просмотров: 149423

comisora
Отправлено: 20.11.2022 - 22:21





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


Цитата(ИНО @ 20.11.2022 - 18:36) *
Цель предложенных Вами методов неясна...

Чтобы зависимые переменные были одинаковые. В таком случае сравниваться будет только результат влияния параметра family.

Цитата(ИНО @ 20.11.2022 - 18:36) *
Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским.


Раз в естественном процессе сбора есть нули в результате отбраковки, почему бы не использовать hurdle-модель?
  Форум: Медицинская статистика · Просмотр сообщения: #27869 · Ответов: 81 · Просмотров: 149423

comisora
Отправлено: 20.11.2022 - 12:47





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


Цитата(ИНО @ 20.11.2022 - 12:07) *
Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать именно то распределение, данные из которого я генерирвал unknw.gif


Код
sim<-rpois(2000, lambda_i)


2 ИНО

Может попробуете погенерировать VGAMdata::rpospois()?

Мне представляется некорректным Ваш подход исходного вычитания единицы в одном случае, так как в другом случае Вы этого не делаете. Считаю, что проверять модели надо в одинаковых условиях, а именно после генерации (два варианта):
1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.
Также попробуйте сравнить модели по другим характеристикам (AIC, MSE и т.п.).
  Форум: Медицинская статистика · Просмотр сообщения: #27862 · Ответов: 81 · Просмотров: 149423

comisora
Отправлено: 19.11.2022 - 13:58





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


2 ИНО

Ещё немного материала.

"In this chapter, we discuss models for zero-truncated and zero-inflated count data. Zero truncated means the response variable cannot have a value of 0. A typical example from the medical literature is the duration patients are in hospital. For ecological data, think of response variables like the time a whale is at the surface before re-submerging, counts of fin rays on fish (e.g. used for stock identification), dolphin group size, age of an animal in years or months, or the number of days that carcasses of road-killed animals (amphibians, owls, birds, snakes, carnivores, small mammals, etc.) remain on the road. These are all examples for which the response variable cannot take a value of 0.

On their own, zero-truncated data are not necessarily a problem. It is the underlying assumption of Poisson and negative binomial distributions that may cause a problem as these distributions allow zeros within their range of possible values. If the mean is small, and the response variable does not contain zeros, then the estimated parameters and standard errors obtained by GLM may be biased. In Section 11.2, we introduce zero-truncated Poisson and zero-truncated negative binomial models as a solution for this problem. If the mean of the response variable is relatively large, ignoring the truncation problem, then applying a Poisson or negative binomial (NB) generalised linear model (GLM), is unlikely to cause a problem. In such cases, the estimated parameters and standard errors obtained by Poisson GLM and truncated Poisson GLM tend to be similar (the same holds for the negative binomial models).

In ecological research, you need to search very hard to find zero-truncated data. Most count data are zero inflated. This means that the response variable contains more zeros than expected, based on the Poisson or negative binomial distribution. A simple histogram or frequency plot with a large spike at zero gives and early warning of possible zero inflation. This is illustrated by the graph in Fig. 11.1, which shows the numbers of parasites for the cod dataset that was used in Chapter 10 to illustrate logistic regression. In addition to presence and absence of parasites in cod, Hemmingsen et al. (2005) also counted the number of parasites, expressed as intensity" (Zuur AF, Ieno EN, Walker N, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R [Internet]. New York, NY: Springer New York; 2009 [cited 2022 May 30]. (Statistics for Biology and Health). Available from: http://link.springer.com/10.1007/978-0-387-87458-6).

Дальнейший поиск следует вести по следующим ключевым словам: "the Shifted Poisson distribution", "the Positive Poisson distribution". Так по второму ключевому слову нашёл краткое упоминание в работе автора библиотеки VGAM (Yee TW. Vector Generalized Linear and Additive Models [Internet]. New York, NY: Springer New York; 2015 [cited 2020 Oct 13]. (Springer Series in Statistics). Available from: http://link.springer.com/10.1007/978-1-4939-2818-7).

Далее приведу пример, демонстрирующий разницу между распределениями:

CODE
set.seed(32167)
table(VGAMdata::rpospois(1000,0.5))

1 2 3 4
763 194 39 4

set.seed(32167)
table(rpois(1000, 0.5))

0 1 2 3
601 303 78 18


Рекомендую к просмотру презентацию, в которой объяснены некоторые идеи.
  Форум: Медицинская статистика · Просмотр сообщения: #27855 · Ответов: 81 · Просмотров: 149423

comisora
Отправлено: 18.11.2022 - 21:57





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


2 ИНО

У меня не хватает компетенции ответить на Ваш исходный вопрос (есть более квалифицированные форумчане). Однако посчитал возможным дать ссылки на материалы, которые мне попались при поиске информации по вопросу:

1. Difference between shifted distribution and zero-truncated distribution - по ссылке демонстрация разницы;
2. Simulation of the Shifted Poisson Distribution with an Application to the CEV Model - теоретическая статья, в которой обсуждается генерация такого распределения.
3. Will adding a constant to a random variable change its distribution? - обсуждение того, как меняется распределение при добавлении константы.

  Форум: Медицинская статистика · Просмотр сообщения: #27854 · Ответов: 81 · Просмотров: 149423

comisora
Отправлено: 19.07.2022 - 00:19





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


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

1. Описание критерия Манна-Уитни, что он умеет и чего не умеет. Рассмотрены случаи с медианами.
2. Вариант использования proportional odds regression для непрерывной зависимой переменной.



Прикрепленные файлы
Прикрепленный файл  Divine_et_al_2018_The_Wilcoxon_Mann_Whitney_Procedure_Fails_as_a_Test_of_Medians.pdf ( 1,15 мегабайт ) Кол-во скачиваний: 136
Прикрепленный файл  Liu_et_al_2017_Modeling_continuous_response_variables_using_ordinal_regression.pdf ( 1,64 мегабайт ) Кол-во скачиваний: 113
 
  Форум: Медицинская статистика · Просмотр сообщения: #27485 · Ответов: 25 · Просмотров: 22923

comisora
Отправлено: 30.06.2022 - 23:43





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


Цитата(salm @ 30.06.2022 - 21:11) *
я вижу следующее, как клиницист: для меня величина коронарного резерва практически всегда определяется числом пораженных сосудов (без направления связи), а возраст, пол, несколько функциональных параметров значимо не влияют на величину этого коронарного резерва.


Величина резерва коронарного кровотока должна быть обратно связана с количеством поражённых сосудов. Если пораженных сосудов нет, то этот резерв максимальный, если все сосуды поражены, то резерв минимальный. Или это не так?
  Форум: Медицинская статистика · Просмотр сообщения: #27423 · Ответов: 20 · Просмотров: 21286

comisora
Отправлено: 22.06.2022 - 10:42





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


Цитата(100$ @ 21.06.2022 - 12:21) *
Мы тут в конце 2017 года развлекались двумя выборками по 4 наблюдения в каждой.
Это тоже было очень весело
Но жизнь не стоит на месте: глядишь, кто-то скоро выборки объемом 2 подтащит...


Цитата(ИНО @ 21.06.2022 - 21:00) *
С каких n начинается статистика?


В приложении книга, в которой обсуждается ситуация n=1.

Ссылка на обсуждение менее экстремальной ситуации
.

Эскизы прикрепленных изображений
Прикрепленное изображение
 

Прикрепленные файлы
Прикрепленный файл  van_de_Schoot_et_al_2020_Small_Sample_Size_Solutions.pdf ( 6,72 мегабайт ) Кол-во скачиваний: 4593
 
  Форум: Медицинская статистика · Просмотр сообщения: #27397 · Ответов: 61 · Просмотров: 33487

comisora
Отправлено: 20.06.2022 - 00:42





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


2 100$

У меня опыта создания пакетов нет (да и не гуру я), хотя наслышан про функцию package.skeleton(). В моей области я редко видел, чтобы делали пакеты к выполненной задаче. Обычно народ колхозит и выкладывает творения в docx файлах. Более продвинутые пользуются ресурсами вроде https://osf.io/. Обученные компьютерной грамоте (т.е. не из моей области) ваяют пакеты так.
На текущем этапе развития форума я вижу возможность создания файлов *.R или постов с кодом, которые бы содержали конкретные функции с источниками, пояснением их работы и т.п. Тот же AtteStat достоин того, чтобы для него была сделана обёртка из R/Python (я размышлял, как его подключить, но квалификации не хватило). Учитывая, что не всё можно найти на CRAN/Github/MATLAB/Python, а на форуме хорошая работа с источниками по различным критериям - есть шанс, что какой-нибудь справочник по статистике будет "оцифрован" и появится пакет forum.disser.ru.misc. Правда тогда придётся писать проверки на тип переменных, на минимальную длину и прочие вещи. И, разумеется, источники. Как пример, было обсуждение на форуме, что критерия история с критерием Крамера-Уэлча не такая простая.
  Форум: Медицинская статистика · Просмотр сообщения: #27376 · Ответов: 61 · Просмотров: 33487

comisora
Отправлено: 18.06.2022 - 23:53





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


2 100$

Конечно. Это я усложнил код, так как можно было легко вызвать функцию напрямую (отредактировал код, вставил функции, в т.ч. ту, которую вызывал через do.call).
do.call применяет функцию однократно к списку, в котором находятся аргументы, что равносильно вызову функции. Функция lapply - это цикл, в ходе которого функция вызывается для каждого элемента списка, просто реализованный в компактом виде. Примеры:

CODE

## Сумма
do.call(sum, list(c(1, 2, 4, 1, 2), na.rm = TRUE)) #Правильно
lapply(c(1, 2, 4, 1, 2), sum) #Неправильно

##Структура операций
L <- list(c(1,2,3), c(4,5,6))
L
[[1]]
[1] 1 2 3

[[2]]
[1] 4 5 6

lapply(L, sum)
[[1]]
[1] 6

[[2]]
[1] 15

list( sum(L[[1]]) , sum(L[[2]]))
[[1]]
[1] 6

[[2]]
[1] 15

do.call(sum, L)
[1] 21
sum(L[[1]], L[[2]])
[1] 21

##Передача аргументов в качестве листа
args <- list(n = 10, mean = 50, sd = 10) #Т-баллы
do.call(rnorm, args) #эквивалентно rnorm(args$n, args$mean, args$sd)

replicate(10, rnorm(10)) #В каждом столбце числа разные
do.call(replicate, list(10, rnorm(10))) #В каждом столбце числа одинаковые

##Вызов функции по её названию
do.call('colnames', list(iris))

##Комбинация списков
rbind(L) #do.call(rbind, list(L))
Reduce('rbind',x = L) #do.call(rbind, L)

Работа функции не самая очевидная, но в некоторых случаях её использование было очень кстати.
  Форум: Медицинская статистика · Просмотр сообщения: #27371 · Ответов: 61 · Просмотров: 33487

comisora
Отправлено: 18.06.2022 - 19:48





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


Я конкретно ответить на Ваш вопрос не могу, но бегло поискал литературу и странички. Возможно, это Вам поможет.

https://stats.stackexchange.com/questions/376273/assumptions-for-lmer-models
https://www.r-bloggers.com/2017/12/linear-mixed-effect-models-in-r/
https://stats.stackexchange.com/questions/77891/checking-assumptions-lmer-lme-mixed-models-in-r
https://www.statalist.org/forums/forum/general-stata-discussion/general/1505099-testing-mixed-model-assumptions
https://vasishth.github.io/Freq_CogSci/checking-model-assumptions.html
Прикрепленные файлы
Прикрепленный файл  mixedeffectsknir.pdf ( 143,03 килобайт ) Кол-во скачиваний: 83
Прикрепленный файл  Schielzeth_et_al_2020_Robustness_of_linear_mixed_effects_models_to_violations_of_distributional.pdf ( 1,67 мегабайт ) Кол-во скачиваний: 84
 
  Форум: Медицинская статистика · Просмотр сообщения: #27369 · Ответов: 2 · Просмотров: 4970

comisora
Отправлено: 18.06.2022 - 19:12





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


2ИНО
Ошибки исправил, результат тоже, спасибо.

Поэтому и написал, что это из серии "продуктивных" методов. Предполагаю, что можно использовать другие максимум и минимумы (например, из данных литературы), zero-one-inflated (gamlss.dist::BEINF) и ещё что-нибудь экзотическое. В моём представлении это не сильно хуже допущения о нормальности распределения или использования уже предложенных подходов к анализу в ходе обсуждения.

Несомненный плюс столь "достоверного" результата - возможность отчитаться за проделанную работу и обосновать необходимость дальнейшего финансирования, чтобы численность обеих выборок довести до 4 crazy.gif .

Касательно других обсуждаемых методов - в приложении статья про применимость теста Стьюдента, модификации Уэлча и ранговых преобразований. Может быть кого-то заинтересует.
Прикрепленные файлы
Прикрепленный файл  Using_the_Students_t_test_with_extremely_small_sample_sizes.pdf ( 333,44 килобайт ) Кол-во скачиваний: 156
 
  Форум: Медицинская статистика · Просмотр сообщения: #27368 · Ответов: 61 · Просмотров: 33487

comisora
Отправлено: 18.06.2022 - 17:49





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


Добрый день.
Предлагаю рассмотреть на жизнеспособность подход из серии "продуктивных методов анализа". Так как данные количественные, минимальное значение концентрации белка 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

  Форум: Медицинская статистика · Просмотр сообщения: #27365 · Ответов: 61 · Просмотров: 33487

comisora
Отправлено: 5.10.2021 - 23:31





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


2 all

Сугубо для поддержания беседы. Вот что нашёл в своей библиотеке:

"The Durbin?Watson test can be used to test the null hypothesis of no residual autocorrelation. More precisely, the null hypothesis of the Durbin?Watson test is that the first p autocorrelation coefficients are all 0, where p can be selected by the user. The p-value for a Durbin?Watson test is not trivial to compute, and different implementations use different computational methods. In the R function durbinWatsonTest() in the car package, p is called max.lag and has a default value of 1. The p-value is computed by durbinWatsonTest() using bootstrapping. The lmtes tpackage of R has another function, dwtest(), that computes the Durbin?Watson test, but only with p = 1. The function dwtest() uses either a normal approximation (default) or an exact algorithm to calculate the p-value." (doi: 10.1007/978-1-4939-2614-5)

Изображение взято из книги https://people.stern.nyu.edu/jsimonof/RegressionHandbook/
Эскизы прикрепленных изображений
Прикрепленное изображение
 
  Форум: Медицинская статистика · Просмотр сообщения: #27008 · Ответов: 7 · Просмотров: 10488

comisora
Отправлено: 5.10.2021 - 17:14





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


2 passant

При крайнем беглом поиске можно встретить упоминания о "нормальном распределении".

https://hal.archives-ouvertes.fr/hal-00642634/document
https://www.tandfonline.com/doi/abs/10.1080...rnalCode=lsta20
https://www.tandfonline.com/doi/abs/10.1080...474939508800333
https://hal.archives-ouvertes.fr/hal-00915951/document
https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2279561
https://www.math.nsysu.edu.tw/~lomn/homepag...nWatsonTest.pdf
https://doi.org/10.2307/1391900
https://doi.org/10.2307/1914122
https://doi.org/10.2307/1914122
https://www.econstor.eu/bitstream/10419/373...010_pid_512.pdf
https://www.inderscience.com/info/inarticle.php?artid=73370
https://www.york.ac.uk/depts/maths/histstat...les/welcome.htm
  Форум: Медицинская статистика · Просмотр сообщения: #27003 · Ответов: 7 · Просмотров: 10488

4 страниц V   1 2 3 > » 

Открытая тема (есть новые ответы)  Открытая тема (есть новые ответы)
Открытая тема (нет новых ответов)  Открытая тема (нет новых ответов)
Горячая тема (есть новые ответы)  Горячая тема (есть новые ответы)
Горячая тема (нет новых ответов)  Горячая тема (нет новых ответов)
Опрос (есть новые голоса)  Опрос (есть новые голоса)
Опрос (нет новых голосов)  Опрос (нет новых голосов)
Закрытая тема  Закрытая тема
Тема перемещена  Тема перемещена