tl;dr: I get warnings about the pareto k diagnostics when the log-likelihood is constant for some data points. Should I bother about these warnings then?
I am fitting a model that describes how participants select iteratively between different slot machines. The model assumes participants start with a random choice between all k slot machines (let’s say, k = 2) and then update their reward expectations, leading to preferences for slot machines with higher payoffs in the following trials. Crucially, the log-likelihood in the very first trial is independent of the model parameters (it is simply log(1/k), assuming all slot machines are picked with the same probability in the first trial) and therefore it has the identical value across all MCMC-iterations.
When using loo(), I get the warning about too high Pareto k diagnostic values and with pareto_k_values() I can trace this down to k=Inf for the very first trial were the likelihood is constant (and from my rather limited understanding, the Pareto distribution probably cannot be estimated). Can I use elpd_loo then safely despite these warnings (since I know that they are not due to bad model fit)? Alternatively, can I simply omit the log-likelihood for this first data point and only include the loglikelihood of the remaining data points when using loo?
Some example code:
# First, simulate some data
winchance <- c(0.7, 0.3) # chances to win for two slot machines
ntrials <- 100 # number of trials
Q <- c(0,0) # initial reward expectations
alpha <- 0.2 # learning rate
theta <- 5 # inverse softmax temperature
choices <- c() # save all choices here
rewards <- c() # save all rewards here
for (t in 1:ntrials) { # iterate over trials
prob <- exp(theta*Q)/sum(exp(theta*Q)) # use softmax to compute choice probabilities
choice <- sample(c(1,2), 1, prob=prob) # sample choice
reward <- rbinom(1, 1, winchance[choice]) # sample reward
Q[choice] <- Q[choice] + alpha * (reward - Q[choice]) # update reward expectation based on temporal-difference learning
choices <- c(choices, choice) # append choice
rewards <- c(rewards, reward) # append reward
}
code <- "
data {
int ntrials;
int choices[ntrials];
int rewards[ntrials];
}
parameters {
real<lower=0, upper=1> alpha; // learning rate
real<lower=0> theta; // inverse softmax temperature
}
transformed parameters {
matrix[2,ntrials] prob; // choice probabilities
vector[2] Q = [0, 0]'; // reward expectations
for (t in 1:ntrials) { // iterate over trials
prob[1,t] = exp(theta*Q[1])/sum(exp(theta*Q)); // probability choosing slot machine 1
prob[2,t] = exp(theta*Q[2])/sum(exp(theta*Q)); // probability choosing slot machine 2
Q[choices[t]] = Q[choices[t]] + alpha * (rewards[t]-Q[choices[t]]); // update reward expectation
}
}
model {
alpha ~ beta(1,1);
theta ~ gamma(7,2);
for (t in 1:ntrials) {
choices[t] ~ categorical(prob[,t]);
}
}
generated quantities {
vector[ntrials] log_lik;
for (t in 1:ntrials) {
log_lik[t] = categorical_lpmf(choices[t] | prob[,t]);
}
}
"
library(rstan)
fit <- stan(model_code = code, data=list(ntrials=ntrials, choices=choices, rewards=rewards), pars=c("alpha", "theta", "log_lik"))
library(loo)
log_lik <- extract_log_lik(fit, merge_chains = FALSE)
r_eff <- relative_eff(exp(log_lik), cores = 2)
loo <- loo(log_lik, r_eff = r_eff, cores = 2)
> loo
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -24.7 5.6
p_loo 1.6 0.5
looic 49.3 11.2
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 99 99.0% 2122
(0.5, 0.7] (ok) 0 0.0% <NA>
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 1 1.0% 2000
See help('pareto-k-diagnostic') for details.
> pareto_k_values(loo)
[1] Inf -0.250820915 -0.231890100 -0.194362530 -0.123999775 -0.082366031 -0.019329143 -0.006330616 0.037577266
[10] 0.091040580 0.095732899 0.094355175 0.067537553 0.094631436 0.001246631 0.104696377 0.064945659 -0.021309253
[19] 0.115981847 0.078319484 -0.058639911 0.009135191 0.061737323 0.074738349 0.094201735 0.116374869 0.009274860
[28] 0.045051165 -0.067315539 0.079272260 0.027958081 0.076988415 0.110051724 -0.012523528 0.094422197 -0.081883085
[37] 0.069715824 -0.059791978 -0.131969375 -0.024692136 0.092425508 0.077963691 0.078872533 0.107749585 0.130824833
[46] 0.150694845 0.043327046 -0.037994837 -0.101710611 -0.066409222 0.119235100 0.063381336 0.075570001 -0.026975655
[55] -0.055317108 -0.143377748 -0.027658296 0.094034307 0.052289974 -0.044402115 0.079836689 -0.091743372 -0.116907409
[64] -0.167224110 0.092147119 0.307332044 -0.127042861 -0.113747312 -0.046923850 0.031718152 0.038639365 0.054683168
[73] -0.103481193 -0.046339586 0.029732864 0.059025017 0.051632538 0.058639093 0.065313094 -0.036766287 -0.044776130
[82] -0.065643602 0.188474906 -0.081540664 -0.072533059 0.155653163 0.385161563 0.041677759 0.056500547 -0.096923764
[91] -0.095223578 -0.102916686 -0.092051979 0.141964451 0.150588115 0.111743144 0.087024258 0.109849397 0.123994752
[100] 0.020656058