Books on data visualization

Here is a compilation of new and classic books on data visualization:

 

Scott Murray (2017) Interactive Data Visualization for the Web 

Elijah Meeks (2017) D3.Js in Action: Data Visualization with JavaScript 

Alberto Cairo (2016) The Truthful Art: Data, Charts, and Maps for Communication 

Andy Kirk (2016) Data Visualization 

David McCandless (2014) Knowledge is Beautiful 

 

Edward Tufte (2006) Beautiful Evidence  

Edward Tufte (2001) The Visual Display of Quantitative Information 

Edward Tufte (1997) Visual Explanations: Images and Quantities, Evidence and Narrative 

Edward Tufte (1990) Envisioning Information 

Olympic medals, economic power and population size

The 2016 Rio Olympic games being officially over, we can obsess as much as we like with the final medal table, without the distraction of having to actually watch any sports. One of the basic questions to ponder about the medal table is to what extent Olympic glory is determined by the wealth, economic power and population size of the countries.

Many news outlets quickly calculated the ratios of the 2016 medal count with economic power and population size per country and presented the rankings of ‘medals won per billion of GDP’ and ‘medals won per million of population’ (for example here and here). But while these rankings are fun, they give us little idea about the relationships between economic power and population size, on the one hand, and Olympic success, on the other. Obviously, there are no deterministic links, but there could still be systematic relationships. So let’s see.

Data

I pulled from the Internet the total number of medals won at the 2016 Olympic games and assigned each country a score in the following way: each country got 5 points for a gold medal, 3 points for silver, and 1 point for bronze. (Different transformations of medals into points are of course possible.) To measure wealth and economic power, I got the GDP (at purchasing power parity) estimates for 2015 provided by the International Monetary Fund, complemented by data from the CIA Factbook (both sets of numbers available here). For population size, I used the Wikipedia list available here.

Olympic medals and economic power

The plot below shows how the total medal points (Y-axis) vary with GDP (X-axis). Each country is represented by a dot (ok, by a snowflake), and some countries are labeled. Clearly, and not very surprisingly, countries with higher GDP have won more medals in Rio. What is surprising however, is that the relationship is not too far from linear: the red line added to the plot is the OLS regression line, and it turns out that this line summarizes the relationship as well (or as badly) as other, more flexible alternatives (like the loess line shown on the plot in grey). The estimated linear positive relationship implies that, on average, each 1,000 billion of GDP bring about 16 more medal points (so ~315 billion earns you another gold medal).olymp1

The other thing to note from the plot is that the relationship is between medal points and total GDP, thus not GDP per capita. In fact, GDP per capita, which measures the relative wealth of a country, has a much weaker relationship with Olympic success with a number of very wealthy, and mostly very small, countries getting zero medals. The correlation of Olympic medal points with GDP is 0.80, while with GDP per capita is 0.21. So it is absolute and not relative wealth that matters more for Olympic glory. This would seem to make sense as it is not money but people who compete at the games, and you need a large pool of contenders to have a chance. But let’s examine more closely whether and how does population size matter.

Olympic medals and population size

The following plot shows how the number of 2016 Rio medal points earned by each country varies with population size. Overall, the relationship is positive, but it is not quite linear, and it is not very consistent (the correlation is 0.40). Some very populous countries, like India, Indonesia, and Pakistan have won very few medals, and some very small ones have won at least one. The implied effect of population size is also small in substantive terms: each 10 million people are associated with 1 more medal point (so, a bronze); for reference three quarters of the countries in the dataset have less than 25 million inhabitants.

olymp2

Putting everything together

Now, we can put both GDP and population size in the same statistical model with the aim of summarizing the observed distribution of medal points as best as we can. In addition to these two predictors, we can add an interaction between the two, as well as different non-linear transformations of the individual predictors. In fact, the possibilities for modeling are quite a few even with only two predictors, so we have to pick a standard for selecting the best model. As the goal is to describe the distribution of medal points, it makes sense to use the sum of the errors (the absolute values of the differences between the actual and predicted medal score for each country) that the models make as a benchmark.

