На днях завершился очередной чемпионат мира по хоккею.

За просмотром матчей родилась идея. Когда в перерывах телевизионная камера показывает уходящих в раздевалку игроков, трудно не заметить, насколько они огромные. На фоне тренеров, функционеров команд, сотрудников ледовой арены, журналистов или просто фанатов они, как правило, выглядят очень внушительно.

Вот, к примеру, восходящие звезды финского хоккея, Патрик Лайне и Александр Барков, вместе с преданными поклонниками

Источник

И я задался вопросами. Действительно ли хоккеисты выше обычных людей? Как изменяется рост хоккеистов со временем в сравнении с обычными людьми? Есть ли устойчивые межстрановые различия?

Данные

IIHF, организация, проводящая чемпионаты мира по хоккею, каждый год публикует составы участвующих команд с информацией о росте и весе каждого игрока. Архив этих данных тут.

Я собрал вместе данные всех чемпионатов мира с 2001 по 2016 годы. От года к году формат предоставления данных слегка меняется, что требует некоторых усилий по их очистке. Не представляя, как грамотно автоматизировать процесс, все данные копировал вручную, что заняло чуть больше 3 часов. Объединенный датасет выложил в открытый доступ.

# load required packages
require(dplyr) # data manipulation
require(lubridate) # easy manipulations with dates
require(ggplot2) # visualization
require(ggthemes) # themes for ggplot2
require(cowplot) # nice alignment of the ggplots
require(RColorBrewer) # generate color palettes
require(texreg) # easy export of regression tables
require(xtable) # export a data frame into an html table

# download the IIHF data set; if there are some problems, you can download manually
# using the stable URL (https://dx.doi.org/10.6084/m9.figshare.3394735.v2)
df <- read.csv('https://ndownloader.figshare.com/files/5303173')

# color palette
brbg11 <- brewer.pal(11,'BrBG')

Растут ли хоккеисты? Грубое (периодное) сравнение

Для начала сравним средний рост игроков на всех 16 чемпионатах мира.

# mean height by championship
df_per <- df %>% group_by(year) %>%
        summarise(height=mean(height))

gg_period_mean <- ggplot(df_per, aes(x=year,y=height))+
        geom_point(size=3,color=brbg11[9])+
        stat_smooth(method='lm',size=1,color=brbg11[11])+
        ylab('height, cm')+
        xlab('year of competition')+
        scale_x_continuous(breaks=seq(2005,2015,5),labels=seq(2005,2015,5))+
        theme_few(base_size = 15)+
        theme(panel.grid=element_line(colour = 'grey75',size=.25))

gg_period_jitter <- ggplot(df, aes(x=year,y=height))+
        geom_jitter(size=2,color=brbg11[9],alpha=.25,width = .75)+
        stat_smooth(method='lm',size=1,se=F,color=brbg11[11])+
        ylab('height, cm')+
        xlab('year of competition')+
        scale_x_continuous(breaks=seq(2005,2015,5),labels=seq(2005,2015,5))+
        theme_few(base_size = 15)+
        theme(panel.grid=element_line(colour = 'grey75',size=.25))

gg_period <- plot_grid(gg_period_mean,gg_period_jitter)

Положительный тренд очевиден. За полтора десятилетия средний рост хоккеиста на чемпионате мира увеличился почти на 2 сантиметра (левая панель). Как будто бы незначительный прирост на фоне довольно большой вариации (правая панель). Много это или мало? Чтобы ответить на вопрос, надо корректно сравнить с населением (но об этом ближе к концу статьи).

Когортный анализ

Более корректный способ изучения изменения в росте подразумевает сравнение по когортам рождения. Тут мы сталкиваемся с любопытным нюансом — некоторые хоккеисты участвовали не в одном чемпионате мира. Вопрос: вычищать ли повторные записи для одних и тех же людей? Если нам интересен средний рост хоккеиста на чемпионате (как на картинке выше), пожалуй, не имеет смысла зачищать. Но если мы хотим проследить изменение роста хоккеистов как таковое, на мой взгляд, было бы неправильно присваивать больший вес тем игрокам, которые регулярнее попадали на чемпионаты мира. Поэтому для дальнейшего анализа я очистил данные от повторных записей одних и тех же игроков.

# remove double counts
dfu_h <- df %>% select(year,name,country,position,birth,cohort,height) %>%
        spread(year,height)
dfu_h$av.height <- apply(dfu_h[,6:21],1,mean,na.rm=T)
dfu_h$times_participated <- apply(!is.na(dfu_h[,6:21]),1,sum)

