@paul.buerkner, I hate to keep pestering you over this family of model, but I wanted to see if I could get your opinion on something before I throw in the towel. I have implemented a variant of the model described in this thread and have been preparing to use it for a paper I am working on. It fits the data in question quite well and the results are reasonable. However, being a little uptight, I usually run parameter recovery simulations for any exotic analyses I might use to ensure there are no biases. Long story short, my model is producing biased estimates and I canât seem to sort out why.
Simulations
I have appended an R script with my current code and simulations, but I quote the relevant portion below:
# This function generates a single subject with a specific d' and mixing proportion (rho)
gen_2afc_sdt_dat = function(d=1, rho=.4, crits=seq(-1,1,.5), ntrials=100, scale=1)
{
expand.grid(trial=1:ntrials, is_right= c(0, 1)) %>%
mutate(signal = rnorm(n(), d, scale),
noise = rnorm(n(), 0, scale),
difference = signal - noise,
is_recollect = rbinom(n(), 1, rho),
rating = case_when(
is_right == 1 & is_recollect == 1 ~ 6,
is_right == 1 & difference < crits[1] ~ 1,
is_right == 1 & difference < crits[2] ~ 2,
is_right == 1 & difference < crits[3] ~ 3,
is_right == 1 & difference < crits[4] ~ 4,
is_right == 1 & difference < crits[5] ~ 5,
is_right == 1 ~ 6,
is_right == 0 & is_recollect == 1 ~ 1,
is_right == 0 & difference < crits[1] ~ 6,
is_right == 0 & difference < crits[2] ~ 5,
is_right == 0 & difference < crits[3] ~ 4,
is_right == 0 & difference < crits[4] ~ 3,
is_right == 0 & difference < crits[5] ~ 2,
is_right == 0 ~ 1,
TRUE ~ 999
)) -> sim2afc
return(sim2afc)
}
# This allows simulation of multilevel data wherein participants vary
gen_2afc_sdt_dat_multi = function(n, m_d=1, sd_d=.3, m_r=.4, sd_r=.1, crits=seq(-.5,.5,.25), ntrials=100, scale=1)
{
subj_d = rnorm(n, m_d, sd_d)
subj_r = rnorm(n, m_r, sd_r)
all_dat = NULL
for(i in 1:n)
{
temp = gen_2afc_sdt_dat(subj_d[i], subj_r[i], crits, ntrials, scale)
temp$id = i
all_dat = rbind(all_dat, temp)
}
return(all_dat)
}
Basically, participants are given two words on each trial and must indicate which they had studied previously using a confidence rating (1 = Very Sure Left and 6 = Very Sure Right). I assume the response is driven by a mixture of two processes: (a) a process modelled by an ordinal regression based on the difference in activation of either response candidate (the difference of which is called dâ); or, (b) a process that always returns a high confidence judgment in favour of the target (which would be 1 if it is on the left and 6 if it is on the right). That is what I simulate above.
The Model
The model I fit in BRMS reflects a mixture model between (1) an ordinal regression with a scale of 2**.5 (the difference between two normal distributions with SD = 1 will be 2**.5), and, (2) a distribution that always returns a high confidence judgment in favour of the correct item (if the item is on the left, this would be 1; if the item is on the right, this would be 6). Theta, reflects the mixture proportion and dâ is calculated by arbitrarily contrasting trials where the item is on the right to trials where the item is on the left. Estimates of dâ are converted post-hoc into its original scale by dividing by 2**.5 (reversing the impact of our disc parameter). The code is here:
# Simulate some data
sim_dat = gen_2afc_sdt_dat_multi(30)
# Define the custom family
rec2_dist <- custom_family(
"rec2_dist",
type = "int",
special = 'ordinal',
vars = 'is_right[n]'
)
rec2_dist$threshold='equidistant'
# Define the corresponding Stan density function
stan_funs <- "
real rec2_dist_lpmf(int y, real mu, int is_right) {
if(is_right == 1) {
if(y == 6)
{
return(log(.9995));
} else {
return(log(.0001));
}
} else {
if(y == 1)
{
return(log(.9995));
} else {
return(log(.0001));
}
}
}
"
# Put them together for inclusion in our model; we also need to provide a vector
# indicating which items are right (used in the above function); I can take this from
# the earlier simulated data so long as the analyzed data have similar structure
stanvars = stanvar(scode = stan_funs, block = "functions") +
stanvar(sim_dat$is_right, "is_right", scode = " int is_right[N];")
# Create the mixture
mix <- mixture(rec2_dist, cumulative(link="probit"))
# Set some priors
Prior <- c(prior(normal(1, 1), class = "b", dpar='mu2'),
prior(normal(0, 1.5), class = "Intercept", dpar='mu2'),
prior(logistic(0, 1), class = 'Intercept', dpar='theta1'),
prior(normal(0, 1), class = "sd", dpar='mu2'),
prior(normal(0, 1), class = 'sd', dpar='theta1'),
prior(lkj(4), class = "cor")
)
# Fit the model
mod_2afc <- brms(bf(rating ~ 1,
mu1 ~ 0,
mu2 ~ is_right + (is_right|id),
theta1 ~ 1 + (1|id),
disc2=(1/2**.5)),
family = mix,
data = sim_dat,
stanvars = stanvars,
prior = Prior,
cores=4, iter=1000, init_r=.5,
control = list(adapt_delta=.9))
Bottom line is that the model fits fine â both to simulated data and to my actual data â but produces biased estimates in the former. Under some circumstances it overestimates theta/rho (I observed this when dâ is small and criteria are tightly clustered, producing an ROC curve with most responses at 1 or 6) whereas under other circumstances it overestimates dâ (this is observed in the simulation provided wherein corrected dâ is estimated closer to 1.3 than to 1 after correction). I have been trying to sort this out for weeks but am at my witâs end. Is there something obvious I am missing?
PS: I have not incorporated a threshold parameter or converted the distribution into a self-encapsulated form yet; I figure the first step is to develop a working model.
2AFC_DPSDT_Sims.R (9.4 KB)