Regression and Other Stories: KidIQ Bayes-R2 and LOO-R2

Linear regression and Bayes-R2 and LOO-R2. See Chapter 11 in Regression and Other Stories.

See also Gelman, Goodrich, Gabry, and Vehtari (2018). R-squared for Bayesian regression models. The American Statistician 73:307-309, doi:10.1080/00031305.2018.1549100

---

Load packages

library("rprojroot")
root<-has_file(".ROS-Examples-root")$make_fix_file()
library("rstanarm")
library("loo")
library("foreign")
library("ggplot2")
theme_set(bayesplot::theme_default(base_family = "sans"))
library(reshape2)
# for reproducability
SEED<-1507

Load children's test scores data

kidiq <- read.csv(root("KidIQ/data","kidiq.csv"))
head(kidiq)

  kid_score mom_hs    mom_iq mom_work mom_age
1        65      1 121.11753        4      27
2        98      1  89.36188        4      25
3        85      1 115.44316        4      27
4        83      1  99.44964        3      25
5       115      1  92.74571        4      27
6        98      0 107.90184        1      18

Compare different models

A single binary predictor

fit_1 <- stan_glm(kid_score ~ mom_hs, data=kidiq,
                  seed=SEED, refresh=0)
loo_1 <- loo(fit_1)
print(loo_1)

Computed from 4000 by 434 log-likelihood matrix

         Estimate   SE
elpd_loo  -1914.8 13.8
p_loo         3.1  0.3
looic      3829.6 27.6
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Two predictors

fit_3 <- stan_glm(kid_score ~ mom_hs + mom_iq, data=kidiq,
                  seed=SEED, refresh=0)
loo_3 <- loo(fit_3)
print(loo_3)

Computed from 4000 by 434 log-likelihood matrix

         Estimate   SE
elpd_loo  -1876.1 14.2
p_loo         4.0  0.4
looic      3752.1 28.5
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

Compare models based on expected log predictive density

loo_compare(loo_1, loo_3)

      elpd_diff se_diff
fit_3   0.0       0.0  
fit_1 -38.7       8.3  

Bayes-R2

Median Bayes-R2 increases

br2_1<-bayes_R2(fit_1)
br2_3<-bayes_R2(fit_3)
round(median(br2_1),3)

[1] 0.056

round(median(br2_3),3)

[1] 0.215

Plot Bayes-R2 posteriors

Increase in R2 is clear

df <- melt(data.frame(fit_3=br2_1,fit_3n=br2_3))
ggplot(df, aes(x=value, linetype=variable)) +
    geom_density(alpha=0.25, show.legend=FALSE) +
    labs(x="Bayes-R^2", y="") +
    scale_y_continuous(breaks=NULL) + 
    annotate("text", x = 0.107, y = 16.2, label = "kid_score ~ mom_hs") +
    annotate("text", x = 0.28, y = 13.2, label = "kid_score ~ mom_hs + mom_iq")

R2 with LOO-CV

LOO-R2 are smaller, but the difference is still large.

round(median(loo_R2(fit_1)),3)

[1] 0.049

round(median(loo_R2(fit_3)),3)

[1] 0.205

Add five pure noise predictors to the data

set.seed(SEED)
n=nrow(kidiq)
kidiqr <- kidiq
kidiqr$noise <- array(rnorm(5*n), c(n,5))

Linear regression with additional noise predictors

fit_3n <- stan_glm(kid_score ~ mom_hs + mom_iq + noise, data=kidiqr,
                   seed=SEED, refresh=0)

Linear regression with interaction

fit_4 <- stan_glm(kid_score ~ mom_hs + mom_iq + mom_iq:mom_hs, data=kidiq,
                  seed=SEED, refresh=0)

Bayes-R2

Median Bayes-R2 increases, but...

br2_3<-bayes_R2(fit_3)
br2_3n<-bayes_R2(fit_3n)
round(median(br2_3),3)

[1] 0.215

round(median(br2_3n),3)

[1] 0.224

Plot Bayes-R2 posteriors

Median Bayes-R2 increases, but that increase is negligible compared to the uncertainty

df <- melt(data.frame(fit_3=br2_3,fit_3n=br2_3n))
ggplot(df, aes(x=value, linetype=variable)) +
    geom_density(alpha=0.25, show.legend=FALSE) +
    labs(x="Bayes-R^2", y="") +
    scale_y_continuous(breaks=NULL) + 
    annotate("text", x = 0.15, y = 11.5, label = "kid_score ~ mom_hs + mom_iq") +
    annotate("text", x = 0.285, y = 12, label = "kid_score ~ mom_hs + mom_iq +\n noise")

R2 with LOO-CV

LOO-R2 decreases when five noise predictors are addded

round(median(loo_R2(fit_3)),3)

[1] 0.205

round(median(loo_R2(fit_3n)),3)

[1] 0.191

LOO-R2 increases when interaction is addded

round(median(loo_R2(fit_4)),2)

[1] 0.22