I find that two models describe the data almost equally well. Both use simple OLS linear regression. The first one features population size, GDP, and GDP squared. In this multivariate model, population size turns out to have a negative relationship with Olympic success, net of economic power. GDP has a positive relationship, but the quadratic term implies that the effect is not truly linear but declines in magnitude with higher values of GDP. The substantive interpretation of this model is something along these lines: Olympic success increases at a slightly declining rate with the economic power of a country, but given a certain level of economic power, less populous countries do better. The sum of errors of Model 1 is 1691 medal points.

The second model is similar, but instead of the squared term for GDP it features an interaction between GDP and population size. The interaction turns out to be negative. This implies that economically powerful but populous countries do less well than their level of GDP alone would suggest. This interpretation is a bit strange as population size is positively associated with GDP and seems to suggest that it is relative wealth (GDP per capita) that matters, but this turns out not to be the case, as any model that features GDP per capita has a bigger sum of errors than either Model 1 or Model 2.

Model 1 Model 2
Population size – 0.20 – 0.09
GDP + 0.04 + 0.03
GDP squared – 0.00000008 /
GDP*Population / -0.0000008
Sum of errors 1691 1678
Adjusted R-squared 0.83 0.81

Both models presented so far are linear which is not entirely appropriate given that the outcome variable – medal points – is constrained to be non-negative and is not normally distributed. The models actually predict that some countries, like Kenya, should get a negative number of medal points, which is clearly impossible. To remedy that, we can use statistical models specifically developed for non-negative (count) data: Poisson, negative binomial, or even hurdle or zero-inflated models that can account for the excess number of countries with no medal points at all. I spend a good deal of time experimenting with these models, but I didn’t find any that improved at all on the simple linear models described above (it is actually quite hard even evaluating the performance of these non-linear models). Let me know if you find a different model that does better than the ones reported here. (But please no geographical dummies or past Olympic performance measures; also, the Olympic delegation size would be a mediator so not a proper predictor).

The one model I can find that outperforms the simple OLS regressions is a generalized additive model (GAM) with a flexible form for the interaction. This model has a sum of errors of 1485, and the interaction surface looks like this:interactionGDPpop

In conclusion, do the population size, economic power and wealth of countries account for their success at the 2016 Olympic games? Yes, to a large extent. It is economic power and not relative wealth that matters more, and population size actually has a negative effect once economic power is taken into account. So the relationships are rather complex and, to remind, far from deterministic.

 

Here is the data (text file): olypm. Let me know if you interested in the R script for the analysis, and I will post it.
Finally, here is a ranking of the countries by the size of the model error (based on Model 2; negative predictions have been replaced with zero). This can be interpreted in the following way: the best way to summarize the distribution of medal points won at the 2016 Rio Olympic games as a function of population size and GDP is the model described above. This model implies a prediction for each country. The ones that outperform their model predictions have achieved more than their level of GDP and economic size imply. The ones with negative errors underperform in the sense that they have achieved less than their level of GDP and economic size imply.

