Введение

Традиционно мы выбираем модель и выполняем вывод и предсказание в зависимости от этой модели. Наши результаты, как правило, не учитывают неопределенность при выборе модели и поэтому могут быть излишне оптимистичными. Они могут быть даже неверными, если выбранная нами модель существенно отличается от истинной модели генерации данных (МГД). В некоторых приложениях мы можем располагать убедительными теоретическими или эмпирическими данными о DGM. В других приложениях, обычно имеющих сложную и нестабильную природу, таких как экономика, психология и эпидемиология, выбор одной надежной модели может быть затруднен. Вместо того чтобы полагаться только на одну модель, усреднение моделей приводит к усреднению результатов по нескольким правдоподобным моделям, основанным на наблюдаемых данных. В BMA «правдоподобность» модели описывается апостериорной вероятностью модели (PMP), которая определяется с использованием фундаментальных принципов Байеса — теоремы Байеса — и применяется универсально для всех анализов данных. BMA может использоваться для учета неопределенности модели при оценке параметров модели и прогнозировании новых наблюдений, что позволяет избежать слишком оптимистичных выводов. Она особенно полезна в приложениях с несколькими правдоподобными моделями, когда нет однозначной причины для выбора конкретной модели, а не других. Но даже если выбор одной модели является конечной целью, BMA может оказаться полезным. Она предоставляет принципиальный способ выявления важных моделей и предикторов в рамках рассматриваемых классов моделей. Его структура позволяет изучать взаимосвязи между различными предикторами с точки зрения их склонности появляться в модели вместе, отдельно или независимо друг от друга. С его помощью можно оценить чувствительность конечных результатов к различным предположениям о важности различных моделей и предикторов. И обеспечивает оптимальные прогнозы в смысле log-score.

Комплект bma

В регрессионных условиях неопределенность модели сводится к неопределенности того, какие предикторы должны быть включены в модель. Для учета выбора предикторов в линейной регрессии можно использовать bmaregress. Он исследует пространство моделей либо исчерпывающим образом с помощью опции перечисления, когда это возможно, либо с помощью алгоритма составления моделей с марковской цепью Монте-Карло (MCMC) с опцией выборки (MC3). При этом выдаются различные сводки по посещенным моделям и включенным в них предикторам, а также постерные распределения параметров моделей. Позволяет задавать группы предикторов, включаемых или исключаемых из модели вместе, а также включаемых во все модели. В опции mprior() предусмотрены различные предшествующие распределения для моделей, а в опции gprior() — для параметра g, управляющего уменьшением коэффициентов регрессии к нулю. Поддерживаются факторные переменные и операторы временных рядов, а также несколько способов обработки взаимодействий в процессе оценки с помощью опции heredity().

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

Давайте посмотрим, как это работает

-> Перечисление моделей -> Достоверные интервалы -> Влиятельные модели -> Важные предикторы ->Распределение моделей по размерам -> Апостериорные распределения коэффициентов -> Прогноз BMA -> Информативность модельных приматов -> Прогнозные характеристики с использованием LPS -> Очистка

Ниже мы приводим описание некоторых функций bma на примере игрушечного набора данных. Этот набор содержит n=200

�=200 наблюдений, p=10 �=10 ортогональных предикторов, а результат y формируется как

y=0,5+1,2×x2+5×x10+ϵ

�=0.5+1.2�2+5�10+�

где ϵ∼N(0,1)

�∼�(0,1) — стандартный нормальный член ошибки.

. webuse bmaintro
(Simulated data for BMA example)

. summarize

Перечисление моделей

С помощью bmaregress мы настраиваем BMA-линейную регрессию y на x1 — x10.

. bmaregress y x1-x10

Enumerating models ...
Computing model probabilities ...

Bayesian model averaging                            No. of obs         =    200
Linear regression                                   No. of predictors  =     10
Model enumeration                                               Groups =     10
                                                                Always =      0
Priors:                                             No. of models      =  1,024
  Models: Beta-binomial(1, 1)                           For CPMP >= .9 =      9
   Cons.: Noninformative                            Mean model size    =  2.479
   Coef.: Zellner's g
       g: Benchmark, g = 200                        Shrinkage, g/(1+g) = 0.9950
  sigma2: Noninformative                            Mean sigma2        =  1.272

Note: Coefficient posterior means and std. dev. estimated from 1,024 models.
Note: Default priors are used for models and parameter g.

При 10 предикторах и фиксированном (эталонном) значении g-prior по умолчанию bmaregress использует перечисление моделей и исследует все 210=1,024

