Zero-truncated Poisson regression Опции

[comisora]

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]

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] ) # число случаи, в которых моя модель оказалась лучше расово правильной

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

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

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

С чем сравнивались обе модели?

[ИНО]

Друг с другом.

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

И какая из них считается правильной?

[ИНО]

Этот вопрос лучше адресовать $100, поскольку мне , в отличие от него, неведомы критерии истинно правильности моделей. Могу лишь сказать, что по результатам тестирования на синтетических данных у моей модели оказалась чуть ниже ошибка предсказания. Но это очень малое отличие, на практике совершенно несущественное. Непонятно иное: почему вопреки теоретическим соображениям, моя модель не оказалась значительно хуже, чем та, которая, по идее, должна идеально описывать именно то распределение, данные из которого я генерирвал

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

[comisora]

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

Код

sim<-rpois(2000, lambda_i)

2 ИНО

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

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

Сообщение отредактировал comisora - 20.11.2022 - 12:49

[100$]

Я на том стою, что базовое распределение Пуассона, из коего командой

Код

sim<-sim[sim>0]

тупо повыбрасывали все нули, не имеет никакого отношения к Shift'у и не имитирует вероятностных свойств случайной величины,
"которая по условиям эксперимента не может принимать значение 0" (с).
Кроме того, такая команда равносильна заявлению, что автор пакета VGAM - дурак, не понимающий простых вещей: он-то зачем-то заморочился генерацией ZTP-распределения (rpospois()), а тут всего-то и делов: из обычного распределения Пуассона знай себе нули выбрасывай - и дело с концом.

Вам comisora вроде бы ясно указал на разницу этих распределений, ан не в коня корм.

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

Изначально речь шла о том, что для ZTP распределения мы сначала прогоняем ZPT-модель, получаем "заковыристо интерпретируемые" коэффициенты, а потом на этом же распределении пытаемся "обмануть" обычную регрессию Пуассона.

Сообщение отредактировал 100$ - 20.11.2022 - 14:03

[ИНО]

100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.

comisora, вечерком попробую этот rpospois, хотя покамест слабо представляю, чего он такого может делать, принципиально отличного от моего ручного усекновения нулей, и почему это что-то должно лучше описывать реальный процесс порождения данных. И, хотя у меня есть некоторые сомнения относительно преимущества квадратичных ошибок при оценке счетной величины, на всякий случай посчитаю также и MSE, и RMSE. А вот AIC тут уж точно лишний - некого штрафовать.

[100$]

100$, а как по-вашему получается этот "позитивный Пуассон" в реальности? По-моему, именно отбрасыванием всех нулевых значений (они просто не попадают в нашу, как однажды выразился 100$, неслучайную выборку) - ровно то же, что делает мой генератор. Конечно, на самом деле реальное изначальное распределение может и не быть пуассоновским (скорее всего, в большинстве реальных случаев это будет "zero-inflated") , но для простоты пока что условимся считать его таковым. В чем согласен: к "шифту" мой генератор действительно не имеет отношения, но чудесным образом на нем шифт-модель работает лучше, чем ZTP.

Вы как-то очень по-своему понимаете разницу между "с.в. не может принимать нулевые значения" и "может, но мы не будем считать это исходом эксперимента и все равно потом повыкидаем".

Поэтому галиматья получается уже на этапе генерации этой вашей "синтетики": я прошу вас сгенерировать мне выборку объемом 2000 из ZTP. Вы честно генерируете выборку из распределения Пуассона (нуачотакова?) объемом 2000, в которой оказывается порядка 1200 нулей. Вы их выбрасываете, и с честными глазами предъявляете мне выборку объемом 800. А я просил 2000...
Словом сделаете так раз-другой... "А потом ваши рыжие кудри примелькаются, и вас просто начнут бить".

P.S. Ваша клоунада несколько затянулась.

Сообщение отредактировал 100$ - 20.11.2022 - 15:21

[ИНО]

Слышал, что люди применительно к уровню знаний о конкретном предмете делятся на категории:

1) знает, что знает;
2) думает, что знает;
3) думает, что не знает;
4) знает, что не знает;
5) не знает, что знает;
6) не знает, что не знает

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