country 2016 medals 2016 medal points predicted medal points model error
Great Britain 67 221 68 153
Russia 56 168 87 81
Australia 29 83 30 53
France 42 118 68 50
Kenya 13 49 0 49
New Zealand 18 52 4 48
Hungary 15 53 6 47
Netherlands 19 65 22 43
Jamaica 11 41 0 41
Croatia 10 36 2 34
Cuba 11 35 2 33
Azerbaijan 18 36 4 32
Germany 42 130 98 32
Uzbekistan 13 33 2 31
Italy 28 84 54 30
Kazakhstan 17 39 10 29
Denmark 15 35 7 28
Ukraine 11 29 5 24
Serbia 8 24 2 22
North Korea 7 21 0 21
Sweden 11 31 12 19
Belarus 9 21 4 17
Ethiopia 8 16 0 16
Georgia 7 17 1 16
South Korea 21 63 47 16
China 70 210 195 15
South Africa 10 30 15 15
Armenia 4 14 0 14
Greece 6 20 7 13
Slovakia 4 16 4 12
Spain 17 53 41 12
Colombia 8 24 14 10
Czech Republic 10 18 8 10
Slovenia 4 12 2 10
Switzerland 7 23 13 10
Bahamas 2 6 0 6
Bahrain 2 8 2 6
Ivory Coast 2 6 0 6
Belgium 6 18 13 5
Fiji 1 5 0 5
Kosovo 1 5 0 5
Tajikistan 1 5 0 5
Lithuania 4 6 2 4
Burundi 1 3 0 3
Grenada 1 3 0 3
Jordan 1 5 2 3
Mongolia 2 4 1 3
Niger 1 3 0 3
Puerto Rico 1 5 2 3
Bulgaria 3 5 3 2
Canada 22 44 43 1
Moldova 1 1 0 1
Romania 5 11 10 1
Vietnam 2 8 7 1
Afghanistan 0 0 0 0
American Samoa 0 0 0 0
Andorra 0 0 0 0
Antigua and Barbuda 0 0 0 0
Aruba 0 0 0 0
Barbados 0 0 0 0
Belize 0 0 0 0
Benin 0 0 0 0
Bermuda 0 0 0 0
Bhutan 0 0 0 0
British Virgin Islands 0 0 0 0
Burkina Faso 0 0 0 0
Cambodia 0 0 0 0
Cameroon 0 0 0 0
Cape Verde 0 0 0 0
Cayman slands 0 0 0 0
Central African Republic 0 0 0 0
Chad 0 0 0 0
Comoros 0 0 0 0
Congo 0 0 0 0
Cook Islands 0 0 0 0
Djibouti 0 0 0 0
Dominica 0 0 0 0
DR Congo 0 0 0 0
Eritrea 0 0 0 0
Estonia 1 1 1 0
Gambia 0 0 0 0
Guam 0 0 0 0
Guinea 0 0 0 0
Guinea-Bissau 0 0 0 0
Guyana 0 0 0 0
Haiti 0 0 0 0
Honduras 0 0 0 0
Iceland 0 0 0 0
Kiribati 0 0 0 0
Kyrgyzstan 0 0 0 0
Laos 0 0 0 0
Lesotho 0 0 0 0
Liberia 0 0 0 0
Liechtenstein 0 0 0 0
Madagascar 0 0 0 0
Malawi 0 0 0 0
Maldives 0 0 0 0
Mali 0 0 0 0
Malta 0 0 0 0
Marshall Islands 0 0 0 0
Mauritania 0 0 0 0
Micronesia 0 0 0 0
Monaco 0 0 0 0
Montenegro 0 0 0 0
Mozambique 0 0 0 0
Nauru 0 0 0 0
Nepal 0 0 0 0
Nicaragua 0 0 0 0
Palau 0 0 0 0
Palestine 0 0 0 0
Papua New Guinea 0 0 0 0
Poland 11 25 25 0
Rwanda 0 0 0 0
Saint Kitts and Nevis 0 0 0 0
Saint Lucia 0 0 0 0
Samoa 0 0 0 0
San Marino 0 0 0 0
Sao Tome and Principe 0 0 0 0
Senegal 0 0 0 0
Seychelles 0 0 0 0
Sierra Leone 0 0 0 0
Solomon Islands 0 0 0 0
Somalia 0 0 0 0
South Sudan 0 0 0 0
St Vincent and the Grenadines 0 0 0 0
Suriname 0 0 0 0
Swaziland 0 0 0 0
Tanzania 0 0 0 0
Timor-Leste 0 0 0 0
Togo 0 0 0 0
Tonga 0 0 0 0
Trinidad and Tobago 1 1 1 0
Tunisia 3 3 3 0
Tuvalu 0 0 0 0
Uganda 0 0 0 0
US Virgin Islands 0 0 0 0
Vanuatu 0 0 0 0
Yemen 0 0 0 0
Zambia 0 0 0 0
Zimbabwe 0 0 0 0
Albania 0 0 1 -1
Bangladesh 0 0 1 -1
Bolivia 0 0 1 -1
Bosnia and Herzegovina 0 0 1 -1
Botswana 0 0 1 -1
Brunei 0 0 1 -1
Cyprus 0 0 1 -1
El Salvador 0 0 1 -1
Equatorial Guinea 0 0 1 -1
FYR Macedonia 0 0 1 -1
Gabon 0 0 1 -1
Ghana 0 0 1 -1
Ireland 2 6 7 -1
Latvia 0 0 1 -1
Mauritius 0 0 1 -1
Namibia 0 0 1 -1
Paraguay 0 0 1 -1
Sudan 0 0 1 -1
Syria 0 0 1 -1
Costa Rica 0 0 2 -2
Dominican Rep. 1 1 3 -2
Guatemala 0 0 2 -2
Libya 0 0 2 -2
Luxembourg 0 0 2 -2
Panama 0 0 2 -2
Turkmenistan 0 0 2 -2
Uruguay 0 0 2 -2
Angola 0 0 3 -3
Lebanon 0 0 3 -3
Myanmar 0 0 3 -3
Ecuador 0 0 4 -4
Morocco 1 1 5 -4
Sri Lanka 0 0 4 -4
Argentina 4 18 23 -5
Finland 1 1 6 -5
Israel 2 2 7 -5
Oman 0 0 5 -5
Qatar 1 3 8 -5
Thailand 6 18 23 -5
Norway 4 4 10 -6
Portugal 1 1 7 -6
Algeria 2 6 13 -7
Brazil 19 59 66 -7
Malaysia 5 13 20 -7
Venezuela 3 5 12 -7
Iran 8 22 30 -8
Pakistan 0 0 8 -8
Peru 0 0 8 -8
Philippines 1 3 11 -8
Singapore 1 5 13 -8
Austria 1 1 11 -10
Chile 0 0 10 -10
Hong Kong 0 0 11 -11
Nigeria 1 1 13 -12
India 2 4 17 -13
Iraq 0 0 13 -13
Japan 41 105 119 -14
U.A.E. 1 1 18 -17
Egypt 3 3 21 -18
Turkey 8 18 37 -19
Chinese Taipei 3 7 29 -22
Mexico 5 11 49 -38
Indonesia 3 11 51 -40
Saudi Arabia 0 0 44 -44
United States 121 379 431 -52

