Brightness perception:
Metaregression model example

Introduction

Data from Ratcliff and Rouder (1998), experiment 1, participant N.H.:

’‘…subjects were asked to decide whether the overall brightness of pixel arrays displayed on a computer monitor was “high” or “low” (…) The brightness of a display was controlled by the proportion of the pixels that were white…’’

The experimental design considered two main factors:

Instructions: Accuracy vs Speed.
Stimulus configuration: Majority black vs Majority white pixels.

For a total of 66 cells in the design.

Cells 1-16 had more black pixels and instructions favored accurate responses.
Cells 18-33 had more white pixels and instructions favored accurate responses.
Cells 34-49 had more black pixels and instructions favored speedy responses.
Cells 51-66 had more white pixels and instructions favored speedy responses.

Cells 17 and 50 had 50/50 black-and-white pixels. These cells are discarded because they can’t provide accuracy measures.

# Load necessary libraries/packages
library(R2jags)
library(here)
seed <- 15

Loading and cleaning the data

Load the data

# Load the data from one participant
data_raw <- read.csv(here("demos", "applications", "brightness_perception", "data", "nh.tsv"), sep = "")
colnames(data_raw) <- c("index","cond","response","RT")
head(data_raw)

##   index cond response  RT
## 1     1    1        1 457
## 2     1    1        1 569
## 3     1    1        1 390
## 4     1    1        1 534
## 5     1    1        1 501
## 6     1    1        1 488

# No. of observations
nrow(data_raw)

## [1] 7889

Clean the data

# Create a copy of the raw data file 
data <- data_raw

# Get 'accuracy' binary coding based on the condition and response in raw data file
accuracy <- as.integer(data_raw$cond > 0  & data_raw$cond < 17 & data_raw$response==1|
                       data_raw$cond > 17 & data_raw$cond < 34 & data_raw$response==2|
                       data_raw$cond > 33 & data_raw$cond < 50 & data_raw$response==1|
                       data_raw$cond > 50 & data_raw$cond < 67 & data_raw$response==2)

# Update the 'response' column of the data copied so it reflects accuracy
data$response <- accuracy

# Remove rows where RT > 3000ms
data <- data[which(data$RT<=3000),]

Get summary statistics

Note: Since the three-parameters considered by the EZDDM formulation don’t predict any difference in shape between the correct and incorrect RT distributions, we use all trials to compute the mean and variance of the RT.

Write custom function `ez_summaries`

# Define a function to compute the summary statistics used by EZ-DDM
ez_summaries <- function(subset){
  # Identify condition ID
  cond <- unique(subset$cond)
  # Return relevant summary statistics
  return(data.frame("nTrials" = nrow(subset),
                    "score" = sum(subset$response),
                    "meanRT"  = mean(subset$RT/1000),
                    "varRT"   = var(subset$RT/1000),
                  # Index variable: Accuracy (-0.5) vs Speed (0.5) condition
                    "Xi"  = as.integer(cond>33)-0.5,
                  # Arbitrary scale of stimulus configuration | 0 is 50/50 black and white 
                    "Xs"  = ((cond-1) %% 33 - 16)/5))
}

Compute summary statistics from data

# Initialize an empty output data frame (df)
tmp <- matrix(0,nrow = max(data$cond),ncol = 6)
df <- as.data.frame(tmp)
colnames(df) <- c("nTrials", "sum_accuracy", "mean_rt",
                  "variance_rt", "Xi", "Xs")

# Populate the df output using the ez_summaries function
for(i in 1:max(data$cond)){
  df[i,] <- ez_summaries(data[which(data$cond==i),])
}

# Remove the two ambiguous conditions (17 and 50, with 50/50 black and white)
df <- df[-which(df$Xs==0),]
head(df,3)

##   nTrials sum_accuracy   mean_rt variance_rt   Xi   Xs
## 1      30           30 0.4917333  0.03115648 -0.5 -3.2
## 2      24           24 0.4888750  0.03118168 -0.5 -3.0
## 3      36           36 0.4994722  0.02872146 -0.5 -2.8

