Tidied up code, added correlation model

Bayes_metad_group.txt

@@ -85,7 +85,7 @@ model {
 	cS2[s,1:nratings-1] <- sort(cS2_raw[s, ])
 
     delta[s] ~ dnorm(0, lambda_delta)
-    logMratio[s] <- mu_Mratio + epsilon_Mratio*delta[s]
+    logMratio[s] <- mu_logMratio + epsilon_logMratio*delta[s]
     Mratio[s] <- exp(logMratio[s])
 	}
 
@@ -101,10 +101,10 @@ model {
 	sigma_c2 ~ dnorm(0, 0.01) I(0, )
 	lambda_c2 <- pow(sigma_c2, -2)
 
-    mu_Mratio ~ dnorm(0, 0.5)
-    sigma_delta ~ dnorm(0, 0.5) I(0,)
+    mu_logMratio ~ dnorm(0, 1)
+    sigma_delta ~ dnorm(0, 1) I(0,)
     lambda_delta <- pow(sigma_delta, -2)
-    epsilon_Mratio ~ dbeta(1,1)
-    sigma_Mratio <- abs(epsilon_Mratio)*sigma_delta
+    epsilon_logMratio ~ dbeta(1,1)
+    sigma_logMratio <- abs(epsilon_logMratio)*sigma_delta
 
 }

Bayes_metad_group_corr.txt

