Hi, Thank you a lot for your reply! Sorry I didn’t upload the code as it seems a bit long.

Here is the script calculating the log_liks. It is written in JAGS.

```
data {
for (s in 1:Nsamps) {
for (i in 1:N) { # trial level
loglik[s, i] <- logdensity.wieners(y[i], 1 , ndt[s,i], bias[s,i], w[s,i], alpha.p[s,idxP[i]])
# generate trial-by-trial nDT
ndt[s,i] <- ndt.p[s,idxP[i]]
dT[s,i] <- rt[i] - ndt[s,i]
# generate trial-by-trial Bias
bias[s,i] <- bias.p[s,idxP[i]]
# The actual drift rate
w[s,i] <- li.hat[s,i]
# rt minus timeHin
timeIn[s,i] <- abs(time.p[s,idxP[i]])
timeStep[s,i] <- (dT[s,i] - timeIn[s,i])/dT[s,i] #time is a weighted variable depending when the health(b1) or taste(b2) comes first
timeInStep[s,i] <- timeIn[s,i]/dT[s,i]
smaller[s,i] <- step(timeStep[s,i]) #if smaller==0 either health or taste doesnt get in, the timeIn > (rt+ndt)
tInsmaller[s,i] <- step(timeInStep[s,i])
li.hat_h[s,i] <- tInsmaller[s,i] * timeInStep[s,i] * vd_h2[s,i] + ( smaller[s,i] * timeStep[s,i] * vd_h[s,i] )
li.hat_t[s,i] <- tInsmaller[s,i] * timeInStep[s,i] * vd_t2[s,i] + ( smaller[s,i] * timeStep[s,i] * vd_t[s,i] )
li.hat[s,i] <- li.hat_h[s,i] + li.hat_t[s,i]
# The linear regression of the value taste and value health
vd_h[s,i] <- b1.p[s,idxP[i]] * wd_S[i]
vd_h2[s,i] <- b12.p[s,idxP[i]] * wd_S[i]
vd_t[s,i] <- b2.p[s,idxP[i]] * td_S[i]
vd_t2[s,i] <- b22.p[s,idxP[i]] * td_S[i]
}
}
}
model{
fake <- 0
}
```

I save the output in Rdata and then call loo in R script. I suppose this should be correct because it was used multiple times for different model comparison:

```
Nsamps = 6000
conds = c(1)
for (cond in conds){
waic_all <- vector(mode="list", length=length(all_models))
names(waic_all) <- all_models
loo_all <- vector(mode="list", length=length(all_models))
names(loo_all) <- all_models
ll_alls <- vector(mode="list", length=length(all_models))
for (i_model in 1:length(all_models)){
model <- all_models[i_model]
folder.results <- file.path(pathToFolder, model)
subG = 1 # initialize with subject 1 and append
resultsFileName <- paste("post_ll_",model,"_", cond, "_", subG, ".RData", sep="")
load(file.path(folder.results, resultsFileName))
ll_alls[[model]] = matrix(post_ll$mcmc[[1]], nrow = Nsamps, ncol = length(post_ll$mcmc[[1]])/Nsamps, byrow = FALSE)
for (subG in 2:n_sub){
resultsFileName <- paste("post_ll_",model,"_", cond, "_", subG, ".RData", sep="")
load(file.path(folder.results, resultsFileName))
ll = matrix(post_ll$mcmc[[1]], nrow = Nsamps, ncol = length(post_ll$mcmc[[1]])/Nsamps, byrow = FALSE)
ll_alls[[model]] <- cbind(ll_alls[[model]],ll)
}
}
# remove and replace inf ----
inf_row_inds <- vector(mode="list", length=length(all_models))
for (i_model in 1:length(all_models)){
model <- all_models[i_model]
ll_all <- ll_alls[[model]]
# find rows with Inf or -Inf or ll_all and does not lead to "replacement has length zero" error
inf_row_inds[[model]] <- which(is.infinite(rowSums(ll_all)))
}
# get the union of all inf_row_inds
inf_row_ind_union <- Reduce(union, inf_row_inds)
if (length(inf_row_ind_union)%%3 != 0){
# randomly sample from 1:6000 which doesnt exist in inf_row_ind_union and add to inf_row_ind_union
inf_row_ind_union <- c(inf_row_ind_union,
sample(setdiff(1:6000, inf_row_ind_union),
3 - length(inf_row_ind_union)%%3))
}
# remove inf_row_ind_union rows from all elements of ll_alls
for (i_model in 1:length(all_models)){
model <- all_models[i_model]
ll_all <- ll_alls[[model]]
if (length(inf_row_ind_union) > 0) {
ll_all <- ll_all[-inf_row_ind_union, ] # Exclude rows if inf_row_ind_union is not empty
}
waic_all[[model]] =waic(ll_all)
chainIDs = c(rep(1:3, each = dim(ll_all)[1]/3))
loo_all[[model]] =loo(ll_all, r_eff = relative_eff(ll_all, chain_id = chainIDs))
}
###-----------------------------------------
comp <- loo_compare(loo_all)
print(comp)
}
```

Here is the output of loo:

Warning messages:

1:

204 (4.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

2: Some Pareto k diagnostic values are too high. See help(‘pareto-k-diagnostic’) for details.

3:

398 (7.8%) p_waic estimates greater than 0.4. We recommend trying loo instead.

4: Some Pareto k diagnostic values are too high. See help(‘pareto-k-diagnostic’) for details.

[The second one is the winning model]

I hope these help!