Code for hierarchical bayesian estimation of CPT with two-parameter weighting...

Code for hierarchical bayesian estimation of CPT with two-parameter weighting function Prelec (1998)

JAGS/cpt_hierarchical.txt

 model
 {
 
-  # specify hierarchical priors for CPT parameters: parameters are estimated for each synthetic agent j (agent level) and across agents (group level)
-  ## alpha, gamma, delta: link group- and agent-level via probit transformation (cf. Rouder & Lu, 2005)
-  ## rho: link group- and agent-level via log transformation (cf. Nilsson et al., 2011)
-
-  for (j in n_agents)
-  {
-
-  # agent level
-  ## draw transformed agent-level parameters from group-level normal distributions
-
-	alpha.prob[j] ~ dnorm(mu.alpha.prob, tau.alpha.prob) T(-3, 3)
-	gamma.prob[j] ~ dnorm(mu.gamma.prob, tau.gamma.prob) T(-3, 3)
-	delta.prob[j] ~ dnorm(mu.delta.prob, tau.delta.prob) T(-3, 3)
-	rho.log[j]   ~ dnorm(mu.rho.log, tau.rho.log)
-
-	## retransform parameters to original scale and use link function to extend the range (cf. Scheibehenne & Pachur, 2015)
-
-	alpha[j] <- 2*phi(alpha.prob[j])
-	gamma[j] <- 2*phi(gamma.prob[j])
-	delta[j] <- 5*phi(delta.prob[j])
-  rho[j]   <- exp(rho.log[j])
-  }
+  # in the hierarchical model, agent level parameters are drawn from a group level distribution
 
   # group level
-  ## specify group-level priors on transformed parameters such that they are uninformative on the original scale
+
+  ## we specify priors for the group level distributions on transformed scales to facilitate the hierarchical implementation (cf. Rouder & Lu, 2005; Nilsson, 2011)
+  ## prob is probit scale, log is the logarithmic scale
+  ## the priors on the prob and log scale translate to uninformative priors on the original line
+  ## we use uninformative priors for CPT parameters to facilitate comparison of different generative mechanisms
 
   ## alpha
 
@@ -51,47 +34,63 @@ model
   tau.rho.log <- pow(sigma.rho.log, -2)
   sigma.rho.log ~ dunif(0, 1.13)
 
-  ## retransfrom group level means to original scale
+  ## retransfrom group level means and extend scales to a reasonable range (cf. Scheibehenne & Pachur, 2015)
 
   mu.alpha <- 2*phi(mu.alpha.prob)
   mu.gamma <- 2*phi(mu.gamma.prob)
   mu.delta <- 5*phi(mu.delta.prob)
   mu.rho <- exp(mu.rho.log)
 
-  # define CPT model
+  # agent level
 
   for (j in n_agents)
   {
 
+  ## transformed agent-level parameters are drawn from group-level normal distributions on the respective scale
+
+	alpha.prob[j] ~ dnorm(mu.alpha.prob, tau.alpha.prob) T(-3, 3)
+	gamma.prob[j] ~ dnorm(mu.gamma.prob, tau.gamma.prob) T(-3, 3)
+	delta.prob[j] ~ dnorm(mu.delta.prob, tau.delta.prob) T(-3, 3)
+	rho.log[j]   ~ dnorm(mu.rho.log, tau.rho.log)
+
+	## retransform parameters and extend scales to a reasonable range
+
+	alpha[j] <- 2*phi(alpha.prob[j])
+	gamma[j] <- 2*phi(gamma.prob[j])
+	delta[j] <- 5*phi(delta.prob[j])
+  rho[j]   <- exp(rho.log[j])
+
+    # define CPT model
+
     for (i in start:stop)
 	  {
 
     # value function
 
-	  v.a.o1[j, i] <- pow(a_o1[j, i], alpha[j])
-	  v.a.o2[j, i] <- pow(a_o2[j, i], alpha[j])
-	  v.b.o1[j, i] <- pow(b_o1[j, i], alpha[j])
-  	v.b.o2[j, i] <- pow(b_o2[j, i], alpha[j])
+	  v.a.o1[i, j] <- pow(a_o1[i], alpha[j])
+	  v.a.o2[i, j] <- pow(a_o2[i], alpha[j])
+	  v.b.o1[i, j] <- pow(b_o1[i], alpha[j])
+  	v.b.o2[i, j] <- pow(b_o2[i], alpha[j])
 
 	  # weighting function (cf. Prelec, 1998)
 
-	  w.a.p2[j, i] <- (delta[j] * (pow(a_p2_exp[j, i], gamma[j]))) / (delta[j] * (pow(a_p2_exp[j, i], gamma[j])) + pow(a_p1_exp[j, i], gamma[j]))
-	  w.a.p1[j, i] <- 1-w.a.p2[j, i]
-	  w.b.p1[j, i] <- (delta[j] * (pow(b_p1[j, i], gamma[j]))) / (delta[j] * (pow(b_p1[j, i], gamma)) + pow(b_p2[j, i], gamma[j]))
+	  w.a.p2[i, j] <- exp(-delta[j]*(pow((-1*(log(a_p2_exp[i, j]))), gamma[j])))
+	  w.a.p1[i, j] <- 1-w.a.p2[i, j]
+	  w.b.p1[i, j] <- exp(-delta[j]*(pow((-1*(log(b_p1[i, j]))), gamma[j])))
 	  w.b.p2[j, i] <- 1-w.b.p1[j, i]
 
-	  Vf.a[j, i]  <- w.a.p1[j, i] * v.a.o1[j, i] + w.a.p2[j, i] * v.a.o2[j, i]
-	  Vf.b[j, i]  <- w.b.p1[j, i] * v.b.o1[j, i] + w.b.p2[j, i] * v.b.o2[j, i]
+	  Vf.a[i, j]  <- w.a.p1[i, j] * v.a.o1[i, j] + w.a.p2[i, j] * v.a.o2[i, j]
+	  Vf.b[i, j]  <- w.b.p1[i, j] * v.b.o1[i, j] + w.b.p2[i, j] * v.b.o2[i, j]
 
 	  # rescale subjective values to alleviate possible parameter intercorrelations (cf. Krefeld-Schwalb et al., 2021)
 
-	  Vf.a.re[j, i] <- pow(Vf.a[j, i], (1/alpha[j]))
-	  Vf.b.re[j, i] <- pow(Vf.b[i, i], (1/alpha[j]))
+	  Vf.a.re[i, j] <- pow(Vf.a[i, j], (1/alpha[j]))
+	  Vf.b.re[i, j] <- pow(Vf.b[i, j], (1/alpha[j]))
 
 	  # stochastic choice rule
 
-	  binval[j, i] <- (1)/(1+exp((-1*rho[j])*(Vf.a.re[j, i]-Vf.b.re[j, i])))
-	  choice[j, i] ~ dbern(binval[j, i])
+	  binval[i, j] <- (1)/(1+exp((-1*rho[j])*(Vf.a.re[i, j]-Vf.b.re[i, j])))
+	  choice[i, j] ~ dbern(binval[i, j])
 	  }
   }
 }