Parameter recovery
Being based on RStan, RStanTVA can be used in different ways, i.e. using maximum-likelihood estimation or Bayesian inference of different hierarchical complexity (single-condition, single-participant/fixed-effects, fully-hierarchical/mixed-effects). Hence we can also recover model parameters in different ways.
Maximum-likelihood estimation (MLE)
Using MLE, we can define four different inference models:
RStanTVA-ML: Single participant, single condition
The most straightforward model uses no regularization by priors
(priors = FALSE) and conducts single fits for all model
parameters:
RStanTVA-MLN: Single participant, condition as fixed effect
By defining a regression model for \(\log
C\), which includes an intercept and a treatment effect of “high”
vs. “low” (log(C) ~ 1 + condition), we can model both
experimental conditions in a single model:
RStanTVA-MLR: Single participant, single condition, with regularization (priors)
Omitting the priors = FALSE argument falls back to the
default priors = NULL, which uses vaguely informative
default priors for all parameters. The model is otherwise identical to
the baseline RStanTVA-ML model above:
RStanTVA-MLRN: Single participant, condition as fixed effect, with regularization (priors)
When using regularizing priors, if a parameter is described by an equation such as \(C\), we need to specify priors for intercepts and slopes:
priors <-
prior(normal(0,.5),dpar=C)+
prior(normal(4.5,.6),dpar=C,coef=Intercept)
m_nested_reg <- stantva_model(
formula = list(log(C) ~ 1 + condition),
locations = 6,
task = "pr",
regions = list(left = 1:3, right = 4:6),
w_mode = "regions",
t0_mode = "gaussian",
K_mode = "probit",
save_log_lik = FALSE,
priors = priors
)For brevity, we are only reporting results of the RStanTVA-MLRN fits here, since those have been found to be most reliable (see manuscript). The model can be fitted to the individual data sets (each participant) as follows:
fit_mle_nested_reg <- lapply(1:50, function(i) {
d <- tva_recovery %>% filter(subject == i)
repeat {
p <- optimizing(m_nested_reg, d)
if(p$return_code == 0L) break
}
tibble(subject = i, param = c("italic(C)[low]","italic(C)[high]","alpha","italic(w)[lambda]","italic(t)[0]","sigma[0]","italic(K)"), fitted_value = c(exp(p$par["C_Intercept"]), exp(p$par["C_Intercept"]+p$par["C_conditionhigh"]), p$par[c("alpha","r[2]","mu0","sigma0")], meanprobitK(p$par["mK"],p$par["sK"],6)), converged = p$return_code == 0L)
}) %>% bind_rows() %>% left_join(true_params_narrow2)The critical function above is
optimizing(m_nested_reg, d), which uses MLE to fit
m_nested_reg to the subject-level subset
d.
RStanTVA-NB: Bayesian non-hierarchical inference
We can, in principle, fit the exact same models as above to the data
using Bayesian inference instead of MLE, simply by replacing the
optimizing(...) with a sampling(...) function
call. So the following code will fit the RStanTVA-MLRN
(m_nested_reg) model to subject-level subsets
d using HMC methods in Stan:
fit_bayesian_nested <- lapply(1:50, function(i) {
d <- tva_recovery %>% filter(subject == i)
sf <- sampling(m_nested_reg, d)
p1 <- predict(sf, data.frame(subject = i, condition = "low"), c("C","alpha", "r", "mu0","sigma0","mK","sK"))
p2 <- predict(sf, data.frame(subject = i, condition = "high"), "C")
tibble(
`italic(C)[low]` = t(p1$C),
`italic(C)[high]` = t(p2$C),
`alpha` = t(p1$alpha),
`italic(w)[lambda]` = t(p1$r[,2,1]),
`italic(t)[0]` = t(p1$mu0),
`sigma[0]` = t(p1$sigma0),
`italic(K)` = t(meanprobitK(p1$mK,p1$sK,6))
) %>%
pivot_longer(everything(), names_to = "param", values_to = "samples") %>%
rowwise(param) %>%
reframe(subject = i, posterior_sd = sd(samples), fitted_value = median(samples), cri = t(quantile(samples, c(.025,.975))), converged = Rhat(matrix(samples, ncol = sf@sim$chains)) < 1.1)
}) %>%
bind_rows() %>%
left_join(true_params_narrow2)Note that we must use the method
predict(fit, data, parameters) to predict the TVA parameter
values from the fitted parameters. This will return a posterior
distribution for each row (subject) in data and each
specified parameter. We are then calculating 95% credible intervals
(CrIs), medians, \(\hat R\) and
posteriors SDs as measures of goodness-of-fit.
RStanTVA-HB*: Bayesian hierarchical inference
For hierarchical inference, we define mixed-effects regressions for the model parameters \(C\), \(\log\alpha\), \(\mu_K\), \(\mu_0\), and \(\textrm{logit}(r)\), while keeping \(\sigma_0\) and \(\sigma_K\) fixed across all subjects:
m_hierarchical <- stantva_model(
formula = list(
C ~ 1 + condition + (1 + condition | subject),
log(alpha) ~ 1 + (1 | subject),
mK ~ 1 + (1 | subject),
mu0 ~ 1 + (1 | subject),
log(r) ~ 1 + (1 | subject)
),
locations = 6,
task = "pr",
regions = list(left = 1:3, right = 4:6),
C_mode = "equal",
w_mode = "regions",
t0_mode = "gaussian",
K_mode = "probit",
priors =
prior(normal(90, 20/sqrt2()), coef = Intercept, dpar = C) +
prior(normal(0, 10/sqrt2()), dpar = C) +
prior(gamma(2, 2/20 * sqrt2()), class = sd, coef = Intercept, dpar = C) +
prior(gamma(2, 2/10 * sqrt2()), class = sd, dpar = C) +
prior(normal(-0.4, 0.6/sqrt2()), coef = Intercept, dpar = alpha) +
prior(normal(0, 0.3/sqrt2()), dpar = alpha) +
prior(gamma(2, 2/0.6 * sqrt2()), class = sd, coef = Intercept, dpar = alpha) +
prior(gamma(2, 2/0.3 * sqrt2()), class = sd, dpar = alpha) +
prior(normal(0, 0.1/sqrt2()), coef = Intercept, dpar = r) +
prior(normal(0, 0.05/sqrt2()), dpar = r) +
prior(gamma(2, 2/0.1 * sqrt2()), class = sd, coef = Intercept, dpar = r) +
prior(gamma(2, 2/0.05 * sqrt2()), class = sd, dpar = r) +
prior(normal(3.5, 0.5/sqrt2()), coef = Intercept, dpar = mK) +
prior(normal(0, 0.25/sqrt2()), dpar = mK) +
prior(gamma(2, 2/0.5 * sqrt2()), class = sd, coef = Intercept, dpar = mK) +
prior(gamma(2, 2/0.25 * sqrt2()), class = sd, dpar = mK) +
prior(normal(20, 15/sqrt2()), coef = Intercept, dpar = mu0) +
prior(normal(0, 7.5/sqrt2()), dpar = mu0) +
prior(gamma(2, 2/15 * sqrt2()), class = sd, coef = Intercept, dpar = mu0) +
prior(gamma(2, 2/7.5 * sqrt2()), class = sd, dpar = mu0) +
prior(normal(0.5, 0.05/sqrt2()), coef = Intercept, dpar = sK) +
prior(normal(0, 0.025/sqrt2()), dpar = sK) +
prior(gamma(2, 2/0.05 * sqrt2()), class = sd, coef = Intercept, dpar = sK) +
prior(gamma(2, 2/0.025 * sqrt2()), class = sd, dpar = sK) +
prior(normal(20, 15/sqrt2()), coef = Intercept, dpar = sigma0) +
prior(normal(0, 7.5/sqrt2()), dpar = sigma0) +
prior(gamma(2, 2/15 * sqrt2()), class = sd, coef = Intercept, dpar = sigma0) +
prior(gamma(2, 2/7.5 * sqrt2()), class = sd, dpar = sigma0)
)We can then simply fit the hierarchical model
m_hierarchical to the entire data set
tva_recovery without subsetting:
fg <- sampling(m_hierarchical, tva_recovery)
fit_bayesian_hierarchical_globals <-
{
p1 <- predict(fg, tibble(subject = 1:50, condition = "low"), c("C","alpha", "r", "mu0","sigma0","mK","sK"))
p2 <- predict(fg, tibble(subject = 1:50, condition = "high"), "C")
tibble(
subject = 1:50,
`italic(C)[low]` = t(p1$C),
`italic(C)[high]` = t(p2$C),
`alpha` = t(p1$alpha),
`italic(w)[lambda]` = t(p1$r[,2,]),
`italic(t)[0]` = t(p1$mu0),
`sigma[0]` = t(p1$sigma0),
`italic(K)` = t(matrix(meanprobitK(p1$mK,p1$sK,6),ncol=50))
)
} %>%
pivot_longer(-subject, names_to = "param", values_to = "samples") %>%
rowwise(c(subject,param)) %>%
reframe(posterior_sd = sd(samples), fitted_value = median(samples), cri = t(quantile(samples, c(.025,.975))), converged = Rhat(t(samples)) < 1.1) %>%
left_join(true_params_narrow2)Comparison of MLE and Bayes
figdat <- bind_rows(
fit_mle_nested_reg %>% add_column(method = "MLRN"),
fit_bayesian_hierarchical_globals %>% add_column(method = "HB"),
fit_bayesian_nested %>% add_column(method = "NB")
) %>%
mutate(method = factor(method, levels = c("MLRN","NB","HB")))
fig <- ggplot(figdat) +
theme_minimal() +
theme(text = element_text(family = "sans", size = 10)) +
labs(x = "True value", y = "Recovered value") +
facet_wrap(method~param, ncol = 7, scales = "free", labeller = function(x) label_parsed(list(sprintf("%s~(\"%s\")", x$param, x$method)))) +
geom_abline(linetype = "dashed") +
geom_linerange(aes(x=true_value,ymin=cri[,1],ymax=cri[,2]), color = "gray50", linewidth = 0.2) +
geom_point(aes(x=true_value,y=fitted_value), size = 0.5, color = "blue") +
geom_point(aes(x=vx,y=vy), color = "transparent", figdat %>% group_by(param) %>% reframe(vx = range(true_value), vy = range(c(fitted_value,cri), na.rm = TRUE)))
print(fig)