Dear aseybolt,
Thank you so much for your help and interesting comments. I very much appreciate.
You are right, it is a within-subject design. Let me give you a little more information about that project. Previous research in developmental psychology has shown that age-related changes in mean RT in various tasks are well described by an exponential function of the form:
More specifically, it has been shown that (1) the rate of decay \beta remains the same whatever the task, and is equal to .334 (beautiful isn’t it?) (2) asymptotic performance levels are attained at the adult age.
We collected RT data in a decision-making task (within-subject design, two experimental conditions) from children of different ages. I then fit each individual dataset with a popular decision-making model to decompose the data into components of processing (this model has 5 parameters called v, a, ter, max_ampl, tau; each parameter quantifies a specific processing component) and determine whether (1) components of processing also follow an exponential-like function, and (2), whether all these components decay at the same rate as the mean RT data. To remove some noise, I averaged data of children from the same grade.
I thus fit the following model to data:
with pm.Model() as model_exp:
α_v = pm.HalfCauchy('α_v', 2)
α_a = pm.HalfCauchy('α_a', 2)
α_Ter = pm.HalfCauchy('α_Ter', 2)
α_max_ampl = pm.HalfCauchy('α_max_ampl', 2)
α_tau = pm.HalfCauchy('α_tau', 2)
α_c = pm.HalfCauchy('α_c',2)
α_i = pm.HalfCauchy('α_i', 2)
β = pm.Normal('β', 0.334, 0.1)
γ_v= pm.Normal('γ_v', mean_param_per_group_model[-1, 0], .1)
γ_a= pm.Normal('γ_a', mean_param_per_group_model[-1, 1], .1)
γ_Ter= pm.Normal('γ_Ter', mean_param_per_group_model[-1, 2], .1)
γ_max_ampl= pm.Normal('γ_max_ampl', mean_param_per_group_model[-1, 3], .1)
γ_tau= pm.Normal('γ_tau', mean_param_per_group_model[-1, 4], .1)
γ_c = pm.Normal('γ_c', Mean_RT_comp_General[-1], .1)
γ_i = pm.Normal('γ_i', Mean_RT_incomp_General[-1], .1)
σ_v = pm.HalfNormal('σ_v', .1)
σ_a = pm.HalfNormal('σ_a', .1)
σ_Ter = pm.HalfNormal('σ_Ter', .1)
σ_max_ampl = pm.HalfNormal('σ_max_ampl', .1)
σ_tau = pm.HalfNormal('σ_tau', .1)
σ_c = pm.HalfNormal('σ_c', .1)
σ_i = pm.HalfNormal('σ_i', .1)
μ_v = -α_v*pm.math.exp(-β*age_categ) + γ_v
μ_a = α_a*pm.math.exp(-β*age_categ) + γ_a
μ_Ter = α_Ter*pm.math.exp(-β*age_categ) + γ_Ter
μ_max_ampl = α_max_ampl*pm.math.exp(-β*age_categ) + γ_max_ampl
μ_tau = α_tau*pm.math.exp(-β*age_categ) + γ_tau
μ_c = α_c*pm.math.exp(-β*Age_General) + γ_c
μ_i = α_i*pm.math.exp(-β*Age_General) + γ_i
y_v = pm.Normal('y_v', μ_v, σ_v, observed=mean_param[:, 0])
y_a = pm.Normal('y_a', μ_a, σ_a, observed=mean_param[:, 1])
y_Ter = pm.Normal('y_Ter', μ_Ter, σ_Ter, observed=mean_param[:, 2])
y_max_ampl = pm.Normal('y_max_ampl', μ_max_ampl, σ_max_ampl, observed=mean_param[:, 3])
y_tau = pm.Normal('y_tau', μ_tau, σ_tau, observed=mean_param[:, 4])
y_meanRT_comp = pm.Normal('y_meanRT_comp', μ_c, σ_c, observed=Mean_RT_comp)
y_meanRT_incomp = pm.Normal('y_meanRT_incomp', μ_i, σ_i, observed=Mean_RT_incomp)
if __name__ == "__main__":
with model_exp:
trace_exp = pm.sample (4000, chains = 4, tune = 2000, target_accept = .95)
az.plot_trace(trace_exp)
which fits the data well and suggests that every processing component evolves at the same rate:
I wonder what I could do to provide more convincing evidence in favor of this conclusion. I am thinking of comparing this model with a variant in which the decay rate \beta is free to vary. However, I can’t compute a LOO statistic (I get the following error: “Data must include log_likelihood in sample_stats” when doing az.LOO). In addition, I would appreciate some advice about the best way to analyze residuum from my posterior predictive checks (right now I have just plotted an histogram of residuals, they are normally distributed and centered on 0 so everything looks ok).
Many thanks!