Visualizing asylum statistics

Note: of potential interest to R users for the dynamic Google chart generated via googleVis in R and discussed towards the end of the post. Here you can go directly to the graph.

02alessandro-penso
An emergency refugee center, opened in September 2013 in an abandoned school in Sofia, Bulgaria. Photo by Alessandro Penso, Italy, OnOff Picture. First prize at World Press Photo 2013 in the category General News (Single).

The tragic lives of asylum-seekers make for moving stories and powerful photos. When individual tragedies are aggregated into abstract statistics, the message gets harder to sell. Yet, statistics are arguably more relevant for policy and provide for a deeper understanding, if not as much empathy, than individual stories. In this post, I will offer a few graphs that present some of the major trends and patterns in the numbers of asylum applications and asylum recognition rates in Europe over the last twelve years. I focus on two issues: which European countries take the brunt of the asylum flows, and the link between the application share that each country gets and its asylum recognition rate.

Asylum applications and recognition rates
Before delving into the details, let’s look at the big picture first. Each year between 2001 and 2012, 370,000 people on average have applied for asylum protection in one of the member states of the European Union (plus Norway and Switzerland). As can be seen from Figure 1, the number fluctuates between 250,000 and 500,000 per year, and there is no clear trend. Altogether, during this 12-year period, approximately 4.5 million people have applied for asylum, which makes slightly less than one percent of the total EU population. Of course, this figure only tracks people who have actually made it to the asylum centers and filed an application – all potential refugees who have perished on the way, or have arrived but been denied the right of formal application, or have remained clandestine are not counted.

asylum_applications_small

Figure 1 also shows the annual number of persons actually recognized as ‘refugees’ under the terms of the Geneva Convention by the European governments: a status which grants considerable rights and protection. This number is quite lower with an average of around 40.000 per year (in the EU+ as a whole) which makes for less than half-a-million in total for the 12 years between 2001 and 2012. While the overall recognition rate remains between 7% and 14%, there is considerable variation between the different European states both in the share from the asylum flows they receive, and in the national asylum recognition rates.

Who takes the brunt of the asylum burden?
Both the asylum flows and the recognition rates are in fact distributed highly unequally across the continent, and in a way that cannot be completely accounted for by the wealth of destination countries, former (colonial) ties between asylum sources and destinations, nor geographical distance. To compare the shares of the total European pool of asylum applications and recognitions that a destination country gets, I create the so-called ‘burden coefficient’. The ‘burden coefficient’ compares the actual share of asylum applications a country received in a year to its ‘fair’ share which is defined as its relative share of the annual  total EU+ GDP. Simply put, if a country accounts for 10% of the European GDP, it would have been expected to receive 10% of all asylum applications filed in Europe that year. Taking account of GDP adjusts the raw asylum application shares in view of the expectation that richer and more populous countries should bear a proportionally higher share of the total European asylum ‘burden’ than poorer and smaller states.

