Я пытаюсь приспособить обобщенные линейные модели к некоторым наборам данных подсчета, которые могут быть или не быть перераспределены. Здесь применимы два канонических распределения: Пуассон и Отрицательный бином (Негбин) с EV и дисперсией.
которые могут быть установлены в R с использованием glm(..,family=poisson)
и glm.nb(...)
, соответственно. Существует также quasipoisson
семья, которая в моем понимании является настроенным Пуассоном с тем же EV и дисперсией
,
то есть попадание где-то между Пуассоном и Негбиным. Основная проблема с семейством квазипуассонов состоит в том, что для него нет соответствующей вероятности, и поэтому многие чрезвычайно полезные статистические тесты и критерии соответствия (AIC, LR и так далее) недоступны.
Если вы сравните дисперсии QP и Negbin, вы можете заметить, что вы можете приравнять их, поставив . Продолжая эту логику, вы можете попытаться выразить квазипуассонное распределение как частный случай Негбина:
,
то есть негин с линейно зависимым от . Я попытался проверить эту идею, сгенерировав случайную последовательность чисел в соответствии с приведенной выше формулой и подгоняя ее к :μglm
#fix parameters
phi = 3
a = 1/50
b = 3
x = 1:100
#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison
mu = exp(a*x+b)
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator
#fit a generalized linear model y = f(x)
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial
> glmQP
Call: glm(formula = y ~ x, family = quasipoisson)
Coefficients:
(Intercept) x
3.11257 0.01854
(Dispersion parameter for quasipoisson family taken to be 3.613573)
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 2097
Residual Deviance: 356.8 AIC: NA
> glmNB
Call: glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)
Coefficients:
(Intercept) x
3.10182 0.01873
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 578.1
Residual Deviance: 107.8 AIC: 824.7
Оба соответствия воспроизводят параметры, и квазипуассон дает «разумную» оценку для . Теперь мы можем также определить значение AIC для квазипуассона:
df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values
#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329
(Мне пришлось вручную скопировать подходящее значение , поскольку я не смог найти его в объекте)summary(glmQP)
glmQP
Поскольку , это указывает на то, что квазипуассон, что неудивительно, лучше подходит; поэтому, по крайней мере, делает то, что должен, и, следовательно, это может быть разумным определением для AIC (и, соответственно, вероятности) квазипуассона. Большие вопросы, которые у меня остались A I C Q P
- Эта идея имеет смысл? Моя проверка основана на циклическом рассуждении?
- Главный вопрос для любого, кто «изобретает» что-то, что, по-видимому, отсутствует в устоявшейся теме: если эта идея имеет смысл, почему она уже не реализована
glm
?
Редактировать: фигура добавлена
Ответы:
Квази-Пуассон - это не модель полного максимального правдоподобия (ML), а модель квази-ML. Вы просто используете функцию оценки (или функцию оценки) из модели Пуассона для оценки коэффициентов, а затем используете определенную функцию дисперсии для получения подходящих стандартных ошибок (или, скорее, полной ковариационной матрицы) для выполнения вывода. Следовательно,
glm()
не поставляет и /logLik()
илиAIC()
здесь и т. Д.size
Если нет регрессоров (только перехват), то параметризация NB1 и параметризация NB2, используемые в
MASS
s,glm.nb()
совпадают. С регрессорами они отличаются. В статистической литературе чаще используется параметризация NB2, но некоторые программные пакеты также предлагают версию NB1. Например, в R вы можете использоватьgamlss
пакет, чтобы сделатьgamlss(y ~ x, family = NBII)
. Обратите внимание, что несколько сбивает с толкуgamlss
использованиеNBI
для параметризации NB2 иNBII
для NB1. (Но жаргон и терминология не едины для всех сообществ.)Тогда вы могли бы спросить, конечно, зачем использовать квази-Пуассона, если есть доступный NB1? Есть еще небольшая разница: первый использует квази-ML и получает оценку из дисперсии от квадратов отклонений (или Пирсона). Последний использует полный ML. На практике разница часто невелика, но мотивы использования любой модели немного отличаются.
источник
gamlss
и, похоже, это именно то, что мне нужно. Не могли бы вы рассказать о мотивах использования квази-правдоподобия по сравнению с полным ОД?vignette("countreg", package = "pscl")
.