Широко используется буржуями для предсказания счетной зависимой переменной, если она в силу дизайна эксперимента не может принимать нулевое значение (ну, например количество прыщей посчитали только на тех больных, которые обратились с жалобами на прыщи). Интерпретация параметров таких моделей заковыриста:
https://stats.stackexchange.com/questions/4...ts?noredirect=1Вопрос: а почему с аналогичной целью нельзя применить обычную регрессию Пуассона, просто предварительно вычтя из всех значений зависимой переменной единицу? Пусть такая модель предсказывает не общее количество счетных единиц, а количество единиц, добавленных к одной обязательно имеющейся. Потом, если надо, добавляем к предсказанной величине единичку, и дело в шляпе! Но, судя по тому, как извращается народ именно с Zero-truncated Poisson regression, и ищет ее программные реализации, очевидно, предложенный мною альтернативный подход неправомерен. Но сколько я не ломал голову, так и не понял, почему. Прошу более головастых подсказать.
И вдогонку еще один маленький вопросик: а как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?
Диагностик
6.11.2022 - 10:40
Цитата(ИНО @ 6.11.2022 - 13:07)

как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?
Ограниченное (усечённое)распределение Пуассона.
По регрессии Пуассона есть хорошая книга. Свободное скачивание
http://qed.econ.queensu.ca/ETM/. С. 475. 2021 года, но по рассматриваемой теме там ничего нет.
Цитата(ИНО @ 6.11.2022 - 08:07)

Интерпретация параметров таких моделей заковыриста:
https://stats.stackexchange.com/questions/4...ts?noredirect=1... как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?
Яндекс.Браузер перевел заглавие ссылки как "Интерпретация коэффициентов регрессии Пуассона с нулевым усечением". По-моему, звучит неплохо.
Для Пуассона с нулевым усечением можно построить аналогичную вычислительную схему (вывести все соотношения, как для обычного Пуассона в показанной выше ссылке)... и составить ПО. Непонятно, в чем сложность?
Сложность в том, что, вероятно, этот подход неправомерен. Я пробовал - для того же набора данных данных получаются существенно другие p, значит, очевидно, есть какой-то теоретический запрет для этого моего ноу-хау

Понять бы, какой.
А гугль переводить каждый раз по-разному, но все время - как-то криво. "С нулвевым усичением" - не годится, потому как нулевое учечение = отсутствие усечения. А тут фишка в том, что усечение совсем не нулевое, но усекают нули.
Цитата
Ограниченное (усечённое)распределение Пуассона.
Это-то понятно, но важно отразить, что усечено оно не где-нибудь, а именно в нуле.
Диагностик
7.11.2022 - 00:50
Цитата(ИНО @ 7.11.2022 - 04:47)

важно отразить, что усечено оно не где-нибудь, а именно в нуле.
Цитата
В общем случае можно рассматривать лево-усечённое распределение (запрещены первые значения ), право-усечённое распределение (запрещены значения, большие некоторого уровня ) и двусторонне-усечённое распределение.
Это как бы самоочевидно, что усеченный ноль находится на левом конце распределения Пуассона. Но какой вывод относительно темы обсуждения из этого следует по-вашему?
Диагностик
7.11.2022 - 11:00
Я отвечал на вопрос.
Цитата(ИНО @ 6.11.2022 - 13:07)

И вдогонку еще один маленький вопросик: а как сабж будет грамотно обозвать по-русски? Регрессия Пуассона с усеченным нулем?
Ну и? Каким образом данная Вами цитата о том, что распределения бывают усеченные слева, справа и спереди дает ответик на мой конкретный вопросик? Я прекрасно понимаю, что означает "zero-truncated", но не знаю, как это корректно перевести на русский.
Диагностик
7.11.2022 - 13:17
Усеченное до нуля распределение означает что нуль исключается.
Диагностик
8.11.2022 - 15:17
Чем богат.
Хотелось бы высказаться не в ответ, но по теме.
У Johnson с соавт. на 188 описано zero-truncated распределение Пуассона (там же и другие варианты упоминаются). Можно сделать, как у Davidson с соавт. (ссылку см. выше) - рассмотреть схему вычисления zero-truncated регрессии Пуассона: записать ФМП, прологарифмировать, дифференцировать, получить систему нелинейных уравнений, аналитически записать градиент и гессиан, решить систему по Ньютону-Рафсону, оценить качество модели и значимость коэффициентов. Можно пойти и дальше - составить программу на языке Си типа как это реализовано в AtteStat для регрессии Пуассона.
Вопрос, а нужно ли практически. Обосную. Есть зарубежные работы, показывающие принципиально лучшее качество zero-truncated регрессии Пуассона по сравнению с регрессией Пуассона для некоторых наборов данных. Но ситуация похожа на использование ТМФ и критерия Барнарда в анализе таблиц сопряженности - может, и лучше обосновано, но никто не интересуется. Ну еще есть пара примеров из практики, когда эффективные алгоритмы интересуют разве что диссертантов для объема работы. Отсутствие русскоязычных работ (даже официального или хотя бы общепринятого перевода термина не существует, кстати, как и в Барнарде) и скромное количество зарубежных доказывают, что по крайней мере здесь эффективные прикладные алгоритмы анализа данных просто игнорируются научным сообществом, для которого эти алгоритмы по идее должны быть интересны: медики, биологи, экономисты.
Однако так просто оставлять тему нельзя. После Нового года, как обычно, несколько дней нечего делать будет, кроме очистки снега. Сяду на даче в библиотеке, сделаю все сказанное выше хотя бы ради красоты формул и законченности темы. Хоть попрактиковаться. Все-равно преподавание ярких прикладных дисциплин никому здесь не нужно.
Так что спасибо автору темы за идею.
Цитата(Игорь @ 9.11.2022 - 16:49)

