No evidence for periodicity in reaction time histograms


In my last lab we discussed findings on the periodicity of reaction times (e.g. referred in Van Rullen 2003). These studies are usually old (Starting with Harter 1968, Pöppel 1968), with small N and not many trials. There was also a extensive discussion in the Max-Planck Journal “Naturwissenschaften” in the 90s (mostly in German, e.g. Vorberg & Schwarz 1987). A methodological critique is from Vorberg & Schwarz 1987. More discussion in Gregson (Gregson, Vorberg, Schwarz 1988). A new method to analyse periodicity is proposed by Jokeit 1990.

This is the newest research I could find on this topic.


Analysing a large corpus of RT-data

I stumbled upon a large reaction time dataset (816 subjects, á 3370 trials, 2.3 million RTs)  from the English Lexicon Project (Balota et al 2007) and decided to look for these oscillations in reaction times .

After outlier correction (3*mad rule, see below), I applied a fourier-transformation on the histogram (1ms bins, accuracy of RT=1ms). Then I looked for peaks in the spectrum which are consistent over subjects.

Each subject is one line, no effect is visible here (in a log-scaled y-axis also no effect can be seen). The range (above ~7Hz) of the subject variance is roughly between +0 and 100


The following graph summarizes  the above graph (blue-smoother-curve = loess, span=0.1, each dot = mean over 800 subjects)):


There are no peaks in the spectrum which I consider as consistent over subjects. I included higher frequencies (up to 250Hz) to get a visual estimate of the noise level (at such high frequencies, an effect seems utterly implausible). But of course, I’m ignoring  within subject information (i.e. a mixed model of some sort could be appropriate).


In this large dataset, I cannot find periodicities of reaction time.


My approach may be too naive. I’m looking for more powerful ways to analyse these data. If you have an idea please leave a comment! I’m also not suggesting that the effects e.g. in Pöppels data are not real. Maybe there is a mistake in my analysis, I don’t know the data by heart, it might depend on the task employed …


I had results like in Jokeit 1990 (but with 50Hz not with 100Hz), when I was using a bin-width of 5ms to 10ms bins.  The peak (in the figure with 6ms bins => 150hz) shifted depending on bin-size. I’m not perfectly sure, but I think it has to do with how integers are binned. In any case, if the effect is real and not an artefact of bin-width, it has to show up also with higher bin-sizes. Please note that Jokeit 1990 used a different methodology, he calculated the FFT on the histogram of reaction time **differences**.

I tried to use density estimates, but so far failed to get better results.

Outlier plot

Percentage of trials marked as outliers. This is well in the recommended range of 10% (Ratcliff).


Bolota, D.A., Yap, M.J., Cortese, M.J., Hutchison, K.A., Kessler, B., Loftis, B., Neely, J.H., Nelson, D.L., Simpson, G.B., & Treiman, R. (2007). The English Lexicon Project. Behavior Research Methods, 39, 445-459. –

d = rbind(read.table('C:\\Users\\behinger\\cloud\\PhD\\exercises\\reactiontimes\\elexicon_1to100.csv',header=T),

d$Sub_ID = factor(d$Sub_ID)
d$D_RT = as.integer(d$D_RT)

# 3*mad outlier correction
d = d%>%group_by(Sub_ID)%>%mutate(outlier = abs(D_RT-median(D_RT))>3*mad(D_RT,na.rm = T))
d$outlier[d$D_RT<1] = TRUE
d$outlier[$D_RT)] = TRUE

# outlier plot
ggplot(d%>%group_by(Sub_ID)%>%summarise(outlier=mean(outlier)),aes(x=outlier))+geom_histogram(binwidth = 0.001)