# Compute accuracy rate per condition from the summary data
df$acc_rate <- df$sum_accuracy/df$nTrials

The model

EZHBDM

The first component of the model is what we called ‘the proxy’ model: We use the sampling distributions for the three summary statistics computed from the data, (i.e., the number of correct responses \(\dot{T}\) and the mean \(\dot{M}\) and variance \(\dot{V}\) of the response times), which are built from the forward-system of equations in the EZDDM that characterize the proportion of correct responses PC, and the mean and precision of the RTs MRT, PRT.

\[\begin{align*} \dot{T} &\sim \mbox{Binomial}(PC,N)\\ \dot{V} &\sim \mbox{Normal}(\frac{1}{PRT},\frac{N-1}{2(PRT)^2})\\ \dot{M} &\sim \mbox{Normal}(MRT, PRT \times N)\\ PC &= \frac{1}{\exp(-\alpha\nu)+1}\\ MRT &= \tau + \frac{\alpha}{2\nu}\frac{\exp(-\alpha\nu)-1}{\exp(-\alpha\nu)+1}\\ PRT &= \frac{\alpha}{2\nu^3}\{\frac{1-2\alpha\nu\exp(-\alpha\nu)-\exp(-\alpha\nu)^2}{(\exp(-\alpha\nu)+1)^2}\} \end{align*}\]

Metaregression structure with a nonlinear regression on the drift rate

The second component of the model is the metaregression structure we build to explain the DDM parameters as a function of covariates of interest.

We incorporate an effect (\(\beta\)) of instruction (i.e., \(x_i\), Xi) on the bound parameter (\(\alpha\)).

\[\begin{align*} \alpha &\sim \mbox{Normal}(\mu_\alpha+\beta X_i,\sigma_\alpha) \end{align*}\]

We also include a nonlinear regression on the drift rate \(\nu\) using instruction (i.e., \(x_i\), Xi) and stimulus configuration (i.e., \(x_s\), Xs) as predictors. For the later, we used the absolute value (i.e., \(|x_s|\), abs(Xs)) to represent the task difficulty getting easier as the black/white balance departs from \(50\%\).

\[\begin{align*} Z &= \Phi(\beta_1+\beta_2 |X_s|+\beta_3X_i|X_s|)\\ \nu_\mbox{pred} &= \mu_\nu+\beta_0 Z+\beta_4 X_i\\ \nu &\sim \mbox{Normal}(\nu_\mbox{pred},\sigma_\nu) \end{align*}\]

In the present example, we focus on the regression parameters capturing the effects of instruction (\(\beta_3\), Beta3, and \(\beta_4\), Beta4).

Write the model in JAGS

modelFile <- here("output", "BUGS-models", "demo_brightness_model.bug")