У Johnson с соавт. на 188 описано
Помню, нам преподаватель философии для разрядки обстановки в промежутках между надиктовкой лекции приводил пример того, как давали литературные ссылки в старинных трактатах: "У Платона находим..." - с комментарием от себя: "У какого Платона? Где находим?". По-видимому, Вы приверженец именно той древней школы цитирования
Как-то в библиотеке заказал редкую статью, типа года 60-го, в которой рассматривалась эквивалентность ранговых критериев. Надеялся, что запросят по МБА. Через какое-то время пригласили придти и выдали листочек с аннотацией. На вопрос, а где же сам текст, ответили, что я оказался первым среди диссертантов, кого не устроила аннотация, а потребовался весь источник. Кстати, о статье - мне потом ее любезно прислал автор, я с ней ознакомился и подарил библиотеке.
Так что - нет, не стиль цитирования, а жизненный опыт. Неужели еще кто-то ходит в библиотеку, сверяется с первоисточником, да просто по ссылкам?
Печальный у Вас жизненный опыт, но не стоит обобщать. По ссылкам ходят и читают, если источник доступен. Если недоступен, ходят, чтобы убедиться в этом. Но первое необходимое условие для этого - наличие этих-самых ссылок в удобоваримом виде. Интересно, сколько есть в природе книг, написанных некоим Ивановым Джонсоном с соавторами с количеством страниц более 188? Прикажете все проверить? И такие с позволения сказать ссылки я встречаю в Ваших постах уже далеко не первый раз.
Цитата(ИНО @ 6.11.2022 - 08:07)

Вопрос: а почему с аналогичной целью нельзя применить обычную регрессию Пуассона, просто предварительно вычтя из всех значений зависимой переменной единицу? Пусть такая модель предсказывает не общее количество счетных единиц, а количество единиц, добавленных к одной обязательно имеющейся. Потом, если надо, добавляем к предсказанной величине единичку, и дело в шляпе! Но, судя по тому, как извращается народ именно с Zero-truncated Poisson regression, и ищет ее программные реализации, очевидно, предложенный мною альтернативный подход неправомерен. Но сколько я не ломал голову, так и не понял, почему. Прошу более головастых подсказать.
Надо думать патамушта моделирование суммы X1+X2+,...,+Xn - An независимых одинаково распределенных величин требует перехода к т.н.
обобщенному распределению Пуассона (ОПР).
Число n слагаемых в этой сумме является Пуассоновой величиной с параметром lambda.
Как любое безгранично делимое распределение (а ОПР является таковым) каждое обобщенное распределение Пуассона является пределом других ОПР, в т.ч. и "сдвинутых" на некоторую неслучайную алгебраическую величину An. (т.е. характеристическая ф-ция ОПР допускает такой "сдвиг").
Разумеется, проверка сказанного с помощью "обычной регрессии Пуассона" приведет к появлению жалобы типа
Цитата
Я пробовал - для того же набора данных данных получаются существенно другие p, значит, очевидно, есть какой-то теоретический запрет для этого моего ноу-хау.
Понять бы, какой.
с единичной вероятностью.
Благодарю за первый ответ качественный по существу проблемы! Но непонятно: почему мы должны подходить к величине X1 (т. е. первому прыщу из описанного примера) как к случайной, а не как к постоянной? Ведь этот первый прыщ есть примененное условие включения пациента в выборку. Не вижу причин не принимать X1 константой, равной 1, и моделировать только X2.
Цитата(ИНО @ 13.11.2022 - 09:55)