fft_on_hist = function (inp){
  maxT = 4
  minT = 0
  fs = 1000
  h = hist(inp$D_RT,breaks = 1000*seq(minT,maxT,1/fs),plot = F)
  h = h$counts;
  # I tried to use density estimates instead of histograms, but it was difficult
  #h = density(inp$D_RT,from = minT,to=4000,n=4000)
  #h = h$y
  f = fft(h)
  f = abs(f[seq(1,length(f)/2)])
  return(data.frame(power = f, freq = seq(0,fs/2-1/maxT,1/maxT)))

d_power = d%>%subset(outlier==F)%>%group_by(Sub_ID)%>%do(fft_on_hist(.))


ggplot(d_power%>%group_by(freq)%>%summarise(power=mean(power)),aes(x=freq,y=(power)))+geom_point()+stat_smooth(method='loess',span=0.1,se=F,size=2)+xlim(c(10,250))+ ylim(c(47,53))


Interaction and Effect/Sum Coding

Some time ago I wrote a blog post on dummy & effect coding. I made some new plots to visualize why the interaction in sum/effect is coded as it is.

Let’s take a typical 2×2 design. We have two 2-level factors $A$ and $B$ and we also allow for an interaction.
$$y \sim A + B + A:B$$

We code A with -1 / 1 and B with -1 / 1 (depending on the level e.g. On=1, Off=-1)

The interaction is coded as the multiplication of A and B: $A * B$. Therefore if $A$ and $B$ are both in the same level (both “off” or both “on”) we get a $+1$, else a $-1$.

Side remark: This is different in dummy/reference coding, where the interaction only codes what is extra if both A & B are “on” (turns out that the magnitude of the interaction is just double – but this is a story for another time).

In the first figure I added the main effects of $A$ and $B$ as Blue and Purple lines. The main effects in reference coding are relative to the means of the group means.

In order to model the original data points, one needs to add the main effects and the interaction together:

Notice that the way we have to add the interactions and main effects is exactly the multiplication I introduced earlier. That is, if we need to take -1 for $A$ and +1 for $B$, you bet we will need -1 for $A:B$.

One way that I like to think about the interaction in effect coding is to think “What would be my prediction if there would be no interaction?”.

“What if there would be a model without interaction” is marked in black (it’s only using the main effects!). Note that the two black lines are parallel. Adding the red interaction-lines again helps us to move to the original datapoints.

Trains and EEG

I had this blog post as a draft for nearly a year now… and finally here it is: In my last lab we moved buildings. Good scientists we are we did measurements to see if the system works as intended…

The EEG Systems

    • An active 64 electrode system from BrainProducts called ActiCap
    • A passive 128 electrode system from ANT, a tmsi REFA-8 with actively shielded cables
    • A passive 128 electrode system from ANT, eegosport with actively shielded cables

Disclaimer: All data shown here are from different subjects. Their individual alpha-peaks etc. will differ, we also switched methods a bit


We measured 2 subjects with the 64 electrode system. When we checked the frequency spectrum we were very worried.

Besides the typical 50Hz artefact that we know, one can clearly see a $16 \frac{2}{3}$ Hz peak in the beta-band. The artefact is very strong and variable over time.

This can’t be right? So we recorded another subject, just to be sure. We also switched analysis methods a bit

Same thing!

Side remark: The ‘activeShield’ button did not improve the shielding (plots not shown, the support told us that the activeShield function is more thought for passive electrode systems, not our system).


Single Amplifier (64 Electrodes)

We then measured a single subject with the ANT system with a single amplifier.

All looks fine here for now… On the second plot, we introduce the signal-to-noise (SNR) index, which is simply the current frequency-bin divided by neighbouring bins (the average of $n_{-12}$:$n_{-2}$ & $n_{+2}$:$n_{+12}$). This will enhance peaks in the spectrum. As we can see, we do not see any peaks anywhere, all is fine with 64 electrodes.

Two Amplifiers (128 Electrodes)

At first we thought, cool! ANT with active shielded cables is perfect, we will just not record with the brainproducts system. But ANT/tmsi allows to “chain” two amplifiers together to effectively get 127 electrodes. That is, the amplifiers share a ground electrode and a single electrode is bridged, so that the data can be referenced to a amplifier-common electrode. When we did this the following picture emerges:

The 16 Hz is back! :(. The electrode depicted in the spectrum plot is from the second amplifier. In the lower right panel we can see, especially in the second plot, that the second amplifier has a lot of 16Hz artefact but not so much in the first amplifier.

ANT eegosport

We much later received an ANT eegosport for testing purposes. Unfortunately for us (and ANT because we did not upgrade ;)), the artefact was also present with comparable strength. Here an interesting new phenomenon arised: when flipping the amplifiers upside down, the artefact switched from channel 1-32 to 33-64 or vice versa.

Typical artefact behavior

The artefact was non-stationairy over time (sometimes its there, sometimes not, stronger, weaker …). Most often observed in the second amplifier, sometimes in the first. The artefact also was stronger (relatively) around the reference electrode. It does not depend on exact placement in room. We tried a lot of factors, at some point even wrapping things in foil – it did help a bit!

What the hell is going on?

Turns out, the new building is quite close to train tracks. In Northern Germany most trains are powered electrically using overhead lines. Germany uses $16 \frac{2}{3}$ Hz to power these lines. So sometimes when a train comes along we have an artefact. Because we did not want to time our experiments to the train schedule (we a) thought about it, but the DeutscheBahn does not give away information when trains are coming and b) the artefact is not _always_ strong when a train comes, but only sometimes, also not train dependent – but who knows…)

