Tutorial to non-linear statistical analysis

If you happen to be interested in running non-linear statistical analysis, I have written a tutorial for running Generalized Additive Models in R. This type of model allows you to fit a dependent variable with a numeric predictor and shows you where the non-linear relationships are. You can find the tutorial either in my Ressources, or jump directly to my osf repository where you can download both, the tutorial and the necessary data.


Phonetic Characteristics of German a-schwa in Words like “Lehrer”

I am very happy to announce that Michael Ramscar and I got a paper published in Frontiers in Psychology. The topics are a) cue-outcome structure in modelling and b) phonetic characteristics of German “a-schwa” in different morphological functions. We have found that a phones morphological functions — in our case those of aschwa — and should be treated as inputs to a computational model of speech cognition, because they also serve as cognitive cues during speech production.

More over, we investigated a-schwa’s phonetic characteristics. We have found that these differ depending on the word’s class that the a-schwa is in. For example, the a-schwa at the end of the content word “Lehrer” has systematically different characteristics than that in the conjunction “aber”.
Crucially, we have found that the degree of experience between these morphological functions and the gesture of that phone — which speakers accumulate over their life time — modulates a-schwa’s phonetic characteristic.

This link leads you to the open source version of the paper.

Paper on vowel articulation

Articulatory trajectories of the tongue tip (left panel) and tongue body (right panel) during the production of American English [i] (top) and [A] (bottom) in verbs. Gray shading indicates frequency of occurrence within a verbal paradigm, with black = high frequency and light gray = low frequency. The onset of the articulatory trajectory is located at the thick end of the trajectory. Original picture here.

I am so happy to announce that we have successfully published a paper on the articulation of stem vowels in morphologically complex words in Morphology (“Paradigmatic enhancement of stem vowels in regular English inflected verb forms”). It has been a long journey. I wrote the first draft of this paper four years and three months ago. The data was recorded in 2013!

In our study, that I performed together with Michael Ramscar, Benjamin Tucker and Harald Baayen, we find that articulatory movements are partially more extreme in high-frequency words than in low-frequency words, supporting our previous findings in articulation and simultaneously contradicting predictions by the Smooth Signal Redundancy Hypothesis.

In our analysis, we used Generalized Additive Mixed Models, a regression technique that allows to investigate non-linear relationships between a dependent variable and one or multiple predictors. If you want to investigate such relationships on your own, the paper provides you with a good introduction to how assess the relative importance of potentially competing and collinear predictors. The supplementary material also provides you with the necessary data and code.

If you want to learn what these findings tell us about speech production, you can find the paper here. If you want to know more about speech production and articulation, check out my other publications here.

Good Scientific Writing is Writing for Readers

A while ago I gave some tips and tricks how to start a scientific manuscript (click here). The idea was to record yourself while you are talking about your research, transcribe the recording and then edit the transcription.

I have not tackled, however, how to actually write and structure your text, build an argument, tell a story. And you know why? Because it is such a complicated process that I do not feel qualified enough to talk about it. At least not yet.

Yes, a scientific text is a story. Scientists are writers and well written papers are key to letting the world knowing about your world. Thanks to a colleague of mine, I discovered the teaching of Judy Swan – Associate Director for Writing in Science and Engineering. One key stone of her career is to communicate about scientific writing. Her key insight is that texts do not have a fixed interpretation. “Words do not have meanings, they have interpretations”, she states. The interpretation of words, sentences, entire texts depends on the contexts. Therefore, once a texts gets published, interpretations depends on the readers and how much they know and how much they understand what we write.

Good scientific writing manages the readers’ expectations. It begins at the level of the paper’s structure – introduction, methods, results, discussion – but it carries on at the level of paragraphs and sentences sentences. In a well written scientific text, important pieces of information are placed in locations where readers expect emphasis on important pieces of information.

Enjoy the video.

Does morphology affect phonetic characteristics?

Yes, is the short answer. The long answer: it depends on how well you learned the connection between a morphological function and  the word, the sound and the context in which the word is located.

Long story short: the better you learned this connection, and consequently the more sure you are the message you want to convey, the longer you make the phonetic signal.

Together with my colleagues Mirjam Ernestus, Ingo Plag and Harald Baayen we were able to show this relation in word-final [s] or [z] in American English and the morphological function it expresses, e.g. plural in ‘dogs’ or third person singular in ‘takes’.

