** Lab 1 656: NMMAPS **

** for the sake of simplicity, we will ignore the standard error of the betas **
** Goal: estimate the effect of air pollution for each region **

use "C:\Documents and Settings\ejohnson\My Documents\MultilevelModels\NMMAPS.pm.coef.dta", clear

** exploratory analysis **
sort region
scatter beta region, xlabel(1 2 3 4 5 6 7, valuelabel)
xtsum beta, i(region)


** Approach A: calculate each region's observed average coefficient on PM **
mean beta, over(region)
** store these as a variable using 
egen region_mean_beta = mean(beta), by(region)

** Approach B: calculate overall average for all regions **
mean beta
** store these as a variable using
egen mean_beta = mean(beta)

** To decide between A and B **
** ANOVA for differences among the means (according to region) **
oneway beta region

** A compromise...

** Approach C: weighted A and B (Empirical Bayes) **
gllamm beta, i(region) adapt nip(15)

* xtreg beta, re i(region) mle ** gives errors
* xtmixed beta || region:, mle ** gives errors

** empirical bayes estimates **
* re-run the gllamm 
matrix a = e(b)
gllamm beta, i(region) adapt nip(15)


gllapred nationalest, linpred
gllapred EBint, u 

gen EBest = EBintm1 + nationalest

scatter beta region_mean_beta mean_beta EBest region, xlabel(1 2 3 4 5 6 7, valuelabel) msymbol(o d X Oh)