...Но непонятно: почему мы должны подходить к величине X1 (т. е. первому прыщу из описанного примера) как к случайной, а не как к постоянной?
Патамушта в теорвере не определено понятие "постоянная величина".
Ergo ваша единица является случайной. Но, чтобы она постоянно присутствовала в каждом испытании, ее - хочешь, не хочешь - надо объявлять достоверным событием. И всю вероятностную массу распределения сбрасывать в точку 1.
Тогда ваша вероятностная модель получаецца самую малость странной: сначала с единичной вероятностью случается (достоверное) событие - некая "постоянная" величина приняла значение, равное
1, а затем к ней еще с некоторой (пуассоновой) вероятностью выпал 0. Так получается пуассонова 1. Для получения пуассоновой двойки надо
1+0+1, etc. Но, увы вам! Вероятность такого "составного" события равна нулю...
Цитата
Ведь этот первый прыщ есть примененное условие включения пациента в выборку
Ну, если каждый раз отсутствие признака считать пуассоновым нулем... (типа "зуб даю, этот чувак болен прыщами, хотя на морде у него их ровно 0"), то тогда ваша (неслучайная) выборка условна относительно нуля. Это как раз и есть хорошо нам известное zero-truncated Poisson dist. в чистом виде.
Вообще, вся эта комедия объясняется тем, что описанным вами способом никак нельзя модифицировать
носитель распределения (0,1,2,...).
Поэтому ребята вынуждены рассматривать ненужные им нули в рамках условного относительно нуля распределения ("нуль-усеченная модель"), трактовать их как число неудач до первого успеха (отрицательное биномиальное распределение и "negative binomial model"), вводить zero-inflated Poisson regression, вспоминать про изученные Кацем в далеком 1963 году Zero-modified Poisson dist. и нет этому конца.
P.S. А вот, a propos, как бы вы перевели словосочетание Zero-inflated Poisson regression? Настоящий пробный камень начинающего переводчика. Или кошмар.
Цитата
(типа "зуб даю, этот чувак болен прыщами, хотя на морде у него их ровно 0")
Вот-такие не доходят до врача, их еще в дверях заворачивают. Сколько их в популяции, не интересно, да и по количеству прышей у прыщавых едва ли возможно сделать вывод о частоте встречаемости непрыщавых.
Цитата
Для получения пуассоновой двойки надо 1+0+1, etc. Но, увы вам! Вероятность такого "составного" события равна нулю...
Ту что-то совсем не понял. Двойка у нас никакая не пуассоновская, она получается сложением постоянной 1 с пуассоновской 1. Предположением, что мы не заворачиваем никого в дверях, но за первый "прыщик" на лице пациента считаем... его нос. Он есть у всех и каждого, в отличие от второго и последующего прыщиков, которые действительно прыщики. Врач - инопланетянин, причем очень тупой, и не умеет отличать нос от прыщиков, а может только посчитать общее количество выпуклостей на лице. Вполне себе попадутся ему и с одной, и с двумя и с тремя и т. д., никаких нулевых вероятностей, за исключением собственно нулевых знаний, которых гарантировано не будет (ну, разве то герой Гоголя на прием придет

). По-моему вполне себе адекватная вероятностная модель выходит.
Цитата(ИНО @ 14.11.2022 - 00:52)

Ту что-то совсем не понял. Двойка у нас никакая не пуассоновская, она получается сложением постоянной 1 с пуассоновской 1.
Ну вот, началось в колхозе утро: я ему про вероятности, он мне - про школьную алгебру. Так, чего доброго, и до таблицы умножения дойдет.
А какая же она? Детерминированная, чтоле?
Двойка - самая что ни на есть пуассонова:
- во-первых, несомненно случайная, потому как зависит от случайного второго слагаемого;
- во-вторых, много ли вы знаете дискретных распределений, носителем которых является {0,1,2,3...}?
Если она не пуассоновская, то какая? Из распределения Скеллама, что ли?
А, может быть, Накагами? Трейси-Видома? Уишарта? Уточните, когда не лень.
P.S. Вы распределение-то сгенерируйте: result<-rpois(1000,.685). При параметре распределения Lambda<1 0 - самое вероятное значение. Как от него избавляться? Распределение Пуассона - это же все-таки закон редких событий.
Что-то я Вас с каждым постом все меньше понимаю. Зачем избавляться от пуассоновского от 0? Он ничуть не мешает, даже наоборот необходим! 1(константа) + 0 (пуассоновский) =1 - как раз самое часто значение в моих подсчетах "прыщиков". И частоты всех последующих значений (1+1=2, 1+2=3...) вполне согласуются с гипотезой о добавлении постоянной единицы к случайной пуассоновской величине. Этот первый прыщик - непременное условие для осмотра больного, достоверное событие, а вот уже количество дополнительных прыщиков - случайная величина. Наезд на алгебру тоже не понял. Прибавление постоянного свободного члена в регрессионной модели Вас тоже смущает?
Цитата(ИНО @ 14.11.2022 - 10:06)

