Почему байесовский вероятный интервал в этой полиномиальной регрессии смещен, тогда как доверительный интервал правильный?

9

Рассмотрим график ниже, на котором я смоделировал данные следующим образом. Мы смотрим на двоичный результат для которого истинная вероятность быть 1 указана черной линией. Функциональная связь между ковариатой и является полиномом 3-го порядка с логистической связью (поэтому она является нелинейной в двустороннем порядке).yobsxp(yobs=1|x)

Зеленая линия - это логистическая регрессия GLM, где вводится как полином 3-го порядка. Пунктирные зеленые линии - это 95% доверительные интервалы вокруг прогноза , где - подогнанные коэффициенты регрессии. Я использовал и для этого.xp(yobs=1|x,β^)β^R glmpredict.glm

Точно так же линия pruple - это среднее значение апостериорного с 95% вероятным интервалом для байесовской модели логистической регрессии с использованием равномерного априора. Для этого я использовал пакет с функцией (настройка дает единый неинформативный априор).p(yobs=1|x,β)MCMCpackMCMClogitB0=0

Красные точки обозначают наблюдения в наборе данных, для которых , черные точки - наблюдения с . Обратите внимание, что, как обычно в классификации / дискретном анализе, наблюдается но не .yobs=1yobs=0yp(yobs=1|x)

введите описание изображения здесь

Можно увидеть несколько вещей:

  1. Я специально симулировал, что редок на левой руке. Я хочу, чтобы доверие и достоверный интервал стали здесь широкими из-за недостатка информации (наблюдений).x
  2. Оба прогноза смещены вверх слева. Это смещение вызвано четырьмя красными точками, обозначающими наблюдения, что ошибочно предполагает, что истинная функциональная форма будет здесь повышаться. Алгоритм не обладает достаточной информацией, чтобы сделать вывод, что истинная функциональная форма имеет нисходящий изгиб.yobs=1
  3. Доверительный интервал становится шире, чем ожидалось, тогда как доверительный интервал - нет . На самом деле доверительный интервал охватывает все пространство параметров, как и должно быть из-за недостатка информации.

Кажется, вероятный интервал здесь неправильный / слишком оптимистичный для части . Это действительно нежелательное поведение для вероятного интервала сужаться, когда информация становится разреженной или полностью отсутствует. Обычно это не то, как реагирует вероятный интервал. Может кто-нибудь объяснить:x

  1. Каковы причины этого?
  2. Какие шаги я могу предпринять, чтобы прийти к более достоверному интервалу? (то есть тот, который включает в себя, по крайней мере, истинную функциональную форму, или, лучше, достигает ширины доверительного интервала)

Код для получения интервалов прогнозирования на графике напечатан здесь:

fit <- glm(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
x_pred <- seq(0, 1, by=0.01)
pred <- predict(fit, newdata = data.frame(x=x_pred), se.fit = T)
plot(plogis(pred$fit), type='l')
matlines(plogis(pred$fit + pred$se.fit %o% c(-1.96,1.96)), type='l', col='black', lty=2)


library(MCMCpack)
mcmcfit <- MCMClogit(y_obs ~ x + I(x^2) + I(x^3), data=data, family=binomial)
gibbs_samps <- as.mcmc(mcmcfit)
x_pred_dm <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=x_pred))
gibbs_preds <- apply(gibbs_samps, 1, `%*%`, t(x_pred_dm))
gibbs_pis <- plogis(apply(gibbs_preds, 1, quantile, c(0.025, 0.975)))
matlines(t(gibbs_pis), col='red', lty=2)

Доступ к данным : https://pastebin.com/1H2iXiew благодаря @DeltaIV и @AdamO

Томка
источник
Если кто-то может объяснить мне, как поделиться таблицей с данными, я могу это сделать.
Томка
Вы можете использовать dputна фрейме данных, содержащем данные, а затем включить dputвывод в виде кода в своем посте.
DeltaIV
1
@ Tomka о, я вижу. Я не дальтоник, но мне очень трудно увидеть разницу между зеленым и синим!
AdamO
1
@AdamO надеюсь , что это лучше
Томка
1
@Flounderer Проверьте, например, stats.stackexchange.com/questions/26450/… или stats.stackexchange.com/questions/6652/…
Тим

