fit_ertmpt_SBC: Simulation-based calibration for RT-MPT models

In rtmpt: Fitting (Exponential/Diffusion) RT-MPT Models

Simulation-based calibration for RT-MPT models

Description

Simulate data from RT-MPT models using ertmpt_model objects. The difference to sim_ertmpt_data is that here only scalars are allowed. This makes it usable for simulation-based calibration (SBC; Talts et al., 2018). You can specify the random seed, number of subjects, number of trials, and some parameters (same as prior_params from fit_ertmpt).

Usage

fit_ertmpt_SBC(
  model,
  seed,
  n.eff_samples = 99,
  n.chains = 4,
  n.iter = 5000,
  n.burnin = 200,
  n.thin = 1,
  Rhat_max = 1.05,
  Irep = 1000,
  n.subj = 40,
  n.trials = 30,
  prior_params = NULL,
  sim_list = NULL
)

Arguments

model	A list of the class ertmpt_model.
seed	Random seed number.
n.eff_samples	Number of effective samples. Default is 99, leading to 100 possible ranks (from 0 to 99).
n.chains	Number of chains to use. Default is 4. Must be larger than 1 and smaller or equal to 16.
n.iter	Number of samples per chain. Default is 5000. Must be larger or equal to n.eff_samples.
n.burnin	Number of warm-up samples. Default is 200.
n.thin	Thinning factor. Default is 1.
Rhat_max	Maximal Potential scale reduction factor: A lower threshold that needs to be reached before the actual sampling starts. Default is 1.05
Irep	Every Irep samples an interim state with the current maximal potential scale reduction factor is shown. Default is 1000. The following statements must hold true for Irep: n.burnin is smaller than or equal to Irep, Irep is a multiple of n.thin and n.iter is a multiple of Irep / n.thin.
n.subj	Number of subjects. Default is 40.
n.trials	Number of trials per tree. Default is 30.
prior_params	Named list of parameters from which the data will be generated. This must be the same named list as prior_params from fit_ertmpt and has the same defaults. It is not recommended to use the defaults since they lead to many probabilities close or equal to 0 and/or 1 and to RTs close or equal to 0. Allowed parameters are: mean_of_exp_mu_beta: This is the expected exponential rate (E(exp(beta)) = E(lambda)) and 1/mean_of_exp_mu_beta is the expected process time (1/E(exp(beta)) = E(tau)). The default mean is set to 10, such that the expected process time is 0.1 seconds. var_of_exp_mu_beta: The group-specific variance of the exponential rates. Since exp(mu_beta) is Gamma distributed, the rate of the distribution is just mean divided by variance and the shape is the mean times the rate. The default is set to 100. mean_of_mu_gamma: This is the expected mean parameter of the encoding and response execution times, which follow a normal distribution truncated from below at zero, so E(mu_gamma) < E(gamma). The default is 0. var_of_mu_gamma: The group-specific variance of the mean parameter. Its default is 10. mean_of_omega_sqr: This is the expected residual variance (E(omega^2)). The default is 0.005. var_of_omega_sqr: The variance of the residual variance (Var(omega^2)). The default is 0.01. The default of the mean and variance is equivalent to a shape and rate of 0.0025 and 0.5, respectivly. df_of_sigma_sqr: degrees of freedom for the individual variance of the response executions. The individual variance follows a scaled inverse chi-squared distribution with df_of_sigma_sqr degrees of freedom and omega^2 as scale. 2 is the default and it should be an integer. sf_of_scale_matrix_SIGMA: The original scaling matrix (S) of the (scaled) inverse Wishart distribution for the process related parameters is an identity matrix S=I. sf_of_scale_matrix_SIGMA is a scaling factor, that scales this matrix (S=sf_of_scale_matrix_SIGMA*I). Its default is 1. sf_of_scale_matrix_GAMMA: The original scaling matrix (S) of the (scaled) inverse Wishart distribution for the encoding and motor execution parameters is an identity matrix S=I. sf_of_scale_matrix_GAMMA is a scaling factor that scales this matrix (S=sf_of_scale_matrix_GAMMA*I). Its default is 1. prec_epsilon: This is epsilon in the paper. It is the precision of mu_alpha and all xi (scaling parameter in the scaled inverse Wishart distribution). Its default is also 1. add_df_to_invWish: If P is the number of parameters or rather the size of the scale matrix used in the (scaled) inverse Wishart distribution then add_df_to_invWish is the number of degrees of freedom that can be added to it. So DF = P + add_df_to_invWish. The default for add_df_to_invWish is 1, such that the correlations are uniformly distributed within [-1, 1].
sim_list	Object of class ertmpt_sim. This is also an output object. Can be used to re-fit the model if n.eff_samples was not achieved in a previous fitting attempt. It will then use the data stored in this object. Default is NULL and this object will be created anew.

Value

A list of the class ertmpt_sbc containing

• ranks: the rank statistic for all parameters,

• sim_list: an object of the class ertmpt_sim,

• fit_list: an object of the class ertmpt_fit,

• specs: some specifications like the model, seed number, etc.,

Author(s)

Raphael Hartmann

References

Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. arXiv preprint arXiv:1804.06788.

Examples

########################################################################################
# Detect-Guess variant of the Two-High Threshold model.
# The encoding and motor execution times are assumed to be different for each response.
########################################################################################

mdl_2HTM <- "
# targets
d+(1-d)*g     ; 0
(1-d)*(1-g)   ; 1

# lures
(1-d)*g       ; 0
d+(1-d)*(1-g) ; 1

# d: detect; g: guess
"

model <- to_ertmpt_model(mdl_file = mdl_2HTM)

params <- list(mean_of_exp_mu_beta = 10,
               var_of_exp_mu_beta = 10,
               mean_of_mu_gamma = 0.5,
               var_of_mu_gamma = 0.0025,
               mean_of_omega_sqr = 0.005,
               var_of_omega_sqr = 0.000025,
               df_of_sigma_sqr = 10,
               sf_of_scale_matrix_SIGMA = 0.1,
               sf_of_scale_matrix_GAMMA = 0.01,
               prec_epsilon = 10,
               add_df_to_invWish = 5)

R = 2 # typically 2000 with n.eff_samples = 99, but this will run many days
rank_mat <- matrix(NA, ncol = 393, nrow = 2)
for (r in 1:R) {
  SBC_out <- fit_ertmpt_SBC(model, seed = r*123, prior_params = params,
                           n.eff_samples = 99, n.thin = 5,
                           n.iter = 5000, n.burnin = 2000, Irep = 5000)
  rank_mat[r, ] <- SBC_out$ranks
}

rtmpt documentation built on May 29, 2024, 3:01 a.m.