Въведение 

Традиционно избираме един модел и извършваме изводи и прогнози, които зависят от този модел. Нашите резултати обикновено не отчитат несигурността при избора на модел и по този начин могат да бъдат прекалено оптимистични. Те могат дори да бъдат неверни, ако избраният от нас модел се различава съществено от истинския модел, генериращ данни (DGM). В някои приложения може да разполагаме с убедителни теоретични или емпирични доказателства за DGM. В други приложения, обикновено със сложен и нестабилен характер, като тези в икономиката, психологията и епидемиологията, изборът на един надежден модел може да бъде труден.

Вместо да се разчита само на един модел, усредняването на моделите осреднява резултатите от множество правдоподобни модели въз основа на наблюдаваните данни. При БМА „достоверността“ на модела се описва от апостериорната вероятност на модела (АВП), която се определя с помощта на основните Байесови принципи – теоремата на Байес – и се прилага универсално за всички анализи на данни.

BMA може да се използва за отчитане на несигурността на модела при оценяване на параметрите на модела и прогнозиране на нови наблюдения, за да се избегнат прекалено оптимистични заключения. Той е особено полезен в приложения с няколко правдоподобни модела, при които няма една категорична причина да се избере конкретен модел пред останалите. Но дори ако крайната цел е избор на един-единствен модел, BMA може да ви бъде от полза. Той осигурява принципен начин за идентифициране на важни модели и предиктори в рамките на разглежданите класове модели. Нейната рамка ви позволява да научите за взаимовръзките между различните предиктори по отношение на тенденцията им да се появяват в модела заедно, поотделно или независимо. Тя може да се използва за оценка на чувствителността на крайните резултати към различни допускания за важността на различните модели и предиктори. И осигурява оптимални прогнози в смисъл на логаритмичен резултат. 

Пакетът bma 

В регресионна среда несигурността на модела е свързана с несигурността на това кои предиктори трябва да бъдат включени в модела. Можем да използваме bmaregress, за да отчетем избора на предиктори в линейна регресия. Той изследва пространството на модела или изчерпателно с опцията за изброяване, когато това е възможно, или чрез използване на алгоритъма за съставяне на модели по веригата на Марков Монте Карло (MCMC) (MC3) с опцията за вземане на проби. Той отчита различни обобщения на посетените модели и включените предиктори и на апостериорните разпределения на параметрите на модела. Позволява да се задават групи от предиктори, които да се включват или изключват заедно от даден модел, както и такива, които се включват във всички модели. Той предоставя различни предварителни разпределения за моделите в опцията mprior() и за параметъра g, който контролира свиването на регресионните коефициенти към нула, в опцията gprior(). Поддържа също така факторни променливи и оператори за времеви редове и предоставя няколко начина за обработка на взаимодействията по време на оценяването с помощта на опцията 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,024210=1,024 възможни модела. Моделът по подразбиране е бета-биномиален с параметри на формата 1, който е еднакъв за размера на модела. A priori нашият модел на 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 не съобщава достоверни интервали (CrIs), тъй като изисква потенциално отнемаща време симулация на апостериорната извадка от параметри на модела. Но ние можем да ги изчислим след оценяването.  

Първо използваме bmacoefsample, за да генерираме MCMC извадка от апостериорното разпределение на параметрите на модела, които включват регресионни коефициенти. След това използваме съществуващата команда bayesstats summary, за да изчислим апостериорните обобщения на генерираната MCMC извадка. 

. bmacoefsample, rseed(18)

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

. bayesstats summary

Posterior summary statistics                      MCMC sample size =    10,000

95% CrI за коефициента x2 е [1,05, 1,34], а този за x10 е [4,9, 5,26], което е в съответствие с нашия DGM. 

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

Оценките на коефициентите на BMA са осреднени за 1024 модела. Важно е да се проучи кои модели са влиятелни. Можем да използваме моделите на bmastats, за да изследваме PMP. 

. 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.

Можем също така да изчислим прогнозни CrI, за да оценим несигурността на нашите прогнози. (Това не е възможно при много от традиционните методи за избор на модел.) Имайте предвид, че тези прогнозни CrI включват и несигурността на модела. За да изчислим прогнозни CrI, първо трябва да запазим апостериорната 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 модел. 

. 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] bmastats lps за това как да сравнявате представянето на моделите на BMA, като използвате кръстосано валидиране. 

Почистване 

По време на нашия анализ генерирахме няколко набора от данни. Вече не се нуждаем от тях, затова ги изтриваме накрая. Но вие може да решите да ги запазите, особено ако те съответстват на потенциално отнемащ време финален анализ. 

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