3.3 interaction

reference:
- ARM chapter 03, github

library(tidyverse)
library(knitr)
# 3.3 interactions
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(shinystan)

interaction

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_interaction2.stan


data {
  int N;
  vector[N] kid_score;
  vector[N] mom_iq;
  vector[N] mom_hs;}
transformed data {           // interaction
  vector[N] inter;
  inter = mom_hs .* mom_iq;
}
parameters {
  vector[4] beta;
  real sigma;
}
model {
  sigma ~ cauchy(0, 2.5);
  kid_score ~ normal(beta[1] + beta[2] * mom_hs + beta[3] * mom_iq
                     + beta[4] * inter, sigma);
}

fit

‘kidiq_interaction.stan’ has an assignment operator in it, <-, in the definition of inter, which provokes a non-fatal warning

‘kidiq_interaction2.stan’ replaces <- with = in the definition of inter

### Model: lm(kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
data.list <- c("N", "kid_score", "mom_hs", "mom_iq")
kidiq_interaction <- stan(file = 'kidiq_interaction2.stan',
                          data = data.list,
                          iter = 1000, chains = 4)
kidiq_interaction
## Inference for Stan model: kidiq_interaction2.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##             mean se_mean    sd     2.5%      25%      50%      75%
## beta[1]   -11.23    0.62 13.35   -37.35   -20.15   -11.28    -2.05
## beta[2]    50.95    0.67 15.16    22.07    40.55    50.90    61.14
## beta[3]     0.97    0.01  0.14     0.69     0.87     0.96     1.06
## beta[4]    -0.48    0.01  0.16    -0.80    -0.59    -0.48    -0.37
## sigma      17.99    0.02  0.64    16.80    17.56    17.97    18.36
## lp__    -1472.30    0.06  1.60 -1476.19 -1473.22 -1471.91 -1471.11
##            97.5% n_eff Rhat
## beta[1]    14.12   469 1.01
## beta[2]    81.95   510 1.00
## beta[3]     1.24   470 1.01
## beta[4]    -0.17   510 1.00
## sigma      19.39   801 1.00
## lp__    -1470.25   608 1.00
## 
## Samples were drawn using NUTS(diag_e) at Wed Aug 23 06:35:21 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_interaction)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)