asylum_applications_burden

Figure 2 shows the (logged) burden coefficient for asylum application shares for each EU+ country, averaged over the period 2010-2012. The solid line at zero indicates an asylum applications share perfectly proportional to a  country’s GDP share (a ‘fair’ burden). Countries with positive values receive a higher share of all applications than implied by their GDP level, and countries with negative values receive a lower than their implied share. (The dotted lines show where a country that is doing twice as much / twice as little as expected would be). Clearly, Spain, Portugal, Italy and many (but not all) of the East European countries underdeliver while Cyprus, Malta, Greece, and several West European states (notably Sweden, Belgium, and Norway) take a disproportionately high  share of the total pool of asylum applications filed in Europe over the last few years. Note that these comparisons already take into account (correct for) the fact that most of the Southern and Eastern European countries are poorer (have lower GDP) than the ones in the Western and Northern parts of the continent.

asylum_recognitions_burden

The picture does not change much when we focus on actual asylum recognitions (under the terms of the Geneva Convention) instead of applications. Figure 3 shows the burden coefficient (again averaged over 2010-2012) for full status refugee recognitions in Europe. The country ranking is similar with a few important exception – Greece grants much fewer asylum recognitions than expected even after we account for the state of its economy; Austria and Switzerland join the ranks of states which do much more than their implied share; and, sadly, many more countries in fact underdeliver when it comes to full refugee status grants. (Note that some states offer alternative protection to those denied the full ‘Geneva Convention’ status but the forms and level of this protection differs significantly across the continent).

Are asylum application shares responsive to the recognition rate?
Given these rather significant discrepancies across Europe in how many asylum applications countries get, and how much protection they offer, it is natural to ask whether the applications shares and the recognition rates are in fact related. Do asylum seekers flock at the gates of the European states which are most generous in their recognition policy? Do low recognition rates deter potential refugees from applying in certain countries? Can the strictness of asylum policy be an effective policy tool shaping future application flows? A comprehensive statistical analysis shows that while application shares and recognition rates are associated, their responsiveness to each other is rather weak. Simply put, manipulating the recognition rates is unlikely to have big practical effects on the asylum application share a country receives, and changes in the applications rates only weakly affect state recognition rates. The details of the analysis are rather technical and can be found here, but a dynamic visualization can help illustrate the patterns.

The dynamic interactive chart linked here shows the relationship between asylum applications and asylum recognition rates for each EU+ country over the last 12 years (the chart cannot be embedded in this post due to WordPress policy, but there is a screenshot below). When you press ‘Play’ each dot traces the experience of one country over time. You can choose to observe all, select a single state to focus upon, or tick a couple to compare their experiences.

dynamic-asylum-1

A movement of a dot (and the trace in leaves) in a horizontal direction means that the number of asylum applications received by a country increases while the recognition rates remains the same. Similarly, a vertical move implies a change in the recognition rate but a stable asylum application flow. A trajectory that follows a diagonal suggests a link between applications and recognition rates.

When paused, the state of the chart at each year shows the cross-sectional association between applications and recognition rates: it is easy to see that there is a (rather stable) weakly-strong positive relationship. But the trajectories of individual countries over time do not suggest that there is a temporal link between the two aspects of asylum policy for particular countries. For example, in the UK between 2001 and 2004 both the recognition rates and the applications fall, which would suggest strong responsiveness, but then the recognition rate moves up from 4% to almost 30% without any significant increase in applications. The trajectory of Denmark (try it out) exhibits something close to a dynamic link with rates depressing applications initially but then when they rise again, applications seem to pick up as well. Of course, asylum flows are driven by many other factors as well, so while suggestive, the patterns in the chart should be interpreted with care.

dynamic-asylum-2

More comprehensive analyses of asylum policy in Europe addressing these questions and more are available in my published articles accessible here and here. The original data comes from the UNHCR annual reports. The dynamic chart is generated using Google Chart Tools through the googleVis library in R, you can find the code here. I found it useful to generate a simple version, adjust the settings manually, and then copy the final settings via the Google Chart’s Advanced Panel back to R.