model <- write("
model {
        ##### Priors for hierarchical DDM parameters
        betaweight ~ dnorm(0.00, 1.00)
        beta0 ~ dnorm(0.00, 1.00)
        beta1 ~ dnorm(0.00, 1.00)
        beta2 ~ dnorm(0.00, 1.00)
        beta3 ~ dnorm(0.00, 1.00)
        beta4 ~ dnorm(0.00, 1.00)
        bound_mean ~ dnorm(1.50, (0.20^-2))T( 0.10, 3.00)
        drift_mean ~ dnorm(0.50, (0.50^-2))
        nondt_mean ~ dnorm(0.30, (0.06^-2))T( 0, )
        bound_sdev ~ dunif(0.01, 1.00)
        drift_sdev ~ dunif(0.01, 3.00)
        nondt_sdev ~ dunif(0.01, 0.50)
        
        # Hierarchical distributions of individual DDM parameters.        
        for (p in 1:length(meanRT)) {
            # Here, we focus on the drift rate
            drift_pred[p] = beta0*phi(beta1 + beta2*abs(Xs[p]) + beta3*Xi[p]*abs(Xs[p])) + beta4 * Xi[p] + drift_mean
            drift[p] ~ dnorm(drift_pred[p], (drift_sdev^-2))
            bound_pred[p] = bound_mean + betaweight * Xi[p]
            bound[p] ~ dnorm(bound_pred[p],(bound_sdev^-2))T( 0.10, 3.00)
            nondt[p] ~ dnorm(nondt_mean, (nondt_sdev^-2))T( 0.05, )
        
            # Forward equations from EZ DDM
            ey[p]  = exp(-bound[p] * drift[p])
            Pc[p]  = 1 / (1 + ey[p])
            PRT[p] = 2 * pow(drift[p], 3) / bound[p] * 
                     pow(ey[p] + 1, 2) / (2 * -bound[p] * 
                     drift[p] * ey[p] - ey[p] * ey[p] + 1)
            MDT[p] = (bound[p] / (2 * drift[p])) * 
                     (1 - ey[p]) / (1 + ey[p])
            MRT[p] = MDT[p] + nondt[p]
            
            # Noiseless predictions from forward EZ DDM
            ey_pred[p]  = exp(-bound_pred[p] * drift_pred[p])
            Pc_pred[p]  = 1 / (1 + ey_pred[p])
            PRT_pred[p] = 2 * pow(drift_pred[p], 3) / bound_pred[p] * 
                     pow(ey_pred[p] + 1, 2) / (2 * -bound_pred[p] * 
                     drift_pred[p] * ey_pred[p] - ey_pred[p] * ey_pred[p] + 1)
            MDT_pred[p] = (bound_pred[p] / (2 * drift_pred[p])) * 
                     (1 - ey_pred[p]) / (1 + ey_pred[p])
            MRT_pred[p] = MDT_pred[p] + nondt_mean
        
            # Sampling distributions for summary statistics
            correct[p] ~ dbin(Pc[p], nTrials[p])
            varRT[p]   ~ dnorm(1/PRT[p], 0.5*(nTrials[p]-1) 
                                         * PRT[p] * PRT[p])
            meanRT[p]  ~ dnorm(MRT[p], PRT[p] * nTrials[p])
      }
}", modelFile)

JAGS

Specify JAGS setup

set.seed(seed)

# General setup
n.chains  <- 4;      n.iter    <- 5000
n.burnin  <- 250;    n.thin    <- 1

# Pass data to JAGS
data_toJAGS <- list("nTrials"  =  df$nTrials,
                    "meanRT"   =  df$mean_rt,
                    "varRT"    =  df$variance_rt,
                    "correct"  =  df$sum_accuracy,
                    "Xi"   =  df$Xi, "Xs"   =  df$Xs)

# Specify parameters to keep track of
parameters <- c('beta3', 'beta4', 'drift', 'drift_pred',
                "Pc_pred", "MRT_pred", "PRT_pred", "Pc", "PRT", "MRT")

# Prepare initial values
myinits <- rep(list(list()), n.chains)
for(i in 1:n.chains){
    myinits[[i]] <- list(drift = rnorm(length(data_toJAGS$nTrials),0,0.1))
}

Run JAGS

set.seed(seed)

start <- Sys.time()
samples <- jags(data=data_toJAGS,
                parameters.to.save=parameters,
                model=modelFile,
                n.chains=n.chains,  n.iter=1000,
                n.burnin=100,  n.thin=n.thin,
                DIC=T, inits=myinits)

## module glm loaded

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 192
##    Unobserved stochastic nodes: 204
##    Total graph size: 3150
## 
## Initializing model

end <- Sys.time()

## [1] "No Rhat greater than 1.05"

Extract posterior samples

##### Drift rate parameters
# Recovered drift rates
drift <- samples$BUGSoutput$sims.list$drift
# Effects of instruction
beta3 <- as.vector(samples$BUGSoutput$sims.list$beta3) # Main
beta4 <- samples$BUGSoutput$sims.list$beta4 # Interaction
# Fitted values / Predicted drift rates
drift_pred <- samples$BUGSoutput$sims.list$drift_pred

##### Summary statistics predicted from the predicted drift and boundary
accRate_hat   <- samples$BUGSoutput$sims.list$Pc_pred
rtMean_hat <- samples$BUGSoutput$sims.list$MRT_pred
rtVar_hat  <- 1/samples$BUGSoutput$sims.list$PRT_pred

##### Summary statistics computed from the recovered drift and boundary
Pc   <- samples$BUGSoutput$sims.list$Pc
PRT  <- samples$BUGSoutput$sims.list$PRT
MRT  <- samples$BUGSoutput$sims.list$MRT

Obtain posterior predicted summary statistics

set.seed(seed)

# Number of trials per condition
nTrials <- df$nTrials
# Number of posterior samples
n <- nrow(Pc) 
# Number of conditions
J <- ncol(Pc)
# Empty matrices to store posterior predictions
pp_accRate <- matrix(NA, nrow=n, ncol=J)
pp_meanRT  <- matrix(NA, nrow=n, ncol=J)
pp_varRT   <- matrix(NA, nrow=n, ncol=J)
# Obtain posterior predictions using sampling distributions
#        and the summary statistics derived from the recovered
#        drift and boundary parameters
for(i in 1:J){
  correct  <-  rbinom(n,nTrials[i],Pc[,i])
  pp_accRate[,i] <- correct/nTrials[i]
  pp_varRT[,i]   <- rnorm(n,1/PRT[,i], sqrt(2/((nTrials[i]-1) * PRT[,i] * PRT[,i])))
  pp_meanRT[,i]  <- rnorm(n,MRT[,i],sqrt(1/(PRT[,i]*nTrials[i])))
}

Results

Instruction effect on drift rate slope

line.color  <- "#A94442"
nBreaks <- 75
hist(beta3, freq = FALSE, breaks = nBreaks, col="#9EA3A8", border = "#9EA3A8", ann=F, axes = T)
lines(density(beta3), lwd=4, col=line.color)
legend("topleft",c("Histogram","KDE"),col=c("#9EA3A8",line.color),lwd=4, cex=1.1, bty = "n")
mtext("Instruction effect on drift rate slope", cex=1.2, f=2, side=3, line=0.8)
mtext("Density",side=2,line=2.15)
mtext("Beta 3",side=1,line=2)
box(col="black")

Instruction main effect on drift rate

hist(beta4, freq = FALSE, breaks = nBreaks, col="#9EA3A8", border = "#9EA3A8", ann=F, axes = T,xaxs = "i", yaxs = "i", xlim=c(0,1.6), ylim = c(0,2.2))
lines(density(beta4), lwd=4, col=line.color)
legend("topleft",c("Histogram","KDE"),col=c("#9EA3A8",line.color),lwd=4, cex=1.1, bty = "n")
mtext("Instruction main effect on drift rate", cex=1.2, f=2, side=3, line=0.8)
mtext("Density",side=2,line=2.15)
mtext("Beta 4",side=1,line=2)
box(col="black")

# Identify conditions to plot on x axis
x_values <- unique(df$Xs)
# Insert a 'jump' in between these values
fit_x <- c(x_values[1:16],NA,x_values[17:32])
full_x <- c(fit_x,NA,fit_x+6.8)

Drift rates predicted and recovered

#### Get mean posterior estimates and percentiles for plotting
# Drift recovered
means <- apply(drift, 2, mean)
percentiles <- apply(drift,2, quantile, probs=c(0.025,0.975))
lower_percentiles <- percentiles[1,]  #  2.5%
upper_percentiles <- percentiles[2,]  # 97.5%
# Drift predictions (i.e., fitted values)
preds <- apply(drift_pred,2,mean)

Posterior predictive checks

## Get mean posteriors and key percentiles
# Accuracy rate
pp.accRate <- apply(pp_accRate, 2, mean)
pp.perc.accRate <- apply(pp_accRate,2, quantile, probs=c(0.025,0.975))
# Mean RT-correct
pp.rtMean <- apply(pp_meanRT, 2, mean)
pp.perc.rtMean <- apply(pp_meanRT,2, quantile, probs=c(0.025,0.975))
# RT-correct Variance
pp.rtVar <- apply(pp_varRT, 2, mean)
pp.perc.rtVar <- apply(pp_varRT,2, quantile, probs=c(0.025,0.975))

Brightness perception: Metaregression model example