I’m working on a task where participants make risky choices on a continuous scale (0 to 50) across 20 trials. My model includes standard risk-related parameters and allows them to drift over trials so each parameter has a baseline value and a corresponding drift term.
The issue I’m facing is with high Pareto k̂ values, but only for some subjects, and even then, only for a subset of their trials, sometimes just a single trial.
I’m also running around 15 variants of the same template model for model comparison purposes. This includes models that are intentionally less predictive, to establish a benchmark against the best-performing (winning) model.
I’ve tried mitigating this by:
-
Adding priors to all base and drift parameters,
-
Increasing randomness
-
And inspecting the problem trials more closely.
Still, I’m unsure how best to handle these local issues without discarding data. Any advice or insights would be greatly appreciated! I assume in this case it would not be acceptable to exclude outliers from data analysis? As this might be related to low number of trials I have per subjects which gives more weight to extreme trials. I would not want to remove these datapoints if I don’t have to as it might be related to exploration strategy by the participant, and there is very few data already (20 trials only per sub).
Thanks!
My stan code:
MODEL_TEMPLATE = “”"
data {{
int<lower=1> N; // number of observations
int<lower=1> B; // number of bet options
vector[N] xmult; // multiplier values
array[N] int<lower=0, upper=B-1> bet; // chosen bet (0-indexed to match your data)
vector[N] trial_scale; // scaled trial numbers for drift
}}
parameters {{
{params}
}}
transformed parameters {{
vector[N] theta = {theta_block};
vector[N] rho = {rho_block};
vector[N] beta = {beta_block};
{clamp}
}}
model {{
{priors}
// Likelihood
for (n in 1:N) {{
vector[B] U;
for (b in 0:(B-1)) {{ // 0-based loop to match data
real k = fmax(50 - b, 1e-6);
real g = fmax(0.5 * (xmult[n] * b), 1e-6);
if (abs(rho[n] - 1.0) < 1e-5) {{
U[b+1] = pow(k, theta[n]) * pow(g, 1 - theta[n]); // Store in 1-based array
}} else {{
real temp = theta[n] * pow(k, rho[n]) + (1-theta[n]) * pow(g, rho[n]);
U[b+1] = pow(fmax(temp, 1e-10), 1.0 / rho[n]);
}}
}}
target += beta[n] * U[bet[n]+1] - log_sum_exp(beta[n] * U); // bet[n] is 0-based, U is 1-based
}}
}}
generated quantities {{
vector[N] log_lik;
for (n in 1:N) {{
vector[B] U;
for (b in 0:(B-1)) {{ // 0-based loop to match data
real k = fmax(50 - b, 1e-6);
real g = fmax(0.5 * (xmult[n] * b), 1e-6);
if (abs(rho[n] - 1.0) < 1e-5) {{
U[b+1] = pow(k, theta[n]) * pow(g, 1 - theta[n]);
}} else {{
real temp = theta[n] * pow(k, rho[n]) + (1-theta[n]) * pow(g, rho[n]);
U[b+1] = pow(fmax(temp, 1e-10), 1.0 / rho[n]);
}}
}}
log_lik[n] = beta[n] * U[bet[n]+1] - log_sum_exp(beta[n] * U); // bet[n] is 0-based, U is 1-based
}}
}}
“”"
model_variants = [
{
“name”: “6_theta_rho_noise_drift_theta_drift_rho_drift_noise_corrected_2”,
“params”: “real<lower=0.01, upper=0.99> theta0;\n real drift_theta;\n real<lower=0.01, upper=0.99> rho0;\n real drift_rho;\n real<lower=0.01, upper=5> beta0;\n real drift_beta;”,
“theta_block”: “theta0 + drift_theta * trial_scale”,
“rho_block”: “rho0 + drift_rho * trial_scale”,
“beta_block”: “beta0 + drift_beta * trial_scale”,
“clamp”: “for (n in 1:N) {{ theta[n] = fmin(fmax(theta[n], 0.01), 0.99); rho[n] = fmin(fmax(rho[n], 0.01), 0.99); beta[n] = fmax(beta[n], 0.01); }}”,
“priors”: “theta0 ~ beta(2, 2);\n rho0 ~ beta(2, 2);\n beta0 ~ gamma(0.5, 2);\n drift_theta ~ normal(0, 0.09);\n drift_rho ~ normal(0, 0.09);\n drift_beta ~ normal(0, 1.0);”
},
]