A common problem when modeling response time data with the wiener-diffusion model are contaminant response times. Especially fast outliers (sometimes called anticipatory responses) can essentially completely determine the non-decision time parameter (\tau). A common approach for dealing with this problem is assuming that the data comes from a mixture between the diffusion process and a uniform distribution bounded at the minimum and maximum observed response time (e.g., Ratcliff & Tuerlinkcx, 2002). It is usually assumed that the vast majority of trials come from the Wiener process (e.g., over 95%) and only a small faction comes from the uniform distribution.
I have generated such data and tried to recover the parameters in Stan. This works for the very simple case of having only one condition (i.e., intercept only for all parameters of the wiener model) after some tinkering with the priors, but not for the case with two conditions which differ in their drift rate. Not working here means that all post-warm up samples are divergent and the chains appear to be stuck.
I first generate the data in R.
library("rtdists")
set.seed(555)
rts <- rbind(
cbind(rdiffusion(500, a=1, v=2, t0=0.5), condition = "A"),
cbind(rdiffusion(500, a=1, v=-1.5, t0=0.5), condition = "B"))
rts$response2 <- as.numeric(rts$response)-1
prop_outlier <- 0.02
n_out <- round(prop_outlier*nrow(rts))
rts[sample(nrow(rts), n_out),"rt"] <- runif(n_out, 0.3, 1.5)
rts[sample(nrow(rts), n_out),"response2"] <- sample(0:1, size = n_out, replace = TRUE)
min(rts$rt) # [1] 0.3520957
It is straight forward to fit the data with the Wiener model alone using brms. However, thanks to the contaminants, the parameters do not recover very well (without the contaminants, the recovery would be basically perfect).
library("brms")
formula <- bf(rt | dec(response2) ~ 0+condition,
bs ~ 1, ndt ~ 1, bias ~ 1)
prior <- c(
prior("normal(0, 2)", class = "b"),
set_prior("normal(1, 1)", class = "Intercept", dpar = "bs"),
set_prior("normal(0.2, 0.2)", class = "Intercept", dpar = "ndt"),
set_prior("normal(0.5, 0.1)", class = "Intercept", dpar = "bias")
)
tmp_dat <- make_standata(formula,
family = wiener(link_bs = "identity",
link_ndt = "identity",
link_bias = "identity"),
data = rts, prior = prior)
initfun <- function() {
list(
b = rnorm(tmp_dat$K),
temp_bs_Intercept = runif(tmp_dat$K_bs, 1, 2),
temp_ndt_Intercept = runif(tmp_dat$K_ndt, 0.1, 0.2),
temp_bias_Intercept = rnorm(tmp_dat$K_bias, 0.5, 0.1),
lambda = runif(1, 0, 0.05)
)
}
fit_orig <- brm(formula,
data = rts,
family = wiener(link_bs = "identity",
link_ndt = "identity",
link_bias = "identity"),
prior = prior, inits = initfun,
iter = 1000, warmup = 500,
chains = 4, cores = 4,
control = list(max_treedepth = 15))
summary(fit_orig)
# Population-Level Effects:
# Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# bs_Intercept 1.53 0.03 1.48 1.58 1225 1.00
# ndt_Intercept 0.34 0.00 0.34 0.34 2000 1.00
# bias_Intercept 0.45 0.01 0.43 0.48 1024 1.00
# conditionA 1.84 0.10 1.64 2.04 960 1.00
# conditionB -1.01 0.09 -1.19 -0.83 1138 1.00
Next, I adapt the brms code by changing the likelihood to a mixture. So the likelihood becomes:
for (n in 1:N) {
target += log_mix(lambda,
uniform_lpdf(Y[n] | min_Y, max_Y),
wiener_diffusion_lpdf(Y[n] | dec[n], bs[n], ndt[n], bias[n], mu[n])
);
}
Where lambda is defined in the parameters block as real<lower=0,upper=1> lambda; and max_Y in the transformed data block as real max_Y = max(Y);.
The prior for lambda is target += normal_lpdf(lambda | 0.02, 0.01);
The full model code is posted below.
If fitting the model then I get only divergences and cannot find a way to deal with them (e.g., changing priors or control options does not appear to work).
library("rstan")
fit_mix <- stan(model_code = mixture_code,
data = tmp_dat,
chains = 4,
warmup = 500, iter = 1000,
init = initfun)
Warning messages:
1: There were 2000 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
3: Examine the pairs() plot to diagnose sampling problems
Any help for getting this model to work would be very much appreciated.
The full Stan model (which needs to be in a string mixture_code for the above code to work) is below.
// adapted from brms 2.4.0 generated code
functions {
real wiener_diffusion_lpdf(real y, int dec, real alpha,
real tau, real beta, real delta) {
if (dec == 1) {
return wiener_lpdf(y | alpha, tau, beta, delta);
} else {
return wiener_lpdf(y | alpha, tau, 1 - beta, - delta);
}
}
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=0,upper=1> dec[N]; // decisions
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int prior_only; // should the likelihood be ignored?
}
transformed data {
real min_Y = min(Y);
real max_Y = max(Y);
}
parameters {
vector[K] b; // population-level effects
real temp_bs_Intercept; // temporary intercept
real temp_ndt_Intercept; // temporary intercept
real temp_bias_Intercept; // temporary intercept
real<lower=0,upper=1> lambda;
}
transformed parameters {
}
model {
vector[N] mu = X * b;
vector[N] bs = temp_bs_Intercept + rep_vector(0, N);
vector[N] ndt = temp_ndt_Intercept + rep_vector(0, N);
vector[N] bias = temp_bias_Intercept + rep_vector(0, N);
// priors including all constants
target += normal_lpdf(b | 0, 2);
target += normal_lpdf(temp_bs_Intercept | 1, 1);
target += normal_lpdf(temp_ndt_Intercept | 0.2, 0.2);
target += normal_lpdf(temp_bias_Intercept | 0.5, 0.1);
target += normal_lpdf(lambda | 0.02, 0.01);
// likelihood including all constants
if (!prior_only) {
for (n in 1:N) {
target += log_mix(lambda,
uniform_lpdf(Y[n] | min_Y, max_Y),
wiener_diffusion_lpdf(Y[n] | dec[n], bs[n], ndt[n], bias[n], mu[n])
);
}
}
}
generated quantities {
// actual population-level intercept
real b_bs_Intercept = temp_bs_Intercept;
// actual population-level intercept
real b_ndt_Intercept = temp_ndt_Intercept;
// actual population-level intercept
real b_bias_Intercept = temp_bias_Intercept;
}
