Zero-truncated Poisson regression Опции

[comisora]

Цель предложенных Вами методов неясна...

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

Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским.

Раз в естественном процессе сбора есть нули в результате отбраковки, почему бы не использовать hurdle-модель?

[ИНО]

100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.

comisora, эта гибрид из логичстической модели и ZPT, применяется. наоборот, когда в выборке нулей овер дофига. Наверное, можно расковырять готовую функцию, и оторвать половину гибридной модели, но получим все тот же pospoisson.

Зависимые переменные в моем коде и так одинаковые, столь огорчившие теоретика арифметические фокусы с единицей - это уже часть регрессионной модели. Кстати, собственно к регрессии мы покамест еще перешли, этот шаг запланирован на завтра. Может, с появлением независимых переменных модели наконец-то поведут себя по-разному, и неправильные облажаются. Хотя в теории загвоздка должна была заключаться именно в той части процесса, которая уже проделана, - в моделировании распределения остатков - однако покамест теория что-то плохо контачит с практикой, так что я уже не удивлюсь любым сюрпризам.

Сообщение отредактировал ИНО - 20.11.2022 - 22:44

[comisora]

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)

[100$]

100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.

Ба! Затянувшаяся клоунада сменилась затянувшейся истерикой.

Так, все-таки, кто это там мои "промахи" себе на что-то наматывает? От чьего имени вы тут вещаете? На гельминтоз давно проверялись?

Вам, видимо, невдомек, что гарантированно ошибку I рода удерживают только перестановочные критерии, вычисляющие точное распределение. У всех остальных (с аппроксимациями) она будет либо чумовая, либо чуть менее чем. Но это общее место в статистике. Зачем проверять теорему Пифагора линейкой?
Так что никакого "конфуза".

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

И, кстати, что такого нерабочего в моем коде?
Что вы разумеете под "работоспособностью"?
Что это у вас там так подгорело?
Здесь же не форум программистов: здесь бесконечное число юзеров способны написать что угодно на чем угодно быстрее и лучше меня. И чо?

Зачем оно нужно вообще и мне в частности? Интересно сначала решить задачу, а потом кропать программы.

...И только это не имеющее аналогов в мире изобретение отличлось!
...(а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков)...

А вот идиотничать я вам не позволю.


Сообщение отредактировал 100$ - 21.11.2022 - 00:45

[ИНО]

100$, с каждым новым постом Вы все менее походите на интеллектуала, и все более - на недалекого гопника из подворотни. Иногда лучше промолчать - за умного сойдешь (с)народная мудрость.

Что касается KWTeest(), приведу Вашу же цитату:

Я тут чуток помонтекарлил для равных выборочных объемов (n1=n2) от 3 до 100 способность CWTest() удерживать номинальный уровень значимости (alpha=,05).
Получается, что сколько-нибудь осмысленное его применение возможно только при n1=n2>60, п.ч. сходимость к предельному распределению очень медленная. Ребята, это же ужасно.

Пусть каждый решает для себя, нужен ли ему такой ущербный критерий при наличие давно известных намного лучших альтернатив. Как по мне, это просто мусор. Но не надо втирать, будто иные критерии ведут себя схожим образом, ибо это наглое вранье, CWTest() - воистину аналоговнет. Ладно, предположим, что великий 100$ не стал создавать пакет имени себя не потому, что осознал ущербность этого своего критерия, а чисто в силу неспособности освоить техническую сторону процесса (вернее "сомнения в собственных силах", ведь, как известно, истинный мэтр может все, но имеет право иногда подвергать эту всемогучесть легкому сомнению ). Но почему же он в этой теме не написал ни строчки кода R в подтверждение своих обвинительных заявлений? Неужто и в этой способности через полгода возникли сомнения? Или просто лень одолела?

Сообщение отредактировал ИНО - 21.11.2022 - 10:44

[ИНО]

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

Пожалуйста, проиллюстрируйте примером. А то я сосем не понимаю, что вы имеете в виду.

Я пока не придумал как адаптировать Ваш код для доказательства/опровержения утверждения, что Ваш генератор и rpospois() даёт одно и тоже.

Это просто:

# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
x<-rep(NA, 200)
y<-rep(NA, 200)
for(i in 1:200){
lambda_i<-lambda[i]
x[i]<-mean(rpospois(1000, lambda_i))
sim<-rpois(1000000, lambda_i)
sim<-sim[sim>0]
y[i]<-mean(sample(sim, size=1000))
}

