Исследование данных о транспортных инцидентах на железной дороге

Аннотация

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

Это то, что я собираюсь изучить: тренды и география ЖД инцидентов в Америке.

Введение

В этом отчете анализируется набор данных о практически 230 000 транспортных инцидентах, связанных с железнодорожными переездами и произошедших в США с 1975 по 2016 гг.

Постановка задачи:

Для полноты картины нужно:

провести разведочный анализ даныых на предмет зависимостей между переменными;
посмотреть историческое изменение ситуации, проверить на наличие сезонных ухудшений ситуации;
определить самые “ненадежные” железно-дорожные компании;
проанализировать географическое распределение инцидентов.

Описание набора данных:

Набор данных incidents.zip загружен с сайта Министерства транспорта США и содержит информацию о транспортных инцидентах, связанных с железнодорожными переездами и произошедших в США с 1975 по 2016 гг.

Набор данных содержит порядка 230 000 строк. В исходном наборе данных имеются следующие столбцы:

Reporting Railroad - название конкретной железной дороги;
Reporting Railroad Code - код железной дороги;
Track Owner - владелец пути;
Track Owner Code - код владельца пути;
Highway User Speed - скорость транспортного средства на шоссе;
Calendar Year - календарный год;
Fiscal Year - отчетный год;
CrossingID - идентификационный номер ж-д переезда;
County - округ;
Highway - шоссе;
DisplayDate - дата происшествия;
Highway User Type - тип транспортного средства на шоссе;
Rail Equipment Type - тип транспортного средства на железной дороге;
Carrying Hazmat - перевозились ли опасные грузы одним из участников ДТП;
Releasing Hazmat - произошел ли выброс опасных грузов в результате ДТП;
Display Time - время происшествия;
id - идентификационный номер инцидента;
nonsuicidefatality - количество погибших;
nonsuicideinjury - количество пострадавших;
City - город;
State Name - штат;
Longitude - географическая долгота;
Latitude - географическая широта;
GeoLocation - местоположение.

Подготовка данных к анализу

Загрузка данных

# Код для загрузки данных
incidents <- read_csv ("incidents.csv") 

# Чтобы названия штатов не различались в написании между таблицами
incidents$`State Name` <- tolower(incidents$`State Name`) 

# Создание новой таблицы: отбираем нужные столбцы 
incidents_new <- select(incidents, `Reporting Railroad`:`Reporting Railroad Code`, `Highway User Speed`:`Calendar Year`, `DisplayDate`:`Carrying Hazmat`, `Display Time`:nonsuicideinjury, `State Name`:Longitude) 

# Переименовывание столбцов для удобства
names(incidents_new) <- c("railroad", "railroad_code", "speed", "year", "date", "highway_type", "rail_type", "dangerous_goods", "time", "id", "fatality", "injury", "region", "lat", "long") 

summary (incidents_new)

##    railroad         railroad_code          speed            year     
##  Length:229665      Length:229665      Min.   : 0.00   Min.   :1975  
##  Class :character   Class :character   1st Qu.: 0.00   1st Qu.:1979  
##  Mode  :character   Mode  :character   Median : 5.00   Median :1985  
##                                        Mean   :11.66   Mean   :1988  
##                                        3rd Qu.:20.00   3rd Qu.:1996  
##                                        Max.   :99.00   Max.   :2016  
##                                        NA's   :26548                 
##      date           highway_type        rail_type        
##  Length:229665      Length:229665      Length:229665     
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##                                                          
##  dangerous_goods        time                id            fatality       
##  Length:229665      Length:229665     Min.   :     1   Min.   : 0.00000  
##  Class :character   Class1:hms        1st Qu.: 57417   1st Qu.: 0.00000  
##  Mode  :character   Class2:difftime   Median :114833   Median : 0.00000  
##                     Mode  :numeric    Mean   :114833   Mean   : 0.09726  
##                                       3rd Qu.:172249   3rd Qu.: 0.00000  
##                                       Max.   :229665   Max.   :12.00000  
##                                                                          
##      injury            region               lat             long        
##  Min.   :  0.0000   Length:229665      Min.   :20.89   Min.   :-891.87  
##  1st Qu.:  0.0000   Class :character   1st Qu.:33.55   1st Qu.: -95.89  
##  Median :  0.0000   Mode  :character   Median :37.98   Median : -89.07  
##  Mean   :  0.3719                      Mean   :37.53   Mean   : -91.59  
##  3rd Qu.:  1.0000                      3rd Qu.:41.43   3rd Qu.: -83.83  
##  Max.   :101.0000                      Max.   :98.30   Max.   : -10.47  
##                                        NA's   :21664   NA's   :21667

