Bayesian Data Analysis 2 - Bayesian Hierarchical Models

Survey weights, the right way!

This is an example of how to use survey design weights with linear mixed models. Here is an example of using weights that sum to the strata-specific sample size. First, I make an id that is the strata, then I follow the first method described by Carle (2009).

eclsk.sub$sampleid<-eclsk.sub$w1c0str
#within each sampling unit, sum the weights
wts<-tapply(eclsk.sub$w1c0,eclsk.sub$sampleid,sum)
#make a data frame from this
wts<-data.frame(id=names(unlist(wts)), wt=unlist(wts))
#get the unique sampling location ids'
t1<-as.data.frame(table(eclsk.sub$sampleid))
#put all of this into a data set
wts2<-data.frame(ids=wts$id, sumwt=wts$wt, jn=t1$Freq)
head(wts2)

##   ids    sumwt  jn
## 1   1 63224.28 293
## 2   2 15747.17  67
## 3   3 23367.05 114
## 4   4 28346.22 129
## 5   5 36348.00 121
## 6   6 34575.17 169

#merge all of this back to the original data file
eclsk.sub<-merge(eclsk.sub, wts2, by.x="sampleid", by.y="ids", all.x=T)

eclsk.sub$swts<-eclsk.sub$w1c0*(eclsk.sub$jn/eclsk.sub$sumwt)

#compare standardized weights to old person weights
hist(eclsk.sub$swts); summary(eclsk.sub$swts)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0561  0.7025  0.9494  1.0000  1.1910  4.1870      28