dfu_w <- df %>% select(year,name,country,position,birth,cohort,weight) %>%
        spread(year,weight)
dfu_w$av.weight <- apply(dfu_w[,6:21],1,mean,na.rm=T)

dfu <- left_join(dfu_h %>% select(name,country,position,birth,cohort,av.height,times_participated),
                 dfu_w %>% select(name,country,position,birth,cohort,av.weight),
                 by = c('name','country','position','birth','cohort')) %>%
        mutate(bmi = av.weight/(av.height/100)^2)

Общее количество наблюдений сократилось с 6292 до 3333. Если хоккеист участвовал более чем в одном чемпионате мира, данные о росте и весе я усреднял, поскольку рост и (в особенности) вес отдельно взятого хоккеиста мог меняться со временем. Сколько же раз хоккеисты удостаиваются чести сыграть за национальные сборные на чемпионатах мира? В среднем чуть менее 2 раз.

# frequencies of participation in world championships

mean(dfu$times_participated)

df_part <- as.data.frame(table(dfu$times_participated))

gg_times_part <- ggplot(df_part,aes(y=Freq,x=Var1))+
        geom_bar(stat='identity',fill=brbg11[9])+
        ylab('# of players')+
        xlab('times participated (out of 16 possible)')+
        theme_few(base_size = 15)

Но есть и уникумы. Посмотрим, кто из игроков принял участие как минимум в 10 чемпионатах мира. Таких игроков оказалось 14.

# the leaders of participation in world championships
# save the table to html
leaders <- dfu %>% filter(times_participated > 9)
View(leaders)
print(xtable(leaders), type="html", file="table_leaders.html")

Александр Овечкин, 11 раз! Но тут надо отметить, что не для всех хоккеистов в принципе возможно было поучаствовать во всех 16 чемпионатах: зависит когорты рождения (насколько игровая карьера пересеклась именно с этим периодом наблюдения), от того, участвовала ли сборная игрока во всех чемпионатах мира (см. рисунок 3) и попадал ли игрок стабильно в сборную; наконец есть еще НХЛ, стабильно отвлекающий лучших из лучших от участия в чемпионатах мира.

# countries times participated
df_cnt_part <- df %>% select(year,country,no) %>%
        mutate(country=factor(paste(country))) %>%
        group_by(country,year) %>%
        summarise(value=sum(as.numeric(no))) %>%
        mutate(value=1) %>%
        ungroup() %>%
        mutate(country=factor(country, levels = rev(levels(country))),
               year=factor(year))

d_cnt_n <- df_cnt_part %>% group_by(country) %>%
        summarise(n=sum(value))

gg_cnt_part <- ggplot(data = df_cnt_part, aes(x=year,y=country))+
        geom_point(color=brbg11[11],size=7)+
        geom_text(data=d_cnt_n,aes(y=country,x=17.5,label=n,color=n),size=7,fontface=2)+
        geom_text(data=d_cnt_n,aes(y=country,x=18.5,label=' '),size=7)+
        scale_color_gradientn(colours = brbg11[7:11])+
        xlab(NULL)+
        ylab(NULL)+
        theme_bw(base_size = 25)+
        theme(legend.position='none',
              axis.text.x = element_text(angle = 90, hjust = 1,vjust=0.5))

Растут ли хоккеисты? Регрессионный анализ

Регрессионный анализ позволяет более корректно ответить на вопрос об изменении роста игроков. В данном случаем с помощью мультиноминальной линейной регрессии предсказыватся рост хоккеиста в зависимость от когорты рождения. Включая в спецификацию регресиионной модели различные дополнительные (контрольные) переменные, мы получаем значение наиболее интересующего нас коэффициента "при прочих равных". Например, добавляя к объясняющим переменным помимо когорты рождения позицию игрока на поле, мы получаем взаимосвязь роста и когорты, очищенную от эффекта различий в зависимости от позиции; добавляя в контрольны переменные страны, получаем результат, очищенный от межстрановых различий. Разумеется, если контрольные переменные сами оказываются значимыми, на это тоже стоит обратить внимание.
Регрессионные модели (особенно линейные регрессии) очень чувствительны к выбросам (см., например, эту статью). Не вдаваясь глубоко в эту обширную тему, я лишь убрал из анализа когорты, для которых мы имеем слишком небольшое количество представителей.

# remove small cohorts
table(dfu$cohort)
dfuc <- dfu %>% filter(cohort<1997,cohort>1963)

Не желая резать данные сильно, я убрал только когорты 1963, 1997 и 1998 годов рождения, для которых у нас есть менее 10 игроков.