Что-то я Вас с каждым постом все меньше понимаю. Зачем избавляться от пуассоновского от 0? Он ничуть не мешает, даже наоборот необходим! 1(константа) + 0 (пуассоновский) =1 - как раз самое часто значение в моих подсчетах "прыщиков". И частоты всех последующих значений (1+1=2, 1+2=3...) вполне согласуются с гипотезой о добавлении постоянной единицы к случайной пуассоновской величине. Этот первый прыщик - непременное условие для осмотра больного, достоверное событие, а вот уже количество дополнительных прыщиков - случайная величина.
Ну, а что тут непонятного? Вы сконструировали некоторое дискретное распределение. Для дискретного распределения можно явно указать вероятность любого события. Вот и укажите вероятность появления суммы вида
1(константа) + 0 (пуассоновский) =1.
Пока этого не сделано, я буду считать это школьной математикой.
Цитата
Наезд на алгебру тоже не понял. Прибавление постоянного свободного члена в регрессионной модели Вас тоже смущает?
Нет. Потому что в регрессионной модели константа - это такой же (оцениваемый) регрессионный коэффициент, и так же точно является случайной величиной. Для нее так же точно строятся доверительные интервалы, тестируется отличие от нуля.
Так что смущает меня только то, что вы этого не знаете.
Но это не приговор. Люди с этим живут и неплохо себя чувствуют. Вы, например.
Что такое "школьная математика" я не в курсе, равно как и о причине Вашего уничижительного отношения к алгебре.
Цитата
Вот и укажите вероятность появления суммы вида 1(константа) + 0 (пуассоновский) =1
Та же, что и вероятность 0 в обычной пуассоновской регрессионной модели, построенной по той же выборке, из каждой варианты которой вычтена 1. Формулой это не напишу, увы.
Цитата
Нет. Потому что в регрессионной модели константа - это такой же (оцениваемый) регрессионный коэффициент, и так же точно является случайной величиной
Может быть оценен, а может быть и взят
с неба из теоретических предпосылок, и коэффициенты, кстати, тоже, так что оцениваться будет только ошибка. И это все равно будет статистическая модель.
И да тут все ж форум о прикладной статистике, а не о математической (это разные дисциплины, хотя и связанные). О последней пусть д. ф.-м. н. семиэтажными формулами дискутируют, теоретически обосновывая очередной элегантный чудо-метод имени себя (то, что на практике он может толком не работать, потому что "забыли про овраги" - это уже дело десятое

) . Так что объяснение "на пальцах" для простых смертных типа меня всяко приветствуется.
Цитата(ИНО @ 15.11.2022 - 09:07)

Та же, что и вероятность 0 в обычной пуассоновской регрессионной модели, построенной по той же выборке, из каждой варианты которой вычтена 1. Формулой это не напишу, увы.
Ну, если на пальцах, то вы сконструировали некое условное относительно неслучайной величины (единицы) распределение и утверждаете, что в нем условная вероятность равна безусловной.
Хотя это ниоткуда не следует.
Вот, например, если мы захотим сдвинуть стандартное нормальное распределение, мы же не будем городить такую словесную конструкцию "а вот мы дождемся с вероятностью 1 единицу, а потом посмотрим, как реализуется ст.норм.расп." Мы этот сдвиг просто добавим к параметру распределения. И оно будет симметрично теперь уже относительно 1.
Вот бы и с Пуассоном также точно.
Как обычно, я попробовал развеять теоретический сумрак моделированием. Решил: какая из двух моделей будет иметь более точный прогноз на смоделированных данных, той буду пользоваться на практике. Все ж я биолог, а не математик. Но возник конфуз: прогнозы у сабжа и "моего безграмотного ноу-хау" и их ошибки оказались практически неотличимыми. Более того, при лямбда<1 (а это больше похоже на реальные данные) моя модель работает даже чуточку лучше! Хотя возможно тут дело не непосредственно в лямбде, а в том, что мой алгоритм моделирования усекает объемы выборок, сгенерированных при малых лямбда, но тогда это означает, что расово правильная модель хуже моей ведет себя на малых выборках, так что хрен редьки не слаще.
Уменьшение n на описанную картину не влияет (правда, при n меньше 100 и очень малых лямбда, алгоритм сталкивается с невозможностью оценки модели из-за полного отсутствия ненулевых значений и стопорится, так что в этом диапазоне условий исследование провести не удалось. Наверное, по-хорошему надо применять вместо цикла сопли-функции

или как их там, а я в этом не силен. Но для понимания проблемы и имеющееся сойдет.
Итак, прошу указать на мою ошибку, не в теории, которую как тут выяснилось, я постичь не в силах, а именно в моделировании!
Код
# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
simmat<-matrix(NA, 200, 1000)
for(i in 1:200){
lambda_i<-lambda[i]
simmat[i,]<-rpois(1000, lambda_i)}
# предложенная мной раскритикованная выше модель:
prediction<-rep(NA, 200)
MAE<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df<-subset(df, x>0)
df$modx<-df$x-1
mod<-glm(modx~1, data=df, family="poisson")
pred<-predict(mod, type="response")+1
prediction[i]<-pred
MAE[i]<-mean(abs(pred-df$x))}
# модель c расово правильно усеченным нулем:
library(VGAM)
prediction2<-rep(NA, 200)
MAE2<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df<-subset(df, x>0)
mod<-vglm(x~1, data=df, family="pospoisson")
pred<-predict(mod, type="response")
prediction2[i]<-pred
MAE2[i]<-mean(abs(pred-df$x))}
# пргнозы обеих моделей (моя - красненькая):
plot(lambda, prediction, type="n")
lines(lambda, prediction, col="red")
lines(lambda, prediction2, col="blue")
# средние абсолютные ошибки обеих моделей:
plot(lambda, MAE, type="n")
lines(lambda, MAE, col="red")
lines(lambda, MAE2, col="blue")
# резульаты их сравнения при различных лямбда:
plot(lambda, MAE>MAE2)
Диагностик
18.11.2022 - 11:15
Цитата(ИНО @ 18.11.2022 - 14:57)