Выведенная сводка по новой таблице, включающей в себя 15 переменных, демонстрирует общую информацию об исходных данных. Стоит отметить, что данные находятся не в идеальном состоянии. Наблюдаются пропущенные значения в таких переменных, как Latitude, Longitude и Highway User Speed. На следующем шаге необходимо детальнее оценить данные и подготовить их к анализу.

В дальнейшем будет удобно использовать карты для визуализации географических данных. С этой целью была подгружена библиотека maps с географическими координатами штатов (таблица map), их центральными координатами и краткими названиями (таблица states):

map <- map_data("state")
states <- data.frame(state.center, state.abb)

ggplot() + 
  geom_polygon(data = map, aes(x = long, y = lat, group = group),col = "#333300", fill= "#CCCCCC") +
  geom_text(data = states, aes(x = x, y = y, label = state.abb))

# Создание новых переменных
incidents_new <- mutate(incidents_new, lost = ifelse(fatality > 0, 1, 0)) # Показывает, есть ли погибшие 
incidents_new <- mutate(incidents_new, victim = fatality + injury) # Общее количество жертв

Проверка технического качества

Исследование на предмет пропущенных значений

С помощью функции is.na() проверяем, является ли значение неизвестным. Функция is.na() для каждого значения вернет true или false. Далее с помощью функции sum() мы просуммируем значения true, выполнив эту операцию для каждого из столбцов набора данных avocado с помощью функции lapply()

Анализ пропущенных значений показал, что пропущенные значения присутствуют в 11 столбцах из 24: speed (26548); highway_type (2); time (11810); Latitude и Longitude (21664).

По некоторым переменным потери незначительные. Может получиться так, что часть пропущенных значений в одной переменной не влияет на выводы при анализе другой. На основании этого, некоторые потери можно устранить непосредственно перед анализом данных с применением тех или иных переменных.

Географическая карта

Как было сказано ранее, данные содержат значения широты и долготы инцидентов. С их помощью можно изучить географическое распереление инцидентов и попробовать выявить некие закономерности. При этом строки с ‘пропущенными значениями’ придется удалить. Замена их на “типичные” приведет к искажению информации.

ggplot(incidents_new, aes(x = long, y = lat, xlab = "Долгота", ylab = "Широта")) + 
  geom_point(col = "#333300", size = 1)

Необходимо избавиться от ‘выбросов’.

Крайние точки Америки: - Остров Сент-Круа, Виргинские острова США (17°45′19″ с. ш., 64°33′54″ з. д.) - самая восточная точка территории США; - Алеутские острова, штат Аляска (52°55′14″ с. ш., 172°26′16″ в. д.) - самая западная точка пятидесяти штатов США; - Мыс Барроу, штат Аляска (71°23’20" с. ш., 150°28’45" з. д.) - самая северная точка территории США (и пятидесяти штатов США); - Атолл Розе, Американское Самоа (14°34’11" ю. ш., 152°9’10" з. д.) - самая южная точка территории США;

Можно выделить следующие границы координат: Долгота (от -64 до -173) и Широта (от 14 до 72). Для дальнейшей работы избавлямся от выбросов:

incidents_new <- filter(incidents_new, (long > -152 & long < -64) & (lat > 14 & lat < 72)) # Убрали выбросы

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

outliers <- subset(incidents_new, long < -125)
unique(outliers$region)

## [1] "alaska"

# Используются при построении карты
incidents_new_2 <- subset(incidents_new, ! region %in% c("hawaii", "alaska"))
incidents_new_2 <- filter(incidents_new, lat < 50)
#ggplot(incidents_new_2, aes(x=long, y=lat)) +
#  geom_point(size = 1)

Разведочный анализ данных

Определение дисперсии и стандартного отклонения

Первым делом в разведочном анализе необходимо посмотреть распределение числовых переменных:

# Дисперсия и станд. отклонение для переменной speed
var(incidents_new$speed, na.rm = TRUE); sd(incidents_new$speed, na.rm = TRUE)

## [1] 190.0709

## [1] 13.78662

a <- ggplot (incidents_new) + 
  geom_histogram (aes (x = incidents_new$speed, y=..density..), binwidth = 25) + 
  stat_function(fun = dnorm, sd = sd(incidents_new$speed), color="blue", size=1)

# Дисперсия и станд. отклонение для переменной fatality
var(incidents_new$fatality, na.rm = TRUE); sd(incidents_new$fatality, na.rm = TRUE)

## [1] 0.1415148

## [1] 0.3761846

b <- ggplot (incidents_new) + 
  geom_histogram (aes (x = incidents_new$fatality, y=..density..), binwidth = 0.5) + 
  stat_function(fun = dnorm, sd = sd(incidents_new$fatality), color="blue", size=1)

# Дисперсия и станд. отклонение для переменной injury
var(incidents_new$injury, na.rm = TRUE); sd(incidents_new$injury, na.rm = TRUE)