Predicting movie ratings with IMDb data and R

It’s Oscars season again so why not explore how predictable (my) movie tastes are. This has literally been a million dollar problem and obviously I am not gonna solve it here, but it’s fun and slightly educational to do some number crunching, so why not. Below, I will proceed from a simple linear regression to a generalized additive model to an ordered logistic regression analysis. And I will illustrate the results with nice plots along the way. Of course, all done in R (you can get the script here).

Data
The data for this little project comes from the IMDb website and, in particular, from my personal ratings of 442 titles recorded there. IMDb keeps the movies you have rated in a nice little table which includes information on the movie title, director, duration, year of release, genre, IMDb rating, and a few other less interesting variables. Conveniently, you can export the data directly as a csv file.

Outcome variable
The outcome variable that I want to predict is my personal movie rating. IMDb lets you score movies with one to ten stars. Half-points and other fractions are not allowed. It is a tricky variable to work with. It is obviously not a continuous one; at the same time ten ordered categories are a bit too many to treat as a regular categorical variable. Figure 1 plots the frequency distribution (black bars) and density (red area) of my ratings and the density of the IMDb scores (in blue) for the 442 observations in the data.

figure1

The mean of my ratings is a good 0.9 points lower than the IMDb scores, which are also less dispersed and have a higher peak (can you say ‘kurtosis’).

Data-generating process
Some reflection on how the data is generated can highlight its potential shortcomings. First, life is short and I try not to waste my time watching bad movies. Second, even if I get fooled to start watching a bad movie, usually I would not bother rating it on IMDb.There are occasional two- and three-star scores, but these are usually movies that were terrible and annoyed me for some reason or another (like, for example, getting a Cannes award or featuring Bill Murray). The data-generating process leads to a selection bias with two important implications. First, the effective range of variation of both the outcome and the main predictor variables is restricted, giving the models less information to work with. Second, because movies with a decent IMDb ratings which I disliked have a lower chance of being recorded in the dataset, the relationship we find in the sample will overestimate the real link between my ratings and the IMDb ones.

Take one: linear regression
Enough preliminaries, let’s get to business. An ordinary linear regression model is a common starting point for analysis and its results can serve as a baseline. Here are the estimates that lm provides for regressing my ratings on IMDb scores:

summary(lm(mine~imdb, data=d))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.6387     0.6669  -0.958    0.339    
imdb          0.9686     0.0884  10.957   ***
---
Residual standard error: 1.254 on 420 degrees of freedom
Multiple R-squared: 0.2223,	Adjusted R-squared: 0.2205

The intercept indicates that on average my ratings are more than half a point lower. The positive coefficient of IMDb score is positive and very close to one which implies that one point higher (lower) IMDb rating would predict, on average, one point higher (lower) personal rating. Figure 2 plots the relationship between the two variables (for an interactive version of the scatter plot, click here):

figure2

The solid black line is the regression fit, the blue one shows a non-parametric loess smoothing which suggests some non-linearity in the relationship that we will explore later.

Although the IMDb score coefficient is highly statistically significant that should not fool us that we have gained much predictive capacity. The model fit is rather poor. The root mean squared error is 1.25 which is large given the variation in the data. But the inadequate fit is most clearly visible if we plot the actual data versus the predictions. Figure 3 below does just that. The grey bars show the prediction plus/minus two predictive standard errors. If the predictions derived from the model were good, the dots (observations) would be very close to the diagonal (indicated by the dotted line). In this case, they are not. The model does a particularly bad job in predicting very low and very high ratings.

figure3

We can also see how little information IMDb scores contain about (my) personal scores by going back to the raw data. Figure 4 plots to density of my ratings for two sets of values of IMDb scores – from 6.5 to 7.5 (blue) and from 7.5- to 8.5 (red). The means for the two sets differ somewhat, but the overlap in the density is great.

figure4

In sum, knowing the IMDb rating provides some information but on its own doesn’t get us very far in predicting what my score would be.