Итак, результаты рагрессионного анализа. В каждой следующей модели я добавляю одну переменную.
Зависимая переменная: рост хоккеиста.
Объясняющие перемеенные: 1) когорта рождения; 2) + позиция на поле (сравнение с защитниками); 3) + страна (сравнение с Россией).

# relevel counrty variable to compare with Russia
dfuc$country <- relevel(dfuc$country,ref = 'RUS')

# regression models
m1 <- lm(data = dfuc,av.height~cohort)
m2 <- lm(data = dfuc,av.height~cohort+position)
m3 <- lm(data = dfuc,av.height~cohort+position+country)

# export the models to html
htmlreg(list(m1,m2,m3),file = 'models_height.html',single.row = T)

Интерпретация моделей

Модель 1. Увеличение когорты на один год соответсвует увеличению роста хоккеистов на 0.1 см. Коэффициент статистически значим, но при этом модель объясняет лишь 1% вариации зависимой переменной. В принципе это не проблема, поскольку моделирование носит объясняющий характер, задача предсказания не ставится. Тем не менее, низкий коэффициент детерминации показывает, что должны быть другие переменные, гораздо лучше объясняющие различия между хоккеистами в росте.

Модель 2. Защитники — самые высокие игроки в хоккее. Вратари ниже на 2 см, нападающие — на 2.6 см. Все коэффициенты статистически значимы. Объясненная вариация зависимой переменной возрастает до 6%. При этом коэффициент при переменной когорта рождения не изменяется.

Модель 3. Добавление контрольных переменных для стран любопытно по двум причинам. Во-первых, некоторые различия статистически значимы и интересны сами по себе. Так, например, шведы, словаки и канадцы статистически значимо выше наших игроков. Большинство же наций значительно ниже нас, японцы аж на 5.2 см, итальянцы — на 3.6 см, французы — на 2.2 см (см. также рисунок 4). Во-вторых, введение контрольных переменных для стран значительно уменьшает коэффициент при переменной когорта рождения — до 0.08. Это значит, что межстрановые различия объясняют часть различий по когортам рождения. Коэффициент детерминации модели возрастает до 13%.

R code. Рисунок 4. Рост хоккеистов по странам

# players' height by country
gg_av.h_country <- ggplot(dfuc ,aes(x=factor(cohort),y=av.height))+
        geom_point(color='grey50',alpha=.25)+
        stat_summary(aes(group=country),geom='line',fun.y = mean,size=.5,color='grey50')+
        stat_smooth(aes(group=country,color=country),geom='line',size=1)+
        #geom_hline(yintercept = mean(height),color='red',size=.5)+
        facet_wrap(~country,ncol=4)+
        coord_cartesian(ylim = c(170,195))+
        scale_x_discrete(labels=paste(seq(1965,1995,10)),breaks=paste(seq(1965,1995,10)))+
        theme_few(base_size = 15)+
        theme(legend.position='none',
              panel.grid=element_line(colour = 'grey75',size=.25))

Наиболее полная модель показывает, что увеличение роста хоккеистов происходит со скоростью 0.08 см в год. Это означает прирост 0.8 см за десятилетие или на 2.56 см за 32 года с 1964 по 1996. Обратите внимание, что при учете контрольных переменных скорость увеличения роста хоккеистов оказывается примерно в полтора раза ниже, чем при более грубом анализе средних значений (рисунок 1): 0.8 см за десятилетие против примерно 1.2 см.

Прежде чем мы, наконец, постараемся понять, насколько значительным оказывается увеличение роста, хочу обратить внимание еще на один любопытный момент. Введение контрольных переменных подразумевает фиксацию различий между категориями при едином наклоне регрессионной линии (единый коэффициент при главной объясняющей переменной). Это не всегда хорошо и может замаскировать значительные различия в тесноте связи между исследуемыми переменными в подвыборках. Так, например, раздельное моделирование зависимости роста игроков от амплуа (рисунок 5) показывает, что взаимосвязь наиболее ярко выражена для вратарей и наименее заметна для защитников.

dfuc_pos <- dfuc
levels(dfuc_pos$position) <- c('Defenders','Forwards','Goalkeeprs')

gg_pos <- ggplot(dfuc_pos ,aes(x=cohort,y=av.height))+
        geom_jitter(aes(color=position),alpha=.5)+
        stat_smooth(method = 'lm', se = T,color=brbg11[11],size=1)+
        scale_x_continuous(labels=seq(1965,1995,5),breaks=seq(1965,1995,5))+
        scale_color_manual(values = brbg11[c(8,4,10)])+
        facet_wrap(~position,ncol=3)+
        xlab('birth cohort')+
        ylab('height, cm')+
        theme_few(base_size = 20)+
        theme(legend.position='none',
              panel.grid=element_line(colour = 'grey75',size=.25))

