3.5 graphical displays

reference:
- ARM chapter 03, github
- saveRDS, readRDS, rdocumentation
- saveRDS, readRDS, R-Bloggers

library(tidyverse)
library(knitr)
# 3.5 graph displays
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(grid)
library(shinystan)

ARM 3.5

one predictor

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

kidscore_momiq.stan


data {
  int N;
  vector[N] kid_score;
  vector[N] mom_iq;
}
parameters {
  vector[2] beta;
  real sigma;
}
model {
  sigma ~ cauchy(0, 2.5);
  kid_score ~ normal(beta[1] + beta[2] * mom_iq, sigma);
}

fit

### Regression line as a function of one input variable
data.list.2 <- c("N", "kid_score", "mom_iq")
stanfit.2 <- stan(file='kidscore_momiq.stan', data=data.list.2,
                  iter=500, chains=4)
print(stanfit.2, pars = c("beta", "sigma", "lp__"))
## Inference for Stan model: kidscore_momiq.
## 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.43    0.33 6.39    13.18    21.20    25.41    29.68    38.23
## beta[2]     0.61    0.00 0.06     0.49     0.57     0.61     0.66     0.73
## sigma      18.27    0.03 0.61    17.14    17.83    18.26    18.69    19.50
## lp__    -1479.43    0.07 1.29 -1482.96 -1480.02 -1479.11 -1478.50 -1477.98
##         n_eff Rhat
## beta[1]   383 1.01
## beta[2]   383 1.01
## sigma     341 1.01
## lp__      329 1.02
## 
## Samples were drawn using NUTS(diag_e) at Mon Aug 21 04:13:10 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(stanfit.2)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)