## [1] 0.8663124

## [1] 0.930759

c <- ggplot (incidents_new) + 
  geom_histogram (aes (x = incidents_new$injury, y=..density..), binwidth = 5) + 
  stat_function(fun = dnorm, sd = sd(incidents_new$injury), color="blue", size=1)

# Графики распределения
ggarrange(a, b, c + rremove("x.text"), 
          labels = c("A", "B", "C"),
          ncol = 2, nrow = 2)

Определение коэффициента корреляции Спирмена

Так как переменные не имеют нормального закона распределения, то стоит использовать коэффициент корреляции Спирмена. Коэффициент корреляции позволит сказать, есть ли зависимость между переменными и насколько она сильна.

cor.test(incidents_new$speed, incidents_new$victim, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  incidents_new$speed and incidents_new$victim
## S = 7.8123e+14, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2569196

Коэффициент корреляции составляет примерно 0,26. На основании этого можно сделать вывод, что зависимость очень слабая.

Анализ исторического изменения

Можно предположить, что ситуация могла меняться с течением времени. С одной стороны, население планеты растет; люди становятся более мобильными, а поехать путешествовать по стране уже не является чем-то особенным; компании стремятся к минимизации издержак и выносят заводы и склад в места с максимально дешевой стоимостью земли, за счет чего возрастает потребность в доставке товара. Все это ведет к тому, что возрастает само количество перевозок и может стать причиной увеличения числа инцидентов. С другой стороны, мир не стоит на месте, технологии совершенствуются. Это может вести, наоборот, к сокращению числа инцидентов.

Более того, сравнивать инциденты нужно не только по их числу, но и по количеству жертв.

lost = 0 - жертв нет

lost = 1 - жертвы есть (пострадавшие или погибшие)

# График "Количество инцидентов в год"
d <- ggplot(incidents_new, aes(x=fct_rev(year), fill = factor(lost))) +
  geom_bar(position = 'dodge') +
  coord_flip()

# График "Количество жертв в год"
e <- ggplot(incidents_new, aes(x=fct_rev(year), fill = factor(lost), y = victim)) +
  stat_summary(fun.y = 'sum', geom = 'bar', position = 'dodge') +
  coord_flip()

ggarrange(d, e + rremove("x.text"), 
          labels = c("Число инцидентов", "Количество жертв"),
          ncol = 2, nrow = 1)

# Сравнение значимых лет (по предыдущему графику)
incidents_new %>%
  filter((year %in% c("1978", "1989", "2016")) & victim > 0) %>%
  ggplot(aes(x = year, y = victim)) +
  geom_boxplot()

Действительно, к концу 70-х и к концу 80-х годов наблюдалось уведичение инцидентов и их жертв, но ближе к настоящему времени ситуация стала значительно лучше.

Выявление сезонности

Можно предположить, что количество инцидентов в течение года может изменяться в зависимости от погодных условий или времени суток. Ведь в плохую погоду хуже сцепление в поверхностью, а ночью или в снегопад - снижается видимость.

# График: Сезонность по месяцам 
f <- ggplot(incidents_new, na.rm = TRUE, aes(x=month(date), fill = factor(lost), y = victim)) +
  stat_summary(fun.y = 'sum', geom = 'density', position = 'dodge', alpha = 0.6) 

# График: Сезонность по часам
g <- ggplot(incidents_new, na.rm = TRUE, aes(x=hour(time), fill = factor(lost), y = victim)) +
  stat_summary(fun.y = 'sum', geom = 'density', position = 'dodge', alpha = 0.6) 

# График: Распределение по часам
h <- ggplot(incidents_new) +
  geom_histogram(aes(x = hour(time), fill = "#FF0000"))

ggarrange(f, g, h + rremove("x.text"), 
          labels = c("Месяцы", "Часы", "Гистограмма по часам"),
          ncol = 2, nrow = 2)

Предположения оправдались лишь частично: на зимние месяцы приходится большее количество тяжелых (с жертвами) инцидентов, но вот по времени наблюдаются пики днем. Скорее всего это связано с тем, что днем совершается большее количество поездок. Правильнее было бы смотреть относительные показатели.

Самая ненадежная и опасная дорога

С точки зрения потребителя мы можем заинтересоваться, какая дорога является самой неспокойной из представленных.

# Дорога, на которой в среднем наблюдается большее количество жертв
incidents_new %>% 
  group_by(railroad) %>% 
  summarise(VICTIM = mean(victim)) %>%
  arrange(desc(VICTIM)) %>%
  top_n(10, wt = VICTIM) %>%
  ggplot(aes(x = fct_rev(fct_inorder(railroad)), y = VICTIM)) +
  geom_point(alpha = 0.6) +
  coord_flip()

## # A tibble: 1 x 17
##   railroad railroad_code speed year  date       highway_type rail_type
##   <chr>    <chr>         <int> <fct> <date>     <chr>        <chr>    
## 1 West Vi… WVC              55 2013  2013-10-11 Truck        Psgr Tra…
## # ... with 10 more variables: dangerous_goods <chr>, time <time>,
## #   id <int>, fatality <int>, injury <int>, region <chr>, lat <dbl>,
## #   long <dbl>, lost <fct>, victim <int>

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

Попробуем посмотреть по количеству инцидентов на той или иной дороге:

incidents_new %>% 
  group_by(railroad) %>% 
  summarise(times = length(unique(id))) %>%
  arrange(desc(times)) %>%
  top_n(10, wt = times) %>%
  ggplot(aes(x = fct_rev(fct_inorder(railroad)), y = times)) +
  geom_point(alpha = 0.6) +
  coord_flip()

ggplot() +
  geom_polygon(data = map, aes(x = long, y = lat, group = group),col = "#333300", fill= "#CCCCCC", alpha = 0.8) +
  geom_text(data = states, aes(x = x, y = y, label = state.abb)) +
  geom_point(data = (subset(incidents_new_2, railroad_code %in% c("UP"))), aes(x = long, y = lat, group = NULL), colour = "#FF0000", alpha = 0.5, size = 1)

Union Pacific Railroad (UP) — американская компания, владеющая самой большой сетью железных дорог в США. Основана в 1862 году.

Большое количество инцидентов объясняется долгим существованием самой дороги. Сравнение лишь инцидентов без данных обо всех поездках на дороге не даст нам полноценных выводов.

Географическое распределение

Первым делом можно спроицировать инциденты на карту.

ggplot() +
  geom_polygon(data = map, aes(x = long, y = lat, group = group),col = "#333300", fill= "#CCCCCC") +
  geom_text(data = states, aes(x = x, y = y, label = state.abb)) +
  geom_point(data = subset(incidents_new_2, year %in% c("2006", "2007", "2008", "2009", "2010", "2011", "2012","2013", "2014", "2015","2016")), aes(x = long, y = lat, group = NULL), colour = "#FF0000", alpha = 0.5, size = 1)

Количество инцидентов в восточной части страны больше скорее всего из-за большей плотности населения здесь.

Следует проследить местоположения крупных инцидентов, на которые стоит обратить внимание государства.

ggplot() +
  geom_polygon(data = map, aes(x = long, y = lat, group = group),col = "#333300", fill= "#CCCCCC") +
  geom_text(data = states, aes(x = x, y = y, label = state.abb)) +
  geom_point(data = subset(incidents_new_2, victim > 10), aes(x = long, y = lat, group = NULL), colour = "red", alpha = 0.9, 
             size = 1)

Итоговые результаты и выводы

Итоговые визуализации и результаты

Целью исследования был анализ железнодорожных инцидентов. Первым делом были проверены переменные на наличие зависимостей между ними. Сильной корреляции обнаружено не было.

Ниже приведены основные визуализации и выводы к ним.

Географическое распределение крупных инцидентов

Проецирование инцидентов на карту дает возможность наглядно увидеть распределение инцидентов. В данном случае были отобраны инциденты с количеством жертв превышающим 10 человек. По карте можно выделить два “центра” с такими инцидентами: районы штатов IL и CA. Во-первых, это улучшает восприятие информации и позволяет обратить внимание государства на отдельные участки дорог. Стоит проверить состояние этих участков, а также проверить, качественно ли компании, в собственности которых находятся эти ЖД, произволят ремонт. Во-вторых, это позволяет показать пользователям места, в которых следует быть предельно внимательными. В-третьих, призвать машинистов отслеживать всю ситуацию на дороге и незамедлительно реагировать на любые угрозы.

Самая “ненадежная” дорога

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

incidents_new %>% 
  group_by(railroad) %>% 
  summarise(times = length(unique(id))) %>%
  arrange(desc(times)) %>%
  top_n(10, wt = times) %>%
  ggplot(aes(x = fct_rev(fct_inorder(railroad)), y = times)) +
  geom_point(alpha = 0.6) +
  coord_flip()

Анализ исторического изменения

Тем не менее, ситуация с годами меняется. На сегодняшний день и самих происшествий, и их жертв в разы меньше, чем в 70-х и 80-х годах 20 столетия. Скорее всего, большую роль играют технологии, которые сейчас появляются и улучшаются большими темпами: отследивание состояния дорог в реальном времени, возможность быстрого реагирования на ЧП, быстрая проверка транспортного средства на предмет поломок… Стоит отметить, что на сегодняшний день происшествия уже не так сильно зависят от погодных условий или условий видимости, нежели ранее.