# separate models for positions
m3d <- lm(data = dfuc %>% filter(position=='D'),av.height~cohort+country)
m3f <- lm(data = dfuc %>% filter(position=='F'),av.height~cohort+country)
m3g <- lm(data = dfuc %>% filter(position=='G'),av.height~cohort+country)
htmlreg(list(m3d,m3f,m3g),file = '2016/160500 Hockey players/models_height_pos.html',single.row = T,
        custom.model.names = c('Model 3 D','Model 3 F','Model 3 G'))

Раздельное моделирование показывает, что в когортах 1964-1996 годов рождения, средний рост хоккеистов, участвовавших в чемпионатах мира в 2001-2016 годах, увеличивался со скоростью 0.4 см за десятиление для защитников, 0.7 см — для нападающих и (!) 2.4 см — для вратарей. За три десятиления средний рост вратарей увеличился на 7 см!

Пришло время сравнить эти изменения со средними значениями для населения.

Сравнение с населением

Результаты регрессионного анализа фиксируют значительные межстрановые различия. Поэтому сравнивать имеет смысл по странам: хоккеистов определенной страны с мужским населением этой же страны.

Для сравнения роста хоккеистов со средними показателями мужского населения я использовал данные из релевантной научной статьи (PDF). Данные я скопировал из статьи (использовав замечательную программку tabula) и тоже разместил в открытом доступе.

# download the data from Hatton, T. J., & Bray, B. E. (2010). 
# Long run trends in the heights of European men, 19th–20th centuries. 
# Economics & Human Biology, 8(3), 405–413. 
# http://doi.org/10.1016/j.ehb.2010.03.001
# stable URL, copied data (https://dx.doi.org/10.6084/m9.figshare.3394795.v1)
df_hb <- read.csv('https://ndownloader.figshare.com/files/5303878')

df_hb <- df_hb %>%
        gather('country','h_pop',2:16) %>%
        mutate(period=paste(period)) %>%
        separate(period,c('t1','t2'),sep = '/')%>%
        transmute(cohort=(as.numeric(t1)+as.numeric(t2))/2,country,h_pop)

# calculate hockey players' cohort height averages for each country
df_hoc <- dfu %>% group_by(country,cohort) %>%
        summarise(h_hp=mean(av.height)) %>%
        ungroup()

К сожалению, данные о динамике роста населения пересекаются лишь с 8 странами из моего хоккейного датасета: Австрия, Дания, Финляндия, Франция, Германия, Италия, Норвегия, Швеция.

# countries in both data sets
both_cnt <- levels(factor(df_hb$country))[which(levels(factor(df_hb$country)) %in% levels(df_hoc$country))]
both_cnt

gg_hoc_vs_pop <- ggplot()+
        geom_path(data = df_hb %>% filter(country %in% both_cnt), aes(x=cohort,y=h_pop),
                  color=brbg11[9],size=1)+
        geom_point(data = df_hb %>% filter(country %in% both_cnt), aes(x=cohort,y=h_pop),
                   color=brbg11[9],size=2)+
        geom_point(data = df_hb %>% filter(country %in% both_cnt), aes(x=cohort,y=h_pop),
                   color='white',size=1.5)+
        geom_point(data = df_hoc %>% filter(country %in% both_cnt), aes(x=cohort,y=h_hp),
                   color=brbg11[3],size=2,pch=18)+
        stat_smooth(data = df_hoc %>% filter(country %in% both_cnt), aes(x=cohort,y=h_hp),
                    method='lm',se=F,color=brbg11[1],size=1)+
        facet_wrap(~country,ncol=2)+
        ylab('height, cm')+
        xlab('birth cohort')+
        theme_few(base_size = 15)+
        theme(panel.grid=element_line(colour = 'grey75',size=.25))

Во всех проанализировнных странах хоккеисты выше стеднестатистических мужчин на 2-5 см. Но это не удивительно — в спорте значительная селекция.
Примечательно другое. В развитых странах мира особенно бурное увеличение роста мужского населения происходило в первой середине 20 века. В когортах примерно 1960-х годов рождения рост мужчин приблизился к плато и пеерстал бурно увеличиваться. Тренд среднего роста хоккеистов во всех странах (кроме почему-то Дании) как будто бы продолжил приостановившийся многолетний тренд всего мужского населения.
Для когорт европейцев, родившихся в первой половине 20 века, темпы увеличения среднего роста варьировались от 1.18 до 1.74 см за десятилетие в зависимости от страны (рисунок 7). Начиная с 1960-х годов этот показатель опустился до уровня 0.15-0.80 за 10 лет.

