Stejně jako u jiných příkazů Staty můžete samozřejmě místo zadávání příkazů použít rozhraní „ukaž a klikni“.

Bayesovské modely panelových dat nejsou určeny pouze pro spojité výsledky. Stejně snadno je můžete zadat i pro binární výsledky

. bayes: xtprobit y x1 x2

nebo pro výsledky počítání

. bayes: xtpoisson y x1 x2

Nebo použijte některý z osmi modelů panelových dat, které podporují předponu bayes, včetně nového modelu multinomického logitu pro panelová data.

Podívejme se, jak to funguje

• Odhad
• Konvergence MCMC
• Vlastní priory
• Posteriorní rozdělení panelových nebo skupinových efektů
• Bayesovské předpovědi
• Posteriorní prediktivní kontroly – shoda modelu
• Úklid

Odhad

Vezměme si podskupinu údajů z Národního longitudinálního výzkumu mladých žen ve věku 14 až 24 let v roce 1968, které žily na jihu. Budeme modelovat logaritmus mezd jako funkci vzdělání jedince, grade, jeho pracovní zkušenosti, ttl_exp, která vstupuje do modelu kvadraticky, a toho, zda žije ve standardní metropolitní oblasti, not_smsa. Na jednoho jedince připadá více pozorování, označených id.

Pro zohlednění individuálních efektů jsme použili bayesovský lineární model panelových dat. Chceme také vypočítat bayesovské předpovědi logaritmických mezd a porovnat individuální (panelové) efekty.

. webuse nlswork6
(Subsample of 1986 National Longitudinal Survey of Young Women)

Abychom ušetřili čas, provedeme namísto výchozího vzorku 10 000 vzorků metodou Markovova řetězce Monte Carlo (MCMC) vzorek 1 000 vzorků. Kvůli reprodukovatelnosti zadáváme také semeno náhodného čísla.

. bayes, mcmcsize(1000) rseed(17): xtreg ln_wage grade c.ttl_exp##c.ttl_exp i.not_smsa
note: Gibbs sampling is used for regression coefficients and variance
      components.

Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 1000 .........1000 done

Model summary

Likelihood: 
  ln_wage ~ normal(xb_ln_wage,{sigma2})

