Hi, the current model seems to work fine if I wanted psychometric thresholds for one time point. However, given that there have been some issues with methods (there were issues with adaptive algorithm), it now seems best to model all thresholds for all timepoints (4) in one model. There are ~ 1000 participants with 14 conditions, and I want to add on an interaction with time point as well. However, I’ve tried a couple different ways of adding the interaction of time point, and I never even get the model to being compiling, so I must be doing something wrong.
This is the code for the first model (without the interaction of time point):
data{
int<lower=4> n_total; // Total number of trials to analyze
int<lower=2> n_subjects; // Number of unique observers
int<lower=2> n_levels; // Number of levels of Factor
real intensity[n_total]; // Intensity of the stimulus on each trial
int<lower=0> subject[n_total]; // Observer on each trial
int<lower=0> level[n_total]; // Level of Factor on each trial
int<lower=0,upper=1> correct[n_total]; // Whether the response was correct (1) on each trial
real<lower=0,upper=1> chance_performance; // Chance performance for experiment (e.g., 1/n for n-AFC)
}
transformed data{
int df_levels_sA; // Number of degrees of freedom in the interaction
int n_levels_sA; // Number of levels in the interaction
int df_levels_sB; // Number of degrees of freedom in the interaction
int n_levels_sB;
df_levels_sA = (n_levels - 1) * (n_subjects - 1);
n_levels_sA = n_levels * n_subjects;
df_levels_sB = (n_levels - 1) * (n_subjects - 1);
n_levels_sB = n_levels * n_subjects;
}
parameters{
vector<lower=0,upper=1>[n_subjects] lapse; // Observer's lapse rate
real mum;
real muw;
vector[n_levels-1] fA;
vector[n_levels-1] fB;
vector[n_subjects-1] sA;
vector[n_subjects-1] sB;
//matrix[n_subjects-1,n_levels-1] fsA;
vector[(n_subjects - 1) * (n_levels-1)] fsA;
vector[(n_subjects - 1) * (n_levels-1)] fsB;
}
transformed parameters {
real<lower=0,upper=1> mean_beta; // Mean of beta prior for individual lapse rate
real<lower=0> betaEta; // Precision parameter of beta prior for individual lapse rate
vector<lower=0,upper=1>[n_total] psi;
real<lower=0> lapse_alpha;
real<lower=0> lapse_beta;
real m[n_subjects,n_levels];
real w[n_subjects,n_levels];
real threshold[n_total];
real width[n_total];
vector[n_levels] factor_alpha = append_row(fA, -sum(fA));
vector[n_levels] factor_beta = append_row(fB, sum(fB));
vector[n_subjects] subject_alpha = append_row(sA, -sum(sA));
vector[n_subjects] subject_beta = append_row(sB, sum(sB));
matrix[n_subjects,n_levels] interaction_alpha;
matrix[n_subjects,n_levels] interaction_beta;
interaction_alpha = rep_matrix(0, n_subjects, n_levels);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_alpha[sj, l] = fsA[(sj - 1) * (n_levels - 1) + l];
}
interaction_alpha[sj, n_levels] = -1 * sum(interaction_alpha[sj,]);
}
for (l in 1:n_levels) {
interaction_alpha[n_subjects,l] = -1 * sum(interaction_alpha[,l]);
}
interaction_beta = rep_matrix(0, n_subjects, n_levels);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_beta[sj, l] = fsB[(sj - 1) * (n_levels - 1) + l];
}
interaction_beta[sj, n_levels] = sum(interaction_beta[sj,]);
}
for (l in 1:n_levels) {
interaction_beta[n_subjects,l] = sum(interaction_beta[,l]);
}
for (sj in 1:n_subjects) {
for (l in 1:n_levels) {
m[sj,l] = mum + factor_alpha[l] + subject_alpha[sj] + interaction_alpha[sj, l];
w[sj,l] = muw + factor_beta[l] + subject_beta[sj] + interaction_beta[sj, l];
}
}
for (tr in 1:n_total) {
threshold[tr] = m[subject[tr],level[tr]];
width[tr] = w[subject[tr],level[tr]];
psi[tr] = chance_performance + (1-lapse[ subject[tr]] - chance_performance) * inv_logit((intensity[tr]-threshold[tr])/width[tr]);
//(1-lapse[ subject[tr] ])*( (1-chance_performance)*inv_logit((4.4/width[tr]) * (intensity[tr]-threshold[tr])) + chance_performance) + chance_performance*lapse[ subject[tr] ];
if (is_nan(psi[tr])) {
print("lapse = ", lapse[subject[tr]]);
print("width = ", width[tr]);
print("threshold = ", threshold[tr]);
print("intensity = ", intensity[tr])
}
}
mean_beta = 0.01;
betaEta = 100;
lapse_alpha = mean_beta * betaEta;
lapse_beta = (1-mean_beta) * betaEta ;
}
model {
//mean_beta ~ beta(1,60);
//betaEta ~ gamma(1,0.01);
fsA ~ normal(0, inv(sqrt(1 - inv(n_levels_sA)))); // Produces a standard normal on interaction_alpha
fA ~ normal(0, inv(sqrt(1 - inv(n_levels))));
sA ~ normal(0, inv(sqrt(1 - inv(n_subjects))));
sB ~ normal(0, inv(sqrt(1 - inv(n_subjects))));
fB ~ normal(0, inv(sqrt(1 - inv(n_levels))));
fsB ~ normal(0, inv(sqrt(1 - inv(n_levels))));
muw ~ gamma(2,.5);
lapse ~ beta(lapse_alpha,lapse_beta);
correct ~ bernoulli(psi);
}
This is what I have for the interaction with timepoint:
data{
int<lower=4> n_total; // Total number of trials to analyze
int<lower=2> n_subjects; // Number of unique observers
int<lower=2> n_levels; // Number of levels of Factor
int<lower=2> n_timepoints;
real intensity[n_total]; // Intensity of the stimulus on each trial
int<lower=0> subject[n_total]; // Observer on each trial
int<lower=0> level[n_total]; // Level of Factor on each trial
int<lower=0> timepoint[n_total];
int<lower=0,upper=1> correct[n_total]; // Whether the response was correct (1) on each trial
real<lower=0,upper=1> chance_performance; // Chance performance for experiment (e.g., 1/n for n-AFC)
}
transformed data{
int df_levels_sA; // Number of degrees of freedom in the interaction
int n_levels_sA; // Number of levels in the interaction
int n_timepoints_sA;
int df_levels_sB; // Number of degrees of freedom in the interaction
int n_levels_sB;
df_levels_sA = (n_levels - 1) * (n_subjects - 1) * (n_timepoints -1);
n_levels_sA = n_levels * n_subjects * n_timepoints;
df_levels_sB = (n_levels - 1) * (n_subjects - 1);
n_levels_sB = n_levels * n_subjects;
}
parameters{
vector<lower=0,upper=1>[n_subjects] lapse; // Observer's lapse rate
real mum;
real muw;
vector[n_levels-1] fA;
vector[n_levels-1] fB;
vector[n_subjects-1] sA;
vector[n_subjects-1] sB;
vector[n_timepoints-1] tA;
//matrix[n_subjects-1,n_levels-1] fsA;
vector[(n_subjects - 1) * (n_levels-1)] fsA;
vector[(n_subjects - 1) * (n_timepoints-1)] fsA1;
vector[(n_levels - 1) * (n_timepoints-1)] fsA2;
vector[(n_subjects - 1) * (n_levels-1)] fsB;
}
transformed parameters {
real<lower=0,upper=1> mean_beta; // Mean of beta prior for individual lapse rate
real<lower=0> betaEta; // Precision parameter of beta prior for individual lapse rate
vector<lower=0,upper=1>[n_total] psi;
real<lower=0> lapse_alpha;
real<lower=0> lapse_beta;
real m[n_subjects,n_levels, n_timepoints];
real w[n_subjects,n_levels];
real threshold[n_total];
real width[n_total];
vector[n_levels] factor_alpha = append_row(fA, -sum(fA));
vector[n_levels] factor_beta = append_row(fB, sum(fB));
vector[n_subjects] subject_alpha = append_row(sA, -sum(sA));
vector[n_subjects] subject_beta = append_row(sB, sum(sB));
vector[n_timepoints] timepoint_alpha = append_row(tA, -sum(tA));
matrix[n_subjects, n_levels] interaction_alpha;
matrix[n_subjects, n_timepoints] interaction_alpha1;
matrix[n_levels, n_timepoints] interaction_alpha2;
matrix[n_subjects, n_levels] interaction_beta;
interaction_alpha = rep_matrix(0, n_subjects, n_levels);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_alpha[sj, l] = fsA[(sj - 1) * (n_levels - 1) + l];
}
interaction_alpha[sj, n_levels] = -1 * sum(interaction_alpha[sj,]);
}
for (l in 1:n_levels) {
interaction_alpha[n_subjects,l] = -1 * sum(interaction_alpha[,l]);
}
}
interaction_alpha1 = rep_matrix(0, n_subjects, n_timepoints);
for(sj in 1:(n_subjects - 1)) {
for(y in 1:(n_timepoints - 1)) {
interaction_alpha1[sj, y] = fsA1[(sj - 1) * (n_timepoints - 1) + y];
}
interaction_alpha1[sj, n_timepoints] = -1 * sum(interaction_alpha1[sj,]);
}
for (y in 1:n_timepoints) {
interaction_alpha1[n_subjects,y] = -1 * sum(interaction_alpha1[,y]);
}
}
interaction_alpha2 = rep_matrix(0, n_levels, n_timepoints);
for(l in 1:(n_levels - 1)) {
for(y in 1:(n_timepoints - 1)) {
interaction_alpha2[l, y] = fsA2[(l - 1) * (n_timepoints - 1) + y];
}
interaction_alpha2[l, n_timepoints] = -1 * sum(interaction_alpha2[l,]);
}
for (l in 1:n_timepoints) {
interaction_alpha2[n_levels,y] = -1 * sum(interaction_alpha2[,y]);
}
}
interaction_beta = rep_matrix(0, n_subjects, n_levels);
for(sj in 1:(n_subjects - 1)) {
for(l in 1:(n_levels - 1)) {
interaction_beta[sj, l] = fsB[(sj - 1) * (n_levels - 1) + l];
}
interaction_beta[sj, n_levels] = sum(interaction_beta[sj,]);
}
for (l in 1:n_levels) {
interaction_beta[n_subjects,l] = sum(interaction_beta[,l]);
}
for (sj in 1:n_subjects) {
for (l in 1:n_levels) {
m[sj,l] = mum + factor_alpha[l] + subject_alpha[sj] + interaction_alpha[sj, l];
w[sj,l] = muw + factor_beta[l] + subject_beta[sj] + interaction_beta[sj, l];
}
}
for (tr in 1:n_total) {
threshold[tr] = m[subject[tr],level[tr]];
width[tr] = w[subject[tr],level[tr]];
psi[tr] = chance_performance + (1-lapse[ subject[tr]] - chance_performance) * inv_logit((intensity[tr]-threshold[tr])/width[tr]);
//(1-lapse[ subject[tr] ])*( (1-chance_performance)*inv_logit((4.4/width[tr]) * (intensity[tr]-threshold[tr])) + chance_performance) + chance_performance*lapse[ subject[tr] ];
if (is_nan(psi[tr])) {
print("lapse = ", lapse[subject[tr]]);
print("width = ", width[tr]);
print("threshold = ", threshold[tr]);
print("intensity = ", intensity[tr])
}
}
mean_beta = 0.01;
betaEta = 100;
lapse_alpha = mean_beta * betaEta;
lapse_beta = (1-mean_beta) * betaEta ;
}
model {
//mean_beta ~ beta(1,60);
//betaEta ~ gamma(1,0.01);
fsA ~ normal(0, inv(sqrt(1 - inv(n_levels_sA)))); // Produces a standard normal on interaction_alpha
fA ~ normal(0, inv(sqrt(1 - inv(n_levels))));
sA ~ normal(0, inv(sqrt(1 - inv(n_subjects))));
tA ~ normal(0, inv(sqrt(1 - inv(n_timepoints))));
sB ~ normal(0, inv(sqrt(1 - inv(n_subjects))));
fB ~ normal(0, inv(sqrt(1 - inv(n_levels))));
fsB ~ normal(0, inv(sqrt(1 - inv(n_levels))));
muw ~ gamma(2,.5);
lapse ~ beta(lapse_alpha,lapse_beta);
correct ~ bernoulli(psi);
}