Estimating errors in predictive accuracy using K-fold cross validation

[James Camac]

Hi all,

I've built a stan model used k-fold cross validation as done in Aki's paper. The only difference from the paper is that I've monitored the sum of the log likelihoods as opposed to the individual log likelihoods. I did this because of memory issues as I'm dealing with a large dataset. 

In my case, the point of doing cross validation was to compare predictive accuracy between a suite of models and ultimately choose a single 'best' model. I was hoping to do this using the mean and credible intervals of the monitored log likelihoods pooled across kfolds.

Specifically I'd:

Pool the samples of the sum log likelihood across all K-fold fits and then estimate the mean and the lower and upper 95% quantiles. My thinking was that this would capture the full uncertainty associated with estimating the log likelihood of the heldout data and allow people to determine whether there was a significant difference between model predictive accuracy. 

However, I'm beginning to think my approach isn't technically correct. Reading Aki's paper again I see that the standard error is estimated as: SQRT(N) x SD(log likelihoods).

Another approach I've seen people use is a point estimate approach whereby they extract the mean log likelihood for each kfold and then use these values to estimate the average log likelihood and confidence interval (i.e. across kfolds).

Does anyone have any guidance on how I should best approach this?

[James Camac]

I guess I should have also stated that I could use a log likelihood ratio between two models (i.e. LL of complex model/ LL of Simpler model) and compare the it to the chi-square distribution... But again I feel like this approach doesn't take into consideration the uncertainty in the estimates as it relies on a point estimate.

[Aki Vehtari]

Hi,

If you need save memory, in K-fold-CV you can do like this

generated quantities {

  real log_likt;

  log_likt <- 0;

  for (n in 1:Nt)

    log_likt <- log_likt + bernoulli_logit_log(yt[n], beta0 + xt[n]*beta);

}

Then outside of Stan compute log of the expected summed densities (summing was made in Stan)

 elpdhats(k)=log(mean(exp(log_likt)) 

Then you will have K elpdhats values, instead of N values as in the paper you mention.You can then compute

standard error and pairwise model comparison using these values and K, instead of elpdhat_i and LOO as in the paper.

This approach is mentioned, e.g., in http://becs.aalto.fi/en/research/bayes/publications/VehtariLampinen_NC2002_preprint.pdf

Aki

[Bob Carpenter]

I'm surprised this was Aki's recommendation since he's the one
who suggested I calculate it more stably using log_sum_exp:

elpdhats[k] <- log(mean(exp(log_likt)));

should be:

elpdhats[k] <- log_sum_exp(log_likt) - log(N);

where N is the number of elements in log_likt and log_sum_exp is
the usual robust implementation as

log_sum_exp(u) = max(u) + log(sum(exp(u - max(u)))

which won't underflow.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[James Camac]

Thanks Aki and Bob,

I have one more question. Why the need for log(mean(exp(log_likt)) ?

Isn't just calculating the mean and standard error log_likt across kfolds sufficient?

As in wouldn't that give you the mean expected log likelihood?

Isn't a lpd the same as a log likelihood?

[Bob Carpenter]

I just went through this while writing a tutorial :-)

The problem is that log doesn't distribute through integrals
or sums; it's easy to see in the sum case for a vector x of
positive values:

mean(log(x)) != log(mean(x))

- Bob

[James Camac]

Ah obvious! Silly me.

Thanks

[Aki Vehtari]

On Monday, February 1, 2016 at 8:33:20 PM UTC+2, Bob Carpenter wrote:

I'm surprised this was Aki's recommendation since he's the one
who suggested I calculate it more stably using log_sum_exp:

  elpdhats[k] <- log(mean(exp(log_likt)));

should be:

  elpdhats[k] <- log_sum_exp(log_likt) - log(N);

Heh, the latter one is exactly what I use :) 

I was going to add it in my answer, too, but I was in a hurry and forgot it...

Aki 

[James Camac]

Thanks again guys. This forum has really taught me a lot about better coding and general queries about probability theory etc.

I really appreciate the time you invest in answering peoples questions.