3.4 stat inference

reference:
- ARM chapter 03, github

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

ARM 3.4

multiple predictors

data

### 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 (kidiq_multi_preds.stan): 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.70    0.29 5.77    14.20    21.69    25.66    29.56    37.30
## beta[2]     5.89    0.09 2.26     1.51     4.37     5.86     7.43    10.23
## beta[3]     0.56    0.00 0.06     0.45     0.53     0.56     0.61     0.68
## sigma      18.14    0.02 0.62    16.94    17.70    18.13    18.58    19.30
## lp__    -1476.26    0.07 1.36 -1479.75 -1476.92 -1475.99 -1475.28 -1474.51
##         n_eff Rhat
## beta[1]   392 1.01
## beta[2]   610 1.00
## beta[3]   392 1.01
## sigma     731 1.00
## lp__      393 1.00
## 
## Samples were drawn using NUTS(diag_e) at Wed Aug 23 06:39:00 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)