@@ -0,0 +1,122 @@
+# Bayesian estimation of meta-d/d for group correlation between domains
+
+data {
+    for (s in 1:nsubj) {
+        # Type 1 counts for task 1
+        N[s,1] <- sum(counts1[s,1:nratings*2])
+        S[s,1] <- sum(counts1[s,(nratings*2)+1:nratings*4])
+        H[s,1] <- sum(counts1[s,(nratings*3)+1:nratings*4])
+        FA[s,1] <- sum(counts1[s,(nratings*2)+1:nratings*3])
+        M[s,1] <- sum(counts1[s,nratings+1:nratings*2])
+        CR[s,1] <- sum(counts1[s,1:nratings])
+
+        # Type 1 counts for task 2
+        N[s,2] <- sum(counts2[s,1:nratings*2])
+        S[s,2] <- sum(counts2[s,(nratings*2)+1:nratings*4])
+        H[s,2] <- sum(counts2[s,(nratings*3)+1:nratings*4])
+        FA[s,2] <- sum(counts2[s,(nratings*2)+1:nratings*3])
+        M[s,2] <- sum(counts2[s,nratings+1:nratings*2])
+        CR[s,2] <- sum(counts2[s,1:nratings])
+    }
+}
+
+model {
+  for (s in 1:nsubj) {
+
+    ## TYPE 2 SDT MODEL (META-D)
+	# Multinomial likelihood for response counts ordered as c(nR_S1,nR_S2)
+    counts1[s,1:nratings] ~ dmulti(prT[s,1:nratings,1],CR[s,1])
+    counts1[s,nratings+1:nratings*2] ~ dmulti(prT[s,nratings+1:nratings*2,1],M[s,1])
+    counts1[s,(nratings*2)+1:nratings*3] ~ dmulti(prT[s,(nratings*2)+1:nratings*3,1],FA[s,1])
+    counts1[s,(nratings*3)+1:nratings*4] ~ dmulti(prT[s,(nratings*3)+1:nratings*4,1],H[s,1])
+
+    counts2[s,1:nratings] ~ dmulti(prT[s,1:nratings,2],CR[s,2])
+    counts2[s,nratings+1:nratings*2] ~ dmulti(prT[s,nratings+1:nratings*2,2],M[s,2])
+    counts2[s,(nratings*2)+1:nratings*3] ~ dmulti(prT[s,(nratings*2)+1:nratings*3,2],FA[s,2])
+    counts2[s,(nratings*3)+1:nratings*4] ~ dmulti(prT[s,(nratings*3)+1:nratings*4,2],H[s,2])
+
+    for (task in 1:2) {
+
+		# Means of SDT distributions]
+		mu[s,task] <- Mratio[s,task]*d1[s,task]
+		S2mu[s,task] <- mu[s,task]/2
+		S1mu[s,task] <- -mu[s,task]/2
+
+		# Calculate normalisation constants
+		C_area_rS1[s,task] <- phi(c1[s,task] - S1mu[s,task])
+		I_area_rS1[s,task] <- phi(c1[s,task] - S2mu[s,task])
+		C_area_rS2[s,task] <- 1-phi(c1[s,task] - S2mu[s,task])
+		I_area_rS2[s,task] <- 1-phi(c1[s,task] - S1mu[s,task])
+
+		# Get nC_rS1 probs
+		pr[s,1,task] <- phi(cS1[s,1,task] - S1mu[s,task])/C_area_rS1[s,task]  
+		for (k in 1:nratings-2) {                
+			pr[s,k+1,task] <- (phi(cS1[s,k+1,task] - S1mu[s,task])-phi(cS1[s,k,task] - S1mu[s,task]))/C_area_rS1[s,task]
+		}
+		pr[s,nratings,task] <- (phi(c1[s,task] - S1mu[s,task])-phi(cS1[s,nratings-1,task] - S1mu[s,task]))/C_area_rS1[s,task]   
+
+		# Get nI_rS2 probs
+		pr[s,nratings+1,task] <- ((1-phi(c1[s,task] - S1mu[s,task]))-(1-phi(cS2[s,1,task] - S1mu[s,task])))/I_area_rS2[s,task]
+		for (k in 1:nratings-2) {                
+			pr[s,nratings+1+k,task] <- ((1-phi(cS2[s,k,task] - S1mu[s,task]))-(1-phi(cS2[s,k+1,task] - S1mu[s,task])))/I_area_rS2[s,task]
+		}
+		pr[s,nratings*2,task] <- (1-phi(cS2[s,nratings-1,task] - S1mu[s,task]))/I_area_rS2[s,task]
+
+		# Get nI_rS1 probs
+		pr[s,(nratings*2)+1, task] <- phi(cS1[s,1,task] - S2mu[s,task])/I_area_rS1[s,task]
+		for (k in 1:nratings-2) {
+			pr[s,(nratings*2)+1+k,task] <- (phi(cS1[s,k+1,task] - S2mu[s,task])-phi(cS1[s,k,task] - S2mu[s,task]))/I_area_rS1[s,task]
+		}
+		pr[s,nratings*3,task] <- (phi(c1[s,task] - S2mu[s,task])-phi(cS1[s,nratings-1,task] - S2mu[s,task]))/I_area_rS1[s,task]
+
+		# Get nC_rS2 probs
+		pr[s,(nratings*3)+1,task] <- ((1-phi(c1[s,task] - S2mu[s,task]))-(1-phi(cS2[s,1,task] - S2mu[s,task])))/C_area_rS2[s,task]
+		for (k in 1:nratings-2) {                
+			pr[s,(nratings*3)+1+k,task] <- ((1-phi(cS2[s,k,task] - S2mu[s,task]))-(1-phi(cS2[s,k+1,task] - S2mu[s,task])))/C_area_rS2[s,task]
+		}
+		pr[s,nratings*4,task] <- (1-phi(cS2[s,nratings-1,task] - S2mu[s,task]))/C_area_rS2[s,task]
+
+		# Avoid underflow of probabilities
+		for (i in 1:nratings*4) {
+			prT[s,i,task] <- ifelse(pr[s,i,task] < Tol, Tol, pr[s,i,task])
+		}
+
+		# Specify ordered prior on criteria (bounded above and below by Type 1 c) 
+		for (j in 1:(nratings-1)) {
+			cS1_raw[s,j,task] ~ dnorm(-mu_c2[task], lambda_c2[task])
+			cS2_raw[s,j,task] ~ dnorm(mu_c2[task], lambda_c2[task])
+		}
+		cS1[s,1:nratings-1,task] <- sort(cS1_raw[s,1:nratings-1,task])
+		cS2[s,1:nratings-1,task] <- sort(cS2_raw[s,1:nratings-1,task])
+
+		Mratio[s,task] <- exp(logMratio[s,task])
+
+    	} 
+
+    	# Draw log(M)'s from bivariate Gaussian
+    	logMratio[s,1:2] ~ dmnorm(mu_logMratio[], TI[,])
+
+	}
+
+	mu_c2[1] ~ dnorm(0, 0.01)
+	mu_c2[2] ~ dnorm(0, 0.01)
+	sigma_c2[1] ~ dnorm(0, 0.01) I(0, )
+	sigma_c2[2] ~ dnorm(0, 0.01) I(0, )
+	lambda_c2[1] <- pow(sigma_c2[1], -2)
+	lambda_c2[2] <- pow(sigma_c2[2], -2)
+
+    mu_logMratio[1] ~ dnorm(0, 1)
+    mu_logMratio[2] ~ dnorm(0, 1)
+    lambda_logMratio[1] ~ dgamma(0.001,0.001)
+    lambda_logMratio[2] ~ dgamma(0.001,0.001)
+    sigma_logMratio[1] <- 1/sqrt(lambda_logMratio[1])
+    sigma_logMratio[2] <- 1/sqrt(lambda_logMratio[2])
+    rho ~ dunif(-1,1)
+
+    T[1,1] <- 1/lambda_logMratio[1]
+    T[1,2] <- rho*sigma_logMratio[1]*sigma_logMratio[2]
+    T[2,2] <- 1/lambda_logMratio[2]
+    T[2,1] <- rho*sigma_logMratio[1]*sigma_logMratio[2]
+    TI[1:2,1:2] <- inverse(T[1:2,1:2])
+
+}