To simulate learning with a computer, we used Naive Discriminative Learning (NDL), a computational formalization of Discriminative Learning. Discriminative learning assumes that speakers learn to discriminate their environment by language. Listeners on the other hand learn to use the perceived signal to predict the message the speaker intended. The theory predicts that the amount of experience speakers and listeners have with the relation between message and signal will affect their behavior. In speakers, it is the behavior is the production of the speech signal which encodes the message; in listeners, the behavior is how they respond to the signal. Here, we studied the speakers.

We operationalized the intended message as the morphological function a word has to convey and the cues as the intended acoustic signal. NDL thus learned to discriminate the morphological function on the basis of the acoustic cues. Crucially, our model used cues from the target word such as ‘dogs’ or ‘takes’, but also from words in the context of those words. The reason for this is that we assume that a word’s morphological function is not only encoded by the acoustics of the word itself, but by the acoustics of the words surrounding it.

NDL further allows to calculate measures about 1) how well the morphological function is supported by a set of cues and 2) how uncertain the model is about the lexical information. Take a driving on the highway as an example. Let’s say, you are waiting for an exit to leave the highway and finally a sign occurs. The larger an exit sign on the street, the more support you have that an exit will occur. However, the size of the sign does not say anything about how certain (or uncertain) you actually are that it is actually the exit you should take. Both will affect how you will behave on the street.

To sum it up, NDL allows you to simulate how morphological function and the phonetic signal are connected and therefore to investigate how the processes guiding speech production affect the phonetic signal

The Figure shows how [s] duration is affected by the interaction between bottom-up support (x-axis) and uncertainty (y-axis) about a word’s morphological function. Blue regions represent short durations, yellow regions represent long durations.

The entire study will be published under the title:

Phonetic effects of morphology and context: Modeling the duration of word-final S in English with NaÏve Discriminative Learning

It has been accepted for publication by the “Journal of Linguistics” and is being prepared for publication. A pre-publication version can be downloaded from psyarxive here.

The paper also provides a very good introduction to Discriminative Learning, how it is performed and how it can be used for predicting phonetic characteristics. If you want to perform an NDL learning simulation on your own, you can find an introduction to this technique in my R introduction here.

How do appearances of female characters in TV shows affect viewer ratings?



Of course, your first reaction is: Come on, that show is like a hundred years old. It is from the last century, people were different back then in the dark ages. My response to anyone who really thinks this at this moment: Really? Are we really that different in light of the current political situation (sorry, different topic) and in light of the reactions to a woman taking over as the main character of a Science-Fiction/Phantasy/Time-travel/The main character can change into every shape possible? Are we really that different? Then prove me wrong. I will continue here with Star Trek: The Next Generation (TNG). And by the way, this happens to be the dataset I have available. If you can provide any other, more “modern” dataset, then I will analyze it (but only if you preprocess the data).

The current material

I collected the data in 2015 (yes, it’s been a while, but science takes time). I downloaded the transcripts for all the TNG episodes from http://www.chakoteya.net/ (many thanks to whoever transcribed the episodes), which, I’m just realizing, also has transcripts for Doctor Who (any volunteers?). I could have taken the average ratings for each episode, but I wanted to have a more fine grained data. This is why I contacted imbd.com who provided me with information on the sex of the rater and their age for all the episodes of all Star Trek franchises (thanks again!).

I decided to calculate the average rating per episode depending on the sex of the rater, as I was interested whether women would rate episodes with a higher proportion of women dialogues better than man (any guesses?). This gives me two ratings per episode: one for female raters and one for male raters. For the analysis, I transformed imdb’s 1 to 10 scale to range between 0 and 1.

I am focusing here on the dialogues of the main characters, which are Wesley, Worf, Dr. Crusher, Dr. Pulaski, Laforge, Troi, Riker, Data and Picard. I ignored all the guest characters, side kicks, Klingons, Vulcans and who ever appears on the show, mainly because it was too complicated to tag each of these characters according to their sex.

I calculated a simple measure: the main characters’ women-to-men ratio in terms of the number of words spoken in each episode. I am ignoring the name of the character (everyone loves Data) and, most importantly, I am ignoring the plot. This is important, because only in this way we can assess whether viewers have a certain preference to the sex of the characters. I would expect that female raters favor episodes in which more female characters appear.

This gives us 176 values for 176 episodes. The larger the value, the more text (measured by the number of words) the main female characters have in that episode. A ratio of 1 would mean, men and woman have an equal amount of text. Above 1, women have more text, below one, men have more text.