100$, с каждой вашей новой нападкой на меня Вы становитесь все менее похожим на умного человека. Поясняю: умный человек на вашем месте, прежде чем искрометно блеснуть своей гениальностью и посрамить оппонента, на всякий пожарный проверил бы, а действительно ли по результатам работы мой генератор "распределения-уродца", существенно отличается от расово правильного генератора ZTP из VGАM (автор которого, как выяснилось, точно не дурак), и, если да, то действительно ли эти различия влияют на эффективность обсуждаемых моделей. А то вдруг окажется, что генератор вообще не причем? Всего-то пара примитивных строчек кода и меньше минуты вычислений - а какая замечательная страховка от конфуза! Но Ваше ЧСВ, очевидно, настолько велико, что всякий раз дает железобетонную априорную уверенность в собственной правоте, не нуждающуюся в каких-либо плебейских проверках. Оттого раз за разом Вы попадаете в, скажем так, неловкие ситуации (хотя, возможно, с высоты Вашего персонального олимпа они и не кажутся таковыми, но окружающие, уж поверьте, каждый промах Акелы замечают и мотают на ус). Предыдущая была с некой функцией CWtest() Вашего авторства. Поэтому, если претендуете на высокое звание человека, принадлежащего к первой из вышеозначенных категорий, придется такие досадные промахи исключить полностью путем дотошных перепроверок любого собственного тезиса, позиционируемого как непреложная истина, до его выноса на публику. Моя-то позиция в этом плане несравненно проще - это категория ?4 Так что могу себе позволить в высказываниях заблуждаться безо всяких ограничений, желая быть поправленным более знающими людьми. Увы, покамест в данной теме это мое желание удовлетворено не было.

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

Только диаграммку приложу, что такое x и y - сами догадайтесь, я в Вас верю!

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

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

 

[ИНО]

1) Без фильтрования нулей добавить ко всему ряду +1;
2) После фильтрации нулей работать с тем, что получилось без добавлений и вычитаний.

Можно много чего посчитать, но нужно понимание целей процесса. Цель предложенных Вами методов неясна, какие реальные ситуации Вы предлогаете при помощи них смоделировать? Первый предложенный Вами подход как раз и есть генерация shifted Poisson distribution, на статью о которой вы дали ссылку в предыдущем посте. Кстати, я так и не понял, зачем его автор предлагал аж три разных метода генерации и почему не было достаточно просто прибавить или отнять константу, т. е. собственно произвести этот самый shift. Но на практике мне не встречались явления, для которых такое распределение можно было бы допустить в качестве процесса порождения данных. Реальные данные, с которыми я работаю как раз и есть результат "фильтрации нулей", таков естественный процесс их сбора. Именно его я попытался восседать в своем генераторе, приняв допущение о том, что до отбраковки распределение было пуассоновским. И именно это фактически делает ZTP (пусть и немножко иным путем).

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

Проверил и снова удивился: обычная пуассоновская модель на данных, сгенерированных rpospois(), вопреки теоретическим соображениям показала себя тоже неплохо, только на больших выборках (n=1000) она систематически отстает от остальных двух (моей сдвинутой и расово правильной ZTP), на малых же (n=20) положение трех моделей на пьедестале почета меняется случным образом в зависимости от конкретной выдачи датчика случайных чисел, но моя неправильная лидирует чаще других. В любом случае разница мизерная.

Читал я, что разница между прогнозами обычной модели Пуассона и модели ZTP тем больше, чем меньше лямбда. Мол при больших лямбда можно смело применять обычную пуассоновскую модель, забив на отсутствие нулей, а вот при малых - ни-ни. По теоретическим соображениям логично. Но проверка на синтетических данных показала странное. На приведенных ниже диаграммах попарно сравниваются средние абсолютные ошибки трех моделей: моей сдвинутой (MAE), ZTP (MAE2) и обычной пуассоновской (MAE3). Данные сгенерированы расово правильной функцией rpospois(), n=1000. На удивление обычная пуассоновская модель обошла остальные именно на малых лямбда, где по идее должна быть неприменима от слова "совсем". Но обратите внимание на порядок величин разностей, этот мизер может быть в большей степени связан с программной реализацией, чем реальными свойствами распределений. Графики прогнозов и ошибок практически идентичны, так что не привожу.

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

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

 

[100$]

А что не так с CWTest()'ом?
Там, вроде бы, перестановочный тест Фишера-Питмана не мог отвергнуть гипотезу менее чем на 5%-ном уровне.

Все остальные тесты отвергали бы ее за счет чумовой ошибки I рода.

P.S. Кстати, а что это за Большой Брат, который тут за мной наблюдает и чего-то там себе на что-то наматывает? Вы там всей палатой сгрудились, штоле?