при лямбда<1
А что такое лямбда?
Параметр распределения Пуассона. В моем случае имеется в виду распределение, из которого были сгенерированы синтетические данные. Следует отметить, что варианты с нулевым значением исключены из последующего анализа. По идее, именно исходное распределение, их включающее, пытается реконструировать zero-truncated Poisson regression. Мой же подход работает по-другому: в нем оценивается обычное распределение Пуассона, по выборочным значениям, из которых была вычтена единица. У этого распределения уже иная лямбда, которая здесь не приводится. Да, для простоты в обоих случаях построена регрессия к 1, т, е. просто произведена оценка параметров распределения. Регрессионный анализ как таковой - это уже следующий шаг.
Диагностик
18.11.2022 - 13:25
Цитата(ИНО @ 18.11.2022 - 16:35)

Параметр распределения Пуассона.
То есть матожидание. А как здесь ТЕХ использовать?
Наконец-то наш естествоиспытатель сам признался, что переливать из пустого в порожнее ему сподручнее безо всякой теории )
Однако теория гласит, что ежели моделировать Пуассоном счетные данные, то с необходимостью надо принять, что эти данные измерены в абсолютной шкале. А в ней, как это знает не только взрослый, но даже карапуз, допустимы только тождественные преобразования, т.е. преобразования вида f(x)=x, переводящие каждый элемент шкалы в самого себя. Ни тебе сложения, ни умножения, ни возведения в степень, ни извлечение квадратного корня. Как из таких данных можно вычесть 1 непонятно.
Далее. В Пуассоново распределение мутирует биномиальное распределение при малой вероятности успеха и неограниченном возрастании N. А что есть биномиальное распределение?
Вот выпало в результате эксперимента из 100 бросков монеты 49 орлов и 51 решка. Как из этих данных вычесть единицу? Это ж фальсификацыя ).
А модель... Ну, что модель: у нее же нет "защиты от дурака". Если она не столкнулась с вычислительными проблемами, то она чего-то там посчитает. Ну и что?
comisora
18.11.2022 - 21:57
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? - обсуждение того, как меняется распределение при добавлении константы.
comisora
19.11.2022 - 13:58
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
Рекомендую к
просмотру презентацию, в которой объяснены некоторые идеи.
100$, я уже видел Ваши 100500 теоретических аргументов в пользу того, что шмель летать моя модель работать не может, спасибо. 100501-й и 100502-й уже погоды не сделают. Теперь я отошел от изначальной формулировки вопроса и интересуюсь иным: почему результаты моего моделирования все же показывает, что шмель летает моя модель работает не хуже (а местами даже самую малость лучше лучше), чем аэроплан zero-truncated Пуассон. Прошу заметить, что смоделированные мною данные являются классическим случаем zero-truncated Пуассона, который обязан адекватно обрабатываться соотвествующей функцией из пакета VGAM (если, конечно, ее не дурак писал, но, вроде, не жалуются) и поэтому Ваши отсылки к условиям применимости идут лесом. Имеется факт: на синтетическом примере моя "неправильная" модель справляется не хуже общепризнанной для смоделированной ситуации расово правильной. Но это лишь в том случае, если я не ошибся в моделировании. в чем совсем не уверен. Поэтому прошу проверить мой код и указать на ошибку, если таковая имеется. Теоретико-философских же аргументов довольно, я их услышал, а далее уже становится непонятно, кто из нас на самом деле тут из пустого в порожнее переливает. Вот я, например, код написал и даже анимацию сотряпал, а Вы?
comisоra, в очередной раз спасибо за полезные ссылки! Судя по ним, то что я здесь предложил, - это именно shifted Poisson distribution. И я не первый, кто задавался вопросом о возможности его применения для "положительной регрессии":
https://stats.stackexchange.com/questions/2...hifting-vs-trun. Но, к сожалению, в этой теме обсуждение сместилось в другую область и заданные автором вопросы так и остались без ответа. Увы, больше ничего полезного по запросу "shifted Poisson distribution" найти не смог. Отличия zero-truncated от обычного Poisson distribution я, в принципе представляю, и это не то, о мне хотелось бы сейчас почитать.
Итак, решил я еще немного "из пустого в порожнее попереливать", а именно немного модифицировать код из своего позапрошлого поста, так чтобы на каждом шаге модели строились по одинаковому числу "наблюдений". В данном примере оно установлено на 10, что приближает условия учений к боевым

