Úvod

Tradičně vybíráme model a provádíme inferenci a predikci podmíněnou tímto modelem. Naše výsledky obvykle nezohledňují nejistotu při výběru modelu, a mohou být proto příliš optimistické. Mohou být dokonce nesprávné, pokud se námi zvolený model podstatně liší od skutečného modelu generujícího data (DGM). V některých aplikacích můžeme mít o DGM silné teoretické nebo empirické důkazy. V jiných aplikacích, obvykle složité a nestabilní povahy, jako jsou aplikace v ekonomii, psychologii a epidemiologii, může být výběr jednoho spolehlivého modelu obtížný. Namísto spoléhání se pouze na jeden model se při průměrování modelů zprůměrují výsledky z více věrohodných modelů založených na pozorovaných datech. V BMA je „věrohodnost“ modelu popsána posteriorní modelovou pravděpodobností (PMP), která je určena pomocí základních bayesovských principů – Bayesovy věty – a je univerzálně aplikována na všechny analýzy dat. BMA lze použít k zohlednění nejistoty modelu při odhadu parametrů modelu a předpovědi nových pozorování, aby se předešlo příliš optimistickým závěrům. Je zvláště užitečná v aplikacích s několika pravděpodobnými modely, kde neexistuje jediný definitivní důvod pro výběr určitého modelu na úkor ostatních. Ale i v případě, že je konečným cílem výběr jediného modelu, může být BMA přínosná. Poskytuje principiální způsob, jak identifikovat důležité modely a prediktory v rámci uvažovaných tříd modelů. Její rámec umožňuje poznat vzájemné vztahy mezi různými prediktory z hlediska jejich tendence objevovat se v modelu společně, samostatně nebo nezávisle. Lze jej použít k vyhodnocení citlivosti konečných výsledků na různé předpoklady o důležitosti různých modelů a prediktorů. A poskytuje optimální predikce ve smyslu logaritmického skóre.

Sada bma

V regresním prostředí se nejistota modelu rovná nejistotě, které prediktory by měly být do modelu zahrnuty. K zohlednění výběru prediktorů v lineární regresi můžeme použít bmaregress. Prozkoumá modelový prostor buď vyčerpávajícím způsobem s možností výčtu, kdykoli je to proveditelné, nebo pomocí algoritmu složení modelu Markovovým řetězcem Monte Carlo (MCMC) (MC3) s možností vzorkování. Uvádí různé přehledy navštívených modelů a zahrnutých prediktorů a posteriorních rozdělení parametrů modelu. Umožňuje určit skupiny prediktorů, které mají být zahrnuty nebo vyloučeny z modelu společně, a prediktory, které jsou zahrnuty ve všech modelech. Poskytuje různá apriorní rozdělení pro modely v parametru mprior() a pro parametr g, který řídí smršťování regresních koeficientů směrem k nule, v parametru gprior(). Podporuje také faktorové proměnné a operátory časových řad a pomocí volby heredity() poskytuje několik způsobů, jak při odhadu zpracovávat interakce.

Podporováno je mnoho funkcí postestimace, které zahrnují i některé standardní funkce bayesovské postestimace.

Podívejme se, jak to funguje

-> Výčet modelů -> Věrohodné intervaly -> Vlivné modely -> Důležité prediktory ->Rozdělení velikosti modelu -> Posteriorní rozdělení koeficientů -> Předpověď BMA -> Informativní modelové priory -> Prediktivní výkonnost pomocí LPS -> Čištění

Níže uvádíme prohlídku některých funkcí bma s využitím souboru dat. Tato datová sada obsahuje n=200

�=200  pozorování, p=10 �=10 ortogonální prediktory a výsledek y generovaný jako

y=0.5+1.2×x2+5×x10+ϵ

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

kde ϵ∼N(0,1)

�∼�(0,1)  je standardní normální chybový člen.

. webuse bmaintro
(Simulovaná data pro příklad BMA)

