Lets say we have hypothetical participant that is presented with 3 stimuli conditions: a flashing dot ($C_x$), a sound blip ($C_y$) and combination of both ($C_z$). We ask this participant to respond as quickly as possible to the onset of the stimuli and we measure the reaction time (RT). It's been observed that responses are often observed to be substantially faster when both stimuli are presented simultaneously ($C_z$), compared to condition when they are presented separately ($C_x$ and $C_y$), and this effect has been labeled redundant signals effect (e.g. Hershenson, 1962). Different models account for this effect, one of those models being race model (Miller, 1982). Quoting Colonius et al (2006) definition:

Raab (1962) was the first to propose a race model for simple RT such that (a) each individual stimulus elicits a detection process performed in parallel to the others and (b) the winner’s time determines the observable RT. This model suggests that RSE is generated by statistical facilitation: If detection latencies are interpreted as (nonnegative) random variables, the time to detect the first of several redundant signals is faster, on average, than the detection time for any single signal.

Urlich et al (2007) continues that:

According to race models, the observed RT distributions should satisfy, for every value of $t$, the so-called race model inequality (Miller, 1982):

$F_z(t) \le F_x(t) + F_y(t), t > 0$

where $F_x$ and $F_y$ are the cumulative density functions (CDFs) of RT in the two single-stimulus conditions $C_x$ and $C_y$, respectively, and $F_z$ is the CDF of RT in the redundant-stimulus condition $C_z$.

I wanted to test this model, but I found it hard to translate available Matlab code into R. Matlab code is available in the appendix of this paper by Urlich et al (2007), which can can be downloaded here, in case you don't have access to Behavioural Research Methods journal. To summarize Urlich et al (2007) approach, his algorithm goes through the follwing steps:

Step 1: Estimating empirical cumulative density functions (CDFs) for every participant (we got just one in the example here) and every stimulus condition ($C_x$, $C_y$, $C_z$). As Urlich et al (2007) points out:

these CDFs will be denoted as the individual CDFs. Specifically, let $G_x$, $G_y$, and $G_z$ be the individual estimates of two single-stimulus CDFs, $F_x$ and $F_y$, and the redundant CDF, $F_z$, respectively.

I thought I managed to compute $G_x$, $G_y$, and $G_z$, but I might be wrong:

#example reaction times for single hypothetical participants for 3 experimental  conditions
#this is the same dat set as one used by Urlich et al 2007
cx <- c(244, 249, 257, 260, 264, 268, 271, 274, 277, 291)
cy <- c(245, 246, 248, 250, 251, 252, 253, 254, 255, 259, 263, 265, 279, 282, 284, 319)
cz <- c(234, 238, 240, 240, 243, 243, 245, 251, 254, 256, 259, 270, 280)

#here I calculate quantiles corresponding to the given probabilities
#I use the default quantile algorithm (type 7 see R manual for 'quantile')
#it seems right, but I'm not 100% sure
psq <- seq(0,1,0.1)
gx <- quantile(cx, probs=psq, names=F)
gy <- quantile(cy, probs=psq, names=F)
gz <- quantile(cz, probs=psq, names=F)

#here I visualize the quantiles, looks more or less like in Urlich etal 2007
xyz <- data.frame(rt=c(gx,gy,gz), 
                  cond=rep(c("gx(t)", "gy(t)","gz(t)"), each=length(psq)))
xyplot(perc ~ rt, groups = cond, xyz, col=1, lty=3:1, type=c("g","l","p"), 
       key=list(lines=list(lty=3:1), text=list(labels=levels(xyz$cond))))

enter image description here

Step 2: From $G_x$ and $G_y$, one computes the bounding sum $B(t) = G_x(t) + G_y(t)$ (for each participant), which provides an estimate for each participant of the upper bound on the right side of race model equation above. But then, I'm not sure how to compute the race model bounds, which seems trivial, but doesn't really work. The authors simply show the equation as , and I applied it in R as:

b = gx + gy

but then when I visualize it, it's wrong:

xyz <- data.frame(rt=c(gx,gy,gz,b), 
                  cond=rep(c("gx(t)", "gy(t)","gz(t)","gx(t)+gy(t)"), each=length(psq)))
xyplot(perc ~ rt, groups = cond, xyz, col=1, lty=4:1, type=c("g","l","p"), 
       key=list(lines=list(lty=4:1), text=list(labels=levels(xyz$cond))))

enter image description here

Step 3 and 4: Percentile values are computed from $G_z$ and $B$ for each participant at certain prespecified probabilities that must be the same for all participants. Finally, the percentile values are aggregated across participants, which is omitted in Urlich etal (2007) example.

I hope my explanation is a bit clearer now then it was before.


Miller, J. (1982). Divided attention: evidence for coactivation with redundant signals. Cognitive psychology, 14(2), 247–79.

Ulrich, R., Miller, J., & Schröter, H. (2007). Testing the race model inequality: an algorithm and computer programs. Behavior research methods, 39(2), 291–302.

Colonius, H., & Diederich, A. (2006). The race model inequality: interpreting a geometric measure of the amount of violation. Psychological review, 113(1), 148–54. doi:10.1037/0033-295X.113.1.148

Raab, D. H. (1962). Statistical facilitation of simple reaction times. Transactions of the New York Academy of Sciences, 24, 574-590.

Hershenson, M. (1962). Reaction time as a measure of intersensory facilitation. Journal of Experimental Psychology, 63, 289-293.

  • $\begingroup$ You have not made it easy for us to understand the question here. What, for example, is the "race model"? Unless someone has read the articles you cite and done the actual analysis themselves how could we help? Your original question appears to be mostly copied (without citation) from an article in the behavioural psychology literature. Something that is "well known" there will not necessarily be well known to CV readers. $\endgroup$ Feb 9, 2013 at 1:50
  • $\begingroup$ Ok, you're right @PeterEllis - thank you for your feedback. I've tried to revise the question with much more detail, I hope it's a bit clearer now. $\endgroup$ Feb 11, 2013 at 17:51
  • 2
    $\begingroup$ It may indeed be clearer, but it still doesn't specify well-defined questions. I assume the issue is that somewhere in either step 1 (estimating the CDF's) or step 2 (calculating race model bounds). With that said, looking at step 2, I don't see how that would possibly be a race model bound. It implies that the bound is as high as BOTH reaction times combined? I would think that the bound would be similar to min(X,Y) not X+Y. $\endgroup$
    – Namey
    Feb 11, 2013 at 22:26

2 Answers 2


UPDATE: All the custom versions of functions rewritten in R and used in this answer are added in the bottom for better clarity.

I guess. You probably knew where's gone wrong. As you said

it seems right, but I'm not 100% sure

Yes. This is the start of the differences between your data and Ulrich et al's.

cx <- c(244, 249, 257, 260, 264, 268, 271, 274, 277, 291)
cy <- c(245, 246, 248, 250, 251, 252, 253, 254, 255, 259, 263, 265, 279, 282, 284, 319)
cz <- c(234, 238, 240, 240, 243, 243, 245, 251, 254, 256, 259, 270, 280)   

Here is the correct psq:

psq <- probSpace(10); psq
[1] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95    

Here is your psq:

