3.7 stat prediction

reference:
- ARM chapter 03, github
- Using Stan and ShinyStan for posterior predictive checking
- Bayesian Basics, Michael Clark

library(tidyverse)
library(knitr)
# 3.4 stat inference
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(shinystan)

ARM 3.7

data

### Data
source("kidiq_prediction.data.R", echo = TRUE)
## 
## > N <- 434
## 
## > kid_score <- c(65, 98, 85, 83, 115, 98, 69, 106, 102, 
## +     95, 91, 58, 84, 78, 102, 110, 102, 99, 105, 101, 102, 115, 
## +     100, 87, 99, 96, 72,  .... [TRUNCATED] 
## 
## > mom_hs <- c(1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 
## +     1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 
## +     1, 0, 1, 1, 1, 1, 1, 1, 1, .... [TRUNCATED] 
## 
## > mom_iq <- c(121.117528602603, 89.3618817100663, 115.443164881725, 
## +     99.4496394360723, 92.7457099982118, 107.901837758501, 138.893106071162, 
## +  .... [TRUNCATED] 
## 
## > mom_hs_new <- 1
## 
## > mom_iq_new <- 100

model

kidiq_prediction2.stan


data {
  int N;
  vector[N] kid_score;
  vector[N] mom_iq;
  vector[N] mom_hs;
  real mom_hs_new;           // for prediction
  real mom_iq_new;
}
parameters {
  vector[3] beta;
  real sigma;
}
model {
  sigma ~ cauchy(0, 2.5);
  kid_score ~ normal(beta[1] + beta[2] * mom_hs + beta[3] * mom_iq, sigma);
}
generated quantities {       // prediction
  vector[N] kid_score_pred;
  for (n in 1:N) 
    kid_score_pred[n] = normal_rng(beta[1] + beta[2] * mom_hs_new
                                   + beta[3] * mom_iq_new, sigma);
}

fit

### Model (kidiq_prediction2.stan): kid_score ~ mom_hs + mom_iq
data.list <- c("N", "kid_score", "mom_hs", "mom_iq", "mom_hs_new", "mom_iq_new")
kidiq_prediction <- stan(file='kidiq_prediction2.stan', data=data.list,
                         iter=500, chains=4)

shinystan

To check prediction in shinystan, under Diagnose click on PPCheck, then for \(T\) choose “kid_score” from the global environment and for \(T_{rep}\) choose “kid_score_pred” from the model

kidiq_prediction.sso <- launch_shinystan(kidiq_prediction)
## 
## Launching ShinyStan interface... for large models this  may take some time.
## 
## Listening on http://127.0.0.1:5801

saveRDS

saveRDS(kidiq_prediction, "3.7.kidiq_prediction.rds")

readRDS

# restore saved object
# kidiq_prediction2 <- readRDS("3.7.kidiq_prediction.rds")
# identical(kidiq_prediction, kidiq_prediction2)
# print(kidiq_prediction2, pars = c("beta", "sigma", "lp__"))
# launch_shinystan(kidiq_prediction2)