# Trying to find WAIC for a simple model from Statistical Rethinking

**URL:** <https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991>\
**Category:** Modeling\
**Tags:** loo\
**Created:** [April 25, 2018, 7:41pm UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991 "2018-04-25T19:41:58Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![rosgori](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rosgori/32/3059_2.png) [@rosgori](https://discourse.mc-stan.org/u/rosgori)\
**Post date:** [April 25, 2018, 7:41pm UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991/1 "2018-04-25T19:41:59Z")

</div>

I have been reading _Statistical Rethinking_ and I want to replicate code 10.23 and code 10.24 in Stan. For the moment, I only care about model `m10.6`. So I have:

```stan
data{
    int<lower=1> N;
    int admit[N];
    int applications[N];
    real male[N];
}
parameters{
    real a;
    real bm;
}
model{
    vector[N] p;
    bm ~ normal( 0 , 10 );
    a ~ normal( 0 , 10 );
    for ( i in 1:N ) {
        p[i] = a + bm * male[i];
    }
    admit ~ binomial_logit( applications , p );
}
generated quantities{
    vector[N] p;
    real dev;
    dev = 0;
    for ( i in 1:N ) {
        p[i] = a + bm * male[i];
    }
    dev = dev + (-2)*binomial_logit_lpmf( admit | applications , p );
}

```

I have extracted this code from `rethinking::stancode(m10.6)`. The command I am using is

```stan
fit_admit <- stan(file = 'admit.stan', data = d_copy, iter = 5000, chains = 4)

```

with `d_copy` as

```stan
d_copy <- list( N = 12,
                admit = c(512, 89, 353, 17, 120, 202, 138, 131, 53, 94, 22, 24),
                applications = c(825, 108, 560, 25, 325, 593, 417, 375, 191, 393, 373, 341),
                male = c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0))

```

Since I want WAIC, I use `library(loo)` and I do the following

```stan
library(loo)
loo::waic(fit_admit)

```

and I get this error

```stan
Error in UseMethod("waic") : 
  no applicable method for 'waic' applied to an object of class "stanfit"

```

then I realize that I need `loo::extract_log_lik()`, according to [this](https://cran.r-project.org/web/packages/loo/loo.pdf), p.6. And this is my problem, it should be easy, but it is not. The document gives an example:

> For example, the following is the generated quantities block for computing and saving the  
> log-likelihood for a linear regression model with N data points, outcome y, predictor matrix X,  
> coefficients beta, and standard deviation sigma:  
> `vector[N] log_lik; for (n in 1:N) log_lik[n] = normal_lpdf(y[n] | X[n,] * beta, sigma);`

but I cannot translate that to my model `m10.6` (because I don’t understand it). The examples that are shown [here](http://www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf) should help me, but…

Anyway, my question is simple: What is the expression for `log_lik` in this case?

By the way, if you really want to know, the WAIC score for `m10.6` shown in the book is presumably incorrect, it doesn’t matter if you use `map()` o `map2stan()`.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [April 25, 2018, 9:02pm UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991/2 "2018-04-25T21:02:47Z")

</div>

The log-likelihood is the same as the deviance without the -2 factor, so you could change the Stan code to output log-likelihoods rather than the deviance. Or you can create the matrix of log-likelihoods in R after the fact with

```
post <- as.data.frame(fit_admit)
ll <- sapply(1:d_copy$N, FUN = function(i) {
  dbinom(d_copy$admit[i], size = d_copy$applications[i], 
         prob = plogis(post$a + post$bm * d_copy$male[i]), log = TRUE)
})

```

At which point you can call `loo::waic(ll)` and it will tell you

```
11 (91.7%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

```

The loo package recommends against `waic` because `loo` is more robust and outputs more diagnostics. If you call `loo::loo(ll))`, you get

```
Pareto k diagnostic values:
                         Count Pct. Min. n_eff
(-Inf, 0.5] (good) 4 33.3% 1396      
 (0.5, 0.7] (ok) 2 16.7% 407       
   (0.7, 1] (bad) 2 16.7% 16        
   (1, Inf) (very bad) 4 33.3% 1         
See help('pareto-k-diagnostic') for details.

```

which basically means the assumptions for calculating these estimates are not satisfied and you should not use them. With only 12 observations, that is not too surprising but it would be fairly easy to estimate your model 12 times each time leaving out one observation in order to do cross-validation.

---

<div class="post-metadata">

**Author:** ![rosgori](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rosgori/32/3059_2.png) [@rosgori](https://discourse.mc-stan.org/u/rosgori)\
**Post date:** [April 25, 2018, 11:18pm UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991/3 "2018-04-25T23:18:29Z")

</div>

Thanks! And that function is easy to understand!

---

<div class="post-metadata">

**Author:** ![richard\_mcelreath](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/richard_mcelreath/32/107_2.png) [@richard\_mcelreath](https://discourse.mc-stan.org/u/richard_mcelreath)\
**Post date:** [May 2, 2018, 7:41am UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991/4 "2018-05-02T07:41:21Z")

</div>

David pointed me here. The issue of how to do WAIC/LOO with binomial models is a subtle one. My rethinking code actually splits these 12 observations into Bernoulli 0/1 outcomes and then does the WAIC/LOO calculation. So the results will differ. If you think the right leave-one-out design is by department/gender combination, then the code in this thread look correct. If instead you think each application should be left out, then my code is (at least in that regard) correct.

This would at least explain why it appears to be only the aggregated binomial models in Chapter 10 that yield conflicting results. Assuming that is the case.

I think there may still be subtle numerical issues with how I do my WAIC calculation, and such an issue would more likely show up in small data sets. So there is perhaps something more to figure out.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [May 2, 2018, 4:59pm UTC](https://discourse.mc-stan.org/t/trying-to-find-waic-for-a-simple-model-from-statistical-rethinking/3991/5 "2018-05-02T16:59:52Z")

</div>

Like Richard said, there’s not a right way to do this as it depends on whether you care about leaving out observations or groups of observations. Then there’s also sometimes the complication that groups are unbalanced in your data set but you don’t think they’re so unbalanced in the population. Then even if you leave out one group at a time it can still make sense to apply an additional correction for having different numbers of observations in groups. But it’s totally context dependent whether that makes sense to do.

I thought we had put an example of doing something like this somewhere in the loo package doc or vignettes but I can’t find it now. Definitely an example worth having somewhere though!