Bayes_metad_group_nodp.txt

@@ -75,7 +75,7 @@ model {
 	cS2[s,1:nratings-1] <- sort(cS2_raw[s, ])
 
     delta[s] ~ dnorm(0, lambda_delta)
-    logMratio[s] <- mu_Mratio + epsilon_Mratio*delta[s]
+    logMratio[s] <- mu_logMratio + epsilon_logMratio*delta[s]
     Mratio[s] <- exp(logMratio[s])
 	}
 
@@ -84,10 +84,10 @@ model {
 	sigma_c2 ~ dnorm(0, 0.01) I(0, )
 	lambda_c2 <- pow(sigma_c2, -2)
 
-    mu_Mratio ~ dnorm(0, 0.5)
-    sigma_delta ~ dnorm(0, 0.5) I(0,)
+    mu_logMratio ~ dnorm(0, 1)
+    sigma_delta ~ dnorm(0, 1) I(0,)
     lambda_delta <- pow(sigma_delta, -2)
-    epsilon_Mratio ~ dbeta(1,1)
-    sigma_Mratio <- abs(epsilon_Mratio)*sigma_delta
+    epsilon_logMratio ~ dbeta(1,1)
+    sigma_logMratio <- abs(epsilon_logMratio)*sigma_delta
 
 }

Bayes_metad_rc_group.txt

@@ -88,10 +88,10 @@ model {
 	cS2[s,1:nratings-1] <- sort(cS2_raw[s, ])
 
     delta_rS1[s] ~ dnorm(0, lambda_delta_rS1)
-    logMratio_rS1[s] <- mu_Mratio_rS1 + epsilon_Mratio_rS1*delta_rS1[s]
+    logMratio_rS1[s] <- mu_logMratio_rS1 + epsilon_logMratio_rS1*delta_rS1[s]
     Mratio_rS1[s] <- exp(logMratio_rS1[s])
     delta_rS2[s] ~ dnorm(0, lambda_delta_rS2)
-    logMratio_rS2[s] <- mu_Mratio_rS2 + epsilon_Mratio_rS2*delta_rS2[s]
+    logMratio_rS2[s] <- mu_logMratio_rS2 + epsilon_logMratio_rS2*delta_rS2[s]
     Mratio_rS2[s] <- exp(logMratio_rS2[s])
 
 	}