Ответы:

6

Для модели, частотная дисперсия в предсказания не увеличится пропорционально квадрату расстояния от центроида . Ваш метод расчета интервалов прогнозирования для байесовской GLM использует эмпирические квантили на основе подобранной кривой вероятности, но не учитывает левереджXX

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

Обратите внимание, что любое полиномиальное представление вероятностей логита приводит к предсказаниям риска, которые сходятся к 0 как и 1 как или наоборот, в зависимости от знака члена высшего полиномиального порядка .XX

Для частых прогнозов доминирует эта тенденция в квадрате отклонений (рычагов) пропорционального увеличения дисперсии прогнозов. Вот почему скорость сходимости к интервалам предсказания, приблизительно равная [0, 1], выше, чем полиномиальная логитная сходимость третьего порядка к вероятностям 0 или 1, в частности.

Это не так для байесовских задних квантилей. Нет явного использования квадрата отклонения, поэтому мы полагаемся просто на долю доминирующих 0 или 1 тенденций для построения интервалов долгосрочного прогнозирования.

Это стало очевидным экстраполяцией очень далеко в крайности .X

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

> x_pred_dom <- model.matrix(~ x + I(x^2) + I(x^3), data=data.frame('x'=c(1000)))
> gibbs_preds <- plogis(apply(gibbs_samps[1000:10000, ], 1, `%*%`, t(x_pred_dom))) # a bunch of 0/1s basically past machine precision
> prop.table(table(gibbs_preds))
gibbs_preds
         0          1 
0.97733585 0.02266415 
> 

Таким образом, в 97,75% случаев третий полиномиальный член был отрицательным. Это подтверждается образцами Гиббса:

> prop.table(table(gibbs_samps[, 4]< 0))

 FALSE   TRUE 
0.0225 0.9775 

Следовательно, предсказанная вероятность сходится к 0, когда уходит в бесконечность. Если мы проверяем SE в байесовской модели, мы находим, что оценка третьего полиномиального члена равна -185,25, а se 108,81 означает, что это 0,70 SD от 0, поэтому, используя нормальные законы вероятности, он должен упасть ниже 0 95,5% времени ( не совсем другой прогноз, основанный на 10 000 итераций). Просто еще один способ понять это явление.X

С другой стороны, подгонка частых ударов до 0,1, как и ожидалось:

freq <- predict(fit, newdata = data.frame(x=1000), se.fit=T)
plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)

дает:

> plogis(freq$fit + c(-1.96, 1.96) %o% freq$se.fit)
     [,1]
[1,]    0
[2,]    1
Adamo
источник
Тем не менее: не является ли байесовская модель чрезмерно уверенной в областях данных , из которых она не видела примеров? Я знаю, что байесовские постеры или предиктивные распределения часто ведут себя очень по-разному (то есть больше похоже на интервал конф.). Я подозреваю, что есть какое-то влияние предыдущего. Если вы манипулируете внутри, вы указываете точность обычного априора и можете наблюдать значительное влияние на вероятный интервал. xB0MCMClogit
Томка
@tomka Я не знаю, как точно ответить на этот вопрос, поскольку он кажется касательным к рассматриваемому вопросу. Наиболее важным является указание на то, что эти методы расчета ИП на самом деле не сопоставимы, особенно в том, что касается экстраполяции. Конечно, с помощью байесовского вывода, если вы используете информативный априор, вы получаете эффективность, когда априор прав, и теряете, когда априор неверен.
AdamO
Просто чтобы сообщить вам, что я все еще думаю о вашем ответе. Я все еще чувствую, что странно, что задняя часть не реагирует на разреженность расширением. Я считаю, что для других приоров можно добиться лучшего поведения в редком регионе. Я не могу определить это точно в данный момент; Возможно, я дополню этот вопрос примером, в котором вероятный интервал работает так, как я ожидал, даже в случае экстраполяции (в частности, я имею в виду нормальную линейную байесовскую регрессию). Когда я это сделаю, я дам вам знать.
Томка