What can we do against it?

It took us the most part of the year (with lots of discussion with ANT & TMSI) to figure out, where exactly the artefact is induced into the signal.

Our solution

  • Keep the cables as close together as possible, all in one orientation. Because the ANT REFA system allows the ~10cm cable to go into various direction, the common-noise which influences these channels can be slightly different. If the electrical field / artefact is strong enough, this slight differences will show up in your data. Thus always keep all cable very similar. We now “comb” the connector cables
  • Keep the cables parallel to the field-lines of the electrical field. I.e. don’t put your cables parallel to the elongated electrical source, this will induce stronger effects. We now lead the cables to the top, so inline with the electric field
  • If you use a cascaded system, the reference-electrode that links the two systems, has to follow both strings of electrodes. I.e. you have to guide it along the connector. If you fail to do so, there will be stronger artefact in one than in the other amplifier (depending on how the cable lays relative to the artefact electrical field). We now take care of this.

Combining all these measures helped us enourmously to reduce the artefact in a stationairy setting. For mobile recordings this will be tougher.

Debugging Lessons Learned

It took us a long time to figure out how to measure the artefact reliably and how to manipulate the system in order to study the artefact. Here are some lessons learned:

  • A waterbucket+salt (=> a head without needing a human) is your friend. BUT: We tried this earlier and failed because the system was DC-maxing out (no signal in range). This was because we put in the whole ground electrode (which has a zink-connector) and not only the tip (AgCl). This gives very large DC offsets due to the zink contact and you cannot measure.
  • Making turn around times super fast. This allows for more tests and easier replication of supposedly found effects. A online/live-script that shows us the current 16Hz power was very helpful here but took quite some programming to get smoothly.
  • You need to show a lot of own initiative, support can help you, but not save you. This problem was very specific for us, a general solution did not exist
  • If there is a non-stationairy noise source, it is very helpful to monitor the artefact with one amplifier while modifying things with the other amplifier. Only like that can you be sure that the artefact occured / your modification actually changed something (and you did not just measured in a low-artefact regime). An interesting phenomenon we saw, was that on one day after ~19.00 the artefact vanishes. We are not sure why this is the case, because trains were still running. The trains passing by do not produce a reliable effect in itself. We think this is because the power to the grid is turned on whenever a train is on the whole section, with a section being multiple kilometers long and possibly reaching into the trainstation (but we do not know this before).


We detrended or highpass filtered the data at 1Hz. Welchs estimate does not like DC-drifts. We generally either bandpass filter the data, or use the pwelch algorithm. We also employed the TMSI matlab-toolbox and programmed a live-script to monitor the 16Hz activity while manipulating the EEG. We also used some normalization to compare the signal to its frequency-bin neighbours. This allows for relative noise effects that show up when re-referencing to a single electrode.

I uploaded the tmsi EEG live-script on github

Tuning your plots

I showed a student how I would improve a plot. The result is I think typical for what I would design. The details change for each paper and I do not have a specific “style” I wanted to emulate. I usually enjoy the before & after of other graphics, so here is mine:

PS: The data shown in the plots are preliminary

Thesis Art: Maria Sokotushchenko

I was a supervisor for Maria Sokotushchenko’s Master’s Thesis.

In her thesis-art, I artistically visualized the brain’s surprise response to a unexpected stimulus change. This response is sorted by how fast subjects responded (late on top, fast in the bottom)


The idea of “thesis art” is to inspire discussion with persons who do not have an academic background or work in a different field. The thesis is hidden in the drawer, but the poster is out there at the wall for everyone to see. You can find all past thesis art pieces here

Thesis Art: Edoardo Pinzuti

I was a co-supervisor for Edoardo Pinzuti’s Master’s Thesis. I finally came around to make this artwork with the text from his thesis.