. За одно пришлось увеличить объем изначальной пуассоновской выборки до 2000, так как из 1000 при малых лямбда было, что не набиралось и десятка единиц. Таким образом мы полностью исключили влияние эффекта неравновеликости выборок на поведение моделей и можем утверждать, что моя неправильная модель таки обходит "правильную" именно при малых лямбда, да и при больших не отстает. В целом моя модель почти в три раза чаще имеет меньшую среднюю абсолютную погрешность, чем "правильная".
Код
# генерация синтетических данных:
lambda<-seq(0.01, 2, 0.01)
simmat<-matrix(NA, 200, 10)
for(i in 1:200){
lambda_i<-lambda[i]
sim<-rpois(2000, lambda_i)
sim<-sim[sim>0]
simmat[i,]<-sample(sim, size=10)
}
# предложенная мной раскритикованная выше модель:
prediction<-rep(NA, 200)
MAE<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
df$modx<-df$x-1
mod<-glm(modx~1, data=df, family="poisson")
pred<-predict(mod, type="response")+1
prediction[i]<-pred
MAE[i]<-mean(abs(pred-df$x))}
# модель c расово правильно усеченным нулем:
library(VGAM)
prediction2<-rep(NA, 200)
MAE2<-rep(NA, 200)
for(i in 1:200){
x<-simmat[i,]
df<-as.data.frame(x)
mod<-vglm(x~1, data=df, family="pospoisson")
pred<-predict(mod, type="response")
prediction2[i]<-pred
MAE2[i]<-mean(abs(pred-df$x))}
y<-MAE>MAE2
length(y[y==TRUE]) # число случаев, в которых расово правильная модель оказалась лучше моей
length(y[y==FALSE] ) # число случаи, в которых моя модель оказалась лучше расово правильной
Диагностик
20.11.2022 - 10:33
Цитата(ИНО @ 19.11.2022 - 21:33)

можем утверждать, что моя неправильная модель таки обходит "правильную" именно при малых лямбда, да и при больших не отстает. В целом моя модель почти в три раза чаще имеет меньшую среднюю абсолютную погрешность, чем "правильная".
С чем сравнивались обе модели?
Диагностик
20.11.2022 - 11:27
И какая из них считается правильной?
Этот вопрос лучше адресовать $100, поскольку мне , в отличие от него, неведомы критерии истинно правильности моделей. Могу лишь сказать, что по результатам тестирования на синтетических данных у моей модели оказалась чуть ниже ошибка предсказания. Но это очень малое отличие, на практике совершенно несущественное. Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать именно то распределение, данные из которого я генерирвал
comisora
20.11.2022 - 12:47
Цитата(ИНО @ 20.11.2022 - 12:07)

Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать
именно то распределение, данные из которого я генерирвал 
Код
sim<-rpois(2000, lambda_i)
2 ИНО
Может попробуете погенерировать
VGAMdata::rpospois()?
Мне представляется некорректным Ваш подход исходного вычитания единицы в одном случае, так как в другом случае Вы этого не делаете. Считаю, что проверять модели надо в одинаковых условиях, а именно после генерации (два варианта):
1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.
Также попробуйте сравнить модели по другим характеристикам (AIC, MSE и т.п.).
Я на том стою, что базовое распределение Пуассона, из коего командой
Код
sim<-sim[sim>0]
тупо повыбрасывали все нули, не имеет никакого отношения к Shift'у и не имитирует вероятностных свойств случайной величины,
"которая по условиям эксперимента не может принимать значение 0" (с).
Кроме того, такая команда равносильна заявлению, что автор пакета VGAM - дурак, не понимающий простых вещей: он-то зачем-то заморочился генерацией ZTP-распределения (rpospois()), а тут всего-то и делов: из обычного распределения Пуассона знай себе нули выбрасывай - и дело с концом.
Вам comisora вроде бы ясно указал на разницу этих распределений, ан не в коня корм.
Ваши "синтетические данные" - это распределение-уродец, на котором одинаково лажают обе модели, поскольку ни одна из них не предназначена для работы с таким "мутантом".
Изначально речь шла о том, что для ZTP распределения мы сначала прогоняем ZPT-модель, получаем "заковыристо интерпретируемые" коэффициенты, а потом на этом же распределении пытаемся "обмануть" обычную регрессию Пуассона.
100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.
comisora, вечерком попробую этот rpospois, хотя покамест слабо представляю, чего он такого может делать, принципиально отличного от моего ручного усекновения нулей, и почему это что-то должно лучше описывать реальный процесс порождения данных. И, хотя у меня есть некоторые сомнения относительно преимущества квадратичных ошибок при оценке счетной величины, на всякий случай посчитаю также и MSE, и RMSE. А вот AIC тут уж точно лишний - некого штрафовать.
Цитата(ИНО @ 20.11.2022 - 14:03)