210=1 024 возможных моделей. По умолчанию используется бета-биномиальный приоритет модели с параметрами формы, равными 1, который равномерно распределяется по размеру модели. Априори в нашей модели BMA предполагалось незначительное уменьшение коэффициентов регрессии к нулю, так как g/(1+g)=0,9950 �/(1+�)=0,9950 близко к 1.

bmaregress определил x2 и x10 как очень важные предикторы — их апостериорные вероятности включения (PIP) практически равны 1. Апостериорные средние оценки их коэффициентов регрессии составляют, соответственно, 1,2 (округленно) и 5,1 и очень близки к истинным значениям 1,2 и 5. Все остальные предикторы имеют гораздо более низкие PIP и оценки коэффициентов, близкие к нулю. Полученные нами результаты BMA согласуются с истинным значением DGM.

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

. bmaregress, saving(bmareg)
note: file bmareg.dta saved.

. estimates store bmareg

Достоверные интервалы

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

Сначала мы используем bmacoefsample для генерации MCMC-выборки из апостериорного распределения параметров модели, включающих коэффициенты регрессии. Затем мы используем существующую команду bayesstats summary для вычисления апостериорных сумм сгенерированной MCMC-выборки.

. bmacoefsample, rseed(18)

Simulation (10000): ....5000....10000 done

. bayesstats summary

Posterior summary statistics                      MCMC sample size =    10,000

95%-ный критерий для коэффициента x2 составляет [1,05, 1,34], а для x10 — [4,9, 5,26], что согласуется с нашей ДГМ.

Влиятельные модели

Оценки коэффициентов BMA усреднены по 1024 моделям. Важно выяснить, какие модели являются влиятельными. Для этого мы можем использовать модели bmastats, чтобы исследовать ПМП.

. bmastats models

Model summary         Number of models:
                               Visited = 1,024
                              Reported =     5

Analytical PMP Model size .6292 2 .1444 3 .0258 3 .0246 3 .01996 3

Variable-inclusion summary

Legend: 
x - estimated

Неудивительно, что модель, содержащая и x2, и x10, имеет самый высокий PMP — около 63%. Фактически, две лучшие модели соответствуют суммарному PMP около 77%:

. bmastats models, cumulative(0.75)

Computing model probabilities ...

Model summary           Number of models:
                                 Visited = 1,024
                                Reported =     2

Variable-inclusion summary

Legend: 
x - estimated

Мы можем использовать bmagraph pmp для построения графика кумулятивного PMP.

. bmagraph pmp, cumulative

По умолчанию команда выводит на экран первые 100 моделей, но можно указать опцию top(), если вы хотите увидеть больше моделей.

Важные предикторы

Мы можем наглядно изучить важность предикторов, используя bmagraph varmap для создания карты включения переменных.

. bmagraph varmap

Computing model probabilities ...

x2 и x10 входят во все 100 лучших моделей, ранжированных по PMP. Их коэффициенты положительны во всех моделях.

Распределение моделей по размерам

Мы можем исследовать сложность нашей модели BMA с помощью bmastats msize и bmagraph msize для изучения предшествующего и последующего распределений размеров модели.

. bmastats msize

Model-size summary

Number of models = 1,024
Model size:
  Minimum =  0
  Maximum = 10

Note: Frequency summaries not available.

. bmagraph msize

Предварительное распределение размера модели равномерно по размеру модели. Апостериорное распределение размера модели имеет перекос влево с модой около 2. Средний размер модели в априорном случае равен 5, а в апостериорном — 2,48. Таким образом, исходя из наблюдаемых данных, BMA склоняется к меньшим моделям с двумя предикторами в среднем по сравнению с нашим априорным предположением.

Апостериорные распределения коэффициентов

С помощью bmagraph coefdensity можно построить апостериорные распределения коэффициентов регрессии.

. bmagraph coefdensity {x2} {x3}, combine

Эти распределения представляют собой смесь точечной массы в нуле, соответствующей вероятности невключения предиктора, и непрерывной плотности, обусловленной включением. Для коэффициента x2 вероятность невключения пренебрежимо мала, поэтому его апостериорное распределение представляет собой, по сути, непрерывную, достаточно симметричную плотность с центром в точке 1,2. Для коэффициента x3, помимо условной непрерывной плотности, в нулевой точке имеется красная вертикальная линия, соответствующая апостериорной вероятности невключения, равной примерно 1-,2 = 0,8.

Прогноз BMA

bmapredict выдает различные прогнозы BMA. Например, вычислим апостериорные прогностические средства.

. bmapredict pmean, mean
note: computing analytical posterior predictive means.

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

. bmacoefsample, saving(bmacoef)
note: saving existing MCMC simulation results without resampling; specify
      option simulate to force resampling in this case.