. summarize

Výčet modelů

K dosazení lineární regrese BMA pro y na x1 až x10 použijeme bmaregress.

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

S 10 prediktory a výchozí pevnou (referenční) hodnotou g-prior použije bmaregress výčet modelů a prozkoumá všech 210=1,024

210=1,024  možných modelů. Výchozí priorita modelu je beta-binomická s tvarovými parametry 1, která je rovnoměrná v celé velikosti modelu. A priori náš model BMA předpokládal malé smrštění regresních koeficientů směrem k nule, protože g/(1+g)=0.9950 �/(1+�)=0.9950  je blízké 1.

bmaregress identifikoval x2 a x10 jako velmi důležité prediktory – jejich posteriorní pravděpodobnosti zahrnutí (PIP) jsou v podstatě 1. Posteriorní střední odhady jejich regresních koeficientů jsou 1,2 (zaokrouhleno), resp. 5,1 a jsou velmi blízké skutečným hodnotám 1,2 a 5. Všechny ostatní prediktory mají mnohem nižší PIP a odhady koeficientů blízké nule. Naše zjištění BMA jsou v souladu se skutečným DGM.

Uložíme si aktuální výsledky odhadu pro pozdější použití. Stejně jako u ostatních bayesovských příkazů musíme nejprve uložit výsledky simulace, abychom mohli použít příkaz estimates store k uložení výsledků odhadu.

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

. estimates store bmareg

Věrohodné intervaly

bmaregress ve výchozím nastavení neuvádí věrohodné intervaly (CrI), protože vyžaduje potenciálně časově náročnou simulaci posteriorního vzorku parametrů modelu. Můžeme je však vypočítat po odhadu.

Nejprve použijeme bmacoefsample k vygenerování MCMC vzorku z posteriorního rozdělení parametrů modelu, které zahrnuje regresní koeficienty. Poté použijeme existující příkaz bayesstats summary k výpočtu posteriorních souhrnů vygenerovaného vzorku MCMC.

. bmacoefsample, rseed(18)

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

. bayesstats summary

Posterior summary statistics                      MCMC sample size =    10,000

95% CrI pro koeficient x2 je [1,05, 1,34] a pro x10 je [4,9, 5,26], což je v souladu s naším DGM.

Vlivné modely

Odhady koeficientů BMA jsou zprůměrovány pro 1 024 modelů. Je důležité prozkoumat, které modely mají vliv. Ke zkoumání PMP můžeme použít modely 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

Není překvapením, že model, který obsahuje x2 i x10, má nejvyšší PMP, přibližně 63 %. Ve skutečnosti dva nejlepší modely odpovídají kumulativnímu PMP přibližně 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

K vykreslení kumulativního PMP můžeme použít bmagraph pmp.

. bmagraph pmp, cumulative

Příkaz ve výchozím nastavení vykreslí prvních 100 modelů, ale pokud chcete zobrazit více modelů, můžete zadat volbu top().

Důležité prediktory

Důležitost prediktorů můžeme vizuálně prozkoumat pomocí bmagraph varmap, který vytvoří mapu zahrnutí proměnných.

. bmagraph varmap

Computing model probabilities ...

Modely x2 a x10 jsou zahrnuty ve všech 100 nejlepších modelech podle PMP. Jejich koeficienty jsou ve všech modelech kladné.

Rozdělení velikosti modelu

Složitost našeho modelu BMA můžeme prozkoumat pomocí bmastats msize a bmagraph msize k prozkoumání předběžného a následného rozdělení velikosti modelu.

. bmastats msize

Model-size summary

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

Note: Frequency summaries not available.

. bmagraph msize

Předběžné rozdělení velikosti modelu je rovnoměrné pro velikost modelu. Posteriorní rozdělení velikosti modelu je zkosené doleva s modem přibližně 2. Prioritní průměrná velikost modelu je 5 a posteriorní 2,48. Na základě pozorovaných údajů má tedy BMA tendenci upřednostňovat menší modely se zhruba dvěma prediktory v průměru ve srovnání s naším apriorním předpokladem.