plot(x, y)
abline(0, 1, col="red")

Здесь x - выдача rpospois(), а y - моего генератора. Фактически он делает ровно то же, но способом, который в более явном виде имитирует реальный процесс порождения "позитивнопуассоновских" данных. Конечно, полного соответствия не будет, т. к. у меня датчик псевдослучайных чисел работает в двух разных функциях и я не знаю как в такой ситуации корректно задать одинаковое семя обоим генераторам. Но и полученной диаграммы вполне достаточно для осознания их эквивалентности. Поскольку мы имеем дело с однопараметрическим распределением, достаточно оценки одного параметра, но для пущей убедительности можно проделать то же самое с любой другой, скажем заменить mean() на sd() - картина не изменится.

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

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

[Диагностик]

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

Для какого лямбда?

[ИНО]

Для любого в интервале от 0,001 до 2 (на больших не проверял, ибо неинтересно).

[ИНО]

Итак, переходим к собственно регрессии. Пока что пусть будет один бинарный предиктор. В этот раз для генерации данных взяты лямбды, похожие на оцененные в одной из моих реальных задач, и объем выборки того же порядка. В этот раз мы будем сравнивать модели на 1000 разных реализаций распределений с одними и теми же параметрами. Для разнообразия вместо MAE взята MSE, но и с MAE картина ровно та же - модели практически не различаются по точности прогнозов. Традиционно под номером 1 идет моя оригинальная шифт-модель, под номером 2 - ZTP, под номером 3 - обычная пуассоновская

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg<-function(pr)
{
resp1<-rpospois(50, 0.6)
resp2<-rpospois(50, 1.35)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
MSE1<-mean((prediction1-df$resp)^2)
mod2<-vglm(resp~pr, data=df, family="pospoisson")
prediction2<-predict(mod2, type="response")
MSE2<-mean((prediction2-df$resp)^2)
mod3<-glm(resp~pr, data=df, family="poisson")
prediction3<-predict(mod3, type="response")
MSE3<-mean((prediction1-df$resp)^2)
return(c(MSE1, MSE2, MSE3))
}
x<-replicate(1000,  TestPoisReg(pr=pr))
MSE1<-x[1,]
MSE2<-x[2,]
MSE3<-x[3,]

boxplot(MSE1, MSE2, MSE3)

Сообщение отредактировал ИНО - 21.11.2022 - 16:38

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

 

[ИНО]

А вот в p разница существенная:

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg2<-function(pr)
{
resp1<-rpospois(50, 0.6)
resp2<-rpospois(50, 1.35)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
p1<-coef(summary(mod1))[2,4]
mod2<-vglm(resp~pr, data=df, family="pospoisson")
p2<-coef(summary(mod2))[2,4]
mod3<-glm(resp~pr, data=df, family="poisson")
p3<-coef(summary(mod3))[2,4]
return(c(p1, p2, p3))
}
x<-replicate(1000,  TestPoisReg2(pr=pr))
p1<-x[1,]
p2<-x[2,]
p3<-x[3,]

boxplot(p1, p2, p3)

Ну, третью модель (обычный Пуассон) можно отправлять сразу в мусорное ведро, потому что она пытается предсказать нули, которых нет, построенные по ней ДИ будут невозможными, соответственно и p тоже должно быть неправильное. А вот как сделать выбор между ?1 и ?2, понять какая из них больше врет в данном случае, я покамест не сообразил. Буду благодарен за подсказку. И все ж странно, почему при таких серьезных различиях в объясняющей составляющей (которая была ожидаема), в составляющей прогностической все три модели оказались практически одинаковыми.

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

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

 

[ИНО]

За отсутствием альтернативных предложений решил попробовать ответить на последний поставленный мною вопрос самостоятельно, а именно при помощи моего оригинального метода, предложенного здесь: http://forum.disser.ru/index.php?showtopic...ost&p=27378

Теорема имени меня : чем ближе распределение p-значений какого-либо статистического критерия при справедливости нулевой гипотезы к равномерному, тем ближе номинальный уровень значbмости к реальному. При этом смещение моды распределения влево от значения 0,5 сигнализирует о завышенной ошибке первого рода, а вправо - о завышенной ошибке второго рода. Доказывает (или опровергает) пусть $100 - у него голова большая он шибко умный

Исследуем при помощи предложенного метода три вышеупомянутые модели:

Код

pr1<-rep(1, 50)
pr2<-rep(2, 50)
pr<-as.factor(c(pr1, pr2))

TestPoisReg2<-function(pr)
{
resp1<-rpospois(50, 1)
resp2<-rpospois(50, 1)
resp<-c(resp1, resp2)
df<-data.frame(pr, resp)
df$modresp<-df$resp-1
mod1<-glm(modresp~pr, data=df, family="poisson")
prediction1<-predict(mod1, type="response")+1
p1<-coef(summary(mod1))[2,4]
mod2<-vglm(resp~pr, data=df, family="pospoisson")
p2<-coef(summary(mod2))[2,4]
mod3<-glm(resp~pr, data=df, family="poisson")
p3<-coef(summary(mod3))[2,4]
return(c(p1, p2, p3))
}
x<-replicate(1000,  TestPoisReg2(pr=pr))
p1<-x[1,]
p2<-x[2,]
p3<-x[3,]

hist(p1)
abline(v=0.5, col="red")

hist(p2)
abline(v=0.5, col="red")

hist(p3)
abline(v=0.5, col="red")

Ожидаемо у модели ?3 (обычный Пуассон) с достигаемым уровнем значимости для эффекта предиктора все плохо - он сильно завышен. В этот раз, когда предиктор является мусорным, оно даже хорошо, однако и для по-настоящему влиятельного предиктора он тоже будет завышен, что уже совсем нехорошо. Так что использовать обычную регрессию Пуассона для данных с отсутствующими нулями с целью оценки влияния независимых переменных однозначно недопустимо (хотя с задачей прогноза в такой ситуации она вполне справляется, что весьма странно, но все результаты моделирования покамест говорят именно об этом). А вот с моделями ?1 (моя шифтонутая) и ?2 (ZTP) все не так однозначно: в зависимости от n и лямбда, задаваемых генератору, распределение p красивее то у одной, то у другой из них, но разница всегда небольшая. В случае, представленном на прикрепленных гистограммах лидирует расово правильная усеченная модель, а моя дает малость заниженное p, но в некоторых иных, мною смоделированных, напротив, моя модель имеет почти идеальное распределение, а расово правильная демонстрирует признаки небольшого завышения. Может поэкспериментировать самостоятельно.

Сообщение отредактировал ИНО - 23.11.2022 - 17:52

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

 

[100$]

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

[ИНО]

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

Почему Вы упомянуты? Так Вы тут единственный (оставшийся), кто полагает (ИМХО необоснованно), что в совершенстве понимает теорию математической статистики. Так что можете считать это обращение знаком особого уважения к вашему интеллекту. Но нет, так нет.

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





Сообщение отредактировал ИНО - 23.11.2022 - 20:31

[Игорь]

Врачи-аспиранты, для которых создан форум, ознакомившись с показанными выше результатами исследований, могут составить впечатление о непостижимости темы для медика или биолога. Сугубо личное мнение, просто такие исследования не нужно обсуждать здесь. А нужно [в рассматриваемом случае] представлять совет в виде:
1. Есть данные, являющиеся результатами практических измерений в эксперименте.
2. Смотрим, какая статистическая модель (модели) подходит (подходят) для их описания (с точностью до коэффициентов).
3. Идентифицируем модель - вычисляем оптимальные значения коэффициентов, при необходимости значимость параметров, включенных в модель.
4. Проверяем соответствие модели экспериментальным данным. Часто применяется критерий AIC. Если подходит несколько моделей, выбираем лучшую.
5. Используем построенную модель для прогнозирования и в публикациях. Всё красиво, все довольны.

Есть еще одна идея. Ранее, когда компании-разработчики статистического программного обеспечения (ПО) рассылали демонстрационные материалы (диски, брошюры) по почте, получил от StaXact очень интересные картонные карточки формата меньше A4 (наверное, B5) с голубой каемкой. На лицевой стороне карточки была описана проблема - реальный случай из медицинской практики, включая ссылки. На оборотной стороне - результаты статистического исследования с интерпретацией результатов. Идея понравилась - найди пример и сделай так же.

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

Не нужно пугать аспирантов.

[Диагностик]

Для любого в интервале от 0,001 до 2 (на больших не проверял, ибо неинтересно).

Распределения близки, можно выбирать любое. На графиках теоретические вероятности.



Сообщение отредактировал Диагностик - 24.11.2022 - 10:42