Priors: 
                {ln_wage:grade} ~ normal(0,10000)                          (1)
              {ln_wage:ttl_exp} ~ normal(0,10000)                          (1)
  {ln_wage:c.ttl_exp#c.ttl_exp} ~ normal(0,10000)                          (1)
           {ln_wage:1.not_smsa} ~ normal(0,10000)                          (1)
                {ln_wage:_cons} ~ normal(0,10000)                          (1)
                        {U[id]} ~ normal(0,{var_U})                        (1)
                       {sigma2} ~ igamma(0.01,0.01)

Hyperprior: 
  {var_U} ~ igamma(0.01,0.01)

(1) Parameters are elements of the linear form xb_ln_wage.

Bayesian RE normal regression                    MCMC iterations  =      3,500
Metropolis–Hastings and Gibbs sampling           Burn-in          =      2,500
                                                 MCMC sample size =      1,000
Group variable: id                               Number of groups =        831
                                                 Obs per group:
                                                              min =          1
                                                              avg =        1.4
                                                              max =          5
                                                 Number of obs    =      1,174
                                                 Acceptance rate  =       .796
                                                 Efficiency:  min =     .02411
                                                              avg =     .07145
Log marginal-likelihood                                       max =      .1676
 


         
 
                                                               
         
 
                                               Equal-tailed
         
 
      Mean   Std. dev.     MSCE     Median  [95% cred. interval]

ln_wage
 
                                                                     
          grade
 
  .0696248   .0052371    .00082    .069638   .0591421     .07999
        ttl_exp
 
  .0404646   .0070371   .000544   .0404901    .026454    .053855
               
 
                                                                
     c.ttl_exp#
 
                                                               
      c.ttl_exp
 
 -.0005053   .0003978   .000039  -.0005132  -.0012609   .0003009
               
 
                                                                
     1.not_smsa
 
 -.1502314   .0251494   .002755   -.150213   -.200979  -.1014476
          _cons
 
  .5646936   .0658804     .0099   .5680976   .4377888   .6957415

          var_U
 
  .0801363   .0073138   .001258   .0795163   .0674855   .0959143
         sigma2
 
  .0673308   .0046724   .000952   .0671238   .0590742   .0771401

Posteriorní střední hodnota koeficientu grade je kladná a dosahuje 7 %. Teorie naznačuje, že mzdy se zvyšují se zkušenostmi, ale jejich nárůst se s časem snižuje. To znamená, že koeficient ttl_exp by měl být kladný a koeficient c.ttl_exp#c.ttl_exp by měl být záporný. To pozorujeme v našich datech. A konečně, bydlení mimo velká městská centra má na mzdy negativní vliv.

Bayesovské modelování vyžaduje zadání priorů pro všechny parametry modelu. Stejně jako u jiných příkazů bayes i u příkazu bayes: xtreg jsou pro větší pohodlí k dispozici výchozí priory. Měli byste si prohlédnout specifikace priorit a zadat vlastní.

Konvergence MCMC

bayes: xtreg používá k získání výsledků metodu MCMC. Před další analýzou je třeba ověřit jeho konvergenci. To můžete provést graficky nebo můžete vypočítat Gelmanovu-Rubinovu statistiku konvergence pomocí více řetězců.

Posuďme konvergenci MCMC vizuálně například pro koeficient grade.

. bayesgraph diagnostics {ln_wage:grade}

V grafu stop není patrný žádný trend. Autokorelace s časem klesá. Nemáme důvod podezřívat nekonvergenci pro grade, ale konvergenci je třeba ověřit pro všechny parametry modelu.

Výsledky MCMC a odhadu si uložme pro pozdější použití.

. bayes, saving(bxtregsim)
note: file bxtregsim.dta saved.

. estimates store bxtreg

Vlastní priory

U bayesovských modelů můžete chtít použít vlastní priory. Tyto priory často pocházejí z historických dat. Například na základě předchozích studií může být rozumné předpokládat, že koeficient pro grade se pohybuje v rozmezí 0 až 25. Můžeme tedy pro něj uvažovat rovnoměrnou prioritu na (0,25). K zadání vlastních priorů můžeme použít bayesovu volbu prior().

. bayes, mcmcsize(1000) rseed(17) prior({ln_wage:grade}, uniform(0,25)):
> xtreg ln_wage grade c.ttl_exp##c.ttl_exp i.not_smsa
note: Gibbs sampling is used for variance components.

  
Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 1000 .........1000 done


Model summary

Bayesian RE normal regression                    MCMC iterations  =      3,500
Metropolis–Hastings and Gibbs sampling           Burn-in          =      2,500
                                                 MCMC sample size =      1,000
Group variable: id                               Number of groups =        831
                                                 Obs per group:
                                                              min =          1
                                                              avg =        1.4
                                                              max =          5
                                                 Number of obs    =      1,174
                                                 Acceptance rate  =      .6181
                                                 Efficiency:  min =    .009861
                                                              avg =     .01259
Log marginal-likelihood                                       max =     .01751



         
 
                                                               
         
 
                                               Equal-tailed
         
 
      Mean   Std. dev.     MSCE     Median  [95% cred. interval]

ln_wage
 
                                                                     
          grade
 
  .0691558   .0063225   .002013   .0692614   .0567288   .0819581
        ttl_exp
 
  .0389882   .0080916   .002355   .0390646   .0228502   .0530886
               
 
                                                                
     c.ttl_exp#
 
                                                               
      c.ttl_exp
 
 -.0004548   .0004804   .000153  -.0004778  -.0013733   .0005973
               
 
                                                                
     1.not_smsa
 
 -.1452854   .0241483   .006649  -.1450175  -.1934294  -.0947368
          _cons
 
  .5751218   .0727991   .022234   .5795395   .4176046   .7018132

          var_U
 
  .0785815   .0082846    .00198   .0784331   .0615956   .0956798
         sigma2
 
  .0688682   .0059815   .001537   .0684258   .0583475   .0823034

Náš vlastní prior výsledky příliš nezměnil.

Posteriorní rozdělení panelových nebo skupinových efektů

Můžeme se zajímat o vyvozování závěrů o (panelových) účincích na jednotlivce. Bayesovské modelování k tomu nabízí přirozený způsob. Na rozdíl od klasických modelů panelových dat s náhodnými efekty odhadují bayesovské modely panelových dat náhodné (panelové) efekty společně se všemi ostatními parametry modelu. Každý náhodný účinek je tak reprezentován celým MCMC vzorkem z jeho posteriorního rozdělení. Tento vzorek lze použít pro inferenci, například pro porovnání panelových nebo skupinových efektů. (Představte si porovnání výkonnosti různých společností nebo nemocnic.)

Vraťme se k našemu původnímu modelu, který používal výchozí neinformativní priory.

. estimates restore bxtreg
(results bxtreg are active now)

Vykresleme posteriorní rozdělení pro prvních devět individuálních efektů.

. bayesgraph histogram {U[1/9]}, byparm normal

Výše uvedené individuální vlivy představují posuny nebo posuny oproti průměrné logaritmické mzdě. Mezi mzdami jednotlivců rozhodně existují rozdíly. Například vidíme, že mzda jedince 8 je vyšší než mzda jedince 5. Mezi platy jednotlivců tedy stále existují rozdíly, které nejsou zohledněny prediktory modelu.

Bayesovské předpovědi

V rámci bayesovské metody lze vypočítat předpovědi a jejich nejistoty bez asymptotických předpokladů. To je možné, protože získané předpovědi jsou vzorky z „přesného“ posteriorního predikčního rozdělení nových dat vzhledem k pozorovaným datům. Nepředpokládá se, že posteriorní predikční rozdělení je asymptoticky normální.

Vypočítejme posteriorní průměry predikovaných logaritmických mezd spolu s jejich 95% věrohodným intervalem

. bayespredict pmean, mean 

Computing predictions ...

. bayespredict cri_l cri_u, cri

Computing predictions ...

Uveďme výsledky prvních deseti pozorování.

. list ln_wage pmean cri_l cri_u in 1/10


     
  ln_wage      pmean       cri_l       cri_u 
  1.
 1.543923    1.55197    .8831815    2.184215
  2.
 1.815738   1.645933    1.006488    2.284888 
  3.
 1.870532   1.718343    1.112855    2.360129 
  4.
 2.340405    2.16819    1.523068    2.747879 
  5.
 1.545327   1.629953    1.092535    2.221188 
  6.
 1.594335   1.687503    1.089151    2.268714 
  7.
 1.848619   1.849651    1.266758    2.408907 
  8.
 1.757637    1.80307    1.140366    2.467546 
  9.
 2.346746   2.224904    1.529079    2.833333 
 10.
 3.539868   2.826903    2.270285    3.511156 

Předpokládané posteriorní střední hodnoty jsou blízké pozorovaným hodnotám, snad s výjimkou pozorování 10. V následující části si ukážeme, jak formálněji ověřit, zda model dobře sedí. K tomu však musíme nejprve uložit všechny předpovědi MCMC.

. bayespredict {_ysim1}, saving(bxtregpred)

Computing predictions ...

file bxtregpred.dta saved.
file bxtregpred.ster saved.

bxtregpred.dta obsahuje vzorek hodnot MCMC pro každé pozorování. Posteriorní průměry (pmean), které jsme předpověděli dříve, jsou pro každé pozorování průměrem ze vzorku MCMC.

Datové soubory predikcí MCMC jsou obvykle velké, takže můžete zvážit jejich uložení pouze v případě potřeby; podrobnosti viz [BAYES] bayespredict.

Posteriorní prediktivní kontroly – shoda modelu

Další výhodou bayesovské analýzy modelů panelových dat jsou formalizované posteriorní prediktivní kontroly pro kontrolu vhodnosti modelu, které u klasické frekvenční analýzy nejsou k dispozici.

K provádění kontrol správnosti shody modelu můžeme použít posteriorní vzorky pro predikovaný výsledek, vzorek predikce MCMC. Porovnejme například minimální a maximální statistiky z predikčních vzorků MCMC se statistikami pozorovanými v datech. Tyto statistiky popisují chvosty rozdělení dat. Místo maxima a minima (nebo jako doplněk k nim) můžete použít jakoukoli jinou statistiku.
bayesstats ppvalues provádí posteriorní prediktivní kontroly výpočtem tzv. posteriorních prediktivních p-hodnot (PPP). PPP popisují, jak často byly statistiky z predikčního vzorku MCMC stejně extrémní nebo extrémnější než statistiky z pozorovaného vzorku; viz [BAYES] bayesstats ppvalues.

. bayesstats ppvalues (pmin:@min({_ysim1})) (pmax:@max({_ysim1})) using bxtregpred

Posterior predictive summary   MCMC sample size =     1,000



              T 
 
      Mean   Std. dev.  E(T_obs)  P(T>=T_obs)

           pmin 
 
  .0299837   .1652143   .0176546         .565
           pmax 
 
  3.515637   .2291669    3.85025         .085

Note: P(T>=T_obs) close to 0 or 1 indicates lack of fit.

PPP pro minimum, 0,565, se neblíží hodnotě 0 ani 1, takže model s ohledem na tuto statistiku dobře vyhovuje. PPP pro maximum, 0,085, se však blíží 0, a proto naznačuje špatnou shodu s touto statistikou.

Prozkoumáme-li data, objevíme shluk pozorování s hodnotami mezi 3,5 a 4, která vypadají jako odlehlá. Možná budeme muset náš model přehodnotit a problém odlehlých hodnot vyřešit.

Úklid

Po analýze můžeme odstranit soubory vytvořené pomocí bayes: var a bayespredict, které již nejsou potřeba.

. erase bxtregsim.dta
. erase bxtregpred.dta
. erase bxtregpred.ster