rjagsr2jags

How to make a variable in the function equal to the sum of that variable in previous trials in jags


In my model, I have 10 options (from 1-10) for each subject and each trial to choose (expectancy). I calculated the value for each option based on the rule in the graph below, so the value for each option updated based on the difference between shock and v in every trial (multiply alpha). Then, I used softmax rule to transform v for each option to a certain probability with the same function in this threat: JAGS errors: "Resolving undeclared variables" and "Invalid vector argument to exp".

I guess the problem here is I can't make jags update the value for the same choice.

data: expectancy = number from 1-10 in each trial. shock=number either 1 or 0 in each trial. (I provided example data below)

The second plot is how this be done in stan with 2 choices/1 subject situation.

RW_model <- function(){
  # data
  for(i in 1:nsubjects) # for each person
  { 
    # initial value for v
    v [i,1,expectancy[i,1]] <- 0
    
    for (j in 2:ntrials) # for each trial
    {
      # expectancy chosen
      expectancy[i,j] ~ dcat(mu[i,j,1:10])
      predk[i,j] ~ dcat(mu[i,j,1:10])
      
      # softmax rule to calculate values of each expectancy for each subject
      # tau is the value sensitivity parameter
      mu[i,j,1:10] <- exp_v[i,j,1:10] / sum(exp_v[i,j,1:10])
      exp_v[i,j,expectancy[i,j-1]] <- exp(v[i,j,expectancy[i,j-1]]/tau[i])
      
      # prediction error: difference between feedback and learned values of the chosen expectancy  
      pe [i,j-1] <- shock [i,j-1] - v [i,j-1,expectancy[i,j-1]]
      # value updating process for expectancy
      v [i,j,expectancy[i,j-1]] <- v [i,j-1,expectancy[i,j-1]] + alpha [i] * pe [i,j-1]
    }
  }
  
  # priors
  for (i in 1:nsubjects){
    tau [i] ~ dunif (0,3)
    alpha [i] ~ dunif (0,1)
  }
  
}


# example data/ initial value/ parameters
nsubjects <- 42
ntrials <- 14
shock <- matrix(c(0,0,1,1,0,0,1,1,0,0,1,0,1,0),nrow=42,ncol = 14,byrow = T)
expectancy <- matrix(c(1,2,3,4,5,6,7,7,8,8,7,10,10,00),nrow=42,ncol = 14,byrow = T)



data <- list('shock','nsubjects','ntrials','expectancy')


myinits <-  list(list(tau = runif (42,0,3),
                     alpha = runif (42,0,1)))
parameters <- c("tau",'alpha','v','predk') 

# jags sampling
samples <- jags(data, inits=myinits, parameters,
                    model.file = RW_model,
                    n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T) 

enter image description here

enter image description here