hist(eclsk.sub$w1c0); summary(eclsk.sub$w1c0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   172.5   240.7   255.6   302.8   958.8

Linear Mixed Model Regression Example

The First model will fit a model that considers the individual and family level variables and a random intercept only. One may state a reserach question such as:

“Do children living in households in poverty have lower height for age/sex z-scores than children living in households that are not in poverty, net of other demographic and socioeconomic factors?”

So the hypothesis of interest here involves us testing whether the effect of household poverty shows a $\beta$ coefficient different from zero.

If we were to write the symbolic form of this model it would be:

$$y_{ij} = \beta_{0j} + \sum {\beta_k x_{ik}} + e_{ij}$$ $$\beta_{0j} = \beta_0 + u_j$$ $$e_{ij} \sim N(0, \sigma^2_e)$$ $$u_j \sim N(0, \sigma^2_u)$$

In terms of priors, we need them for each fixed effect $\beta$ and for $\sigma^2$. We can a Normal prior for $\beta$ parameters, since they can be both positive or negative, but the null hypothesis is that they are 0. Here we will use half-cauchy priors for the residual variance and the between-intercept variance. and priors:

$\beta \sim N(0, 10)$

$\sigma_e \sim \text {half-Cauchy(5)}$

$\sigma_u \sim \text {half-Cauchy(5)}$

prior<-c(set_prior("normal(0,10)", class="b"),#prior for the beta's
         set_prior("cauchy(0,5)", class="sigma"),#prior for the residual std. deviation
         set_prior("cauchy(0,5)", class="sd", group="s1_id"))#prior for the random intercept std. deviation
fit.1<-brm(height_z1|weights(swts)~male+hisp+black+asian+nahn+other+lths+gths+single+notmar+urban+hhsize+pov+(1|s1_id)+(1|w1c0str:w1c0psu), data=eclsk.sub, family=gaussian(),
              prior = prior, 
              warmup=500, #burnin for 500 interations for each chain = 1000 burnin
              iter=1000, chains=2, #2*1000 =2000 - 1000 burnin = 1000 total iterations
              cores=2,seed = 1115, #number of computer cores, 1 per chain is good.
              save_model="~/Google Drive/dem7903_App_Hier/code/code16/fit_lmm_stan.txt")

## Compiling the C++ model

print(fit.1, digits=3)

##  Family: gaussian (identity) 
## Formula: height_z1 | weights(swts) ~ male + hisp + black + asian + nahn + other + lths + gths + single + notmar + urban + hhsize + pov + (1 | s1_id) + (1 | w1c0str:w1c0psu) 
##    Data: eclsk.sub (Number of observations: 10665) 
## Samples: 2 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 1000
##    WAIC: Not computed
##  
## Random Effects: 
## ~s1_id (Number of levels: 853) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)    0.152     0.016    0.121    0.181        298 1.009
## 
## ~w1c0str:w1c0psu (Number of levels: 234) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)    0.111     0.019    0.075    0.148        162 1.012
## 
## Fixed Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept    0.051     0.047   -0.036    0.143        806 1.003
## male         0.155     0.018    0.119    0.189       1000 0.999
## hisp        -0.046     0.031   -0.104    0.013        619 1.000
## black        0.375     0.034    0.309    0.439        768 0.999
## asian       -0.271     0.050   -0.365   -0.172       1000 1.003
## nahn         0.278     0.091    0.105    0.452       1000 1.000
## other        0.079     0.046   -0.005    0.170       1000 0.998
## lths         0.009     0.037   -0.069    0.079        761 1.001
## gths         0.046     0.026   -0.006    0.094        738 0.999
## single      -0.048     0.031   -0.109    0.010        807 0.999
## notmar      -0.015     0.031   -0.075    0.048        781 0.999
## urban        0.001     0.008   -0.013    0.018        758 0.999
## hhsize      -0.027     0.007   -0.042   -0.013       1000 0.998
## pov         -0.087     0.022   -0.131   -0.044        973 1.001
## 
## Family Specific Parameters: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma(height_z1)    0.956     0.007    0.942    0.969       1000 0.999
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

rstan_gg_options(fill = "skyblue", color = "skyblue4", pt_color = "red")
stanplot(fit.1,pars=c("sigma_height_z1","sd_s1_id_Intercept"), type="dens"  )+ggtitle("Posterior Density for Variances")

stanplot(fit.1,pars="sd_s1_id_Intercept", type="trace" )+ggtitle("Traceplot for Between school Intercept Variance")

stanplot(fit.1, pars = "^b", type = "plot", ci_level = 0.95)+ggtitle("95% Credible Intervals for fixed effects")

## ci_level: 0.95 (95% intervals)

## outer_level: 0.95 (95% intervals)

#We can do a hypothesis test for the poverty parameter
h0pov<-hypothesis(fit.1, "hisp<0", class = "b") #larger evidence ratio = more support
h0pov

## Hypothesis Tests for class b:
##          Estimate Est.Error l-95% CI u-95% CI Evid.Ratio 
## hisp < 0    -0.05      0.03     -Inf        0      13.71 
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

#extract the MCMC samples from the model
samps<-fit.1$fit

#calculate the bayesian p value for a parameter, in this case, Pr(b_black > 0)
#Remember, higher values mean you have more support for your hypothesis!
mean( samps@sim$samples[[1]]$b_black  > 0)

## [1] 1

In terms of our research question, we have evidence that the poverty effect is negative, with a credible interval that does not contain 0, so we can conclude that our research hypothesis has support.

Hierarchical Logistic Regression Example

The next model will fit a model that considers a binary response of childhood growth stunting. . A research question corresponding to this could be:

“Do children living in households in poverty have have higher odds of growth stunting than children living in households that are not in poverty, net of other demographic and socioeconomic factors?”

So the hypothesis of interest here involves us testing whether the effec of household poverty shows a $\beta$ coefficient different from zero.

If we were to write the symbolic form of this model it would be the logistic regression model: $$ln \frac {p_{ij}}{1-p{ij}} = \beta_{0j} + \sum {\beta_k x_{ik}}$$ $$\beta_{0j} = \beta_0 + u_j$$ $$u_j \sim N(0, \sigma^2_u)$$

In terms of priors, we need them for each fixed effect $\beta$ and for $\sigma^2_u$. We can a Normal prior for $\beta$ parameters, since they can be both positive or negative, but the null hypothesis is that they are 0. Here we will use half-cauchy priors for the between-intercept variance, and there is no prior for the residual variance, becuase in the logistic model, this is the constant $\frac{\pi^2}{3}$.

$\beta \sim N(0, 10)$

$\sigma_u \sim \text {half-Cauchy(5)}$

eclsk.sub$stunt<-ifelse(eclsk.sub$height_z1 < -1, 1,0)

#bayesian model
prior<-c(set_prior("normal(0,1)", class="b"),
         set_prior("cauchy(0,5)", class="sd", group="s1_id"))
#no residual sigma for the binomial
fit.b<-brm(stunt|weights(swts)~male+hisp+black+asian+nahn+other+lths+gths+single+notmar+urban+hhsize+pov+(1|s1_id)+(1|w1c0str:w1c0psu), data=eclsk.sub, family=bernoulli(),
              prior = prior, 
              warmup=500, 
              iter=1000, chains=2, 
              cores=2,seed = 1115,
              save_model="~/Google Drive/dem7903_App_Hier/code/code16/fit_glmm_b_stan.txt")
print(fit.b, digits=3)

##  Family: bernoulli (logit) 
## Formula: stunt | weights(swts) ~ male + hisp + black + asian + nahn + other + lths + gths + single + notmar + urban + hhsize + pov + (1 | s1_id) + (1 | w1c0str:w1c0psu) 
##    Data: eclsk.sub (Number of observations: 10665) 
## Samples: 2 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 1000
##    WAIC: Not computed
##  
## Random Effects: 
## ~s1_id (Number of levels: 853) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)    0.243     0.064    0.106    0.354        102 1.003
## 
## ~w1c0str:w1c0psu (Number of levels: 234) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)     0.15     0.063    0.019    0.267        105 1.003
## 
## Fixed Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept   -1.879     0.118   -2.120   -1.660       1000 0.999
## male        -0.345     0.053   -0.453   -0.245        841 1.000
## hisp         0.063     0.081   -0.100    0.220       1000 1.000
## black       -0.743     0.114   -0.964   -0.510       1000 1.000
## asian        0.629     0.125    0.378    0.867       1000 1.003
## nahn        -0.121     0.268   -0.681    0.374       1000 0.998
## other        0.115     0.123   -0.129    0.356       1000 0.998
## lths         0.025     0.101   -0.172    0.222        872 0.998
## gths        -0.103     0.074   -0.249    0.043       1000 1.001
## single      -0.057     0.093   -0.229    0.122        842 0.999
## notmar      -0.176     0.101   -0.375    0.017       1000 0.999
## urban       -0.023     0.020   -0.061    0.017        813 1.000
## hhsize       0.070     0.020    0.031    0.109       1000 0.999
## pov          0.191     0.066    0.065    0.322        911 0.999
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

rstan_gg_options(fill = "skyblue", color = "skyblue4", pt_color = "red")
stanplot(fit.b,pars="sd_s1_id_Intercept", type="dens" )+ggtitle("Posterior Density for Between school Intercept Variance")

stanplot(fit.b,pars="sd_s1_id_Intercept", type="trace" )+ggtitle("Traceplot for Between school Intercept Variance")

stanplot(fit.b, pars = "^b", type = "plot", ci_level = 0.95)+ggtitle("95% Credible Intervals for fixed effects")

h0pov<-hypothesis(fit.b, "pov>0")
h0pov

## Hypothesis Tests for class b:
##         Estimate Est.Error l-95% CI u-95% CI Evid.Ratio  
## pov > 0     0.19      0.07     0.08      Inf        499 *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

samps<-fit.b$fit
#Bayesian p value for the poverty effect, combining the samples from chains 1 and 2
mean(c(samps@sim$samples[[1]]$b_pov, samps@sim$samples[[2]]$b_pov) >0)

## [1] 0.9945

#Credible interval for the odds ratio
or_pov<-exp(c(samps@sim$samples[[1]]$b_pov, samps@sim$samples[[2]]$b_pov))
quantile(or_pov, p=c(0, .025, .05, .5, .95, .975, 1))

##        0%      2.5%        5%       50%       95%     97.5%      100% 
## 0.4160436 1.0602785 1.0807511 1.2108747 1.3586282 1.3823291 2.3214061

In terms of our research question, we have evidence that the poverty effect is positive, with a credible interval that does not contain 0, so we can conclude that our research hypothesis has support.

prior<-c(set_prior("normal(0,10)", class="b"),#prior for the beta's
         set_prior("cauchy(0,5)", class="sigma"),#prior for the residual std. deviation
         set_prior("cauchy(0,5)", class="sd", group="s1_id"),#prior for the random intercept std. deviation
         set_prior("lkj(2)", class = "cor", group="s1_id"))
fit.2d<-brm(height_z1|weights(swts)~male+hisp+black+asian+nahn+other+lths+gths+single+notmar+urban+hhsize+hhses1+pov+(pov|s1_id)+(1|w1c0str:w1c0psu), data=eclsk.sub, family=gaussian(),
              prior = prior, 
              warmup=500, #burnin for 500 interations for each chain = 1000 burnin
              iter=1000, chains=2, #2*1000 =2000 - 1000 burnin = 1000 total iterations
              cores=2,seed = 1115, #number of computer cores, 1 per chain is good.
           save_model="~/Google Drive/dem7903_App_Hier/code/code16/fit_lmm_2d_stan.txt"
              )
print(fit.2d, digits=3)

##  Family: gaussian (identity) 
## Formula: height_z1 | weights(swts) ~ male + hisp + black + asian + nahn + other + lths + gths + single + notmar + urban + hhsize + hhses1 + pov + (pov | s1_id) + (1 | w1c0str:w1c0psu) 
##    Data: eclsk.sub (Number of observations: 10665) 
## Samples: 2 chains, each with iter = 1000; warmup = 500; thin = 1; 
##          total post-warmup samples = 1000
##    WAIC: Not computed
##  
## Random Effects: 
## ~s1_id (Number of levels: 853) 
##                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sd(Intercept)         0.173     0.022    0.130    0.215        193 1.005
## sd(pov)               0.196     0.056    0.047    0.283         36 1.067
## cor(Intercept,pov)   -0.492     0.227   -0.774    0.101         54 1.044
## 
## ~w1c0str:w1c0psu (Number of levels: 234) 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept)    0.106     0.019    0.069    0.142        323    1
## 
## Fixed Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept    0.045     0.046   -0.044    0.128       1000 0.999
## male         0.156     0.019    0.119    0.192       1000 0.999
## hisp        -0.044     0.030   -0.106    0.013       1000 1.000
## black        0.381     0.036    0.312    0.449       1000 1.001
## asian       -0.281     0.052   -0.383   -0.178       1000 0.999
## nahn         0.290     0.096    0.104    0.477       1000 1.000
## other        0.079     0.051   -0.015    0.176       1000 0.999
## lths         0.019     0.039   -0.057    0.093       1000 1.000
## gths         0.025     0.028   -0.032    0.077       1000 0.998
## single      -0.042     0.032   -0.103    0.021       1000 0.998
## notmar      -0.015     0.031   -0.075    0.048       1000 0.999
## urban        0.002     0.008   -0.014    0.017       1000 0.998
## hhsize      -0.027     0.008   -0.043   -0.013       1000 0.999
## hhses1       0.034     0.018   -0.002    0.071       1000 0.998
## pov         -0.061     0.028   -0.115   -0.006       1000 0.999
## 
## Family Specific Parameters: 
##                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma(height_z1)    0.952     0.007    0.938    0.964       1000 1.002
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

rstan_gg_options(fill = "skyblue", color = "skyblue4", pt_color = "red")
stanplot(fit.2d,pars="sigma_height_z1", type="dens"  )+ggtitle("Posterior Density for Residual Variance")

stanplot(fit.2d,pars=c("sigma_height_z1","sd_s1_id_Intercept","sd_s1_id_pov"), type="dens" )+ggtitle("Posterior Density for Variances")

stanplot(fit.2d,pars="sd_s1_id_Intercept", type="trace" )+ggtitle("Traceplot for Between school Intercept Variance")

stanplot(fit.2d,pars="sd_s1_id_pov", type="trace" )+ggtitle("Traceplot for Between school poverty effect Variance")

stanplot(fit.2d,pars="cor_s1_id_Intercept_pov", type="dens" )+ggtitle("Posterior Density for correlation between school random slope and Intercept")

stanplot(fit.2d,pars="cor_s1_id_Intercept_pov", type="trace" )+ggtitle("Traceplot for correlation between school random slope and Intercept")

stanplot(fit.2d, pars = "^b", type = "plot", ci_level = 0.95)+ggtitle("95% Credible Intervals for fixed effects")