Take two: adding predictors
Let’s add more variables to see if things improve. Some playing around shows that among the available candidates only the year of release of the movie and dummies for a few genres and directors (selected only from those with more than four movies in the data) give any leverage.

 summary(lm(mine~imdb+d$comedy +d$romance+d$mystery+d$"Stanley Kubrick"+d$"Lars Von Trier"+d$"Darren Aronofsky"+year.c, data=d))

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.074930   0.651223   1.651  .  
imdb                  0.727829   0.087238   8.343  ***
d$comedy             -0.598040   0.133533  -4.479  ***
d$romance            -0.411929   0.141274  -2.916  ** 
d$mystery             0.315991   0.185906   1.700  .  
d$"Stanley Kubrick"   1.066991   0.450826   2.367  *  
d$"Lars Von Trier"    2.117281   0.582790   3.633  ***
d$"Darren Aronofsky"  1.357664   0.584179   2.324  *  
year.c                0.016578   0.003693   4.488  ***
---
Residual standard error: 1.156 on 413 degrees of freedom
Multiple R-squared: 0.3508,	Adjusted R-squared: 0.3382

The fit improves somewhat. The root mean squared error of this model is 1.14. Moreover, looking again at the actual versus predicted ratings, the fit is better, especially for highly rated movies – no surprise given that the director dummies pick these up.

figure5

The last variable in the regression above is the year of release of the movie. It is coded as the difference from 2014, so the positive coefficient implies that older movies get higher ratings. The statistically significant effect, however, has no straightforward predictive interpretation. The reason is again selection bias. I have only watched movies released before the 1990s that have withstood the test of time. So even though in the sample older films have higher scores, it is highly unlikely that if I pick a random film made in the 1970s I would like it more than a random film made after 2010. In any case, Figure 6 below plots the year of release versus the residuals from the regression of my ratings on IMDb scores (for the subset of films after 1960). We can see that the relationship is likely nonlinear (and that I really dislike comedies from the 1980s).

figure6

So far both regressions assumed that the relationship between the predictors and the outcome is linear. Needless to say, there is no compelling reason why this should be the case. Maybe our predictions will improve if we allow the relationships to take any form. This calls for a generalized additive model.

Take three: generalized additive model (GAM)
In R, we can use the mgcv library to fit a  GAM. It doesn’t make sense to hypothesize non-linear effects for binary variables, so we only smooth the effects of IMDb rating and year of release. But why stop there, perhaps the non-linear effects of IMDb rating and release year are not independent, why not allow them to interact!

library(mgcv)
summary(gam(mine ~ te(imdb,year.c)+d$"comedy " +d$"romance "+d$"mystery "+d$"Stanley Kubrick"+d$"Lars Von Trier"+d$"Darren Aronofsky", data = d)) 

PParametric coefficients:
                     Estimate Std. Error t value Pr(|t|)    
(Intercept)           6.80394    0.07541  90.225   ***
d$"comedy "          -0.60742    0.13254  -4.583   ***
d$"romance "         -0.43808    0.14133  -3.100   ** 
d$"mystery "          0.32299    0.18331   1.762   .  
d$"Stanley Kubrick"   0.83139    0.45208   1.839   .  
d$"Lars Von Trier"    2.00522    0.57873   3.465   ***
d$"Darren Aronofsky"  1.26903    0.57525   2.206   *  
---
Approximate significance of smooth terms:
                  edf Ref.df     F p-value    
te(imdb,year.c) 10.85  13.42 11.09

Well, the root mean squared error drops to 1.11 and the jointly smoothed (with a full tensor product smooth) variables are significant, but the added predictive value is minimal in this case. Nevertheless, the plot below shows the smoothed terms are more appropriate than the linear ones, and that there is a complex interaction between the two:

figure7

Take four: models for categorical data
So far we treated personal movie ratings as if they were a continuous variable, but they are not – taking into account that they are essentially an ordered categorical variable might help. But ten categories, while possible to model, would make the analysis rather unwieldy, so we recode the personal ratings into five categories without much loss of information: 5 and less, 6,7,8,9 and more.

We can first see a nonparametric conditional destiny plot of the newly created categorical variable as a function of IMDb scores:
figure8

The plot shows the observed density for each category of the outcome variable along the range of the predictor. For example, for a film with an IMDb rating of ‘6’, about 35% of the personal scores are ‘5’, a further 50% are ‘6’, and the remaining 15% are ‘7’. Remember that the plot is based on the observed conditional frequencies only (with some smoothing), not on the projections of a model. But the small ups and downs seem pretty idiosyncratic. We can also fit an ordered logistic regression model, which would be appropriated for the categorical outcome variable we have, and plot its predicted probabilities given the model.