The following Figure  illustrates the distribution of the women-to-men ratio in all episodes. It becomes very apparent that men have more text than women in TNG. There are only 9 episodes where men have less text. On average, women have roughly 34% of the text that men have (median of 22%). For analysis, I am ignoring those episodes that have more text, because excluding them makes the data better distributed (dont’t worry, I have also performed the analysis including the data and it does not change a thing).

I have used the betareg library that allows us to model values that range between 0 and 1 (the betareg library is actually designed to model ratios by means of beta-regression, that transforms probabilities ranging between 0 and 1 into logits. Hence, all statistical results presented below (i.e. beta = X), are logits. If you want to know more about logits, click here). The summary of beta-regression gives us coefficients (beta estimates), standard errors (sde), z-values (larger than 2 represents significant) and p-values (smaller than 0.05 represents significant).

I fitted the ratings with three predictors in one model.

The first predictor was the episode number in order to investigate how ratings evolved across the seasons. Indeed, the effect was significantly positive (beta = 0.028, standard error = 0.0005, z = 5.7, p < 0.001). Indeed, TNG’s ratings increased with every season. This effect is illustrated in the next plot. Although there is a strong variation within each season, the ratings got better and better. Nicely done, Star Trek producers.


Unsurprisingly, male raters gave TNG higher ratings then female raters: on a scale between 0 and 1, female raters gave on average 0.66 (beta = 0.86, sde = 0.07, z =12.7, p < 0.001), male raters gave on average 0.71 (beta = 0.23, sde = 0.8, z = 3, p = 0.003).

Now, what about the women-to-men-text ratio? Well, the effect is highly negative. The more text woman had in an episode, the worse the episode was rated (beta = -0.86, sde = 0.17, z = -5.1, p < 0.001). The following plots illustrates this effect.

Two insights follow from the figure. First, the distribution of the women-to-men-text-ratio is skewed. This that there are less episodes with more text for women. Second, that the variability in ratings is really large for episodes with a high percentage of text for male characters. One potential interpretation of this finding is that since there are fewer female-oriented eposides in the series, there is also a smaller probability that there will be a good episode episodes with higher women-to-men-text-ratios.

What is more important: female and male raters did not differ in their ratings. This means that female viewers have the same opinions about how often female characters should occur on Television like male viewers. In my opinion this result is devastating. Not only does it mean that viewers do not accept female characters to be present in Television. It also means that this opinion is supported by those who are represented by these characters: The women themselves. Given that today’s television has a strong influence on public opinion and on how each person defines their role in this society, this is absolutely unacceptable. Of course female viewers rate episodes with more female appearances worse then with more male appearances. This is what they have learned to be the status quo by watching Television.

Coming back to Doctor Who (have not forgotten this one). I claim that the reason why the latest season of Doctor Who has significantly worse ratings than all those before is simple: The main character is portrait by a woman, the plots center around female characters. For example, the episode “Rosa” that focuses on Rosa Parks. In the episode “The Tsuranga Conundrum” there is a female General. And so forth.

I totally acknowledge that the current finding has to be regarded with caution because Star Trek targets a male audience. Not only it should be replicated with Doctor Who (that started all this idea) but also with TV shows that target both sexes. Maybe someone will provide this data to me.

Statistical investigation of Mass Shootings in the US in 2018

After the recent Mass Shooting in Thousand Oaks, California, together with some colleagues we were discussion potential reasons for these outbreaks of violence in the US. We had a look at the Wikipedia list on Mass Shootings in 2018. We were shocked to see that there as many as 106 shootings noted. On abc15.com, we found a map illustrating the location of every mass shooting in 2018:

It becomes immediately apparent that most of the shootings happened in the Eastern part of the US. Given that there is already structure in the data, we were wondering whether there is more.
For example, is unemployment one reason for the shooting? Or population density? Concretely, the question arose whether the number of casualties can be predicted by such information. I expected that the number of casualties would increase, the larger the size of the population, the larger the population growth and the larger unemployment.


On the basis of the Wikipedia list on Mass Shootings in 2018, I collected information about the city/town/village the mass shooting occurred. The lists contains information on 106 mass shootings. I used Wikipedia, because it allowed me to obtain easily additional information on the locality of the mass shooting by following the links on the page. I collected the following pieces of information:

  • Population of the locality (city) in 2010 and 2016 (as provided by Wikipedia)
  • Population of the state in 2010 and 2016
  • Unemployment rate of the state

