Data Analysis

serial_reaction_time_mixture_model.R: A script for fitting (in JAGS) a mixture model of log-transformed reaction times from a serial reaction time learning task. The goal in developing this model was to account for bimodal distributions of RTs due to a combination of reactive and anticipatory responses.

multilevel_multinomial_logistic_regression.R: A script for fitting a multilevel multinomial logistic regression model with item-specific intercepts and subject-specific intercepts and other predictors in JAGS. This model uses Gelman & Hill’s scaled-inverse-Wishart approach to model the covariance between the subject-specific parameters.

multilevel_sub_x_item_bernoulli_logistic_regression.R: A script for fitting a multilevel logistic regression model with subject- and stimulus-specific intercepts with JAGS.

MI_ROC_function_w_xval.R: A function for calculating ROC curves, area under ROC curves, hit and false alarm rates, and various other things, for random-intercept logistic regression models and multiply-imputed data sets. It is unlikely that anyone else will need exactly the code I wrote, but it should have chunks that can be readily adapted to similar cases. I will provide a detailed description at some point soon.

mvconfit.R/mvconfit.xval.R: Functions to fit multivariate linear models allowing both mean vectors and covariance matrices to differ across various discrete predictor values (as described in detail in this paper), with AIC and BIC fit statistics (mvconfit.R) and with 10-fold cross-validation and associated log-likelihood values (mvconfit.xval.R). Hyper-specific to the work described in the linked paper, but likely has chunks that can be adapted to similar analyses.

cilim.R: A function for estimating (using the ‘rising water’ method) the limits of the highest density interval for a set of samples from a posterior parameter distribution.

r.hat.R: A function for calculating \hat{R} values for fitted parameters stored in a JAGS output list. Note that \hat{R} assumes that your parameter chains are initially overdispersed, but that, by default, JAGS does not randomize the starting points of different chains. So, you either need to initialize the model yourself to ensure overdispersion, or you need to accept responsibility for not meeting this assumption of \hat{R}. A function for calculating generalized d^\prime between two multivariate normal densities.

grt_bugs_script.R: R script for reading in data and fitting Bayesian multilevel 2×2 Gaussian General Recognition Theory/Multidimensional Signal Detection Theory/Bivariate Probit model in WinBUGS. Appropriate data are collected in a 2 x 2 feature-complete factorial identification experiment (i.e., stimuli and responses consisting of the factorial combination two levels on each of two dimensions; e.g., dimension one = {red, purple}, dimension two = {square, circle} → red square, red circle, purple square, purple circle; more generally, we can refer to A = “level one” and B = “level two” on each dimension, and describe the stimuli and responses in terms of their discrete (x,y) feature specifications: AA, AB, BA, BB). The rows/columns of each individual subject’s confusion matrix are assumed to correspond to the feature specifications AA, AB, BA, BB; the full data set consists of multiple subjects’ confusion matrices stacked vertically, such that S1’s matrix spans rows 1 to 4, S2’s matrix spans rows 5 to 8, …, SN’s matrix spans rows 1+4(N-1) to 4+4(N-1).

BUGS can’t calculate bivariate normal integrals (neither can JAGS), and all of the numerous approximation formulas I tried produced unacceptably large errors, so this model reads in a great big matrix of pre-calculated bivariate normal CDF values and uses trilinear interpolation to get the predicted probabilities for each set of sampled parameters. Hence, this file (13 MB) is required. Be warned that, largely because of this issue, it can take quite a while for the model to initialize in WinBUGS.

plot.grt.R/plot.grt.simple.R: R functions for plotting the fitted model from grt_bugs_script.R described above. Both take the BUGS output, a subject index, and a number of optional input arguments. Both functions produce plots like the panels in Figures 4 and 6 in this paper.

ggplot_white_background.R: I usually code my plotting functions from scratch, but I use ggplot on occasion, too. When I do, I usually want a plain white background. This code produces just such a background.

plot.grt.predobs.R: R function for plotting the observed and predicted identification-confusion probabilities from grt_bugs_script.R described above. Takes the BUGS output and data as input, produces a plot like Figures 3 and 5 in this paper.

grtnewtrap.m: Newton-Raphson algorithm for finding maximum likelihood estimates for a GRT model fit to identification-confusion count data. The function takes a 4 x 4 confusion matrix with frequency counts as input, and it returns maximum likelihood mean and correlation estimates for latent bivariate normal distributions (it assumes that the decision bounds/cutoffs are the coordinate axes).  The output is a structure with seven fields: p: fitted parameters; M: the original data; pr: predicted confusion probabilities; nll: negative log-likelihood; aic: AIC; bic: BIC; icomp: another fit statistic based on the likelihood, the information matrix, and the number of parameters.

grt_plotfit_newtrap.m: Function for plotting fitted GRT model. Takes (part of) the output of grtnewtrap (p; see above), the confusion matrix that the model is fit to, and an optional figure label argument, and plots a bird’s eye view of the fitted model. Here’s an example for a model fit to one of the confusion matrices from a speech perception experiment (see here for lots and lots of details):

Fitted GRT model

The main panel shows equal likelihood contours for fitted perceptual distributions and decision bounds; the rectangular panels to the left and on the bottom show fitted marginal perceptual distributions, and the panel in the bottom left shows the observed and predicted confusion probabilities (the better the fit, the closer the dots in this panel will be to the diagonal dotted line).

To be updated again soon… (3.18.12)

Comments are closed.