He wrote an impressive matlab toolbox to analyze causality directions in time series based on Takens Theorem. The whole idea is about reconstructing embeddings of chaotic systems, with the Lorenz system (the one depicted in this artwork) being a simulation example in his thesis. Please find the DDIFTOOL toolbox here

The idea of “thesis art” is to inspire discussion with persons who do not have an academic background or work in a different field. The thesis is hidden in the drawer, but the poster is out there at the wall for everyone to see. You can find all past thesis art pieces here

Thesis Art: Judith Schepers

I was a supervisor for Judith Scheper’s Bachelor’s Thesis.

In this thesis-art, I visualized the guided-bubble paradigm used in a recent publication in the Journal of Vision. Judith generalized the paradigm to more than five bubbles, therefore, many more bubbles are visible in the thesis-art.

The idea of “thesis art” is to inspire discussion with persons who do not have an academic background or work in a different field. The thesis is hidden in the drawer, but the poster is out there at the wall for everyone to see. You can find all past thesis art pieces here

New Paper: The temporal dynamics of eye movements as an exploitation-exploration dilemma.

We just published a new paper in the Journal of Vision

The temporal dynamics of eye movements as an exploitation-exploration dilemma
Ehinger Kaufhold & König, 2018

The highlights:

  • We put eye movements as a decision process between exploitating the current view and exploring more of the scene
  • We use gaze-contingent eye-tracking to control the When and Where of eye movements
  • We find large effects of how long a subject fixates on their reaction time to continue exploring
  • We find large effects of the number of possible future target locations (Hick’s effect)

Check out the paper at the Journal of Vision

Ubuntu 16+: Recover ctrl+alt+bksp to restart X server

Often when developing with psychotoolbox or psychopy/opensesame your program crashes. And I often then have a full-screen window open and cannot click somewhere else. I then try to Alt+Tab and execute “sca” (screen close all) into the matlab console, with often mixed success. Sometimes restarting the computer is the last option. Instead of restarting, a useful command in older ubuntu versions was: STRG + ALT + Backspace => restart X server (=> restart GUI).

In order to activate this again use:

setxkbmap -option terminate:ctrl_alt_bksp

source on askubuntu



PS: I wrote this blogpost because I looked up this thing multiple times – now I know where to look 😉

Stretching the axes; visualizing non-linear linear regression


From time to time I explain concepts and ideas to my students.


Often this pops up in a statistical context, when one has a non-linear dependency between the to-be-predicted variable and the predictor-variables. By transforming the predictors, relationships can be made linear, i.e. a logarithmic (exponential, quadratic etc.) relationships can be modeled by a **linear** model.

The idea

I have a very visual understanding on basis-functions / non-linear transformation of variables in terms of stretching / condensing the basis (the x-axis here). This can also be applied to the generalized linear model (here for logistic regression).

Imagine that the x-axis of a plot is made of some kind of elastic material, you can stretch and condense it. Of course you do not need to stretch every part equally, one example would be to stretch parts that are far away from zero, exponentially more than parts that are close to zero. If you would have an exponential relationship ($ y = e^x$) then $y$ would now lie on a straight line.


Imagine you have a non-linear relationship, by stretching the x-axis in accordance to that non-linear relationship, you will have a linear relationship.

An exemplary non-linear relationship:

We want to do $y = b_0 + b_1x $ but obviously a linear line does not fit well. We can do something called polynomial expansion, i.e. add more predictors which are simple transformations of the predictor $x$. i.e. $y = b_0 + b_1x + b_2x^2 + b_3x^3$

The trick comes here: We can interpret the new $x^3$ basis function as a stretching of the x-axis. I.e. the further we move out on the x-axis, the longer we need to stretch the parts (exactly by x^3 times)

This can be shown also for other functions:



Note that the logarithm is not defined for negative numbers


Note how the stretching can be negative, i.e. the original negative values are stretching/transformed to positive values

Using the trick on the y-axis

One can interprete **logistic regression** with the same trick:
$$ g{-1}(y) = b_0 + b_1*x <=> y = g(b_0+b_1x)$$
with $g$, the logistic (logit) function and $g^{-1}$ the inverse logistic function (invlogit)
$$ g^{-1} = \ln\frac{p}{1-p} <=> g = \frac{1}{1+e^{-x}}$$

Usually we would have some non-linear dependency on a probabilty of e.g. success. That means, with a low value of x, your success-chance are low. To model this kind of data, one can transform the y-axis using $g$ above.

1 of 4