WAIC for multiple subjects

[glae...@gmail.com]

Dear Bob and other Stan Users,

I am looking at the WAIC example code for the simple regression example in BDA, p. 177

https://github.com/stan-dev/stan/issues/473

and I am working on applying it to my models, which are variants of reinforcement learning using a softmax function on the expected values to compute an action probability (= model prediction of data (subject choices)). This in turn is usually used to compute the log likelihood of the model.

From the example code I gather that in order to compute WAIC I need to get a log likelihood for every data point. In my case, I have several subjects doing an experiment with about 64 - 200 trials depending on the task. Following the example, I think I should get log likelihood distribution for every data point in every subject, which get into the range of 5000 - 120000 posteriors of log likelihoods. Would this be the recommendation of the experts or should I rather calculate a likelihood per subject by summing/averaging over all trials in a participant?

Furthermore, the example code also uses colMeans to average over MCMC samples. Is this also a valid approach, if I have very skewed posteriors, or is this only valid for almost symmetric unimodal posteriors? (The histograms for log_lik in the simple regression example definitely look like symmetric and unimodal.)

Thanks a lot for your insights.
Jan


--
Jan Glaescher
glae...@gmail.com

[Avraham Adler]

Hello, Jan.

Regarding the code, below is a slightly more robust version of the R code that Drs. Gelman and Vehtari used in their paper on WAIC, perhaps you may find it helpful. Just note that the Gelman & Vehtari code does not put loo-cv on the deviance scale. If you want to, just multiply it by -2.

WAIC <- function(stanfit){
  colVars <- function(M){
    vars <- M[1, ]
    for (n in seq_along(vars)) vars[n] <- var(M[, n])
    return(vars)
  }
  log_lik <- extract (stanfit, "log_lik")$log_lik
  dim(log_lik) <- if (length(dim(log_lik)) == 1) c(length(log_lik), 1) else
    c(dim(log_lik)[1], prod(dim(log_lik)[2:length(dim(log_lik))]))
  S <- nrow(log_lik)
  n <- ncol(log_lik)
  lpd <- log(colMeans(exp(log_lik)))
  p_waic <- colVars(log_lik)
  elpd_waic <- lpd - p_waic
  waic <- -2 * elpd_waic
  loo_weights_raw <- 1 / exp(log_lik - max(log_lik))
  loo_weights_normalized <- loo_weights_raw/
    matrix(colMeans(loo_weights_raw), nrow=S, ncol=n, byrow=TRUE)
  loo_weights_regularized <- pmin(loo_weights_normalized, sqrt(S))
  elpd_loo <- log(colMeans(exp(log_lik) * loo_weights_regularized) /
                  colMeans(loo_weights_regularized))
  p_loo <- lpd - elpd_loo
  pointwise <- cbind(waic, lpd, p_waic, elpd_waic, p_loo, elpd_loo)
  total <- colSums(pointwise)
  se <- sqrt(n * colVars(pointwise))
  return(list(waic=total["waic"], elpd_waic=total["elpd_waic"],
              p_waic=total["p_waic"], elpd_loo=total["elpd_loo"], p_loo=total["p_loo"],
              pointwise=pointwise, total=total, se=se))
}

As for the use of colMeans, all of the xAIC measures are based on an expectation--that of the log predictive density for the existing data. Much of the difference is in how these are adjusted for the inherent bias of testing against the training data. The benefit of WAIC vs. DIC {and yes, I'm probably breaking rigor here} is that the calculations are not done conditioned against the "plug in" Bayesian posterior of the parameters, but at the mean of the distributions of each of the points themselves, but it is still a penalized expectation measure, thus the need for for colMeans/colVars. If the distribution of the posterior of each data point is highly skewed, some of that should be captured in colVars. If you haven't seen it already, you'll probably find this helpful:<www.stat.columbia.edu/~gelman/research/unpublished/waic_stan.pdf>

Hope that helps!

Avi

[Aki Vehtari]

You should first decide whether you are interested in the predictive performance for one new trial or for one new subject. For example, in brain research it is usually more interesting to estimate the model predictive performance for new subject, so that it is more justified to say that the model generalizes to other people than just the observed people.

WAIC is an approximation of leave-one-out cross-validation and works if removing one data point changes the posterior only a little. In this case I would assume that removing a single trial observation would change the posterior only a little and WAIC would work well.

If you are interested in leave-one-subject-out you could consider the likelihood per subject, but since it is the total likelihood contribution from that subject, you need to compute as much as in leave-one-trial-out. WAIC approximation for leave-one-subject-out can work if removing one subject changes the posterior only a little, but since you are removing more information it is more likely to have a bigger change, too. In that case, you would need to make the leave-one-subject-out cross-validation by using k-fold-cross-validation where in each fold you leave on subject out.

I don't have a figure showing how WAIC would probably behave in leave-one-subject-out case, but WAIC has a certain similarity to importance sampling cross-validation, and you may look Figure 11 in 

 Aki Vehtari and Jouko Lampinen (2002). Bayesian model assessment and comparison using cross-validation predictive densities. Neural Computation, 14(10):2439-2468.

 http://becs.aalto.fi/en/research/bayes/publications/VehtariLampinen_NC2002_preprint.pdf

In that specific case, leave-one-data point-out worked, but leave-one-group of data-out did not.

Aki

[glae...@gmail.com]

Dear Avi,

thanks a lot for your answer and your code, both of which I find really helpful. In the email digest I have seen that you actually posted a much shorter R code, but in the web interface it is marked as deleted. Should I ignore this shorter code version then?

Thanks a lot,

Jan

[glae...@gmail.com]

Dear Aki,

thanks a lot for your insightful answer. I think you are right in that I am more interested in the leave-one-subject-out case, but you haveconvinced me that I need to check that carefully, although I guess the impact of leave-one-subject-out will also be influenced by the sample size, which in one study is quite big  (~600).

Thanks a lot,

Jan

[Avraham Adler]

Hi, Jan.

What I forgot was that shorter code could not be applied to a stanfit item, but to the extracted matrix of log probabilities. As I'm not sure how your data is set up, it isn't directly applicable. However, realizing that the pointmatrix (or whatever the input was named) is a matrix where the columns are the data points and the rows are the observations of the densities at those point at each iteration of the simulation, you should be able to adjust the shorter code as needed. Otherwise, the longer code should be able to be applied to any stanfit item.

One other possible difference of which you should be aware, is that for the stan models I built, I capture the densities themselves. I beleive Gelman and Vetari capture log-densities. Looking at the code should identify which was used.

Thanks,

Avi

[glae...@gmail.com]

Dear Avi,

thanks for the clarification. I am collecting the matrix of simulations posterior densities as you describe so the short code version should be fine. But in the future, I’ll probably use the other version that is suitable for stanfit objects. Thanks also for mentioning the difference between different versions in terms of densities and log densities. I’ll watch out for that.

Best,

Jan

--

Jan Glaescher

glae...@gmail.com

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/G9VuPvQchFY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.