OLS and Bayesian anti-essentialism

m7.1_ols <- lm(brain_std  ~ mass_std, data = d)
post <- extract.samples(m7.1_ols)
str(post) # no posterior for sigma

## 'data.frame':    10000 obs. of  2 variables:
##  $ Intercept: num  0.61 0.591 0.486 0.444 0.559 ...
##  $ mass_std : num  0.1151 0.155 0.2036 0.0876 0.021 ...

\(R^2\) - proportion of variance explained by model. Computation:

• Compute posterior predictive distributino for each observation (with sim)
• Subtract each observation from its prediction to get residual
• get variance of residuals and actual observations (actual empirical variance not sample variance that R calculates by default)

In principle, Bayesian approach requires we do this for each sample from the posterior. Traditionally R-squared is only computed at the mean prediction which is what we do here.

set.seed(12)
s <- sim(m7.1)                       # predictions
r <- apply(s, 2, mean) - d$brain_std # residuals

resid_var <- var2(r)
outcome_var <- var2(d$brain_std)

1 - resid_var / outcome_var

## [1] 0.4774589

Function to compute R2 again and again

R2_is_bad <- function(quap_fit){
    s <- sim(quap_fit, refresh = 0)
    r <- colMeans(s) - d$brain_std
    1 - var2(r)/var2(d$brain_std)
}

We will now make some models to compare to m7.1. 5 models will be considered each more complex and involving polynomial degrees.