100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.
Вы как-то очень по-своему понимаете разницу между "с.в. не может принимать нулевые значения" и "может, но мы не будем считать это исходом эксперимента и все равно потом повыкидаем".
Поэтому галиматья получается уже на этапе генерации этой вашей "синтетики": я прошу вас сгенерировать мне выборку объемом 2000 из ZTP. Вы честно генерируете выборку из распределения Пуассона (нуачотакова?) объемом 2000, в которой оказывается порядка 1200 нулей. Вы их выбрасываете, и с честными глазами предъявляете мне выборку объемом 800. А я просил 2000...
Словом сделаете так раз-другой... "А потом ваши рыжие кудри примелькаются, и вас просто начнут бить".
P.S. Ваша клоунада несколько затянулась.
Слышал, что люди применительно к уровню знаний о конкретном предмете делятся на категории:
1) знает, что знает;
2) думает, что знает;
3) думает, что не знает;
4) знает, что не знает;
5) не знает, что знает;
6) не знает, что не знает
Это не полный список комбинаций, но дальше все уж совсем печально, и хочется верить, что среди собравшихся здесь таковых нет.
100$, с каждой вашей новой нападкой на меня Вы становитесь все менее похожим на умного человека. Поясняю: умный человек на вашем месте, прежде чем искрометно блеснуть своей гениальностью и посрамить оппонента, на всякий пожарный проверил бы, а действительно ли по результатам работы мой генератор "распределения-уродца", существенно отличается от расово правильного генератора ZTP из VGАM (автор которого, как выяснилось, точно не дурак), и, если да, то действительно ли эти различия влияют на эффективность обсуждаемых моделей. А то вдруг окажется, что генератор вообще не причем? Всего-то пара примитивных строчек кода и меньше минуты вычислений - а какая замечательная страховка от конфуза! Но Ваше ЧСВ, очевидно, настолько велико, что всякий раз дает железобетонную априорную уверенность в собственной правоте, не нуждающуюся в каких-либо плебейских проверках. Оттого раз за разом Вы попадаете в, скажем так, неловкие ситуации (хотя, возможно, с высоты Вашего персонального олимпа они и не кажутся таковыми, но окружающие, уж поверьте, каждый промах Акелы замечают и мотают на ус). Предыдущая была с некой функцией CWtest() Вашего авторства. Поэтому, если претендуете на высокое звание человека, принадлежащего к первой из вышеозначенных категорий, придется такие досадные промахи исключить полностью путем дотошных перепроверок любого собственного тезиса, позиционируемого как непреложная истина, до его выноса на публику. Моя-то позиция в этом плане несравненно проще - это категория ?4

Так что могу себе позволить в высказываниях заблуждаться безо всяких ограничений, желая быть поправленным более знающими людьми. Увы, покамест в данной теме это мое желание удовлетворено не было.
Код, подтверждающий мои слова, будет позже, сначала попробуйте сами до него дойти, а то надоело приносить сюда готовый на блюдечке, а в благодарность получать лишь упреки и насмешки. Право слово, там все просто, как дважды два.
Только диаграммку приложу, что такое x и y - сами догадайтесь, я в Вас верю!
Цитата(comisora @ 20.11.2022 - 12:47)

1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.
Можно много чего посчитать, но нужно понимание целей процесса. Цель предложенных Вами методов неясна, какие реальные ситуации Вы предлогаете при помощи них смоделировать? Первый предложенный Вами подход как раз и есть генерация shifted Poisson distribution, на статью о которой вы дали ссылку в предыдущем посте. Кстати, я так и не понял, зачем его автор предлагал аж три разных метода генерации и почему не было достаточно просто прибавить или отнять константу, т. е. собственно произвести этот самый shift. Но на практике мне не встречались явления, для которых такое распределение можно было бы допустить в качестве процесса порождения данных. Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским. И именно это фактически делает ZTP (пусть и немножко иным путем).
Второй подход однозначно исключает обычное распределение Пуассона и связанные с ним регрессионные модели. Т. е. спецификация модели будет сильно неверной, и она будет давать сильно смещенный прогноз. Конечно, можно попробовать ради спортивного интереса, но едва ли взлетит.
Проверил и снова удивился: обычная пуассоновская модель на данных, сгенерированных rpospois(), вопреки теоретическим соображениям показала себя тоже неплохо, только на больших выборках (n=1000) она систематически отстает от остальных двух (моей сдвинутой и расово правильной ZTP), на малых же (n=20) положение трех моделей на пьедестале почета меняется случным образом в зависимости от конкретной выдачи датчика случайных чисел, но моя неправильная лидирует чаще других. В любом случае разница мизерная.
Читал я, что разница между прогнозами обычной модели Пуассона и модели ZTP тем больше, чем меньше лямбда. Мол при больших лямбда можно смело применять обычную пуассоновскую модель, забив на отсутствие нулей, а вот при малых - ни-ни. По теоретическим соображениям логично. Но проверка на синтетических данных показала странное. На приведенных ниже диаграммах попарно сравниваются средние абсолютные ошибки трех моделей: моей сдвинутой (MAE), ZTP (MAE2) и обычной пуассоновской (MAE3). Данные сгенерированы расово правильной функцией rpospois(), n=1000. На удивление обычная пуассоновская модель обошла остальные именно на малых лямбда, где по идее должна быть неприменима от слова "совсем". Но обратите внимание на порядок величин разностей, этот мизер может быть в большей степени связан с программной реализацией, чем реальными свойствами распределений. Графики прогнозов и ошибок практически идентичны, так что не привожу.
А что не так с CWTest()'ом?
Там, вроде бы, перестановочный тест Фишера-Питмана не мог отвергнуть гипотезу менее чем на 5%-ном уровне.
Все остальные тесты отвергали бы ее за счет чумовой ошибки I рода.
P.S. Кстати, а что это за Большой Брат, который тут за мной наблюдает и чего-то там себе на что-то наматывает? Вы там всей палатой сгрудились, штоле?
comisora
20.11.2022 - 22:21
Цитата(ИНО @ 20.11.2022 - 18:36)

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

Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским.
Раз в естественном процессе сбора есть нули в результате отбраковки, почему бы не использовать
hurdle-модель?
100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.
comisora, эта гибрид из логичстической модели и ZPT, применяется. наоборот, когда в выборке нулей овер дофига. Наверное, можно расковырять готовую функцию, и оторвать половину гибридной модели, но получим все тот же pospoisson.
Зависимые переменные в моем коде и так одинаковые, столь огорчившие теоретика арифметические фокусы с единицей - это уже часть регрессионной модели. Кстати, собственно к регрессии мы покамест еще перешли, этот шаг запланирован на завтра. Может, с появлением независимых переменных модели наконец-то поведут себя по-разному, и неправильные облажаются. Хотя в теории загвоздка должна была заключаться именно в той части процесса, которая уже проделана, - в моделировании распределения остатков - однако покамест теория что-то плохо контачит с практикой, так что я уже не удивлюсь любым сюрпризам.
comisora
20.11.2022 - 23:53
Цитата(ИНО @ 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)
Цитата(ИНО @ 20.11.2022 - 22:29)

