Indices of model fit quality for SSMs (R2)

[DominiqueMakowski]

I've been looking in the literature for some examples or discussions about assessing the fit quality and/or the predictive power of SSMs, but I found nothing.

I think the Bayes factor comparison allowed by Pigeons is one amazing feature, but having an "absolute" index (like an R2) would be a good addition. Especially since a model (especially basic ones with no predictor variables) could in principle be very "good" (by posterior predictive checks; the predicted distribution matches exactly the data distribution), but have an R2 of 0 (each observation is not predicted at all).

The two problems I can foresee is that SSMs are fundamentally distributional, but that's not necessarily an issue in itself. Bayesian R2 exist and does return a distribution as well, from which one can compute the mean & CI.

The second specificity is that SSMs jointly predict two outcomes (RT and choice). But I wonder if it would make sense, in principle, to compute say the R2 for the RTs (within each choice), and then the R2 for the choice (an R2 variant made for categorical outcomes) separately. Of course, by doing that separately we lose some nuance, but is this approach fundamentally misguided?

[itsdfish]

@DominiqueMakowski, I think you are correct on both accounts. Using measures of absolute fit is necessary to determine whether the best fitting model is adequate. Moreover, absolute measures can help diagnose deficiencies in the model. I also think you are right about R2 not being the right measure because we want to know how well the model fits the entire distribution. I think there are a few reasonable approaches. One would be to compare the probability mass between different quantiles of the conditional rt distributions. For example, the probability mass between the .4 and .45 might be .15 for correct responses, but the 95% credible interval of the posterior predictive distribution for this statistic might range between .16 and .19. You could possibly do something with the conditional CDF, but a downside in my opinion is that a misfit in the leading edge of the distribution would potentially carry through (e.g., shift the CDF to the left or right for the values > t).
A similar idea might be to use the quantile probability functions developed by Ratcliff (see Figure 6 in The Diffusion Decision Model: Theory and Data for Two-Choice Decision Tasks).

Yet another idea might be to use a distance metric for the conditional distributions. The challenge would be selecting the right metric. KL divergence is popular, but has the disadvantage of being asymmetric. Wasserstein distance is symmetric, but is more challenging to compute.

[itsdfish]

Building on the idea the quantile idea described above, I like this the qq-plot envelope described here. You can make one for each response option.

[itsdfish]

Here is an example using the Wald. Some type of interpolation might be needed to make it work for the more general case where there are multiple options. The problem is that the number of observations per option vary across posterior predictions. There are some improvements that can be made with using a ribbon in the plot

using Plots
using SequentialSamplingModels
using Turing
using DynamicPPL


function compute_qq(n_obs, ν, α, τ)
    return map(_ -> rand(Wald(ν, α, τ)), 1:n_obs) |> sort!
end

n_obs = 100
dist = Wald(ν=3.0, α=.5, τ=.130)
data = rand(dist, n_obs)
sort!(data)

@model function model(data; min_rt = minimum(data))
    # Priors
    ν ~ truncated(Normal(3, 1), 0.0, Inf)
    α ~ truncated(Normal(.5, .5), 0.0, Inf)
    τ  ~ Uniform(0.0, min_rt)

    # Likelihood
    data ~ Wald(;ν, α, τ)

    return compute_qq(length(data), ν, α, τ)
end

chain = sample(model(data), NUTS(1000, .85), MCMCThreads(), 1000, 1)

post_pred = generated_quantities(model(data), chain)

post_pred = reduce(vcat, transpose(post_pred))

lb = map(c -> quantile(post_pred[:,c], .025), 1:n_obs)
ub = map(c -> quantile(post_pred[:,c], .975), 1:n_obs)

plot(data, data, lims = (0, 1.5), leg=false, grid=false, ribbon=(data .- lb, ub .- data))
plot!(0:.1:2, 0:.1:2, color=:black, linestyle=:dash)


# mispecified model
uniform_data = sort!(rand(n_obs))
plot(uniform_data, uniform_data, lims = (0, 1.5), leg=false, grid=false, ribbon=(uniform_data .- lb, ub .- uniform_data))
plot!(0:.1:2, 0:.1:2, color=:black, linestyle=:dash)