The 2010 data was based on the 2010 United States Census. I would have wanted to have information on the Unemployment rate of the locality. However, this surpassed my abilities. I also calculated the percentage of growth in population size between 2010 and 2016 for both, locality and state. In total, I used five variables to predict the the number of casualties (NumberOfCasualties).

  • PopulationCity (in 2016, ranging from 3365 to 8622698)
  • PopulationState (in 2016, ranging from 601723 to 37253956 )
  • PopulationGrowthCity (between 2010 and 2016, ranging from 0.875 to 1.315)
  • PopulationGrowthState (between 2010 and 2016, ranging from 0.994 to 1.35)
  • UnemploymentRateState (in 2016)

For 9 of the localities, population size was available only for 2016. Those were excluded from the analysis. In pilot analyses, the inclusion of those 9 localities did not change the overall results of the present analysis. The shootings in the Bronx and in Brooklyn were tagged as shootings in New York City.

All of those pieces information were collected automatically from Wikipedia and a data tables using a custom-made script in R. The download of the html-pages was performed with the function getURL() from the RCurl package.

Analysis and Results

City population size in 2016 had to be log transformed in order to obtain normal distribution. All predictors were centered and scaled for analysis (z-scaled). I first performed a standard Spearman-Rank correlation analysis between the predictors. I want to highlight here two results:

  • PopulationGrowthState and PopulationGrowthCity were strongly correlated (R = 0.56). This is not surprising at all, as the state population depends on the city population.
  • PopulationGrowthState was negatively correlated with UnemploymentRateState (R= – 0.48). The same effect, weaker, was for PopulationGrowthCity (R=-0.32). This means that when the population size grew, the unemployment rate went down. While this might not be nothing new for people working in demographics and economics, this surprised me, as I would have expected these two variables to be positively correlated. I would like to see, whether this was observed on a larger scale.

NumberOfCasualties ranged between 0 (N = 40) and 17 (N = 1).

I fit a generalized linear model to predict the NumberOfCasualties (function glm, family = poisson). The analysis is quite trick, because of the high collinearity in the data, i.e. some predictor can be used to predict another predictor. This becomes apparent above. This is why I used a step-by-step inclusion procedure and checked for indicators of collinearity in the model. If you want to read more about how to address collinearity in regression analyses, together with two other colleagues I have published a paper here. I have not tested any interactions between the predictors because of the small sample number in the data.

As it turned out, UnemploymentRateState and PopulationGrowthState were not significantly predictive for NumberOfCasualties (this means that the effect they cause cannot be used to support our initial hypothesis). I found this surprising, as I indeed would have predicted that unemployment drives people to commit horrible things.
The three remaining predictors were significant. The model’s intercept was 0.271 (std = 0.09, z = 2.9, p = 0.003). This is the estimated logit, where 0 equals 50% (in R, you can transform those values back using the function inv.logit() from the boot package. I did not apply this because it changes the visualization)

  • PopulationState (estimate = 0.2, std = 0.07, z = 2.9, p = 0.004)
  • PopulationCity (estimate = -0.4, std = 0.08, z = -4.8, p < 0.001)
  • PopulationGrowthCity (estimate = 0.3, std = 0.07, z = 4.4, p < 0.001)

The effects are illustrated in the following figure (illustrated with the function visreg() from the package visreg. I restricted the y-axis which is why not all data points are illustrated). The y-axis represents partial effect of the estimated logit, i.e. how the intercept (logit 0.271 = probability of 0.57 that a large number of people are killed) has to be changed depending on the predictor. 0 therefore represents no change (horizontal gray dashed line), -2 means that the intercept has to lowered to logit -1.729 = 0.15 probability; +2 means that the intercept has to be increased to logit 2.271 = 0.90 probability of a high number of casualties.


Panel (A) represents the population size of the state the mass shooting occurred. The larger the state, the larger the probability of a high number of casualties. This supports my hypothesis. Surprisingly and contrary to my hypothesis, the effect is reversed for the population size of the village/town/city where the shooting occurred: the smaller the city, the larger the probability of a higher number of casualties; the larger the city, the smaller the probability of a high number of casualties. Finally, the probability of a high number of casualties increases, when the city population has strongly grown between 2010 and 2016. This makes sense as: when cities grow, the total state population increases.

I will avoid any further interpretation of the data and the analysis. I also have not included other, potentially very interesting, variables into this analysis, in order to keep things as simple as possible.