3.2 multiple predictors

reference:
- ARM chapter 03, github

library(tidyverse)
library(knitr)
# 3.2 multiple predictors.R
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(shinystan)

two predictors

data

source("kidiq.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_multi_preds.stan


data {
  int N;
  vector[N] kid_score;
  vector[N] mom_iq;
  vector[N] mom_hs;
}
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);
}

fit

### Model: kid_score ~ mom_hs + mom_iq
data.list <- c("N", "kid_score", "mom_hs", "mom_iq")
kidiq_multi_preds <- stan(file = 'kidiq_multi_preds.stan',
                          data = data.list,
                          iter = 500, chains = 4)
kidiq_multi_preds
## Inference for Stan model: kidiq_multi_preds.
## 4 chains, each with iter=500; warmup=250; thin=1; 
## post-warmup draws per chain=250, total post-warmup draws=1000.
## 
##             mean se_mean   sd     2.5%      25%      50%      75%    97.5%
## beta[1]    25.87    0.26 5.70    15.10    21.84    25.71    29.69    37.08
## beta[2]     6.01    0.10 2.38     1.21     4.37     6.02     7.55    10.79
## beta[3]     0.56    0.00 0.06     0.45     0.52     0.56     0.60     0.68
## sigma      18.10    0.02 0.64    16.98    17.67    18.05    18.50    19.45
## lp__    -1476.38    0.07 1.44 -1479.88 -1477.09 -1476.07 -1475.29 -1474.50
##         n_eff Rhat
## beta[1]   475    1
## beta[2]   585    1
## beta[3]   475    1
## sigma     747    1
## lp__      425    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Aug 23 06:29:43 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
plot(kidiq_multi_preds)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)