psq <- seq(0,1,0.1); psq
[1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Also (or should say therefore), your gx, gy and gz are different from Ulrich et al's. Here I used Ulrich et al's ties function (off course, after rewriting it in R) to generate a data frame that stores (1) the unique response time (RT) data, (2) the number of iteration for each RT datum, and (3) the cumulative frequency. This is in line with Urlich et al's MATLAB function (note that they stored them as a structure, I suppose). These three pieces of information are then served as inputs to Ulrich et al's CDF function, which differs from the R's quantile function. I rewrote their MATLAB CDF function to work with my ties data frame output. They in fact wrote their CDF function, following the equation in Appendix A, not the equation (2) in the main article. Here is probably the place you might not notice, so did not work out the question. Basically, they wanted to handle the problem of ties, so they wrote their variation of CDF function.

dfx <- ties(cx)
dfy <- ties(cy)
dfz <- ties(cz)
tmax <- max(cx,cy,cz)

gx <- cdf.ulrich(data=dfx, maximum=tmax)
gy <- cdf.ulrich(data=dfy, maximum=tmax)
gz <- cdf.ulrich(data=dfz, maximum=tmax)

I then added gx and gy together and stored them to b. Because R does element-by-element calculation, if one uses the operator as usual, I did not bother to use a for loop, as Ulrich et al.'s did. I then used Ulrich et al's GerPercentile function (rewrote it in R) to calculate the percentiles for each individual CDF, including the b. Also note that you need tmax to confine the domain of RT. This was possible the place you had not noticed.

b <- gx + gy
xp <- GetPercentile(psq, gx, tmax);
yp <- GetPercentile(psq, gy, tmax);
zp <- GetPercentile(psq, gz, tmax);
bp <- GetPercentile(psq, b, tmax);

Finally, construct a data frame to draw the figure. Here is the panel F in their Figure 1.

gdf <- data.frame(RT =c(xp,yp,zp,bp), Probability =rep(psq, 4),
                  Condition =rep(c("gx(t)", "gy(t)","gz(t)","gx(t)+gy(t)"), each=length(xp)))
panelf <- ggplot(gdf, aes(x = RT, y = Probability, group=Condition,
          colour=Condition, shape=Condition)) + 
          geom_point() + geom_line() 
panelf + coord_cartesian(xlim = c(230, 330), ylim=c(-.01,1.01)) +
         theme(legend.position= c(.85, .20),  
               legend.title = element_text(size=12),
               legend.text = element_text(size=12))

I hope these would be helpful.

Figure 1 Panel F

Also I should mention that you may need to load the required libraries (here I have used ggplot2 and grid) and source those home-made functions, if you did not put them together in one RaceModel function as Ulrich et al did. I suggest write individual file for each function, which was what I did and source them. Some of the functions you could translate directly from the equations in Ulrich et al's article and others from their MATLAB or Pascal code. MATLAB uses parentheses, (), to index vector, but R uses square bracket, [].


Here is the probSpace function. I wrote it according to Ulrich et al's equation (3) (p. 293) by myself, so you probably won't be able to find it in any R packages.

probSpace <- function(len){ 
# PROBSPACE Determine the interval of percentiles 
#   The function used equation (3) in Ulrich, Miller and Schroter (2007)
  P <- numeric(len);
  for(i in 1:len){
    P[i] <- (i - .5) / len;

Here is Ulrich et al.'s (2007) CDF function (p. 296 and 299). I wrote it according to their MATLAB code and the equation in Appendix A.

cdf.ulrich <- function(data=NULL, maximum=3000){
  # Create a container, whose length is the longest data vector
  # data is the output from `ties` function
  U <- data[,1]
  R <- data[,2]
  C <- data[,3]
  G <- numeric(maximum);

  # The length of the processed data vector, trimming off ties, if there is any.
  k <- length(U);  # U contains data in millisecond, e.g., 320 ms etc. 

  # The last element of the cumulative frequency supposely is the 
  # length of the data vector.
  n <- C[k]

  for(i in 1:k) { U[i] <- round(U[i]) }

  # from 1 ms to the first value of the data set, their probability should be 0.
  for(t in 1:U[1]) { G[t] <- 0 }   

  for(t in U[1]:U[2]){
  G[t] <- ( R[1]/2 + (R[1]+R[2]) / 2*(t-U[1]) / (U[2] - U[1]) ) / n;

  for(i in 2:(k-1)){
    for(t in U[i]:U[i+1]){
      G[t] <- (C[i-1] + R[i] / 2+(R[i] +R[i+1]) / 2*(t-U[i]) / (U[i+1] - U[i])) / n;

  for(t in U[k]:maximum){
    G[t] <- 1;

Below are the ties() and GetPercentile() functions written by Ulrich et al. (2007). I merely rewrote them into R, based on their description in the article, the equations and MATLAB codes they provided.

Firstly, the ties function:

ties <- function(W){
  # Count number k of unique values
  # and store these values in U.
  U <- NULL; W <- sort(W); n = length(W); k = 1; U[1] <- W[1]
  for (i in 2:n) {
    if (W[i] != W[i-1]) {
      k <- k+1;
      U <- cbind(U, W[i])
  U <- U[1,]

 # Determine number of replications R
 # k is the length of the vector, after trimming off the ties
  R <- numeric(k) 
  for (i in 1:k){
    for (j in 1:n){
      if (U[i] == W[j]) R[i] <- R[i] + 1;

  # Determine the cumlative frequency
  C <- numeric(k) 
  C[1] <- R[1]
  for(i in 2:k){
    C[i] <- C[i-1] + R[i];
  res <- list(U, R, C)
  names(res) <- c("U", "R", "C")

Then, the GetPercentile function:

GetPercentile <- function( P, G, tmax ){
# Determine minimum of |G(Tp[i]) - P[i]|
  np <- length(P);
  Tp <- numeric(np)
  for( i in 1:np) {
    cc <- 100;
    for(t in 1:tmax) {
      if ( abs(G[t] - P[i]) < cc ) {
        c <- t;
        cc <- abs(G[t] - P[i]);

    if( P[i] > G[c] ){
      Tp[i] <-  c + (P[i] - G[c]) / (G[c+1] - G[c]);
    } else {
      Tp[i] <- c + (P[i] - G[c]) / (G[c] - G[c-1]);
 return( Tp )
  • $\begingroup$ Thank you, excellent replication or Urlich et al (2007) method but using R! I slightly edited your answer to make it flow a bit better. Cheers mate! $\endgroup$ Apr 18, 2013 at 17:09
  • $\begingroup$ this is such an excellent answer. $\endgroup$ Feb 2, 2021 at 14:52

In short, I think that at the stage of b you find yourself combining RTs when what you should be combining are proportions. This is speculative, but I took the stab below to get things rolling. I do love me some race models and have been meaning to implement one form or another in R for a while, but never got around to it. If you do complete this project, I hope you might share.

With those disclaimers in place, I'll add... I don't know MATLAB code... but I love speculating.

The MATLAB code to calculate the bounding sum is:

for t=1:tmax;

This looks like a loop for all possible values of RT ranging from 1 to the maximum of cx, cy, and cz. The bounding sum therefore is not the sum of the RTs at various quantiles as suggested by your code above. In fact, in the MATLAB code it looks like Gx is the result of the CDF function across the full time-span. For example, from the mathworks website I can see that cdf('Normal',-2:2,0,1) is expected to yield 0.0228 0.1587 0.5000 0.8413 0.9772. But what the code actually calls is CDF, which is defined in Appendix B. I'd try to turn it into an R function, but without understanding MATLAB data structures it is kind of a nightmare. But, the description of what is happening is in the referred to text... "Following the standard procedures for percentile estimation (see, e.g., Gilchrist, 2000), each step function is next used to generate a corresponding cumulative frequency polygon... [etc]" ... "these values are stored, millisecond by millisecond, in the vectors Gx, Gy, and Gz in the MATLAB program shown in Appendix B". Hey, cool, Gx, Gy, and Gz are clearly results from CDF, so... we know how we need to specify CDF (in abstract). Now it is just a coding problem. I bet you can do that, if you don't, I love co-authorship (half - j/k).

At any rate, I think this means that at this stage in your code b should be giving you something like the sum of proportions correct for various levels of RT not the sum of RT values (as you have coded). That seems a little odd to me and I imagine I'm missing a step because the proportion that will have completed the race when .5 are due to have completed race 1 and .5 are due to have completed race 2 is certainly not 1 for the joint race... but I might be wrong because the authors claim that "Because the sum is only compared against the values of Gz, the computation of the sum can be terminated at the first value of t for which the sum exceeds a value of 1.0"... which seems a little scary, but is certainly beyond the scope of the question and would mean presuming I know something the authors don't which seems unlikely.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.