First, here is the output of the model:

library(MASS)
summary(polr(as.factor(mine.c) ~ imdb+year.c,  Hess=TRUE, data = d)
Coefficients:
        Value Std. Error t value
imdb   1.4103   0.149921   9.407
year.c 0.0283   0.006023   4.699

Intercepts:
    Value   Std. Error t value
5|6  9.0487  1.0795     8.3822
6|7 10.6143  1.1075     9.5840
7|8 12.1539  1.1435    10.6289
8|9 14.0234  1.1876    11.8079

Residual Deviance: 1148.665 
AIC: 1160.665

The coefficients of the two predictors are significant. The plot below shows the predicted probability of the outcome variable – personal movie rating – being in each of the five categories as a function of IMDb rating and illustrates the substantive scale of the effect.

figure9

Compared to the non-parametric conditional density plot above, these model-based predictions are much smoother and have ‘disciplined’ the effect of the predictor to follow a systematic pattern.

It is interesting to ponder which of the two would be more useful for out-of-sample predictions. Despite the fact that the non-parametric one is more faithful to the current data, I think I would go for the parametric model projections. After all, is it really plausible that a random film with an IMDb rating of 5 would have lower chance a getting a 5 from me than a film with an IMDb rating of 6, as the non-parametric conditional density plot suggests? I don’t think so. Interestingly, in this case the parametric model has actually corrected for some of the selection bias and made for more plausible out-of-sample predictions.

Conclusion
In sum, whatever the method, it is not very fruitful to try to predict how much a person (or at least, the particular person writing this) would like a movie based on the average rating the movie gets and covariates like the genre or the director. Non-linear regressions and other modeling tricks offer only marginal predictive improvements over a simple linear regression approach, but bring plenty of insight about the data itself.

What is the way ahead? Obviously, one would want to get more relevant predictors, but, unfortunately, IMDb seems to have a policy against web-scrapping from its database, so one would either have to ask for permission or look at a different website with a more liberal policy (like Rotten Tomatoes perhaps). For me, the purpose of this exercise has been mostly in its methodological educational value, so I think I will leave it at that. Finally, don’t forget to check out the interactive scatterplot of the data used here which shows a user’s entire movie rating history at a glance.

Endnote
As you would have noted, the IMDb ratings come at a greater level of precision (like 7.3) than the one available for individual users (like 7). So a user who really thinks that a film is worth 7.5 has to pick 7 or 8, but its average IMDb score could well be 7.5. If the rating categories available to the user are indeed too coarse, this would show up in the relationship with the IMDb score: movies with an average score of 7.5 would be less predictable that movies with an average score of either 7 or 8. To test this conjecture, a rerun the linear regression models on two subsets of the data: one comprising the movies with an average IMDb rating between 5.9 and 6.1, 6.9 an 7.1, etc., and a  second one comprising those with an average IMDb rating between 5.4 and 5.6, 6.4 and 6.6, etc. The fit of the regression for the first group was better than for the second (RMSE of 1.07 vs. 1.11), but, frankly, I expected a more dramatic difference. So maybe ten categories are just enough.

The evolution of EU legislation (graphed with ggplot2 and R)

During the last half century the European Union has adopted more than 100 000 pieces of legislation. In this presentation I look into the patterns of legislative adoption over time. I tried to create clear and engaging graphs that provide some insight into the evolution of law-making activity: not an easy task given the byzantine nature of policy making in the EU and the complex nomenclature of types of legal acts possible.

The main plot showing the number of adopted directives, regulations and decisions since 1967 is pasted below. There is much more in the presentation. The time series data is available here, as well as the R script used to generate the plots (using ggplot2). Some of the graphs are also available as interactive visualizations via ManyEyes here, here, and here (requires Java). Enjoy.

EU laws over time

New data source for political science researchers

Political Data Yearbook Interactive is a new source for data on election results, turnout and government composition for all EU and some non-European countries. It is basically an online version of the yearbooks that ECPR printed as part of the European Journal for Political Research for many years now.

The interactive online tool has some (limited) visualization options and can export data in several formats.