Posteriorní rozdělení koeficientů

Pomocí bmagraph coefdensity můžeme vykreslit posteriorní rozdělení regresních koeficientů.

. bmagraph coefdensity {x2} {x3}, combine

Tato rozdělení jsou směsí bodové hmotnosti v nule odpovídající pravděpodobnosti nezařazení prediktoru a spojité hustoty podmíněné zařazením. Pro koeficient x2 je pravděpodobnost nezařazení zanedbatelná, takže jeho posteriorní rozdělení je v podstatě spojitá, poměrně symetrická hustota se středem v 1,2. Pro koeficient x3 je kromě podmíněné spojité hustoty v bodě nula červená svislá čára odpovídající posteriorní pravděpodobnosti nezařazení zhruba 1-,2 = 0,8.

Předpověď BMA

bmapredict vytváří různé předpovědi BMA. Vypočítejme například posteriorní predikční průměry.

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

Můžeme také vypočítat prediktivní CrI, abychom odhadli nejistotu našich předpovědí. (To není možné u mnoha tradičních metod výběru modelu.) Všimněte si, že tyto prediktivní CrI zahrnují také nejistotu modelu. Abychom mohli vypočítat prediktivní CrI, musíme nejprve uložit posteriorní vzorek parametrů modelu 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 ...

Shrnujeme předpokládané výsledky:

.  summarize y pmean cri*

Uváděné průměry za pozorování pro prediktivní průměry a dolní a horní 95% prediktivní hranice CrI vypadají vzhledem k pozorovanému výsledku y přiměřeně.

Informativní modelové priory

Jednou ze silných stránek metody BMA je možnost zahrnout předchozí informace o modelech a prediktorech. To nám umožňuje zkoumat citlivost výsledků na různé předpoklady o významu modelů a prediktorů. Předpokládejme, že máme spolehlivou informaci, že všechny prediktory kromě x2 a x10 mají menší pravděpodobnost vztahu k výsledku. Můžeme zadat například prioritu binomického modelu s nízkou prioritou pravděpodobnosti zahrnutí těchto prediktorů.

Abychom mohli později porovnat výkonnost modelů BMA mimo vzorek, rozdělili jsme náhodně vzorek na dvě části a model BMA jsme použili v prvním podvzorku. Výsledky z tohoto informativnějšího modelu BMA rovněž uložíme.

. 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

Důsledkem tohoto modelu prior je, že PIP všech nedůležitých prediktorů je nyní menší než 2 %.

Všimněte si, že při výběru jednoho modelu v podstatě fitujeme model BMA s velmi silnou prioritou, že všechny vybrané prediktory musí být zahrnuty s pravděpodobností 1. Můžeme například donutit bmaregress, aby do modelu zahrnul všechny proměnné:

. 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

Prediktivní výkonnost pomocí LPS

Pro porovnání uvažovaných modelů BMA jsme znovu sestavili náš výchozí model BMA s použitím prvního podvzorku a uložili výsledky odhadu.

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

Predikční výkonnost našich modelů BMA mimo vzorek můžeme porovnat pomocí nástroje bmastats lps, který vypočítá logaritmické predikční skóre (LPS) pro pozorování mimo vzorek.

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

Informativnější model bmareg_inf má o něco menší průměrnou hodnotu LPS, ale souhrny LPS jsou u všech modelů velmi podobné. Postup porovnání výkonnosti modelů BMA pomocí křížové validace naleznete v [BMA] bmastats lps.

Čištění

Během analýzy jsme vytvořili několik souborů dat. Již je nepotřebujeme, a proto je na závěr vymažeme. Můžete se však rozhodnout, že si je ponecháte, zejména pokud odpovídají potenciálně časově náročné závěrečné analýze.

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