@@ -108,16 +108,16 @@ model {
 	sigma_c2 ~ dnorm(0, 0.01) I(0, )
 	lambda_c2 <- pow(sigma_c2, -2)
 
-    mu_Mratio_rS1 ~ dnorm(0, 0.5)
+    mu_logMratio_rS1 ~ dnorm(0, 0.5)
     sigma_delta_rS1 ~ dnorm(0, 1) I(0,)
     lambda_delta_rS1 <- pow(sigma_delta_rS1, -2)
-    epsilon_Mratio_rS1 ~ dbeta(1,1)
-    sigma_Mratio_rS1 <- abs(epsilon_Mratio_rS1)*sigma_delta_rS1
+    epsilon_logMratio_rS1 ~ dbeta(1,1)
+    sigma_logMratio_rS1 <- abs(epsilon_logMratio_rS1)*sigma_delta_rS1
 
-    mu_Mratio_rS2 ~ dnorm(0, 0.5)
+    mu_logMratio_rS2 ~ dnorm(0, 0.5)
     sigma_delta_rS2 ~ dnorm(0, 1) I(0,)
     lambda_delta_rS2 <- pow(sigma_delta_rS2, -2)
-    epsilon_Mratio_rS2 ~ dbeta(1,1)
-    sigma_Mratio_rS2 <- abs(epsilon_Mratio_rS2)*sigma_delta_rS2
+    epsilon_logMratio_rS2 ~ dbeta(1,1)
+    sigma_logMratio_rS2 <- abs(epsilon_logMratio_rS2)*sigma_delta_rS2
 
 }

Bayes_metad_rc_group_nodp.txt

@@ -78,10 +78,10 @@ model {
 	cS2[s,1:nratings-1] <- sort(cS2_raw[s, ])
 
     delta_rS1[s] ~ dnorm(0, lambda_delta_rS1)
-    logMratio_rS1[s] <- mu_Mratio_rS1 + epsilon_Mratio_rS1*delta_rS1[s]
+    logMratio_rS1[s] <- mu_logMratio_rS1 + epsilon_logMratio_rS1*delta_rS1[s]
     Mratio_rS1[s] <- exp(logMratio_rS1[s])
     delta_rS2[s] ~ dnorm(0, lambda_delta_rS2)
-    logMratio_rS2[s] <- mu_Mratio_rS2 + epsilon_Mratio_rS2*delta_rS2[s]
+    logMratio_rS2[s] <- mu_logMratio_rS2 + epsilon_logMratio_rS2*delta_rS2[s]
     Mratio_rS2[s] <- exp(logMratio_rS2[s])
 
 	}
@@ -91,16 +91,16 @@ model {
 	sigma_c2 ~ dnorm(0, 0.01) I(0, )
 	lambda_c2 <- pow(sigma_c2, -2)
 
-    mu_Mratio_rS1 ~ dnorm(0, 0.5)
+    mu_logMratio_rS1 ~ dnorm(0, 1)
     sigma_delta_rS1 ~ dnorm(0, 1) I(0,)
     lambda_delta_rS1 <- pow(sigma_delta_rS1, -2)
-    epsilon_Mratio_rS1 ~ dbeta(1,1)
-    sigma_Mratio_rS1 <- abs(epsilon_Mratio_rS1)*sigma_delta_rS1
+    epsilon_logMratio_rS1 ~ dbeta(1,1)
+    sigma_logMratio_rS1 <- abs(epsilon_logMratio_rS1)*sigma_delta_rS1
 
-    mu_Mratio_rS2 ~ dnorm(0, 0.5)
+    mu_logMratio_rS2 ~ dnorm(0, 1)
     sigma_delta_rS2 ~ dnorm(0, 1) I(0,)
     lambda_delta_rS2 <- pow(sigma_delta_rS2, -2)
-    epsilon_Mratio_rS2 ~ dbeta(1,1)
-    sigma_Mratio_rS2 <- abs(epsilon_Mratio_rS2)*sigma_delta_rS2
+    epsilon_logMratio_rS2 ~ dbeta(1,1)
+    sigma_logMratio_rS2 <- abs(epsilon_logMratio_rS2)*sigma_delta_rS2
 
 }

exampleFit.m

@@ -20,5 +20,5 @@
 hdi = calc_HDI(fit.mcmc.samples.meta_d(:));
 fprintf(['\n HDI: ', num2str(hdi) '\n\n'])
 
-% Visualise fits
+% Visualise single-subject fits
 metad_visualise

exampleFit_group.m

@@ -35,6 +35,6 @@
 fit = fit_meta_d_mcmc_group(nR_S1, nR_S2, mcmc_params);
 
 % Call plotSamples to plot posterior of group Mratio
