Organize output of MCMC loop

JAGS/cpt_trial_level.txt

@@ -18,32 +18,22 @@ model
   for (i in start:stop)
 	{
 
-	# value function
+	v.a.o1[i] <- pow(a_o1[i],alpha)
+	v.a.o2[i] <- pow(a_o2[i],alpha)
+	v.b.o1[i] <- pow(b_o1[i],alpha)
+	v.b.o2[i] <- pow(b_o2[i],alpha)
 
-	v.a.o1[i] <- pow(a_o1[i], alpha)
-	v.a.o2[i] <- pow(a_o2[i], alpha)
-	v.b.o1[i] <- pow(b_o1[i], alpha)
-	v.b.o2[i] <- pow(b_o2[i], alpha)
-
-	# weighting function
-
-	w.a.p2[i] <- (delta * (pow(a_p2_exp[i], gamma))) / (delta * (pow(a_p2_exp[i], gamma)) + pow(a_p1_exp[i], gamma))
+	w.a.p2[i] <- (delta * (pow(a_p2_exp[i],gamma))) / (delta * (pow(a_p2_exp[i],gamma)) + pow(a_p1_exp[i],gamma))
 	w.a.p1[i] <- 1-w.a.p2[i]
-	w.b.p1[i] <- (delta * (pow(b_p1[i], gamma))) / (delta * (pow(b_p1[i], gamma)) + pow(b_p2[i], gamma))
+	w.b.p1[i] <- (delta * (pow(b_p1[i],gamma))) / (delta * (pow(b_p1[i],gamma)) + pow(b_p2[i],gamma))
 	w.b.p2[i] <- 1-w.b.p1[i]
 
-	# valuation
-
 	Vf.a[i]  <- w.a.p1[i] * v.a.o1[i] + w.a.p2[i] * v.a.o2[i]
 	Vf.b[i]  <- w.b.p1[i] * v.b.o1[i] + w.b.p2[i] * v.b.o2[i]
 
-	## retransform
-
 	Vf.a.re[i] = pow(Vf.a[i], (1/alpha))
 	Vf.b.re[i] = pow(Vf.b[i], (1/alpha))
 
-	# choice rule
-
 	binval[i] <- (1)/(1+exp((-1*rho)*(Vf.a.re[i]-Vf.b.re[i])))
 	choice[i] ~ dbern(binval[i])
 	}

R/cpt.R

@@ -17,8 +17,7 @@ choices <- read_csv("data/choices/choices.csv", col_types = cols)
 
 choices_cpt <- choices %>%
   filter(!(is.na(a_ev_exp) | is.na(b_ev_exp))) %>% # remove choices where prospects were not attended
-  mutate(choice_A = if_else(choice == "A", 1, 0), # for logit choice rule
-         i = row_number()) # assign trial numbers
+  mutate(choice_A = if_else(choice == "A", 1, 0)) # to apply logit choice rule
 
 ## prepare computation of CPT estimates for each distinct strategy-parameter combination
 
@@ -26,7 +25,8 @@ params_sim <- choices_cpt %>% distinct(strategy, s, boundary, a) # group trials
 choices_grouped <- vector("list", nrow(params_sim))
 for(set in seq_len(nrow(params_sim))){
   choices_grouped[[set]] <- choices_cpt %>%
-    filter(strategy == theta[[set, "strategy"]] & s == theta[[set, "s"]] & boundary == theta[[set, "boundary"]] & a == theta[[set, "a"]])
+    filter(strategy == params_sim[[set, "strategy"]] & s == params_sim[[set, "s"]] & boundary == params_sim[[set, "boundary"]] & a == params_sim[[set, "a"]]) %>%
+    mutate(i = row_number()) # assign trial numbers
 }
 estimates_cpt <- vector("list", nrow(params_sim)) # allocate space for JAGS output
 
@@ -68,8 +68,9 @@ for(set in seq_len(nrow(params_sim))){
                                 DIC = FALSE,
                                 jags.seed = 8362)
 
-  # to get posterior estimates, credibility intervals, and MCMC diagnostics
-
+  # get posterior estimates, credibility intervals, and MCMC diagnostics
 
+  summary <- MCMC_samples$BUGSoutput$summary %>%
+    as_tibble(rownames = "parameter") %>%
+  estimates_cpt[[set]] <- expand_grid(set, summary)
 }
-