Hi,
I’m rather new to working with stan (as you will see) and trying to run the toy covariance model from Palestro et al. 2018 (http://compmem.org/assets/pubs/Palestro.etal.2018a.pdf). I consulted the Stan documentation and the forums, but I fail to make the model run properly.
The following data are generated in R:
# need both logit and logit^{-1} functions
logit=function(x)log(x/(1-x))
invlogit=function(x){1/(1+exp(-x))}
# set up model specification
n <- 500 # total number of trials
# establish the hyperparameters
sig1 <- .5 # std. dev. of single-trial BOLD responses
sig2 <- 1 # std. dev. of item memory strength (logit scale)
rho <- .6 # cor b/n brain activation and memory strength
# set up hyper variance-covariance matrix Sigma
sigma <- matrix(c(sig1^2, # element [1,1]
sig1*sig2*rho, # element [1,2]
sig1*sig2*rho, # element [2,1]
sig2^2 # element [2,2]
),2,2,byrow=TRUE)
# set up hyper mean vector phi
phi <- c(2,0)
# simulate single-trial delta and theta matrix
DeltaTheta <- rmvnorm(n,phi,sigma)
# generate observed variable nodes
ts <- seq(0,4,1) # scan times
sig <- .5 # the std. dev. of BOLD responses
# declare some storage objects
N=matrix(NA,n,length(ts))
B=numeric(n)
# loop over trials
for(i in 1:n){
# N is a normal deviate with mean controlled by delta
N[i,]=rnorm(length(ts),DeltaTheta[i,1]*ts,sig)
# B is a Bernoulli deviate with probability controlled by theta
B[i]=rbinom(1,1,invlogit(DeltaTheta[i,2]))
}
NB = matrix(c(N,B),nrow = n,ncol = 6)
data = list(n=n,
nPar=2,
B=B,
N=N,
ts=ts,
Nt=length(ts),
sig=sig,
I0=diag(2),
n0=2,
phi0=rep(0,2),
s0=diag(2))
Then they fit the model in JAGS using
model {
# convert sig to tau for convenience
tau <- pow(sig, -2)
# loop through trials to define likelihood
for (i in 1:n){
for (t in 1:Nt){
# likelihood for neural data
N[i,t] ~ dnorm(DeltaTheta[i,1]*ts[t],tau);
}
# likelihood for behavioral data
B[i] ~ dbin(1/(1+exp(-DeltaTheta[i,2])),1);
}
# loop through trials to define prior on (delta, theta)
for(i in 1:n){
DeltaTheta[i,1:2] ~ dmnorm(phi,Omega);
}
# priors on hyperparameters
phi ~ dmnorm(phi0,s0);
Omega ~ dwish(I0, n0);
# convert Omega to Sigma for convenience
Sigma <- inverse(Omega);
}
I tried to translate the model to Stan and also estimate the correlation not the covariance matrix.
// joint neural and behavioural model
data{
int<lower=1> n; //number of trials
int B[n]; // Behavioral data matrix
int<lower=1> Nt; // number of time points
vector[Nt] ts; //time points
real N[n,Nt]; // Neural data matrix
real<lower=0,upper=1> sig; // sd of neural data
int<lower=1> nPar; // number of parameters in the joint model
vector[nPar] phi0; // joint mean vector
matrix[nPar, nPar] s0; // joint sigma
}
parameters{
row_vector[nPar] phi; // joint mean vector
corr_matrix[2] omega; // Correlation matrix of N and B data
vector<lower=0>[2] tau; // scale of the correlation matrix
matrix[n, nPar] raw_DeltaTheta;
}
transformed parameters{
matrix[n, nPar] DeltaTheta; //neural and behavioral parameters
DeltaTheta= raw_DeltaTheta*quad_form_diag(omega,tau);
DeltaTheta[,1] = DeltaTheta[,1]+phi[1];
DeltaTheta[,2] = DeltaTheta[,2]+phi[2];
}
model{
// prior for scale of the hyperparameter correlation matrix sigma
tau ~ normal(0,2.5);//cauchy(0,2.5);
// Joint model priors
phi ~ multi_normal(phi0,s0);
// prior on the correlation matrix
omega ~ lkj_corr(6);
// for the non-centered reparametrization
to_vector(raw_DeltaTheta) ~ std_normal();
// likelihood
for(i in 1:n){
for(t in 1:Nt){
N[i,t]~normal(DeltaTheta[i,1]*ts[t],sig);
}
B[i]~binomial(1,inv_logit(DeltaTheta[i,2]));
}
}
This results in the following error messages (4 chains with 6000 iterations)
1: In system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
'-E' not found
2: There were 3021 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
3: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
4: Examine the pairs() plot to diagnose sampling problems
5: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
6: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
omega[1,1] always has an Rhat = NaN.
It’s a fairly simple model, but I can’t make it run properly. Could you point me in the right direction? Thanks!