-plotSamples(exp(fit.mcmc.samples.mu_Mratio))
-hdi = calc_HDI(exp(fit.mcmc.samples.mu_Mratio(:)));
-fprintf(['\n HDI: ', num2str(hdi) '\n\n'])
+plotSamples(exp(fit.mcmc.samples.mu_logMratio))
+hdi = calc_HDI(exp(fit.mcmc.samples.mu_logMratio(:)));
+fprintf(['\n HDI on meta-d/d: ', num2str(hdi) '\n\n'])

exampleFit_group_rc.m

@@ -23,10 +23,6 @@
     % Generate data
     sim = type2_SDT_sim(d(i), noise, c, c1, c2, Ntrials);
 
-    % Avoid zero counts
-    sim.nR_S1(sim.nR_S1 == 0) = 1;
-    sim.nR_S2(sim.nR_S2 == 0) = 1;
-
     nR_S1{i} = sim.nR_S1;
     nR_S2{i} = sim.nR_S2;
 
@@ -41,5 +37,5 @@
 fit = fit_meta_d_mcmc_group(nR_S1, nR_S2, mcmc_params);
 
 % Plot output
-plotSamples(fit.mcmc.samples.mu_Mratio_rS1)
-plotSamples(fit.mcmc.samples.mu_Mratio_rS2)
+plotSamples(exp(fit.mcmc.samples.mu_logMratio_rS1))
+plotSamples(exp(fit.mcmc.samples.mu_logMratio_rS2))

exampleFit_twoGroups.m

@@ -44,11 +44,11 @@
 fit2 = fit_meta_d_mcmc_group(DATA(2).nR_S1, DATA(2).nR_S2);
 
 % Compute HDI of difference
-sampleDiff = fit1.mcmc.samples.mu_Mratio - fit2.mcmc.samples.mu_Mratio;
+sampleDiff = fit1.mcmc.samples.mu_logMratio - fit2.mcmc.samples.mu_logMratio;
 hdi = calc_HDI(sampleDiff(:));
 fprintf(['\n HDI on difference in log(meta-d''/d''): ', num2str(hdi) '\n\n'])
 
 % Plot group Mratio and the difference
-plotSamples(exp(fit1.mcmc.samples.mu_Mratio))
-plotSamples(exp(fit2.mcmc.samples.mu_Mratio))
+plotSamples(exp(fit1.mcmc.samples.mu_logMratio))
+plotSamples(exp(fit2.mcmc.samples.mu_logMratio))
 plotSamples(sampleDiff)

fit_meta_d_mcmc.m

@@ -3,7 +3,7 @@
 %
 % Given data from an experiment where an observer discriminates between two
 % stimulus alternatives on every trial and provides confidence ratings,
-% fits both d' and meta-d' using MCMC implemented in
+% fits meta-d' using MCMC implemented in
 % JAGS. Requires JAGS to be installed
 % (see http://psiexp.ss.uci.edu/research/programs_data/jags/)
 %
@@ -113,11 +113,7 @@
 % Parts of this code are adapted from Brian Maniscalco's meta-d' toolbox
 % which can be found at http://www.columbia.edu/~bsm2105/type2sdt/
 %
-% Updated 12/10/15 to include estimation of type 1 d' within same model
-
-% toy data
-% nR_S1 = [1552  933  954  720  448  220   78   27];
-% nR_S2 = [33   77  213  469  729 1013  975 1559];
+% Updated 12/10/15 to allow estimation of type 1 d' within same model
 
 cwd = pwd;
 
@@ -175,7 +171,7 @@
     mcmc_params.response_conditional = 0;
     mcmc_params.estimate_dprime = 1;    % also estimate dprime in same model?
     mcmc_params.nchains = 3; % How Many Chains?
-    mcmc_params.nburnin = 3000; % How Many Burn-in Samples?
+    mcmc_params.nburnin = 1000; % How Many Burn-in Samples?
     mcmc_params.nsamples = 10000;  %How Many Recorded Samples?
     mcmc_params.nthin = 1; % How Often is a Sample Recorded?
     mcmc_params.doparallel = 0; % Parallel Option