note: file bmacoef.dta saved.

. bmapredict cri_l cri_u, cri rseed(18)
note: computing credible intervals using simulation.

Computing predictions ...

Мы обобщаем прогнозируемые результаты:

.  summarize y pmean cri*

Приведенные средние значения по наблюдениям для прогностических средств, а также нижняя и верхняя 95%-ные прогностические границы CrI выглядят разумными относительно наблюдаемого результата y.

Информативность модельных приматов

Одним из достоинств BMA является возможность учета предварительной информации о моделях и предикторах. Это позволяет исследовать чувствительность результатов к различным предположениям о важности моделей и предикторов. Предположим, что у нас есть достоверная информация о том, что все предикторы, кроме x2 и x10, с меньшей вероятностью связаны с результатом. Мы можем задать, например, предварительную биномиальную модель с низкой вероятностью включения этих предикторов.

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

. splitsample, generate(sample) nsplit(2) rseed(18)

. bmaregress y x1-x10 if sample == 1, mprior(binomial x2 x10 0.5 x1 x3-x9 0.05)
  saving(bmareg_inf)

Enumerating models ...
Computing model probabilities ...

Bayesian model averaging                          No. of obs         =    100
Linear regression                                 No. of predictors  =     10
Model enumeration                                             Groups =     10
                                                              Always =      0
Priors:                                           No. of models      =  1,024
  Models: Binomial, IP varies                         For CPMP >= .9 =      1
   Cons.: Noninformative                          Mean model size    =  2.072
   Coef.: Zellner's g
       g: Benchmark, g = 100                      Shrinkage, g/(1+g) = 0.9901
  sigma2: Noninformative                          Mean sigma2        =  1.268

Note: Coefficient posterior means and std. dev. estimated from 1,024 models.
Note: Default prior is used for parameter g.
Note: 4 predictors with PIP less than .01 not shown.

file bmareg_inf.dta saved.

. estimates store bmareg_inf

Эффект этого приоритета модели заключается в том, что PIP всех неважных предикторов теперь составляет менее 2%.

Обратите внимание, что при выборе одной модели мы, по сути, подгоняем BMA-модель с очень сильным приоритетом, согласно которому все выбранные предикторы должны быть включены с вероятностью 1. Например, мы можем заставить bmaregress включить в модель все переменные:

. bmaregress y (x1-x10, always) if sample == 1, saving(bmareg_all)

Enumerating models ...
Computing model probabilities ...

Bayesian model averaging                           No. of obs         =    100
Linear regression                                  No. of predictors  =     10
Model enumeration                                              Groups =      0
                                                               Always =     10
Priors:                                            No. of models      =      1
  Models: Beta-binomial(1, 1)                          For CPMP >= .9 =      1
   Cons.: Noninformative                           Mean model size    = 10.000
   Coef.: Zellner's g
       g: Benchmark, g = 100                       Shrinkage, g/(1+g) = 0.9901
  sigma2: Noninformative                           Mean sigma2        =  1.192


Mean Std. dev. Group PIP

.1294521 .105395 0 1

1.166679 .1129949 0 1

-.1433074 .1271903 0 1

-.1032189 .1223152 0 1

-.0819008 .1261309 0 1

.0696633 .1057512 0 1

.0222949 .1215404 0 1

-.0252311 .1124352 0 1

.0587412 .1166921 0 1

4.949992 .1276795 0 1

.6006153 .1127032 0 1

Note: Coefficient posterior means and std. dev. estimated from 1 model.
Note: Default priors are used for models and parameter g.

file bmareg_all.dta saved.

. estimates store bmareg_all

Эффективность прогнозирования с использованием LPS

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

. qui bmaregress y x1-x10 if sample == 1, saving(bmareg, replace)
. estimates store bmareg

Мы можем сравнить предсказательную эффективность наших BMA-моделей за пределами выборки, используя bmastats lps для вычисления логарифмической предсказательной оценки (LPS) для наблюдений за пределами выборки.

. bmastats lps bmareg bmareg_inf bmareg_all if sample == 2, compact

Log predictive-score (LPS)

Number of observations = 100

Notes: Using analytical PMPs.
       Result bmareg_inf has the smallest mean LPS.

Более информативная модель bmareg_inf имеет несколько меньшее среднее значение LPS, но итоговые значения LPS для всех моделей очень близки. О том, как сравнить производительность моделей BMA с помощью перекрестной валидации, см. в [BMA] bmastats lps.

Очистка

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

. erase bmareg.dta
. erase bmacoef.dta
. erase bmareg_inf.dta
. erase bmareg_all.dta