Posterior predictive model checks

Bayesian Regression Modeling with rstanarm

Jake Thompson

Psychometrician, ATLAS, University of Kansas

Using posterior distributions

stan_model <- stan_glm(kid_score ~ mom_iq, data = kidiq)

spread_draws(stan_model, `(Intercept)`, mom_iq) %>%
  select(-.draw)

# A tibble: 4,000 x 4
   .chain .iteration `(Intercept)` mom_iq
    <int>      <int>         <dbl>  <dbl>
 1      1          1          19.9  0.654
 2      1          2          20.7  0.643
 3      1          3          27.2  0.604
 4      1          4          24.9  0.613
 5      1          5          26.4  0.610
 6      1          6          25.2  0.619
 7      1          7          17.8  0.702
# ... with 3,993 more rows

Posterior predictions

predictions <- posterior_linpred(stan_model)

predictions[1:10, 1:5]

iterations         1        2        3        4        5
      [1,] 100.18694 79.04791 96.40964 85.76310 81.30045
      [2,] 100.24843 82.00786 96.98905 87.80231 83.95155
      [3,] 100.85608 81.13109 97.33146 87.39709 83.23295
      [4,] 102.31392 80.81881 98.47300 87.64712 83.10930
      [5,]  97.25617 81.18278 94.38404 86.28879 82.89553
      [6,] 100.86263 79.89830 97.11655 86.55800 82.13223
      [7,]  99.36166 81.10329 96.09910 86.90339 83.04887
      [8,] 101.13487 80.97878 97.53321 87.38173 83.12658
      [9,]  98.72686 79.97596 95.37629 85.93252 81.97403
     [10,] 100.22835 81.04603 96.80069 87.13964 83.09007

predictions <- posterior_linpred(stan_model)

# First replication
iter1 <- predictions[1,]
# Second replication
iter2 <- predictions[2,]

# Data summaries
summary(kidiq$kid_score)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   20.0    74.0    90.0    86.8   102.0   144.0

summary(iter1)
summary(iter2)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  68.54   79.86   85.80   87.14   93.74  112.12

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  70.05   80.19   85.51   86.71   92.62  109.08

Comparing single scores

predictions <- posterior_linpred(stan_model)
kidiq$kid_score[24]
summary(predictions[, 24])

87

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  83.34   86.17   86.77   86.75   87.34   90.23

kidiq$kid_score[185]
summary(predictions[, 185])

111

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  82.81   85.65   86.25   86.24   86.83   89.69

Let's practice

Bayesian Regression Modeling with rstanarm