MAX: Trial by Trial Response
Songtao Song
Data
## Loading the data
df <- read.csv("data/min_shift_trial_by_trial_response.csv", header = TRUE, sep = ",")
df$Subj <- factor(df$Subj)
df$Trial <- factor(df$Trial_num)
df$ROI <- factor(df$ROI_name)
df$Rating <- scale(df$Rating)
df <- df %>% select(Subj, Trial, ROI, cond, Rating, Response_early, Response_late)
df %>% head()print(paste0('Number of Trials: ',nlevels(df$Trial)))## [1] "Number of Trials: 16"print(paste0('Number of ROIs: ',nlevels(df$ROI)))## [1] "Number of ROIs: 70"print(paste0('Number of subjects: ',nlevels(df$Subj)))## [1] "Number of subjects: 86"print(paste0("Number of conditions: ", df %>% select(cond) %>% n_distinct()))## [1] "Number of conditions: 2"print(paste0('Length of dataset: ', dim(df)[1]))## [1] "Length of dataset: 192640"dens_plot <- function(data, var, color='#69b3a2',
title="",
xlab='Response at Early Phrase'){
var = enquo(var)
plot = data %>%
ggplot(aes(x = (!!var))) +
geom_histogram(aes(y = ..density..), bins = 50, alpha = 0.4, fill = 'gray') +
geom_density(alpha = 0.7, fill = color, color = color) +
geom_vline(aes(xintercept = mean(!!var)), linetype = "dashed", size = 1) +
ggtitle(title) +
xlab(xlab) +
xlim(c(0, 20)) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 12))
return(plot)
}
scatter_plot <- function(data, y, x, title=""){
y = enquo(y)
x = enquo(x)
plot = data %>%
ggplot(aes(x = !!x, y = !!y)) +
geom_point(shape = 1, size = 3) +
stat_smooth(method = 'lm') +
ggtitle(title) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 12))
return(plot)
}df_safe <- filter(df, cond == -0.5)
df_threat <- filter(df, cond == 0.5)plot.dens.threat.early <- dens_plot(df_threat, Response_early, title = "Trial by Trial Response of \nFake and Threat at Early phrase")
plot.dens.threat.late <- dens_plot(df_threat, Response_late, title = "Trial by Trial Response of \nFake and Threat at Late phrase",
xlab = "Response at Late Phrase")
grid.arrange(plot.dens.threat.early, plot.dens.threat.late, nrow = 1)## Warning: Removed 686 rows containing non-finite values (stat_bin).## Warning: Removed 686 rows containing non-finite values (stat_density).## Warning: Removed 2 rows containing missing values (geom_bar).## Warning: Removed 188 rows containing non-finite values (stat_bin).## Warning: Removed 188 rows containing non-finite values (stat_density).## Warning: Removed 2 rows containing missing values (geom_bar).
mean(df$Response_early)## [1] 3.509909BML
Model
lme format:
y ~ 1 + cond + rating + ( 1 + cond + rating | gr(Subj, dist = “student”)) + ( 1 + cond + rating | gr(ROI, dist = “student”))
mathematics format:
\[ y \sim \text{LogNormal}(\mu, \sigma) \]
where
\[ \mu = \alpha + \alpha_{[\text{subj}]} + \alpha_{[\text{roi}]} + (\beta_{\text{cond}} + \beta_{\text{cond}, [\text{subj}]} + \beta_{\text{cond}, [\text{roi}]}) \times \text{Cond} + (\beta_{\text{rating}} + \beta_{\text{rating}, [\text{subj}]} + \beta_{\text{rating}, [\text{roi}]}) \times \text{Rating} \]
and
\[ \begin{bmatrix} \alpha_{[j]} \\ \beta_{\text{cond}, [j]} \\ \beta_{\text{rating}, [j]} \end{bmatrix} \sim \text{Multivariate t} \begin{pmatrix} \nu_{[j]}, \text{ } \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \text{ } \mathbf{S}_{[j]} \end{pmatrix} \]
where
\[ \begin{aligned} \mathbf{S}_{[j]} & = \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & \rho_{1,2}\sigma_{\alpha_{[j]}}\sigma_{\beta_{\text{cond}, [j]}} & \rho_{1,3}\sigma_{\alpha_{[j]}}\sigma_{\beta_{\text{rating}, [j]}} \\ & \sigma_{\beta_{\text{cond}, [j]}}^2 & \rho_{2,3}\sigma_{\beta_{\text{cond}, [j]}}\sigma_{\beta_{\text{rating}, [j]}} \\ & & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \\ & = \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & & \\ & \sigma_{\beta_{\text{cond}, [j]}}^2 & \\ & & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \mathbf{R}_{[j]} \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & & \\ & \sigma_{\beta_{\text{cond}, [j]}}^2 & \\ & & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \end{aligned} \]
and
\[\alpha \sim \text{Student t}(3, \mu_{y}, 2.5)\] \[\beta_{i} \sim \text{Student t}(3, 0, 2.5)\] \[\sigma \sim \text{Half Student}(3, 0, 2.5)\] \[\nu_{[j]} \sim \text{Gamma}(3.325, 1)\] \[\sigma_{\alpha_{[j]}} \sim \text{Half Student}(3, 0, 2.5)\] \[\sigma_{\beta_{i, [j]}} \sim \text{Half Student}(3, 0, 2.5)\] \[\mathbf{R}_{[j]} \sim \text{LKJcorr}(2)\]
notice that \(\mu_{y}\) is the sample mean of response; \(i=\text{cond, rating}\), \(j=\text{subj, roi}\); and \(\text{subj}=1,2,...,N\), \(\text{roi}=1,2,...,M_{n}\). \(N\) is the number of participants; \(M_{n}\) is the number of ROI for \(nth\) participant.
mod <- '1 + Rating'
modelForm <- paste0('Response_early ~ ', mod, ' + (', mod, '| gr(Subj, dist = "student")) + (', mod, '| gr(ROI, dist = "student"))')
priorRBA <- get_prior(formula = modelForm, data = df, family = lognormal(link_sigma = "identity"))
priorRBA$prior[1] <- "student_t(3, 0, 2.5)"
priorRBA$prior[3] <- "lkj(2)"
priorRBA$prior[6:7] <- "gamma(3.325, 2)"print(priorRBA)## prior class coef group resp dpar nlpar bound
## student_t(3, 0, 2.5) b
## student_t(3, 0, 2.5) b Rating
## lkj(2) cor
## lkj(2) cor ROI
## lkj(2) cor Subj
## gamma(3.325, 2) df ROI
## gamma(3.325, 2) df Subj
## student_t(3, 1.1, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd ROI
## student_t(3, 0, 2.5) sd Intercept ROI
## student_t(3, 0, 2.5) sd Rating ROI
## student_t(3, 0, 2.5) sd Subj
## student_t(3, 0, 2.5) sd Intercept Subj
## student_t(3, 0, 2.5) sd Rating Subj
## gamma(0.01, 0.01) shape
## source
## default
## (vectorized)
## default
## (vectorized)
## (vectorized)
## default
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## defaultModel
lme format:
y ~ 1 + rating + ( 1 + rating | gr(Subj, dist = “student”)) + ( 1 + rating | gr(ROI, dist = “student”))
mathematics format:
\[ y \sim \text{LogNormal}(\mu, \sigma) \]
where
\[ \mu = \alpha + \alpha_{[\text{subj}]} + \alpha_{[\text{roi}]} + (\beta_{\text{rating}} + \beta_{\text{rating}, [\text{subj}]} + \beta_{\text{rating}, [\text{roi}]}) \times \text{Rating} \]
and
\[ \begin{bmatrix} \alpha_{[j]} \\ \beta_{\text{rating}, [j]} \end{bmatrix} \sim \text{Multivariate t} \begin{pmatrix} \nu_{[j]}, \text{ } \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \text{ } \mathbf{S}_{[j]} \end{pmatrix} \]
where
\[ \begin{aligned} \mathbf{S}_{[j]} & = \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & \rho\sigma_{\alpha_{[j]}}\sigma_{\beta_{\text{rating}, [j]}} \\ \rho\sigma_{\alpha_{[j]}}\sigma_{\beta_{\text{rating}, [j]}} & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \\ & = \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & \\ & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \mathbf{R}_{[j]} \begin{pmatrix} \sigma_{\alpha_{[j]}}^2 & \\ & \sigma_{\beta_{\text{rating}, [j]}}^2 \end{pmatrix} \end{aligned} \]
and
\[\alpha \sim \text{Student t}(3, \mu_{y}, 2.5)\] \[\beta_{\text{rating}} \sim \text{Student t}(3, 0, 2.5)\] \[\sigma \sim \text{Student t}(3, 0, 2.5)\] \[\nu_{[j]} \sim \text{Gamma}(3.325, 2)\] \[\sigma_{\alpha_{[j]}} \sim \text{Half Student}(3, 0, 2.5)\] \[\sigma_{\beta_{\text{rating}, [j]}} \sim \text{Half Student}(3, 0, 2.5)\] \[\mathbf{R}_{[j]} \sim \text{LKJcorr}(2)\]
notice that \(\mu_{y}\) is the sample mean of response, \(j=\text{subj, roi}\); and \(\text{subj}=1,2,...,N\), \(\text{roi}=1,2,...,M_{n}\). \(N\) is the number of participants; \(M_{n}\) is the number of ROI for \(nth\) participant.