Solution

  • Since there are no stochastic nodes in the setup, I'm not sure what the MCMC simulation will get you, but here is some code that works. From what I can see, the main problem is that when the table says c(Vmax-Vn) it is not using c() as a concatenation operator. c is a constant defined as 0.3.

    dat <- list(
      v = c(0, rep(NA, 8)), 
      vn = rep(NA, 8), 
      vmax=1, 
      c=0.3, 
      Ntrials = 9
    )
    
    jmod <- "model{
    for(i in 2:Ntrials){
      v[i] <- c*(vmax - v[(i-1)]) + v[(i-1)]
    }
    }"
    
    library(runjags)
    
    out <- run.jags(jmod, data=dat, monitor="v")
    

    Edit: Answer to edited question.

    You'll have to see if this is doing what you want, but it produces a result. As far as I can tell, this is a direct extension of the Stan code to multiple subjects and multiple choices per trial ported to JAGS. First, the model:

    RW_model <- function(){
      for(i in 1:Nsubjects){
      for(j in 1:Ntrials){
        expectancy[i,j] ~ dcat(p[i, 1:Nposs, j])
        mu[i,1:Nposs, j] <- tau[i] * v[i,1:Nposs, j]
        for(k in 1:Nposs){
          q[i,k,j] <- exp(mu[i,k,j])
          p[i,k,j] <- q[i,k,j]/sum(q[i,,j])
        }
        
        pe[i,j] <- shock[i,j] - v[i,expectancy[i,j], j]
        for(k in 1:Nposs){
          v[i,k,j+1] <- v[i,k,j] + ifelse(expectancy[i,j] == k, 
                                          alpha[i] * pe[i,j], 
                                          0)  
        }
      }
      }
      for (i in 1:Nsubjects){
        tau [i] ~ dunif (0,3)
        alpha [i] ~ dunif (0,1)
      }
    }
    

    Next, we could make up some data. The data here are for 20 subjects, 5 possible choices in expectancy and 14 trials.

    set.seed(40120)
    # v has to be a Nsubjects x Nchoices x Ntrials+1 matrix
    v <- array(NA, dim=c(20, 5, 15))
    # The first trial of v is initialized to 0 fo all subjects and choices
    v[,,1] <- matrix(0, nrow=20, ncol=5)
    # expectancy and shock are both Nsubjects x Ntrials matrices
    expectancy <- matrix(sample(1:5, 20*14, replace=TRUE), ncol=14)
    shock <- matrix(sample(c(0,1), 20*14, replace=TRUE), ncol=14)
    
    dlist <- list(
      Nsubjects = 20,
      Ntrials = 14,
      Nposs = 5,
      expectancy = expectancy, 
      shock = shock, 
      v=v
    )
    

    Finally, we can specify initial values and run the model.

    myinits <-  list(list(tau = runif (nrow(expectancy),0,3),
                          alpha = runif (nrow(expectancy)0,1)))
    parameters <- c("tau",'alpha') 
    library(R2jags)
    # jags sampling
    samples <- jags(dlist, inits=myinits, parameters,
                    model.file = RW_model,
                    n.chains=1, n.iter=1000, n.burnin=500, n.thin=1, DIC=T) 
    

    You can see from the output that it generated samples and estimates different tau and alpha parameters by subject.

    samples$BUGSoutput
    Inference for Bugs model at "/var/folders/55/q9y1hbcx13b5g50kks_p0mb00000gn/T//Rtmp32zeG4/model1726bdee0c99.txt", fit using jags,
     1 chains, each with 1000 iterations (first 500 discarded)
     n.sims = 500 iterations saved
               mean  sd  2.5%   25%   50%   75% 97.5%
    alpha[1]    0.6 0.3   0.1   0.4   0.7   0.9   1.0
    alpha[2]    0.4 0.3   0.0   0.1   0.3   0.7   0.9
    alpha[3]    0.3 0.2   0.0   0.1   0.2   0.4   0.9
    alpha[4]    0.4 0.3   0.0   0.1   0.3   0.6   0.9
    alpha[5]    0.4 0.3   0.0   0.1   0.3   0.6   1.0
    alpha[6]    0.3 0.3   0.0   0.1   0.2   0.5   0.9
    alpha[7]    0.4 0.3   0.0   0.2   0.4   0.6   0.9
    alpha[8]    0.3 0.2   0.0   0.1   0.2   0.5   0.9
    alpha[9]    0.3 0.3   0.0   0.1   0.3   0.5   0.9
    alpha[10]   0.4 0.2   0.0   0.2   0.4   0.5   0.9
    alpha[11]   0.3 0.3   0.0   0.1   0.2   0.5   0.9
    alpha[12]   0.3 0.3   0.0   0.1   0.2   0.5   0.9
    alpha[13]   0.3 0.3   0.0   0.1   0.3   0.5   0.9
    alpha[14]   0.4 0.3   0.0   0.2   0.4   0.6   0.9
    alpha[15]   0.5 0.3   0.0   0.2   0.5   0.7   1.0
    alpha[16]   0.4 0.3   0.0   0.1   0.3   0.7   0.9
    alpha[17]   0.4 0.3   0.0   0.2   0.4   0.6   0.9
    alpha[18]   0.4 0.3   0.0   0.2   0.4   0.7   1.0
    alpha[19]   0.4 0.3   0.0   0.2   0.4   0.7   1.0
    alpha[20]   0.3 0.3   0.0   0.0   0.2   0.5   0.9
    deviance  909.3 5.1 901.1 905.8 908.8 912.5 919.9
    tau[1]      1.1 0.6   0.2   0.7   1.1   1.5   2.3
    tau[2]      0.8 0.7   0.0   0.3   0.6   1.2   2.7
    tau[3]      1.0 0.8   0.0   0.3   0.8   1.6   2.8
    tau[4]      1.1 0.7   0.0   0.4   0.9   1.6   2.7
    tau[5]      0.9 0.7   0.0   0.3   0.7   1.3   2.7
    tau[6]      1.0 0.8   0.0   0.3   0.9   1.6   2.8
    tau[7]      1.0 0.7   0.1   0.5   0.9   1.5   2.7
    tau[8]      1.1 0.8   0.0   0.4   0.9   1.7   2.9
    tau[9]      1.0 0.8   0.0   0.3   0.8   1.6   2.7
    tau[10]     1.6 0.8   0.1   0.9   1.7   2.3   2.9
    tau[11]     1.1 0.9   0.0   0.3   0.9   1.8   2.8
    tau[12]     1.0 0.8   0.0   0.4   0.8   1.6   2.9
    tau[13]     0.9 0.8   0.0   0.3   0.6   1.4   2.7
    tau[14]     1.1 0.8   0.0   0.5   0.9   1.6   2.7
    tau[15]     1.2 0.7   0.1   0.7   1.1   1.6   2.7
    tau[16]     1.0 0.8   0.1   0.4   0.9   1.5   2.8
    tau[17]     1.1 0.8   0.1   0.4   0.8   1.6   2.9
    tau[18]     1.2 0.8   0.1   0.5   1.1   1.7   2.7
    tau[19]     1.2 0.8   0.1   0.5   1.0   1.8   2.8
    tau[20]     1.0 0.9   0.0   0.3   0.8   1.5   2.9
    
    DIC info (using the rule, pD = var(deviance)/2)
    pD = 12.9 and DIC = 922.1
    DIC is an estimate of expected predictive error (lower deviance is better).
    

    Since I don't know what to expect from this model, you'll have to verify that it is actually doing what you expect, but it seems like this is a good place to start.