\[ \begin{align} b_i &\sim Normal(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta_1 m_i + \beta_2 m_i^2 \\ \alpha &\sim Normal(0.5, 1) \\ \beta_j &\sim Normal(0, 10) \tag*{for j = 1..2} \\ \sigma &\sim \operatorname{Log-Normal}(0,1) \end{align} \] Same as before, only adds \(\beta_2\) parameter.

To define this model we could use earlier approach of using b1 and b2, we can also use a vector b and define its lengtrh using start list.

m7.2 <- quap(
    alist(
        brain_std ~ dnorm(mu, exp(log_sigma)), 
        mu <- a + b[1] * mass_std + b[2]*mass_std^2, 
        a ~ dnorm(0.5, 1), 
        b ~ dnorm(0, 10), 
        log_sigma ~ dnorm(0, 1)
    ), data = d, start = list(b = rep(0, 2))
)

We also create next four models in similar fashion with 3rd, 4th, 5th and 6th degree polynomials added.

m7.3 <- quap(
    alist(
        brain_std ~ dnorm(mu, exp(log_sigma)), 
        mu <- a + b[1] * mass_std + b[2] * mass_std^2 + b[3] * mass_std^3, 
        a ~ dnorm(0.5, 1), 
        b ~ dnorm(0, 10), 
        log_sigma ~ dnorm(0, 1)
    ), data = d, start = list(b = rep(0, 3))
)
m7.4 <- quap(
    alist(
        brain_std ~ dnorm(mu, exp(log_sigma)), 
        mu <- a + b[1] * mass_std + b[2] * mass_std^2 + b[3] * mass_std^3 + b[4] * mass_std^4, 
        a ~ dnorm(0.5, 1), 
        b ~ dnorm(0, 10), 
        log_sigma ~ dnorm(0, 1)
    ), data = d, start = list(b = rep(0, 4))
)
m7.5 <- quap(
    alist(
        brain_std ~ dnorm(mu, exp(log_sigma)), 
        mu <- a + b[1] * mass_std + b[2] * mass_std^2 + b[3] * mass_std^3 + b[4] * mass_std^4 + b[5] * mass_std^5,
        a ~ dnorm(0.5, 1), 
        b ~ dnorm(0, 10), 
        log_sigma ~ dnorm(0, 1)
    ), data = d, start = list(b = rep(0, 5))
)
# for m7.6 we replace sigma with constant 0.001. Model won't work otherwise. Why not? This will become evident when we plot them
m7.6 <- quap(
    alist(
        brain_std ~ dnorm(mu, 0.001), 
        mu <- a + b[1] * mass_std + b[2] * mass_std^2 + b[3] * mass_std^3 + b[4] * mass_std^4 + b[5] * mass_std^5 + b[6]*mass_std^6, 
        a ~ dnorm(0.5, 1), 
        b ~ dnorm(0, 10), 
        log_sigma ~ dnorm(0, 1)
    ), data = d, start = list(b = rep(0, 6))
)

Now we will plot each model following steps from previous chapters (extract samples from posterior, compute posterior predictive distribution at each of several locations on the horizontal axis, summarise and plot).

plot_model <- function(quap_fit){
    post <- extract.samples(quap_fit)
    mass_seq <- seq(from = min(d$mass_std), to = max(d$mass_std), length.out = 100)
    l <- link(quap_fit, data = list(mass_std = mass_seq))
    mu <- colMeans(l)
    ylim1 <- min(mu, d$brain_std)
    ylim2 <- max(mu, d$brain_std)
    ci <- apply(l, 2, PI)
    r2 <- round(R2_is_bad(quap_fit), 2)
    model_name <- deparse(substitute(quap_fit))
    plot(brain_std ~ mass_std, data = d, col = "deepskyblue4", pch = 20, cex = 1.8 , 
         xlab = "brain volume (cc)", ylab = "body mass (kg)", ylim = c(ylim1, ylim2))
    lines(mass_seq, mu)
    shade(ci, mass_seq)
    mtext(paste0(model_name,": R^2 = ", r2), adj = 0, cex = 0.8)
}

par(mfrow = c(2, 3))
models <- list(m7.1, m7.2, m7.3, m7.4, m7.5, m7.6)
# plot_model(m7.1)
# plot_model(m7.2)
# plot_model(m7.3)
# plot_model(m7.4)
# plot_model(m7.5)
# plot_model(m7.6)

# below function unscales the original axis too
brain_plot(m7.1)
brain_plot(m7.2)
brain_plot(m7.3)
brain_plot(m7.4)
brain_plot(m7.5)
brain_plot(m7.6)

\(R^2\) improves with the degree of the polynomial increasing in model. Sixth degree polynomial passes through each point and has no residual. This is why we had to set the sigma to 0.001, if it were estimated it would shrink to zero. (why sixth degree passes through each point? Because the equation in that model has 7 parameters, the same as the number of data points).

Underfitting

Models that are inaccurate both within and out of sample. Learn too little and fail to recover regular features of the sample.

Another way to conceptualise is to notice that it is insensitive to the sample. Removing any one point from the sample will and we will still get the same regression line.

Below we remove one point at a time and remodel m7.1 and m7.6.

par(mfrow = c(1, 2))
brain_loo_plot(m7.1)
brain_loo_plot(m7.6)

Firing the weather person

• cost-benefit anlaysis (cost of being wrong)
• accuracy in context - even ignoring cost-benefit, need a way to judge ‘accuracy’ that accounts for how much a model could possibly improve prediction.

Predictions of Current Weatherman-

New weatherman predicts sunny everyday.

Hit Rate - average chance of correct prediction.

Above hit rate for old guy is:

\[ \operatorname{Hit Rate} = 3 \times 1 + 7 \times 0.4 = 5.8 = \textrm{ 5.8 hits in 10 days} = 0.58 \textrm{ correct predictions per day} \]

Hit Rate = 3 X 1 + 7 X 0.4 = 5.8 hits in 10 days = 0.58 hits per day For the new guy, Hit Rate = 3 X 0 + 7 X 1 = 7 hits in 10 days = 0.70 hits per day (new guy wins?!)

p <- c(0.3, 0.7)
-sum(p * log(p))

## [1] 0.6108643

p <- c(0.7, 0.15, 0.15)
-sum(p * log(p))

## [1] 0.8188085

p <- c(0.3, 0.7)
q1 <- seq(0.01, 0.99, by = 0.01)
q2 <- 1 - q1

dkl <- apply(cbind(q1, q2), 1, function(q) sum(p * log(p / q))) 

par(mfrow = c(1,1))

plot(q1, dkl, type = 'l', xlab = "q[1]", ylab = "Divergence fo q from p", lwd = 2, col = 'blue')
abline(v = 0.3, lty = 2, lwd = 1, col = 'black', )
abline(h = 0, lty = 2, lwd = 1, col = 'black')
text(x = 0.3, y = 1.5, "q = p", adj = -0.1)

Estimating Divergence

We typically don’t konw true p. So how do we even compute \(D_{KL}\)? There’s a way out. If we have two competing models we only want to find out which one is closer to target. Let’s say two models q and r. computing divergence for both has log(p) in it which subtracts out if we compute distance of q from r.

So we only need a model’s average log-prob : \(E \log(q_i)\)

Thus

\[ S(q) = \sum_i \log (q_i) \]

For Bayesian model, we have to use entire posterior distribution (otherwise we are throwing away info).

So for each observation, we find the log of the avg prob over entire posterior dist. Log Pointwise Predictive Density.

set.seed(1)
lppd(m7.1, n = 1e4)

## [1]  0.6098665  0.6483436  0.5496095  0.6234932  0.4648144  0.4347607 -0.8444635

We get 7 values, each a log-prob score for a specific obs. Summing these can provide the total log-prob score for model and data.

• Larger values are better (indicates larger average accuracy)
• Deviance is like a lppd score, but multiplied by -2, so that smaller values are better.

set.seed(1)
logprob <- sim(m7.1, ll = TRUE, n = 1e4)
n <- ncol(logprob)
ns <- nrow(logprob)
f <- function(i) log_sum_exp( logprob[,i] ) - log(ns)
(lppd <- sapply( 1:n, f))

## [1]  0.6098665  0.6483436  0.5496095  0.6234932  0.4648144  0.4347607 -0.8444635

set.seed(1)
sapply(models, function(m) sum(lppd(m)))

## [1]  2.490390  2.566165  3.707343  5.333750 14.090061 39.448518

# N <- 20 
# kseq <- 1:5 
# dev <- sapply( kseq , function(k) { 
#     print(k); 
#     r <- replicate( 5 , sim_train_test( N=N, k=k, cv.cores = 4 ) );
#     c( mean(r[1,]) , mean(r[2,]) , sd(r[1,]) , sd(r[2,]) )
# } )
# 
# plot( 1:5 , dev[1,] , 
#       ylim = c( min(dev[1:2,])-5 , max(dev[1:2,])+10 ) , 
#       xlim=c(1,5.1) , 
#       xlab="number of parameters" , 
#       ylab="deviance" , 
#       pch=16 , col=rangi2 )
# 
# mtext( concat( "N = ",N ) ) 
# points( (1:5)+0.1 , dev[2,] ) 
# for ( i in kseq ) { 
#     pts_in <- dev[1,i] + c(-1,+1)*dev[3,i] 
#     pts_out <- dev[2,i] + c(-1,+1)*dev[4,i] 
#     lines( c(i,i) , pts_in , col=rangi2 ) 
#     lines( c(i,i)+0.1 , pts_out )
# }

Model mis-selection

Good PSIS or WAIC only indication of predictive accuracy not causality.

Re-running fungus example from previous chapter.

suppressPackageStartupMessages(library(rethinking))
set.seed(71)
# number of plants
N <- 100

# simulate initial heights
h0 <- rnorm(N, 10, 2)
# assign treatments and simulate fungus and growth
treatment <- rep(0:1, each = N / 2)
fungus <- rbinom(N, size = 1, prob = 0.5 - treatment * 0.4)
h1 <- h0 + rnorm(N, 5 - 3 * fungus)

# compose a clean data frame
d <- data.frame(h0 = h0, h1 = h1, treatment = treatment, fungus = fungus)

m6.6 <- quap(
    alist(
        h1    ~ dnorm(mu, sigma), 
        mu    <- h0 * p, 
        p     ~ dlnorm(0, 0.25), 
        sigma ~ dexp(1)
    ), data = d
)

m6.7 <- quap(
    alist(
        h1    ~ dnorm(mu, sigma), 
        mu    <- h0 * p, 
        p     <- a + bt * treatment + bf * fungus, 
        a     ~ dlnorm(0, 0.25), 
        bt    ~ dnorm(0, 0.5), 
        bf    ~ dnorm(0, 0.5), 
        sigma ~ dexp(1)
    ), data = d
)

m6.8 <- quap(
    alist(
        h1    ~ dnorm(mu, sigma), 
        mu    <- h0 * p, 
        p     <- a + bt * treatment, 
        a     ~ dlnorm(0, 0.25), 
        bt    ~ dnorm(0, 0.5), 
        sigma ~ dexp(1)
    ), data = d
)

m6.6 is the model with just an intercept.

m6.7 is the model with both treatment and fungus (post-treatment var)

m6.8 is the model with treatment but not fungus.

set.seed(11)
WAIC(m6.7)

##       WAIC      lppd  penalty  std_err
## 1 361.4371 -177.1597 3.558906 14.17059

WAIC is the guess for the out-of-sample deviance.

rethinking provides convenience function compare to … compare models!

set.seed(77)
compare(m6.6, m6.7, m6.8, func = WAIC)

##          WAIC       SE    dWAIC      dSE    pWAIC       weight
## m6.7 361.8803 14.26260  0.00000       NA 3.845703 1.000000e+00
## m6.8 402.7684 11.29251 40.88804 10.53858 2.649165 1.322129e-09
## m6.6 405.9174 11.66153 44.03705 12.24953 1.582861 2.738272e-10

Cols are - WAIC, std error of WAIC, differnce of each WAIC from the best model, standard error of this difference, prediction penalty (pWAIC), and the Akaike weight.

• WAIC - smaller values are better. Models automatically sorted from best to worst. Model with fungus variable has smallest WAIC.

• pWAIC - close to (but slightly below) the number of dimensions in the posterior of each model (to be expected in lin reg with reg priors)

• dWAIC - difference in WAIC from best model. e.g. m6.7 is 40 units deviance away from m6.8. and intercept is 3 units deviance smaller m6.8. We can’t say directly whether these are big or small differences. Better question, are the models easily distinguished by their expected out of sample accuracy? We need to consider error in the WAIC estimates to answer this. We don’t have a target sample, so these are just guesses.

• SE - approx std error of each WAIC. approximately we expect uncertainty in out of sample accuracy to be normally distributed with mean equal to the reported WAIC value and a std deviation equal to the std error. (this approximation tends to drastically underestimate uncertainty on smaller samples).

Judging whether two models are easy to distinguish.

We use std error of their difference not SE of WAIC. To compute the std error fo the differnece between models m6.7 and m6.8, we need the pointwise breakdown of the waic values:

set.seed(91)
waic_m6.7 <- WAIC(m6.7, pointwise = TRUE)$WAIC
waic_m6.8 <- WAIC(m6.8, pointwise = TRUE)$WAIC
n <- length(waic_m6.7)
diff_m6.7_m6.8 <- waic_m6.7 - waic_m6.8

sqrt(n * var(diff_m6.7_m6.8))

## [1] 10.41851

This is the value in the dSE column. (difference due to simulation variance).

Difference between models is 40.9 and std error is 10.5. The 99% interval of difference will be (z-score of 2.6):

40.0 + c(-1,1) * 10.5 * 2.6

## [1] 12.7 67.3

doesn’t contain zero, so models are verry easy to distinguish by expected out of sample accuracy. m6.7 is a lot better. A plot helps visualise this easily.

plot(compare(m6.6, m6.7, m6.8))

• Filled points are the in-sample deviance values.
• open points are the WAIC values. (each model does better in-sample)
• line segments show the std error of each WAIC (values from col SE)
• lighter line segment shows std error of difference

So WAIC not good for making causal inferences (which it is not supposed to do anyway).

difference between m6.8 and m6.6 (with treatment vs only intercept model)

m6.8 provides good evidence that treatment works

precis(m6.8)

##             mean         sd       5.5%     94.5%
## a     1.38168639 0.02519766 1.34141567 1.4219571
## bt    0.08366284 0.03431330 0.02882356 0.1385021
## sigma 1.74629890 0.12190853 1.55146553 1.9411323

but difference in deviance is only 3 units of deviance.

set.seed(92)
waic_m6.6 <- WAIC( m6.6 , pointwise=TRUE )$WAIC 
diff_m6.6_m6.8 <- waic_m6.6 - waic_m6.8 
sqrt( n*var( diff_m6.6_m6.8 ) )

## [1] 4.78333

compare table doesn’t show this value by default but it did calculate it and we can access using dSE slot.

set.seed(93)
compare(m6.6, m6.7, m6.8)@dSE

##           m6.6     m6.7      m6.8
## m6.6        NA 12.22586  4.860145
## m6.7 12.225864       NA 10.487119
## m6.8  4.860145 10.48712        NA

we get matrix with pariwise dSE for the models compared.

dSE for m6.6 vs m6.8 is larger than the difference itself. So we can’t distinguish these models on the basis of WAIC. This is since even though treatment does work, the impact on final plant height is minor. So it alone does’t do much to improve prediction of plant height.

• weight column - these are a traditional way to summarise relative support for each model. They always sum to 1, within a set of compared models. The weight of a model i is computed as:

\[ w_i = \frac{\exp(-0.5\Delta_i)}{\sum_j \exp(-0.5\Delta_j)} \]

\(\Delta_i\) is the difference between model i’s WAIC value and the best WAIC in the set. (dWAIC) These can be a quick way to see how big the differences are among models. But we still need to inspect the std errors since weights don’t reflect them.

Weights are also used in Model Averaging (combining the predictions of multiple models).

Outliers and other illusions

In the Divorce example, some states (observations in that dataset) were outliers (like Idaho). Let’s see how PSIS and WAIC represen that importantce. redoing divorce example.

library(rethinking) 
data(WaffleDivorce) 
d <- WaffleDivorce 
d$A <- standardize( d$MedianAgeMarriage ) 
d$D <- standardize( d$Divorce ) 
d$M <- standardize( d$Marriage )
m5.1 <- quap( 
    alist( 
        D ~ dnorm( mu , sigma ) , 
        mu <- a + bA * A , 
        a ~ dnorm( 0 , 0.2 ) , 
        bA ~ dnorm( 0 , 0.5 ) , 
        sigma ~ dexp( 1 )
    ) , data = d )
m5.2 <- quap( 
    alist( 
        D ~ dnorm( mu , sigma ) , 
        mu <- a + bM * M , 
        a ~ dnorm( 0 , 0.2 ) ,
        bM ~ dnorm( 0 , 0.5 ) , 
        sigma ~ dexp( 1 )
    ) , data = d )
m5.3 <- quap( 
    alist( 
        D ~ dnorm( mu , sigma ) , 
        mu <- a + bM*M + bA*A ,
        a ~ dnorm( 0 , 0.2 ) ,
        bM ~ dnorm( 0 , 0.5 ) ,
        bA ~ dnorm( 0 , 0.5 ) ,
        sigma ~ dexp( 1 )
    ) , data = d )

M has no influence on D after controlling for A.

set.seed(24071847)
compare(m5.1, m5.2, m5.3, func = PSIS)

##          PSIS       SE     dPSIS       dSE    pPSIS       weight
## m5.1 127.5665 14.69485  0.000000        NA 4.671425 0.8340176286
## m5.3 130.8062 16.15696  3.239685  1.808671 6.578663 0.1650769863
## m5.2 141.2178 11.56557 13.651299 10.923820 4.057204 0.0009053851

Model without M is on top. since M has little relation with D. So model w/o it has slightly better expected out-of-sample performance even though it actually fits teh sample slightly worse than m5.3, the model wiht both predictors.

We also get message that some Pareto k values are very high. This implies that the smoothing approximation that PSIS uses is unreliable for some points. (when a point’s Pareto k value > 0.5, importance weight can be unreliable). These points tend to be outliers with unlikely values according to the model. As a result they are highly influential and make it difficult to estimate out of sample predictive accuracy. Because any new sample is unlikely to contain these same outliers, and since these were highly influential, they could make out-of-sample predictions worse than expected.

Looking at pointwise PSIS

set.seed(24071847)
PSIS_m5.3 <- PSIS(m5.3, pointwise = TRUE)
set.seed(24071847)
WAIC_m5.3 <- WAIC(m5.3, pointwise = TRUE)

plot(PSIS_m5.3$k, WAIC_m5.3$penalty, xlab = "PSIS Pareto k", ylab = "WAIC penalty", col = rangi2, lwd = 2)

#identify(PSIS_m5.3$k, WAIC_m5.3$penalty, label = d$Loc)

State of Idaho has Pareto k and high WAIC penalty.

WAIC penalty is sometimes called the “effective number of parameters” since in ordinary linear regression, sum of all penalty terms from all points tends to be equal to the number of free parameters in the model. However, in this case there are 4 parameters in the model while the penalty is closer to 6. Outlier Idaho is causing this additional overfitting risk.

WAIC(m5.3)

##      WAIC      lppd  penalty  std_err
## 1 128.483 -58.79191 5.449613 14.35521

m5.3t <- quap(
    alist(
        D ~ dstudent(2, mu, sigma), 
        mu <- a + bM * M + bA*A, 
        a ~ dnorm(0, 0.2), 
        c(bM, bA) ~ dnorm(0, 0.5), 
        sigma ~ dexp(1)
    ), data = d
)
PSIS(m5.3t) # no message that some Pareto k values are very high this time

##       PSIS      lppd  penalty  std_err
## 1 133.5325 -66.76626 6.745632 11.87595

plot(coeftab(m5.3, m5.3t))