100$, А Вы тему просмотрите (благо тут из свежих негусто) и уведите, что не так с CWTest()'ом. "Чумовая ошибка" первого рода была только у него, у всех прочих - второго, что и следовало ожидать при столь малых размерах выборок. И только это не имеющее аналогов в мире изобретение отличлось! Причем, пока я не продемонстрировал данный конфуз моделированием, Вы уж чуть было не написали пакет имени себя, даже не удосужившись проверить работоспособность собственного кода (а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков). Теперь, похоже, наблюдается рецидив. Вместо бестолковых рассуждений о Большом брате и палатах, лучше б R запустили. и убедились что мой обруганный Вами генератор ровно то же, что и хваленый rpospois() - подсказка относительно рисунка из поста #43 для тех, кто сам не догадался.
Ба! Затянувшаяся клоунада сменилась затянувшейся истерикой.
Так, все-таки, кто это там мои "промахи" себе на что-то наматывает? От чьего имени вы тут вещаете? На гельминтоз давно проверялись?
Вам, видимо, невдомек, что гарантированно ошибку I рода удерживают только перестановочные критерии, вычисляющие точное распределение. У всех остальных (с аппроксимациями) она будет либо чумовая, либо чуть менее чем. Но это общее место в статистике. Зачем проверять теорему Пифагора линейкой?
Так что никакого "конфуза".
Что касается "пакета имени себя", то мотивация была немного иная.
Вот
ссылка, прочитав которую я усомнился в своей способности проделать то, о чем там написано (н-р, п.п. 6-7). Если бы тема получила продолжение, то объяснил бы, в чем дело. Ну, вот объяснил теперь. Лучше поздно, чем никогда.
И, кстати, что такого нерабочего в моем коде?
Что вы разумеете под "работоспособностью"?
Что это у вас там так подгорело?
Здесь же не форум программистов: здесь бесконечное число юзеров способны написать что угодно на чем угодно быстрее и лучше меня. И чо?
Зачем оно нужно вообще и мне в частности? Интересно сначала решить задачу, а потом кропать программы.
Цитата
...И только это не имеющее аналогов в мире изобретение отличлось!
...(а чо, проверки всякие - для слабаков, истинный же мастер защищен святой ЦПТ и авторитетом гуру во веки веков)...
А вот идиотничать я вам не позволю.
100$, с каждым новым постом Вы все менее походите на интеллектуала, и все более - на недалекого гопника из подворотни. Иногда лучше промолчать - за умного сойдешь (с)народная мудрость.
Что касается KWTeest(), приведу Вашу же цитату:
Цитата
Я тут чуток помонтекарлил для равных выборочных объемов (n1=n2) от 3 до 100 способность CWTest() удерживать номинальный уровень значимости (alpha=,05).
Получается, что сколько-нибудь осмысленное его применение возможно только при n1=n2>60, п.ч. сходимость к предельному распределению очень медленная. Ребята, это же ужасно.
Пусть каждый решает для себя, нужен ли ему такой ущербный критерий при наличие давно известных намного лучших альтернатив. Как по мне, это просто мусор. Но не надо втирать, будто иные критерии ведут себя схожим образом, ибо это наглое вранье, CWTest() - воистину аналоговнет. Ладно, предположим, что великий 100$ не стал создавать пакет имени себя не потому, что осознал ущербность этого своего критерия, а чисто в силу неспособности освоить техническую сторону процесса (вернее "сомнения в собственных силах", ведь, как известно, истинный мэтр может все, но имеет право иногда подвергать эту всемогучесть легкому сомнению

). Но почему же он в этой теме не написал ни строчки кода R в подтверждение своих обвинительных заявлений? Неужто и в этой способности через полгода возникли сомнения? Или просто лень одолела?