Análisis de información de accidentes de tránsito en UK
Jaime Paz
Mayo 06, 2018
1. Descripción
El gobierno de UK registra el numero accidentes de tráfico diariamente y los publica dentro de sus sitios web al público en general. La información se obtiene a través de informes de la policía, hay que tomar en cuenta que la data que se registra no toma en cuenta accidentes menores. Toda la información y los datasets son obtenidos por medio del sitio web del Departamento de Transportes de UK:
https://www.dft.gov.uk/traffic-counts/download.php
Descripción de los datasets:
- ukTrafficAADF.csv: registra que tanto tráfico hubo en las carreteras principales en un tiempo definido.
El Core estadístico se concentra en la llamada AADT, que representa Promedio Anual de Flujo Diario y mide como la actividad en segmentos de carretera específicos, en cuanto a la cantidad de vehículos que transitan en ella
accidents_2005_to_2007.csv, accidents_2009_to_2011.csv, accidents_2012_to_2014: estos archivos en Excel en formato CSV concentran toda la información de los accidentes tales como, numero de ocurrencia, severidad, fecha, sitio en donde fue reportado, entre otros. Los periodos incluyen los años del 2005 al 2014, pero el año 2008 no fue registrado.
Local Authority Districts.csv: este archivo fue generado manualmente, descargando dentro del mismo sitio web mencionado, las relaciones existentes entre los ID de los sitios y country name en donde fue reportado el accidente. Cabe mencionar que esto se realizó para poder relacionar los datasets 1 y 2.
Para el desarrolo de este análisis se ha empleado el software R en su versión 3.4.2 (2017-09-28). Para que el código sea reproducible se debe de instalar las librerías apropiadas como se indicará más adelante.
2. Preguntas de investigación y motivación
El objetivo de este análisis de tránsito en UK, es conocer más a fondo, la naturaleza de los accidentes de tránsito. Como principales preguntas de investigación hemos destacado las siguientes:
Que patrones podemos observar en la información?
Que variables aumentan las probabilidades de accidentes?
3.Dar una recomendación para la reducción de accidentes de tránsito, basado en todos los análisis realizados.
- Hacer un análisis geográfico
3. Obteniendo y preparando las fuentes de datos
Una vez la información ha sido descargada de internet y también observado que la data viene en formato tidy: columnas y filas; procedemos a realizar un EDA para conocer más a detalle la información que se nos proporciona y esto nos pueda ayudar a responder las preguntas de investigación como se vaya desarrollando.
Como primer paso es cargar las librerías adecuadas en R a utilizar:
#library(ggmap)
library(dplyr) #util para realizar transformaciones
library(ggplot2) #paquete para realizar gráficos
library(gridExtra) #adiciones al paquete GGPLOT## Warning: package 'gridExtra' was built under R version 3.4.4library(ggmap)## Warning: package 'ggmap' was built under R version 3.4.4#library(vcd)
#library(rcompanion)
#library(MASS) #paquete util que contiene funciones estadísticas
#library(car) #paquete util que contiene funciones estadísticas
#detach("package:dplyr", character.only = TRUE)Cargando los datasources y data cleaning
Cargando los datasets en R:
accidents_2005_2007 <- read.csv("accidents_2005_to_2007.csv", header = TRUE, stringsAsFactors = FALSE)
accidents_2009_2011 <- read.csv("accidents_2009_to_2011.csv", header = TRUE, stringsAsFactors = FALSE)
accidents_2012_2014 <- read.csv("accidents_2012_to_2014.csv", header = TRUE, stringsAsFactors = FALSE)
AADF <- read.csv("ukTrafficAADF.csv", header = TRUE, stringsAsFactors = FALSE)Ordenando la data:
a <- colnames(accidents_2005_2007)
b <- colnames(accidents_2009_2011)
c <- colnames(accidents_2012_2014)
#luego integramos todas las fuentes de accidentes
data_accidents <- rbind(accidents_2005_2007, accidents_2009_2011, accidents_2012_2014)
cat("Total de accidentes reportados:", nrow(data_accidents))## Total de accidentes reportados: 1504150Procemos a seleccionar variables de interés:
accidents_uk <- data_accidents[, c("Year", "Number_of_Casualties", "Date", "Time",
"Road_Type", "Urban_or_Rural_Area",
"Light_Conditions", "Weather_Conditions",
"Road_Surface_Conditions",
"Special_Conditions_at_Site", "Carriageway_Hazards",
"Did_Police_Officer_Attend_Scene_of_Accident")]El nuevo dataframe obtenido de accidents_uk, contiene información reportada sobre:
El año de ocurrencia del accidente; Número de personas afectadas; Fecha y tiempo de registro; qué tipo de carretera es; si es área urbana o rural; condiciones climáticas; condiciones de la carretera en sí; objetos externos en carretera que provocaron el accidente y la policía atendió la escena del accidente.
Ahora se procede a integrar la información de tráfico en el segundo datasets, que contiene las siguientes variables reportadas por área y región: PedalCycles, Motorcycles, CarsTaxis, BusesCoaches, LightGoodsVehicles, LightGoodsVehicles y AllHGVs. En resumen, estas variables contienen información del promedio de volumen por día registrado para cada uno de los automóviles indicados: motocicletas, taxis, buses de transporte urbano, vehículos normales y tipo HGV.
traffic_uk <- AADF %>%
select(AADFYear, Region, LocalAuthority, Road, LinkLength_miles,
PedalCycles,Motorcycles, CarsTaxis, BusesCoaches, LightGoodsVehicles,AllHGVs)
#extraemos la información del promedio en volúmen de tráfico para cada uno de los automóviles:
traffic_total <- traffic_uk %>%
filter(AADFYear == "2005" | AADFYear == "2006" | AADFYear == "2007" | AADFYear ==
"2008" |
AADFYear == "2009" | AADFYear == "2010" | AADFYear == "2011" | AADFYear ==
"2012" | AADFYear == "2013" | AADFYear == "2014") %>%
select(AADFYear, LocalAuthority, PedalCycles, Motorcycles, CarsTaxis, BusesCoaches, LightGoodsVehicles,AllHGVs) %>%
group_by(AADFYear, LocalAuthority) %>%
summarise_at(c("PedalCycles", "Motorcycles", "CarsTaxis", "BusesCoaches", "LightGoodsVehicles", "AllHGVs"), mean, na.rm=TRUE)
#mandamos a llamar la tabla data_accidents , extrayendo año y Local Authority con el conteo de accidentes
accidents_uk_r <- data_accidents %>%
select(Year, Local_Authority_.Highway.) %>%
group_by(Year, Local_Authority_.Highway.) %>%
summarise(cout_acc = n())
colnames(accidents_uk_r)[2] <- "districid"
#ahora mandamos a llamar el datset con los ID de los distritos (dataset # 3)
districts <- read.csv("Local Authority Districts.csv", header = TRUE)
colnames(districts)[1] <- "districid"
joined <- merge(x=accidents_uk_r, y=districts, by = "districid", all.x = TRUE )
#ahora unimos la tabla joined contra la tabla de traffic total
colnames(joined)[4] <- "Area"
colnames(traffic_total)[2] <- "Area"
colnames(traffic_total)[1] <- "Year"
#este dataset servirá para futuros análisis ya que integra la información de tráfico y del conteo de accidentes de tránsito.
total_merged <- merge(joined,traffic_total,by=c("Area","Year")) 4. Análisis Exploratorio (EDA)
Como primera pregunta de las variables seleccionadas, sería interesante saber el número de accidentes que ha provocado afectación en la integridad de las personas en los años reportados. Esto lo separamos por área rural y urbana:
#utilizando librería dplyr para filtrar y mutar la data:
casualities.urban <- accidents_uk %>%
filter(Urban_or_Rural_Area ==1) %>%
select(Year, Road_Type, Number_of_Casualties) %>%
group_by(Year, Road_Type) %>%
summarise(total = sum(Number_of_Casualties))
casualities.rural <- accidents_uk %>%
filter(Urban_or_Rural_Area ==2) %>%
select(Year, Road_Type, Number_of_Casualties) %>%
group_by(Year, Road_Type) %>%
summarise(total = sum(Number_of_Casualties))
# en este punto procedemos a guardar en variables cada una de las gráficas, en donde se utilizará la librería ggplot:
urban <- qplot(Year, total, data = casualities.urban,
geom ="line", color = Road_Type,
xlab = "Year", ylab = "Total casualties",
main = "TOTAL AFECTATION OF URBAN ACCIDENTS IN THE UK (2005-2014)"
) + theme(plot.title = element_text(hjust=0.5)) + scale_x_continuous(
breaks = c(2005, 2006,2007,2008,2009,2010,2011, 2012, 2013, 2014)) + coord_cartesian(xlim = c(2005, 2014)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic"))
rural <- qplot(Year, total, data = casualities.rural,
geom ="line", color = Road_Type,
xlab = "Year", ylab = "Total casualties",
main = "TOTAL AFECTATION OF RURAL ACCIDENTS IN THE UK (2005-2014)"
) + theme(plot.title = element_text(hjust=0.5)) + scale_x_continuous(
breaks = c(2005, 2006,2007,2008,2009,2010,2011, 2012, 2013, 2014)) + coord_cartesian(xlim = c(2005, 2014)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic"))
#extensión de ggplot2 para unificar las gráficas:
grid.arrange(urban, rural,ncol=2)
Tanto para las zonas urbanas y rurales se observa que la mayor cantidad de personas han sido afectadas por carreteras de tipo de un solo carril. Este comportamiento fue reduciendo desde el año 2005 al 2011, llegando a un total de 98,205 y 53,524 respectivamente. A partir de este año ha habido comportamientos ascendentes y descendentes para ambos casos. Ya para finales del año 2014 las áreas urbanas han sumado 93,853 casos para este tipo de carreteras.
Las carreteras de doble carril también han sido causantes de afectaciones en las personas, durante los años especialmente para especialmente en las áreas urbanas, aunque comparado con el año 2005, los casos han reduciendo pero no tan pronunciado como el anterior caso.
Al parecer los accidentes reportados en los “roundabouts” han permanecido casi constantes a través de los años. Si esto no se ha ido reduciendo, es probable que aún no ha habido intención alguna en el gobierno de poder manejar y reducir estos accidentes.
A que horas del dia se tiene el mayor registro de afectacion en la integridad de las personas?
#Para el caso de zonas urbanas:
casualities.urban.day <- accidents_uk %>%
filter(Urban_or_Rural_Area ==1) %>%
select(Time, Number_of_Casualties) %>%
mutate(daytime = ifelse(
Time >= "06:00" & Time <="12:00", "Morning",
ifelse (Time > "12:00" & Time <="20:00", "Evening", "Night"))) %>%
group_by(daytime) %>%
summarise(sum(Number_of_Casualties))
#Para el caso de zonas rurales
casualities.rural.day <- accidents_uk %>%
filter(Urban_or_Rural_Area ==2) %>%
select(Time, Number_of_Casualties) %>%
mutate(daytime = ifelse(
Time >= "06:00" & Time <="12:00", "Morning",
ifelse (Time > "12:00" & Time <="20:00", "Evening", "Night"))) %>%
group_by(daytime) %>%
summarise(sum(Number_of_Casualties))
#"Casos reportados en áreas urbanas:"
print(casualities.urban.day)## # A tibble: 3 x 2
## daytime `sum(Number_of_Casualties)`
## <chr> <int>
## 1 Evening 699032
## 2 Morning 323532
## 3 Night 221708#"Casos reportados en áreas rurales:"
print(casualities.rural.day)## # A tibble: 3 x 2
## daytime `sum(Number_of_Casualties)`
## <chr> <int>
## 1 Evening 416536
## 2 Morning 233334
## 3 Night 137726Para ambos casos la mayoría de los accidentes se han ocasionado por la tarde y esto ha afectado la integridad de las personas siendo la zona urbana la más afectada. Este tiempo del día ha sido clasificado en horas entre 12:00 y 8:00 PM (que es donde aún permanece el sol en UK).
¿Ahora bien, en qué condiciones físicas se encontraban dichas carreteras a la hora de reportar los accidentes?
#tipos de estado de carretera (iluminación)
table(accidents_uk$Light_Conditions)##
## Darkeness: No street lighting
## 82559
## Darkness: Street lighting unknown
## 16120
## Darkness: Street lights present and lit
## 296340
## Darkness: Street lights present but unlit
## 6909
## Daylight: Street light present
## 1102222ggplot(accidents_uk, aes(x= Light_Conditions, fill = ..count..)) + geom_bar(aes(y = (..count..)/sum(..count..))) + scale_y_continuous(labels=scales::percent) + guides(fill=FALSE) + xlab("Group of Road") + ylab("Total") + theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="TOTAL % AFECTATION OF ACCIDENTS BY ROAD STATE IN UK (2005-2014)") 
La mayoría de los accidentes ocurren a plena luz del día y esto corresponde a casi 75% de todos los casos. Existe una parte que necesita ser controlada, y esta corresponde a las carreteras en donde no se tiene iluminación, representando casi un 5% de los casos reportados.
Casi un 20% de todos los casos reportados que han afectado a las personas ocurren por la noche, destacando que las condiciones de iluminación se reportaron adecuadas.
Ahora, sería interesante verificar como las variables climáticas y condiciones de carretera respecto al pavimento (seco, con nieve, resbaloso, etc.).
#tipos de variables climáticas reportadas:
table(accidents_uk$Weather_Conditions)##
## Fine with high winds
## 126 18355
## Fine without high winds Fog or mist
## 1203943 8190
## Other Raining with high winds
## 33503 20813
## Raining without high winds Snowing with high winds
## 177663 1960
## Snowing without high winds Unknown
## 11301 28296#tipos de condiciones de la superficie de la carretera:
table(accidents_uk$Road_Surface_Conditions)##
## Dry
## 1958 1034670
## Flood (Over 3cm of water) Frost/Ice
## 2143 31405
## Snow Wet/Damp
## 10497 423477accidents_uk$Road_Surface_Conditions[accidents_uk$Road_Surface_Conditions == "Flood(Over 3cm of water)"] <- "Flood"
weather <- accidents_uk %>%
select(Weather_Conditions, Road_Surface_Conditions) %>%
group_by(Weather_Conditions, Road_Surface_Conditions) %>%
summarise(n = n()) %>%
mutate(percent = round(n * 100 / nrow(accidents_uk),2)) %>%
select(Weather_Conditions, Road_Surface_Conditions, percent) %>%
arrange(desc(percent))
head(weather)## # A tibble: 6 x 3
## # Groups: Weather_Conditions [5]
## Weather_Conditions Road_Surface_Conditions percent
## <chr> <chr> <dbl>
## 1 Fine without high winds Dry 66.1
## 2 Fine without high winds Wet/Damp 12.6
## 3 Raining without high winds Wet/Damp 11.6
## 4 Unknown Dry 1.44
## 5 Raining with high winds Wet/Damp 1.32
## 6 Other Wet/Damp 1.31Tenemos un total de casi un 66% de casos en los cuales las condiciones climáticas son normales, pero se ha reportado que las carreteras se encuentran en estado seco o en buenas condiciones. Pero existe una parte de los casos en las que las carreteras se han encontrado en estado húmedo, esto corresponde al 12.6 % de los casos.
Tenemos otro 11.6% de los casos en donde el clima se ha encontrado en estado lluvioso, por lo que esto hace que las carreteras estén susceptibles a humedad.
En conclusión, podemos decir que la mayoría de los accidentes suceden a plena luz del día y con buenas condiciones climáticas.
Existe otro tipo anomalías reportadas en las carreteras y que han resultado en accidentes de tránsito con afectación a las personas. Seria interesante ver si dentro de este 66% de accidentes reportados anteriormente son causa alguna de estas anomalías en las carreteras que usualmente se ven.
table(accidents_uk$Carriageway_Hazards)##
##
## 29
## Any animal (except a ridden horse)
## 8014
## Dislodged vehicle load in carriageway
## 1606
## Involvement with previous accident
## 2282
## None
## 1476871
## Other object in carriageway
## 11762
## Pedestrian in carriageway (not injured)
## 3586table(accidents_uk$Special_Conditions_at_Site)##
##
## 15
## Auto traffic signal partly defective
## 785
## Auto traffic singal out
## 2788
## Mud
## 4610
## None
## 1467553
## Ol or diesel
## 5243
## Permanent sign or marking defective or obscured
## 2269
## Road surface defective
## 3664
## Roadworks
## 17223Procediento a realizar gráficos para destacar un resúmen:
external <- accidents_uk %>%
filter(Weather_Conditions == "Fine without high winds" &
Road_Surface_Conditions == "Dry") %>%
select(Number_of_Casualties, Carriageway_Hazards, Special_Conditions_at_Site)
ggplot(external, aes(x= Carriageway_Hazards, fill = ..count..)) + geom_bar(aes(y = (..count..)/sum(..count..))) + scale_y_continuous(labels=scales::percent) + guides(fill=FALSE) + xlab("Carriageway Hazzards") + ylab("Percent") + theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="TOTAL HAZZARDS - ENDING UP ON ACCIDENTS IN THE UK (2005-2014)") 
#Una gran parte se ha reportado por por ningún tipo detectado, se elimina este valor:
external.none.hz <- subset(external, Carriageway_Hazards != "None" )
ggplot(external.none.hz, aes(x= Carriageway_Hazards, fill = ..count..)) + geom_bar(aes(y = (..count..)/sum(..count..))) + scale_y_continuous(labels=scales::percent) + guides(fill=FALSE) + xlab("Carriageway Hazzards") + ylab("Percent") + theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="TOTAL HAZZARDS - ENDING UP ON ACCIDENTS IN THE UK (2005-2014)") 
Solo un 2% de accidentes han sido reportados en donde se ha encontrado anomalías en los carriles de las carreteras y que han causado afectación. Se este 2%, la mayoría es a consecuencia objetos tirados en la carretera con un 45% y el otro grupo destacado es por estrellarse con algún tipo de animal (exceptuando caballos)
Ahora veamos otras causas externas a las anteriores en donde también vale la pena destacar:
external.none.special <- subset(external, Special_Conditions_at_Site != "None" )
ggplot(external.none.special, aes(x= Special_Conditions_at_Site, fill = ..count..)) + geom_bar(aes(y = (..count..)/sum(..count..))) + scale_y_continuous(labels=scales::percent) + guides(fill=FALSE) + xlab("Special conditions at site") + ylab("Percent") + theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="EXTERNAL FACTORS ON ROADS ENDING UP ON ACCIDENTS IN THE UK (2005-2014)") 
Al parecer como en el anterior caso, condiciones especiales tales como daños o trabajos en la carretera, o incluso fallos en las señalizaciones tampoco afectan solamente a un 2% de la población. De este 2%, un 60% de las personas al parecer han sido afectadas por trabajos en la carretera y dejando un 10% de afectación para los casos en donde hay problemas con señalización y otro 10% de casos en donde la carretera se encuentra en malas condiciones
Pre-Análisis Inferencial
De las variables seleccionadas al inicio surge una pregunta. De todos los casos reportados, cuantos de ellos fueron atendidos efectivamente por la policía? Para responder esta pregunta retomamos la data del inicio, y seleccionamos la variable: “Did Police Officer Attend Escene of Accident”
police <- data_accidents[, c("Time","Accident_Severity", "Number_of_Vehicles",
"Number_of_Casualties","Did_Police_Officer_Attend_Scene_of_Accident")]
nrow(police)## [1] 1504150En total como vimos son 1,504,150 casos de accidentes reportados.
#Casos atendidos
table(police$Did_Police_Officer_Attend_Scene_of_Accident)##
## No Yes
## 2922 282351 1218877#severidad de los casos
table(police$Accident_Severity)##
## 1 2 3
## 19441 204504 1280205Aunque se tuvo casos sin indicar 1,218,877 casos si fueron atendidos por la policía, dejando 282,351 casos que no. La mayoría de casos reportados son de tipo severidad 3, es decir con mayor criticidad.
ggplot(police, aes(Accident_Severity, Number_of_Casualties)) +
geom_bar(stat = "identity", aes(fill = Did_Police_Officer_Attend_Scene_of_Accident)) +
xlab("Accident Severity") + ylab("Number of Casualities") +
theme(plot.title = element_text(hjust=0.5)) +
theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="TOTAL POLICE ATTENDANCY ON ACCIDENTS IN THE UK") 
Del gráfico anterior observamos con mayor detalle esta información. Hay una cierta porción de casos en donde a pesar de ser casos críticos, la policía no se presentó a lugar de los hechos. Ahora surge la duda si en realidad hay algún factor que determina este hecho, y uno de ellos es el horario del día en que ha ocurrido los hechos. Debido a esto realizaremos un análisis inferencial para responder la siguiente pregunta:
Existe alguna relación o asociación entre la hora reportada de los accidentes de tránsito severos, y de los cuales, la policía no se presentó al lugar de los hechos?
Pre-Análisis Predictivo
Procemos a análizar la información de los reportes de tráfico en UK. Para esto hemos seleccionado las variables que determinan el volúmen en las carreteras de automóviles. Cabe mencionar que a este dataframe se le ha añadido la variable de conteo de accidentes y también se ha promediado cada una de las variables de automóviles para poder enlazarla con el número de años y el ID del área en donde se reportó el accidente. Procedemos a realizar un análisis descriptivo de estas variables:
colnames(total_merged)## [1] "Area" "Year" "districid"
## [4] "cout_acc" "Region" "PedalCycles"
## [7] "Motorcycles" "CarsTaxis" "BusesCoaches"
## [10] "LightGoodsVehicles" "AllHGVs"summary(total_merged[c(4,5,7,8,9,10,11)])## cout_acc Region Motorcycles
## Min. : 13.0 London :297 Min. : 13.2
## 1st Qu.: 388.0 Scotland :261 1st Qu.: 103.7
## Median : 635.0 North West :180 Median : 144.6
## Mean : 869.1 South East :171 Mean : 272.0
## 3rd Qu.: 960.0 West Midlands :126 3rd Qu.: 242.7
## Max. :6094.0 Yorkshire and the Humber:126 Max. :3132.8
## (Other) :288
## CarsTaxis BusesCoaches LightGoodsVehicles AllHGVs
## Min. : 1636 Min. : 27.72 Min. : 382.7 Min. : 86.84
## 1st Qu.:13612 1st Qu.: 133.94 1st Qu.:2149.0 1st Qu.: 757.10
## Median :17653 Median : 208.09 Median :2881.6 Median :1097.58
## Mean :17867 Mean : 309.33 Mean :2841.8 Mean :1245.23
## 3rd Qu.:21840 3rd Qu.: 373.88 3rd Qu.:3488.7 3rd Qu.:1560.12
## Max. :46892 Max. :1920.55 Max. :7099.7 Max. :4233.56
## Nos encontramos con informacion de conteo de accidentes con una media total de 896.1 accidentes. Tenemos la variable categorica “Región” en la cual se describe la región geográfica de donde se reportan los accidentes. Tanto los Buses Coaches como las MotorCycles se caracterizan por tener valores máximos altos, esto indica que sus distribuciones se encuentran con skewness. Para esto procedemos a realizar un análisis con un boxplot:
colnames(total_merged)## [1] "Area" "Year" "districid"
## [4] "cout_acc" "Region" "PedalCycles"
## [7] "Motorcycles" "CarsTaxis" "BusesCoaches"
## [10] "LightGoodsVehicles" "AllHGVs"ggplot(total_merged, aes(factor(Region), cout_acc)) + geom_boxplot() + coord_flip() + xlab("Region") + ylab("Median") + theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="BOXPLOT BY GROUP - Total Accidents in UK (2005-2014)")
#datas <- total_merged %>%
# select(Region, cout_acc) %>%
# group_by(Region) %>%
# summarise(median(cout_acc))Como vemos, es evidente que la data tiene outliers representativos, especialmente para las regiiones de North West y South East. No parece existir una igualdad entre grupos y en algunos casos se presenta variablidad. Por ejemplo en el caso de South East y East Midlands con East of England. Quizá este último grupo mencionado vale la pena destacarlo, ya que su mediana se encuentra en 2,288, que en este caso es la medida más robusta a favor de la media.
Conviene también ver como está distribuida aproximadamente nuestro total de accidentes reportados en este dataes. Para esto procedemos a realizar un histograma de frecuencias:
hist(total_merged$cout_acc, breaks = 20, xlab = "No. of Accidents", main="HISTOGRAM OF FREQUENCIES: COUNT NUMBER OF ACCIDENTS IN THE UK (2005-2014)")
De esto podemos concluir que nuestra data se asemeja a una distribución de Poisson, aunque también vale la pena mencionar que podría seleccionarse una distribución binomial negativa por su similitud a su distribución. Conviene también aplicar la función de distribución de Kernel en R, para determinar la misma:
d <- density(total_merged$cout_acc) # retorna de función de densidad de probabilidad
plot(d) 
#Como dato adicional, conviene corroborar que nuestra data esté distribuida equitativamente conforme a cada año.
print(table(total_merged$Year))##
## 2005 2006 2007 2009 2010 2011 2012 2013 2014
## 161 161 161 161 161 161 161 161 161Con esta información analizada procederemos a destacar una análisis de predicción:
Es posible establecer una función de predicción multivariada que permita determinar el número de accidentes de tráfico (variable dependiente) tomando como input las variables de volúmen de vehículos (variables independientes) que circulan en las carreteras de UK?
5. Análisis Inferencial
Procedemos a responder nuestra pregunta realizando un análisis para establecer si existe una relación entre la hora del día y la ausencia de la policía para los accidentes clasificados como severos. Se procede a seleccionar la data:
police.attend <- police %>%
filter(Accident_Severity == 3) %>%
select(Did_Police_Officer_Attend_Scene_of_Accident, Time) %>%
mutate(daytime = ifelse(
Time >= "06:00" & Time <="12:00", "Morning",
ifelse (Time > "12:00" & Time <="20:00", "Evening", "Night"))) Esto suma a 26,4154 casos en los que la policía no se presentó al sitio de los hechos, de hecho hay 1,280,205 de casos en donde la severidad reportada es de 3. Ahora recopilamos la información sobre los casos y su relación con el tiempo atendido
police.hour <- xtabs(~Did_Police_Officer_Attend_Scene_of_Accident + daytime ,police.attend)
prop.table(police.hour, margin = 2)## daytime
## Did_Police_Officer_Attend_Scene_of_Accident Evening Morning
## 0.002116516 0.002389467
## No 0.220919181 0.216023927
## Yes 0.776964303 0.781586606
## daytime
## Did_Police_Officer_Attend_Scene_of_Accident Night
## 0.001740588
## No 0.137861493
## Yes 0.860397919De este resultado se observa de los accidentes que fueron severos y que ocurrieron por el día, casi el 77.77% de ellos fueron atendidos por la policia. Existe un 22% del cual ellos no atienden.
Una parte representativa son los accidentes que ocurren de noche, ya que el 86% de ellos fueron atendidos igualmente. Otra parte y similar a los accidentes por la tarde, es que casi un 21.60% no fueron atendidos.
Para establecer si existe una relación en la pregunta realizada, es hacer un análisis de Chi-Square para determinar si existe una independencia entre las variables mencionadas. De los resultados se obtendrá un valor p-value el cual nos ayudará para poder aceptar o rechazar esta hipótesis planteada.
HO (Hipótesis NULA): No existe alguna asociación entre las horas del día y que la policía no se presente al lugar de los hechos, cuando se tienen accidentes de tránsito de severidad crítica.
HA (Hipóteis Alternativa): No existe alguna asociación entre las horas del día y que la policía no se presente al lugar de los hechos, cuando se tienen accidentes de tránsito de severidad crítica.
Un requerimiento para utilizar la fórmula de Chi-Square, es que por lo menos los valores esperados de los counts en las tablas sea mayor a 5. En R, podemos utilizar la función chisq.test():
chisq.test(police.hour)$expected## daytime
## Did_Police_Officer_Attend_Scene_of_Accident Evening Morning
## 1506.0 795.7319
## No 145453.7 76854.0289
## Yes 557972.3 294818.2392
## daytime
## Did_Police_Officer_Attend_Scene_of_Accident Night
## 433.2679
## No 41846.2293
## Yes 160525.5029Todos los valores esperados son mayores a 5, por lo que la prueba de Chi-Square es una buena aproximación para responder nuestra hipótesis. Procemos a calcular el p-value:
chisq.test(police.hour)##
## Pearson's Chi-squared test
##
## data: police.hour
## X-squared = 6981.3, df = 4, p-value < 2.2e-16Los resultados muestran que se tiene una fuerte evidencia que si existe una asociación entre el tiempo del día y que la policía no proceda a presentarse al lugar de los hechos para los casos severos de accidentes. Esto se demuestra por el resultado del valor p-value < 5%, rechazando la hipótesis nula Ho.
6. Análisis predictivo
Para poder establecer el análisis predictivo procederemos a seleccionar dos tipos de data:
Trainning: Esta data contendrá la información de accidentes y sus variables de volúmen de automóviles entre los años 2005 y 2013.
Test: Esta data contendrá la información para el año 2014
Cada año, se registra data (calculada con forecasting o bien registrada físicamente) por lo que se construirá el análisis predictivo haciendo un ponderación de volúmen de automóviles agrupandola por año y luego por la región.
Seleccionamos nuestra data de entrenamiento:
training <- total_merged %>%
filter(Year %in% c(2005, 2006, 2007, 2009, 2010, 2011, 2012, 2013)) %>%
select(cout_acc, Area, Region, PedalCycles, Motorcycles, CarsTaxis, BusesCoaches, LightGoodsVehicles, AllHGVs ) %>%
group_by(Area, Region) %>%
summarise_at(c("cout_acc","PedalCycles", "Motorcycles", "CarsTaxis", "BusesCoaches", "LightGoodsVehicles", "AllHGVs"), funs(mean))
training$cout_acc <- as.integer(training$cout_acc)
training.test <- total_merged %>%
filter(Year %in% c(2014)) %>%
select(cout_acc, Area, Region, PedalCycles, CarsTaxis, BusesCoaches, LightGoodsVehicles) %>%
group_by(Area, Region) %>%
summarise_at(c("cout_acc","PedalCycles", "CarsTaxis", "BusesCoaches", "LightGoodsVehicles"), funs(mean))Intentando hacer la predicción con modelo de Poisson
La distribución de Poisson se caracteriza por poder aplicarse para data del tipo “count”. Debido a la similitud de la densidad de probabilidad establecida en puntos anteriores, se hará un intento por establecer un modelo predictivo de Poisson:
fit.pois <- glm(cout_acc ~ PedalCycles + Motorcycles + CarsTaxis + BusesCoaches + LightGoodsVehicles + AllHGVs + Region,
data = training, family = poisson())
summary(fit.pois)##
## Call:
## glm(formula = cout_acc ~ PedalCycles + Motorcycles + CarsTaxis +
## BusesCoaches + LightGoodsVehicles + AllHGVs + Region, family = poisson(),
## data = training)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.347 -11.368 -4.974 5.311 84.795
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.143e+00 2.040e-02 301.115 <2e-16 ***
## PedalCycles -6.575e-04 2.701e-05 -24.343 <2e-16 ***
## Motorcycles -5.082e-06 2.888e-05 -0.176 0.860
## CarsTaxis -5.275e-05 1.251e-06 -42.177 <2e-16 ***
## BusesCoaches 8.689e-04 2.697e-05 32.213 <2e-16 ***
## LightGoodsVehicles 4.188e-04 1.021e-05 40.996 <2e-16 ***
## AllHGVs -1.237e-04 6.682e-06 -18.504 <2e-16 ***
## RegionEast Midlands 9.867e-01 1.979e-02 49.869 <2e-16 ***
## RegionEast of England 1.535e+00 1.989e-02 77.170 <2e-16 ***
## RegionLondon 4.685e-03 2.213e-02 0.212 0.832
## RegionNorth East -2.582e-01 2.341e-02 -11.030 <2e-16 ***
## RegionNorth West 3.849e-01 1.945e-02 19.786 <2e-16 ***
## RegionScotland -5.084e-01 2.098e-02 -24.231 <2e-16 ***
## RegionSouth East 9.274e-01 1.925e-02 48.174 <2e-16 ***
## RegionSouth West 1.138e+00 3.077e-02 36.987 <2e-16 ***
## RegionWest Midlands 5.805e-01 1.966e-02 29.534 <2e-16 ***
## RegionYorkshire and the Humber 4.695e-01 1.985e-02 23.650 <2e-16 ***
## RegionYorkshire and The Humber 1.386e+00 2.859e-02 48.469 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 88818 on 160 degrees of freedom
## Residual deviance: 49473 on 143 degrees of freedom
## AIC: 50845
##
## Number of Fisher Scoring iterations: 5Al parecer el modelo parece explicar la predicción. Sin embaro hay otros factores por tomar en cuenta. Un problema que caracteriza a la distribución de Poisson es que sufre de sobre dispersión, y esto ocurre cuando las medias alrededor los grupos no es igual.
with(training, tapply(cout_acc, Region, mean)) #overdispersion## East East Midlands East of England
## 536.3333 1414.6667 2516.0000
## England London North East
## NA 773.0606 459.3000
## North West Northern Ireland Scotland
## 850.1000 NA 325.5862
## South East South West United Kingdom
## 1313.7895 1634.0000 NA
## Wales West Midlands Yorkshire and the Humber
## NA 1064.5714 945.4286
## Yorkshire and The Humber
## 2044.0000Es evidente que este hecho se refleja. Se intentará utilizar la distribución Negativa Binomial, la cual se caracteriza por trabajar también con “count” data y el tema de sobre dispersión no es un problema. Es importa destacar también otro tema, que es la Desviación Residual, para el caso de Poisson se tiene “Residual deviance: 49473”.
library(MASS) ## Warning: package 'MASS' was built under R version 3.4.4##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectfit.bin.neg <- glm.nb(cout_acc ~ PedalCycles + Motorcycles + CarsTaxis + BusesCoaches + LightGoodsVehicles + AllHGVs + Region,
data = training, na.action = na.omit)
summary(fit.bin.neg)##
## Call:
## glm.nb(formula = cout_acc ~ PedalCycles + Motorcycles + CarsTaxis +
## BusesCoaches + LightGoodsVehicles + AllHGVs + Region, data = training,
## na.action = na.omit, init.theta = 3.053856059, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0726 -0.7455 -0.3504 0.3467 4.1180
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.077e+00 2.881e-01 21.095 < 2e-16 ***
## PedalCycles -6.096e-04 4.761e-04 -1.281 0.20036
## Motorcycles -1.638e-04 5.108e-04 -0.321 0.74843
## CarsTaxis -4.592e-05 2.169e-05 -2.117 0.03428 *
## BusesCoaches 1.199e-03 4.786e-04 2.505 0.01223 *
## LightGoodsVehicles 3.758e-04 1.817e-04 2.069 0.03859 *
## AllHGVs -9.731e-05 1.202e-04 -0.810 0.41802
## RegionEast Midlands 9.817e-01 3.031e-01 3.239 0.00120 **
## RegionEast of England 1.528e+00 3.493e-01 4.376 1.21e-05 ***
## RegionLondon -8.270e-02 3.233e-01 -0.256 0.79808
## RegionNorth East -3.144e-01 3.063e-01 -1.026 0.30471
## RegionNorth West 3.655e-01 2.716e-01 1.346 0.17836
## RegionScotland -5.475e-01 2.707e-01 -2.022 0.04313 *
## RegionSouth East 8.849e-01 2.798e-01 3.163 0.00156 **
## RegionSouth West 1.157e+00 6.241e-01 1.854 0.06369 .
## RegionWest Midlands 6.009e-01 2.831e-01 2.123 0.03379 *
## RegionYorkshire and the Humber 4.230e-01 2.840e-01 1.490 0.13634
## RegionYorkshire and The Humber 1.388e+00 6.230e-01 2.227 0.02593 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0539) family taken to be 1)
##
## Null deviance: 314.11 on 160 degrees of freedom
## Residual deviance: 169.74 on 143 degrees of freedom
## AIC: 2413.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.054
## Std. Err.: 0.326
##
## 2 x log-likelihood: -2375.895En este caso, la desviación residual es “Residual deviance: 169.74”. Mucho menor a la de Poisson, por lo que se procederá con este modelo.
Ahora bien, de los resultados vemos que no todas las variables son siginificativas. Por ejemplo, las que contienen “*" parecen explicar bien el modelo a un nivel determinado de significancia. Sin embargo, existen otras como PedalCycles y MotorCycles que no parecen explicar bien el modelo.
Para poder aproximar de mejor manera el modelo se toman en cuenta muchas situaciones, pero dentro de las que más destacan son las siguientes:
- Presencia de Outliers
- Correlación entre las mismas variables
- Heteroscedasticity
- Valores que contienen gran peso que son capaces de mover la curva de predicción de manera significativa.
- Presencia de errores de estimación (RMSE)
Al iniciar vale la pena calcular primeramente el Root Mean Square Errors, que en su manera ideal debería de ser cero o si se compara con otros modelos, debe ser menor. En nuestro caso tenemos:
rmse = function(y, ypred) {
rmse = sqrt(mean((y - ypred)^2))
return(rmse)
}
fit.bin.neg.yhat <- predict(fit.bin.neg)
rmse(training$cout_acc, fit.bin.neg.yhat)## [1] 1184.206Por lo que mas adelante, valdría la pena considerar el poder reducirlos. Antes de poder el assessing de la predicción, vale la pena seleccionar también aquellas variables que son significativas en el modelo. Para esto se utiliza la metodología STEP WISE, la cual permite ir quitando aquellas variables que no aportan al modelo de dos formas: STEP BACK (dejar todas las variables e irlas quitando una a una) y STEP FORWARD (Incluir las variables una a una hasta que se obtenga el mejor modelo).
La clave del modelo STEPWISE es poder reducir los errores AIC (Akaike Information Criterion) el cual es un estimador de la calidad de modelos estadísticos; este se resume a
AIC = 2k - 2 Ln(L)
en donde k es el número de parámetros estimados y L el valor máximo de la función de parámetros estadísticos. El software R permite realizar estos modelos por medio del paquete “MASS”.
# En este caso se selecciona el mejor modelo utilizando la opción BOTH
library(MASS)
step <- stepAIC(fit.bin.neg, direction ="both")## Start: AIC=2411.9
## cout_acc ~ PedalCycles + Motorcycles + CarsTaxis + BusesCoaches +
## LightGoodsVehicles + AllHGVs + Region
##
## Df AIC
## - Motorcycles 1 2410.0
## - AllHGVs 1 2410.4
## - PedalCycles 1 2411.4
## <none> 2411.9
## - CarsTaxis 1 2413.4
## - LightGoodsVehicles 1 2413.4
## - BusesCoaches 1 2416.5
## - Region 11 2481.0
##
## Step: AIC=2409.99
## cout_acc ~ PedalCycles + CarsTaxis + BusesCoaches + LightGoodsVehicles +
## AllHGVs + Region
##
## Df AIC
## - AllHGVs 1 2408.4
## <none> 2410.0
## - CarsTaxis 1 2411.4
## - LightGoodsVehicles 1 2411.8
## + Motorcycles 1 2411.9
## - PedalCycles 1 2414.4
## - BusesCoaches 1 2414.8
## - Region 11 2479.8
##
## Step: AIC=2408.45
## cout_acc ~ PedalCycles + CarsTaxis + BusesCoaches + LightGoodsVehicles +
## Region
##
## Df AIC
## <none> 2408.4
## - CarsTaxis 1 2409.7
## - LightGoodsVehicles 1 2409.9
## + AllHGVs 1 2410.0
## + Motorcycles 1 2410.4
## - PedalCycles 1 2412.9
## - BusesCoaches 1 2415.4
## - Region 11 2478.1print(step$anova)## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## cout_acc ~ PedalCycles + Motorcycles + CarsTaxis + BusesCoaches +
## LightGoodsVehicles + AllHGVs + Region
##
## Final Model:
## cout_acc ~ PedalCycles + CarsTaxis + BusesCoaches + LightGoodsVehicles +
## Region
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 143 169.7446 2411.895
## 2 - Motorcycles 1 0.004400593 144 169.7490 2409.986
## 3 - AllHGVs 1 0.016687126 145 169.7657 2408.447Con este último parámetro se hace una comparación ANOVA (Análisis de Varianza para determinar si los resultados son significativos). Como se representa en el resultado final, el modelo ha determinado que la variable MotorCycles y AllHGVs no aportan al modelo y no son significativas. Se procede a crear el nuevo modelo:
fit.bin.neg.final <- glm.nb(cout_acc ~ PedalCycles + CarsTaxis + BusesCoaches + LightGoodsVehicles +
Region, data = training, na.action = na.omit)
summary(fit.bin.neg)##
## Call:
## glm.nb(formula = cout_acc ~ PedalCycles + Motorcycles + CarsTaxis +
## BusesCoaches + LightGoodsVehicles + AllHGVs + Region, data = training,
## na.action = na.omit, init.theta = 3.053856059, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0726 -0.7455 -0.3504 0.3467 4.1180
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.077e+00 2.881e-01 21.095 < 2e-16 ***
## PedalCycles -6.096e-04 4.761e-04 -1.281 0.20036
## Motorcycles -1.638e-04 5.108e-04 -0.321 0.74843
## CarsTaxis -4.592e-05 2.169e-05 -2.117 0.03428 *
## BusesCoaches 1.199e-03 4.786e-04 2.505 0.01223 *
## LightGoodsVehicles 3.758e-04 1.817e-04 2.069 0.03859 *
## AllHGVs -9.731e-05 1.202e-04 -0.810 0.41802
## RegionEast Midlands 9.817e-01 3.031e-01 3.239 0.00120 **
## RegionEast of England 1.528e+00 3.493e-01 4.376 1.21e-05 ***
## RegionLondon -8.270e-02 3.233e-01 -0.256 0.79808
## RegionNorth East -3.144e-01 3.063e-01 -1.026 0.30471
## RegionNorth West 3.655e-01 2.716e-01 1.346 0.17836
## RegionScotland -5.475e-01 2.707e-01 -2.022 0.04313 *
## RegionSouth East 8.849e-01 2.798e-01 3.163 0.00156 **
## RegionSouth West 1.157e+00 6.241e-01 1.854 0.06369 .
## RegionWest Midlands 6.009e-01 2.831e-01 2.123 0.03379 *
## RegionYorkshire and the Humber 4.230e-01 2.840e-01 1.490 0.13634
## RegionYorkshire and The Humber 1.388e+00 6.230e-01 2.227 0.02593 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.0539) family taken to be 1)
##
## Null deviance: 314.11 on 160 degrees of freedom
## Residual deviance: 169.74 on 143 degrees of freedom
## AIC: 2413.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.054
## Std. Err.: 0.326
##
## 2 x log-likelihood: -2375.895Encontramos que el modelo STEP WISE no fue capaz de retirar algunas variables que incluso dentro del modelo final no parecen ser significativas. Ahora bien, esto no significa que no aporten al modelo, de hecho algunas de ellas pueden ser del tipo “Coufunding” ya que el que estén dentro del modelo hacen que el modelo sea representativo. Continuaremos con los demás análisis para poder llegar a un modelo final
Podemos expresar el valor de cada coeficiente con su valor original (ya que aplicamos una distribución binomial negativa) :
exp(coef(fit.bin.neg.final))## (Intercept) PedalCycles
## 432.1244192 0.9992648
## CarsTaxis BusesCoaches
## 0.9999571 1.0013306
## LightGoodsVehicles RegionEast Midlands
## 1.0002871 2.6834412
## RegionEast of England RegionLondon
## 4.7414457 0.9144142
## RegionNorth East RegionNorth West
## 0.7819008 1.4775199
## RegionScotland RegionSouth East
## 0.5834277 2.5536480
## RegionSouth West RegionWest Midlands
## 3.1789326 1.8662007
## RegionYorkshire and the Humber RegionYorkshire and The Humber
## 1.5755222 3.8471596Por ejemplo, una unidad promedio de volúmen por día en Buses Coaches, está asociada con un incremento de 1.0013306 accidentes, manteniendo las otras variables constantes.
Algo que debemos determinar siempre es la calidad de los Residuals, que son la diferencia existente entre cada punto Y y el valor predictivo. Se trata de minimizar estos errores y podemos obtener una gráfica del comportamiento de los mismos, que en el mejor de los casos deben seguir la línea de predicción.
par(mfrow = c(2,2))
plot(fit.bin.neg.final)## Warning: not plotting observations with leverage one:
## 30, 96
## Warning: not plotting observations with leverage one:
## 30, 96
st.residuals <- rstandard(fit.bin.neg.final)
hist(st.residuals)
Según vemos en las gráficas anteriores y la distribución de frecuencias de los residuales, se puede observar una aproximación a una curva normal lo cual nos beneficia. Por otro lado observamos en la primera gráfica los valores residuales vs los Leverage values (Cooks Distance). Estos últimos como se mencionó anteriormente, son los que tienen un gran peso en la predicción del modelo. Es decir, si se eliminan existiría una gran diferencia a favor. Ahora bien, estos valores no simplemente pueden ser eliminados, debe haber justificaciones adecuadas y conocer el entorno en materia del caso para poder ejecutarlo; por lo que estos puntos deben ser analizados detenidamente.
cooksd <- cooks.distance(fit.bin.neg.final)
plot(cooksd, pch="*", cex=2, main="Influential Obs by Cooks distance") # plot cook's distance
abline(h = 4*mean(cooksd, na.rm=T), col="red") # add cutoff line
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red") # add labels
influential <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
head((training[influential, ]))## # A tibble: 6 x 9
## # Groups: Area [6]
## Area Region cout_acc PedalCycles Motorcycles CarsTaxis BusesCoaches
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Birming~ West M~ 3130 63.7 162 23473 420
## 2 City of~ Scotla~ 1256 141 170 20988 545
## 3 <NA> <NA> NA NA NA NA NA
## 4 Hampshi~ South ~ 3359 45.9 247 23202 132
## 5 Kent South ~ 4471 43.2 205 17635 186
## 6 Lancash~ North ~ 4069 51.7 128 17306 169
## # ... with 2 more variables: LightGoodsVehicles <dbl>, AllHGVs <dbl>Vemos algunos outliers que sobrepasan los límites establecidos por Cooks Distance. Por ejemplo el número 73 nos indica que este puede afectar de gran manera nuestro modelo. Como se dijo anteriormente, estos valores requiren investigación para poder eliminarlos.
Para realizar el análisis de regresión se elegirán dos variables como ejemplo: PedalCyles y Buses Coaches (Manteniendo las demás variables constantes en su media), y luego se ploteará dos gráficas:
- La proyección de accidentes para el año 2014 para cada variable
- La proyección real de accidentes del 2014 para cada variable
Por último se traza una gráfica de efectos, en la cual se relaciona el modelo binomial en general proyectado para los accidentes contra todas las variables que se eligieron para la predicción.
regs <- factor(c("East Midlands", "East of England","London", "North East",
"North West", "Scotland", "South East",
"South West", "West Midlands", "Yorkshire and the Humber", "Yorkshire and the Humber"))
detach("package:dplyr", character.only = TRUE)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniontraining.1 <- total_merged %>%
filter(Year %in% c(2005, 2006, 2007, 2009, 2010, 2011, 2012,
2013)) %>%
select(cout_acc, Area, Region, PedalCycles, CarsTaxis,
BusesCoaches, LightGoodsVehicles) %>%
group_by(Area, Region) %>%
summarise_at(c("cout_acc","PedalCycles", "CarsTaxis",
"BusesCoaches", "LightGoodsVehicles"), funs(mean))
training.1$cout_acc <- as.integer(training.1$cout_acc)
mynb <- glm.nb(cout_acc ~ PedalCycles + CarsTaxis + BusesCoaches + LightGoodsVehicles +
Region,
data = training.1, na.action = na.omit)
#variable Pedal Cycles
newdata1.p <- with(training.1,
data.frame(PedalCycles = mean(PedalCycles), CarsTaxis =
mean(CarsTaxis), BusesCoaches = mean(BusesCoaches), LightGoodsVehicles = mean(LightGoodsVehicles),
Region = regs))
newdata2.p <- with(training.1,
data.frame(PedalCycles = training.test$PedalCycles,
CarsTaxis = mean(CarsTaxis), BusesCoaches =
mean(BusesCoaches),
LightGoodsVehicles = mean(LightGoodsVehicles),
Region = Region))
pred.p = predict(mynb, newdata=newdata2.p, type="response")
newdata2.p$pred.p <- pred.p
newdata2.p$pred.2014 <- training.test$cout_acc
se <- predict(mynb, newdata=newdata2.p, type = "response", se.fit = TRUE)
newdata2.p$se <- se$se.fit
LL <- newdata2.p$pupper_logit <- newdata2.p$pred.p + (1.96 * newdata2.p$se) # 95% CI upper bound
UL <- newdata2.p$plower_logit <- newdata2.p$pred.p - (1.96 * newdata2.p$se) # 95% CI lower bound
newdata2.p$ll <- LL
newdata2.p$ul <- UL
pedalcycles.pred.p <- ggplot(newdata2.p, aes(training.test$PedalCycles, pred.p)) +
geom_ribbon(aes(ymin = LL, ymax = LL, fill = Region), alpha = .25) +
geom_line(aes(colour = Region), size = 2) +
labs(x = "Pedal Cycles", y = "Predicted Values") + xlim(0, 500) +
theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="PREDICTED ACCIDENTS FOR PEDAL CYCLES IN 2014")
pedalcycles.real.p <- ggplot(newdata2.p, aes(training.test$PedalCycles, pred.2014)) +
geom_ribbon(aes(ymin = LL, ymax = UL, fill = Region), alpha = .25) +
geom_line(aes(colour = Region), size = 2) +
labs(x = "Pedal Cycles", y = "Predicted Values") + xlim(0, 500) +
theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="REAL ACCIDENTS FOR PEDAL CYCLES IN 2014")
######
#Proyección de Buses Coaches
newdata2.c <- with(training.1,
data.frame(PedalCycles = mean(PedalCycles),
CarsTaxis = mean(CarsTaxis), BusesCoaches =
training.test$BusesCoaches,
LightGoodsVehicles = mean(LightGoodsVehicles),
Region = Region))
pred.c = predict(mynb, newdata=newdata2.c, type="response")
newdata2.c$pred.c <- pred.c
newdata2.c$pred.2014 <- training.test$cout_acc
se <- predict(mynb, newdata=newdata2.c, type = "response", se.fit = TRUE)
newdata2.c$se <- se$se.fit
LL <- newdata2.c$pupper_logit <- newdata2.c$pred.c + (1.96 * newdata2.c$se) # 95% CI upper bound
UL <- newdata2.c$plower_logit <- newdata2.c$pred.c - (1.96 * newdata2.c$se) # 95% CI lower bound
newdata2.c$ll <- LL
newdata2.c$ul <- UL
buses.pred.C <- ggplot(newdata2.c, aes(training.test$BusesCoaches, pred.c)) +
geom_ribbon(aes(ymin = LL, ymax = UL, fill = Region), alpha = .25) +
geom_line(aes(colour = Region), size = 2) +
labs(x = "Pedal Cycles", y = "Predicted Values") + xlim(0, 500) +
theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="PREDICTED ACCIDENTS FOR BUSES COACHES IN 2014")
buses.real.c <- ggplot(newdata2.p, aes(training.test$BusesCoaches, pred.2014)) +
geom_ribbon(aes(ymin = LL, ymax = UL, fill = Region), alpha = .25) +
geom_line(aes(colour = Region), size = 2) +
labs(x = "Pedal Cycles", y = "Predicted Values") + xlim(0, 500) +
theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic")) + labs(title="REAL ACCIDENTS FOR BUSES COACHES IN 2014")
grid.arrange(pedalcycles.pred.p , pedalcycles.real.p)## Warning: Removed 12 rows containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_path).
grid.arrange(buses.pred.C, buses.real.c)## Warning: Removed 22 rows containing missing values (geom_path).## Warning: Removed 22 rows containing missing values (geom_path).
library(effects)## Warning: package 'effects' was built under R version 3.4.4## Loading required package: carData## Warning: package 'carData' was built under R version 3.4.4## lattice theme set by effectsTheme()
## See ?effectsTheme for details.plot(allEffects(mynb, lattent = FALSE))
7. Dashboard de accidentes en UK
map.acc.1 <- data_accidents %>%
filter(Number_of_Vehicles >= 10) %>%
select(Local_Authority_.Highway., Latitude, Longitude, Number_of_Vehicles )
uk <- as.numeric(geocode("United Kingdom"))## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=United%20Kingdom&sensor=false#-3.435973 55.378051
ukmap <- ggmap(get_googlemap(center=uk, zoom=5), extent="normal")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=55.378051,-3.435973&zoom=5&size=640x640&scale=2&maptype=terrain&sensor=false#ukmap + geom_point(data=map.acc.1, aes(x=Longitude, y=Latitude, size = Number_of_Vehicles, color = Number_of_Vehicles))
ukmap + geom_point(data=map.acc.1, alpha=0.8, aes(x=Longitude, y=Latitude,
colour = cut(Number_of_Vehicles, c(-Inf, 25, 40, Inf))), size = 4) +
scale_color_manual(name = "Number of Accidents",
values = c("(-Inf,25]" = "gold",
"(25,40]" = "darkorange",
"(40, Inf]" = "firebrick1"),
labels = c("10 - 25", "25 - 40", " > 40")
) + labs(title="Total Severe Accidents in the UK (2005-2014)") +
theme(plot.title = element_text(hjust=0.5)) + theme(plot.title = element_text(hjust=0.5, size = 15, face = "italic"))
8. CONCLUSIONES
```
De todo el análisis realizado para el reporte de accidentes de tránsito en UK podemos destacar lo siguiente:
Se ha visto que los roundabouts o redondeles, han sido consecuencia de accidentes de tráfico. Aunque estos están dentro del 2% de todas las causas en general, pudo detectar que este no se ha controlado en todos años reportados
Las horas en donde ha ocurrido mayor afectación en todas las regiones de UK ha sido por horas de la tarde.
El 75% de todos los accidentes ocurren a plena luz del día. El otro 25% se debe a otras razones, pero destaca entre ellas en el que no se cuenta con iluminación adecuada en la carretera, esto representa un 5% de este total.
Con un análisis de Inferencia Estadística, se pudo aceptar la hipótesis de que la ausencia de las autoridades cuando existen accidentes severos, está relacionada dependiendo de las horas del día (mañana, tarde o noche).
Se estableció un modelo de regresión binomial, el cual nos pudo dar una idea del comportamiento de los accidentes conforme al volúmen de automóviles promedio en las carreteras. Se tomó una data de training y una de testing correspondiendo esta última para el año 2014. Según los resultados, se puede observar que para las predicciones de accidentes causados por Pedal Cycles difiere conforme a las regiones observando tendencias en cada una. Se observa que la predicción más óptima es para la región North West y London.
Para el caso de Buses Coaches, las regiones de London y North West siguen siendo los mejores predictores comparado con los resultados del año 2014. Existen fluctuaciones destacas en el año 2014 especialmente para los casos de South East y West Midlands.
Se detectaron outliers y valores leverage durante el análisis. La eliminación de estos outliers pueden mejorar el desempeño del modelo de predicción planteado.
En la gráfica de interacciones entre variables, es evidente que los comportamientos en general para todas las áreas, decrementan conforme el volúmen de Pedal Cycles y Car Taxis aumentan. Para el caso de Light Good Vehicles y Buses Coaches, a mayor volúmen en las carreteras se observa un incremento en accidentes en general. Se observa que en el gráfico de interacciones, las diferencias por área son evidentes conforme a el número de accidentes
Recomenaciones:
De las conclusiones, se puede ver que el número de accidentes puede controlarse en ciertos aspectos, especialmente para el caso de detectar las carreteras que se encuentran en mal estado de iluminación.
Conviene investigar las causas del por qué la policía no ha llegado a atender las causas severas de accidentes. Si en realidad las condiciones del día dependen de ello, probablemente en ciertas áreas no se tiene la mayor disponibilidad de agentes en cada una de las zonas.
Existe una relación que con los Buses Coaches y los LightGoodsVehicles el número de accidentes incremente conforme el volúmen de ellos incremente. Si estos son las causas principales de los accidentes, conviene investigar las rutas y carreteras en donde ellos transitan. Quizá algunas no sean las adecuadas y a través de la ubicación geográfica se puede determinar las mayores afluencias en ellas.