# growth in population

df_hb_w <- df_hb %>% spread(cohort,h_pop) 
names(df_hb_w)[2:26] <- paste('y',names(df_hb_w)[2:26])

diffs <- df_hb_w[,3:26]-df_hb_w[,2:25]

df_hb_gr<- df_hb_w %>%
        transmute(country,
                  gr_1961_1980 = unname(apply(diffs[,22:24],1,mean,na.rm=T))*2,
                  gr_1901_1960 = unname(apply(diffs[,9:21],1,mean,na.rm=T))*2,
                  gr_1856_1900 = unname(apply(diffs[,1:8],1,mean,na.rm=T))*2) %>%
        gather('period','average_growth',2:4) %>%
        filter(country %in% both_cnt) %>%
        mutate(country=factor(country,levels = rev(levels(factor(country)))),
               period=factor(period,labels = c('1856-1900','1901-1960','1961-1980')))

gg_hb_growth <- ggplot(df_hb_gr, aes(x=average_growth,y=country))+
        geom_point(aes(color=period),size=3)+
        scale_color_manual(values = brbg11[c(8,3,10)])+
        scale_x_continuous(limits=c(0,2))+
        facet_wrap(~period)+
        theme_few()+
        xlab("average growth in men's height over 10 years, cm")+
        ylab(NULL)+
        theme_few(base_size = 20)+
        theme(legend.position='none',
              panel.grid=element_line(colour = 'grey75',size=.25))

На фоне стагнирующего тренда в населении увеличение роста хоккеистов выглядит весьма внушительным. А акселерация среди вратарей вообще беспрецедентна.
Не стоит забывать и про селекцию. Расхождение трендов в населении и среди хоккеистов, вероятно, свидетельствует об усиливающейся селекции — хоккей требует все большего роста для успешной карьеры.

Селекция в спорте

Проглядывая научную литературу по теме я наткнулся на примечательный результат. Оказывается, в профессиональном спорте преобладают люди, рожденные в первой половине года. Объясняется это тем, что спортивные секции, как правило, формируют детские команды по когортам рождения. Таким образом, рожденные в начале года, всегда имеют чуть больше прожитого времени за плечами, что зачастую прямо выражается в физическом превосходстве над сверстниками, рожденными под конец года. Нетрудно проверить этот результат на нашем датасете.

# check if there are more players born in earlier months
df_month <-  df %>% mutate(month=month(birth)) %>%
        mutate(month=factor(month,levels = rev(levels(factor(month)))))

gg_month <- ggplot(df_month,aes(x=factor(month)))+
        geom_bar(stat='count',fill=brbg11[8])+
        scale_x_discrete(breaks=1:12,labels=month.name)+
        xlab('month of birth')+
        coord_flip()+
        theme_few(base_size = 20)+
        theme(legend.position='none',
              panel.grid=element_line(colour = 'grey75',size=.25))

Действительно, респределение довольно сильно смещено в сторону ранних месяцев. Если разбить данные по декадам рождения, то невооруженным глазом видно, что эффект усиливается со временем (рисунок 9). Косвенно это свидетельствует о том, что селекция в хоккее становится жестче.

# facet by decades
df_month_dec <- df_month %>%
        mutate(dec=factor(substr(paste(cohort),3,3),labels = paste('born in',c('1960s','1970s','1980s','1990s'))))

gg_month_dec <- ggplot(df_month_dec,aes(x=factor(month)))+
        geom_bar(stat='count',fill=brbg11[8])+
        scale_x_discrete(breaks=1:12,labels=month.abb)+
        xlab('month of birth')+
        facet_wrap(~dec,ncol=2,scales = 'free')+
        theme_few(base_size = 20)+
        theme(legend.position='none',
              panel.grid=element_line(colour = 'grey75',size=.25))

На будущее

Любопытно будет посмотреть, влияют ли физические данные на игровую статистику хоккеистов. Наткнулся на занимательную статью, опубликованную в очень приличном научном журнале, в которой авторы нашли корреляцию между соотношением пропорций лица хоккеиста и средним количеством штрафных минут за игру.

График из указанной статьи

Источник

Reproducibility

Полный R скрипт, воспроизводящий результаты моей статьи, тут.
Использована версия R-3.2.4
Все пакеты по состоянию на 2016-03-14. В случае пакетных несовместимостей, данный код будет гарантированно воспроизведен при использовании пакета checkpoint с указанием соответствующей даты.