No se trata de una pauta para repetir sino diferentes ejemplos para daros ideas e inspiraros.
Se ha seleccionado un conjunto de datos del National Highway Traffic Safety Administration. El sistema de informes de análisis de mortalidad fue creado en los Estados Unidos por la National Highway Traffic Safety Administration para proporcionar una medida global de la seguridad en las carreteras. (Fuente Wikipedia). Los datos pertenecen al año 2020. Se trata de un conjunto de registros de accidentes que recogen datos significativos que los describen. Todos los accidentes tienen alguna víctima mortal como mínimo. El objetivo analítico que tenemos en mente es entender que hace que un accidente sea grave y que quiere decir que sea grave. https://www.nhtsa.gov/crash-data-systems/fatality-analysis-reporting-system
Queremos hacer una primera aproximación al conjunto de datos escogido y responder a las preguntas más básicas: ¿Cuánto registros tiene? ¿Cuántas variables? ¿De qué tipología son? ¿Cómo se distribuyen los valores de las variables? ¿Hay problemas con los datos, por ejemplo, campos vacíos? ¿Puedo intuir ya el valor analítico de los datos? ¿Qué primeras conclusiones puedo extraer?
El primer paso para realizar un análisis exploratorio es cargar el fichero de datos.
path = 'accident_2020.CSV'
accidentData <- read.csv(path, row.names=NULL)Verificamos la estructura del juego de datos principal. Vemos el número de columnas que tenemos y ejemplos de los contenidos de las filas.
structure = str(accidentData)## 'data.frame': 35766 obs. of 81 variables:
## $ STATE : int 1 1 1 1 1 1 1 1 1 1 ...
## $ STATENAME : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ ST_CASE : int 10001 10002 10003 10004 10005 10006 10007 10008 10009 10010 ...
## $ VE_TOTAL : int 1 4 2 1 1 2 1 2 2 2 ...
## $ VE_FORMS : int 1 4 2 1 1 2 1 2 2 2 ...
## $ PVH_INVL : int 0 0 0 0 0 0 0 0 0 0 ...
## $ PEDS : int 0 0 0 0 0 0 1 0 0 0 ...
## $ PERSONS : int 4 6 2 5 1 3 1 2 4 3 ...
## $ PERMVIT : int 4 6 2 5 1 3 1 2 4 3 ...
## $ PERNOTMVIT : int 0 0 0 0 0 0 1 0 0 0 ...
## $ COUNTY : int 51 73 117 15 37 103 73 25 45 95 ...
## $ COUNTYNAME : chr "ELMORE (51)" "JEFFERSON (73)" "SHELBY (117)" "CALHOUN (15)" ...
## $ CITY : int 0 350 0 0 0 0 330 0 0 1500 ...
## $ CITYNAME : chr "NOT APPLICABLE" "BIRMINGHAM" "NOT APPLICABLE" "NOT APPLICABLE" ...
## $ DAY : int 1 2 2 3 4 4 7 8 9 10 ...
## $ DAYNAME : int 1 2 2 3 4 4 7 8 9 10 ...
## $ MONTH : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MONTHNAME : chr "January" "January" "January" "January" ...
## $ YEAR : int 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ DAY_WEEK : int 4 5 5 6 7 7 3 4 5 6 ...
## $ DAY_WEEKNAME: chr "Wednesday" "Thursday" "Thursday" "Friday" ...
## $ HOUR : int 2 17 14 15 0 16 19 7 20 10 ...
## $ HOURNAME : chr "2:00am-2:59am" "5:00pm-5:59pm" "2:00pm-2:59pm" "3:00pm-3:59pm" ...
## $ MINUTE : int 58 18 55 20 45 55 23 15 0 2 ...
## $ MINUTENAME : chr "58" "18" "55" "20" ...
## $ NHS : int 0 0 0 0 0 0 0 0 0 1 ...
## $ NHSNAME : chr "This section IS NOT on the NHS" "This section IS NOT on the NHS" "This section IS NOT on the NHS" "This section IS NOT on the NHS" ...
## $ ROUTE : int 4 6 3 4 4 3 4 4 4 2 ...
## $ ROUTENAME : chr "County Road" "Local Street - Municipality" "State Highway" "County Road" ...
## $ TWAY_ID : chr "cr-4" "martin luther king jr dr" "sr-76" "CR-ALEXANDRIA WELLINGTON RD" ...
## $ TWAY_ID2 : chr "" "" "us-280" "" ...
## $ RUR_URB : int 1 2 1 1 1 1 2 1 1 1 ...
## $ RUR_URBNAME : chr "Rural" "Urban" "Rural" "Rural" ...
## $ FUNC_SYS : int 5 4 4 7 5 4 4 5 5 3 ...
## $ FUNC_SYSNAME: chr "Major Collector" "Minor Arterial" "Minor Arterial" "Local" ...
## $ RD_OWNER : int 2 4 1 2 2 1 4 2 2 1 ...
## $ RD_OWNERNAME: chr "County Highway Agency" "City or Municipal Highway Agency" "State Highway Agency" "County Highway Agency" ...
## $ MILEPT : int 0 0 49 0 0 390 0 0 0 3019 ...
## $ MILEPTNAME : chr "None" "None" "49" "None" ...
## $ LATITUDE : num 32.4 33.5 33.3 33.8 32.8 ...
## $ LATITUDENAME: chr "32.43313333" "33.48465833" "33.29994167" "33.79507222" ...
## $ LONGITUD : num -86.1 -86.8 -86.4 -85.9 -86.1 ...
## $ LONGITUDNAME: chr "-86.09485" "-86.83954444" "-86.36964167" "-85.88348611" ...
## $ SP_JUR : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SP_JURNAME : chr "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" ...
## $ HARM_EV : int 42 12 34 42 42 12 8 12 12 12 ...
## $ HARM_EVNAME : chr "Tree (Standing Only)" "Motor Vehicle In-Transport" "Ditch" "Tree (Standing Only)" ...
## $ MAN_COLL : int 0 6 0 0 0 2 0 1 1 2 ...
## $ MAN_COLLNAME: chr "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" "Angle" "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" ...
## $ RELJCT1 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ RELJCT1NAME : chr "No" "No" "No" "No" ...
## $ RELJCT2 : int 1 1 3 1 1 1 3 1 8 1 ...
## $ RELJCT2NAME : chr "Non-Junction" "Non-Junction" "Intersection-Related" "Non-Junction" ...
## $ TYP_INT : int 1 1 3 1 1 1 2 1 1 1 ...
## $ TYP_INTNAME : chr "Not an Intersection" "Not an Intersection" "T-Intersection" "Not an Intersection" ...
## $ WRK_ZONE : int 0 0 0 0 0 0 0 0 0 0 ...
## $ WRK_ZONENAME: chr "None" "None" "None" "None" ...
## $ REL_ROAD : int 4 1 4 4 4 1 1 1 1 1 ...
## $ REL_ROADNAME: chr "On Roadside" "On Roadway" "On Roadside" "On Roadside" ...
## $ LGT_COND : int 2 3 1 1 2 2 3 1 2 1 ...
## $ LGT_CONDNAME: chr "Dark - Not Lighted" "Dark - Lighted" "Daylight" "Daylight" ...
## $ WEATHER : int 1 2 2 10 2 1 1 1 10 10 ...
## $ WEATHERNAME : chr "Clear" "Rain" "Rain" "Cloudy" ...
## $ SCH_BUS : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SCH_BUSNAME : chr "No" "No" "No" "No" ...
## $ RAIL : chr "0000000" "0000000" "0000000" "0000000" ...
## $ RAILNAME : chr "Not Applicable" "Not Applicable" "Not Applicable" "Not Applicable" ...
## $ NOT_HOUR : int 99 17 14 99 0 17 19 7 20 10 ...
## $ NOT_HOURNAME: chr "Unknown" "5:00pm-5:59pm" "2:00pm-2:59pm" "Unknown" ...
## $ NOT_MIN : int 99 18 58 99 45 0 23 21 0 3 ...
## $ NOT_MINNAME : chr "Unknown" "18" "58" "Unknown" ...
## $ ARR_HOUR : int 3 17 15 99 0 17 19 7 20 10 ...
## $ ARR_HOURNAME: chr "3:00am-3:59am" "5:00pm-5:59pm" "3:00pm-3:59pm" "Unknown EMS Scene Arrival Hour" ...
## $ ARR_MIN : int 10 26 15 99 55 19 29 28 10 7 ...
## $ ARR_MINNAME : chr "10" "26" "15" "Unknown EMS Scene Arrival Minutes" ...
## $ HOSP_HR : int 99 99 99 99 88 18 88 88 99 10 ...
## $ HOSP_HRNAME : chr "Unknown" "Unknown" "Unknown" "Unknown" ...
## $ HOSP_MN : int 99 99 99 99 88 51 88 88 99 29 ...
## $ HOSP_MNNAME : chr "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" ...
## $ FATALS : int 3 1 1 1 1 1 1 1 1 1 ...
## $ DRUNK_DR : int 1 0 0 0 0 0 0 0 0 0 ...Vemos que tenemos 81 variables y 35766 registros
Revisamos la descripción de las variables contenidas en el fichero y si los tipos de variables se corresponden con las que hemos cargado. Las organizamos lógicamente para darles sentido y construimos un pequeño diccionario de datos utilizando la documentación auxiliar.
HECHOS A ESTUDIAR
DIMENSIÓN GEOGRÁFICA
DIMENSIÓN TEMPORAL
DIMENSIÓN CONDICICIONES ACCIDENTE
DIMENSIÓN METEOROLOGIA
OTROS
DIMENSIÓN SERVICIO EMERGENCIAS
DIMENSIÓN FACTORES RELACIONADOS ACCIDENTE
El siguiente paso será la limpieza de datos, mirando si hay valores vacíos o nulos.
print('NA')## [1] "NA"colSums(is.na(accidentData))## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS PVH_INVL
## 0 0 0 0 0 0
## PEDS PERSONS PERMVIT PERNOTMVIT COUNTY COUNTYNAME
## 0 0 0 0 0 0
## CITY CITYNAME DAY DAYNAME MONTH MONTHNAME
## 0 0 0 0 0 0
## YEAR DAY_WEEK DAY_WEEKNAME HOUR HOURNAME MINUTE
## 0 0 0 0 0 0
## MINUTENAME NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 0 0 0 0 0 0
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER
## 0 0 0 0 0 0
## RD_OWNERNAME MILEPT MILEPTNAME LATITUDE LATITUDENAME LONGITUD
## 0 0 0 0 0 0
## LONGITUDNAME SP_JUR SP_JURNAME HARM_EV HARM_EVNAME MAN_COLL
## 0 0 0 0 0 0
## MAN_COLLNAME RELJCT1 RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT
## 0 0 0 0 0 0
## TYP_INTNAME WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 0 0 0 0 0 0
## LGT_CONDNAME WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 0 0 0 0 0 0
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME ARR_HOUR
## 0 0 0 0 0 0
## ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR HOSP_HRNAME HOSP_MN
## 0 0 0 0 0 0
## HOSP_MNNAME FATALS DRUNK_DR
## 0 0 0print('Blancos')## [1] "Blancos"colSums(accidentData=="")## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS PVH_INVL
## 0 0 0 0 0 0
## PEDS PERSONS PERMVIT PERNOTMVIT COUNTY COUNTYNAME
## 0 0 0 0 0 0
## CITY CITYNAME DAY DAYNAME MONTH MONTHNAME
## 0 0 0 0 0 0
## YEAR DAY_WEEK DAY_WEEKNAME HOUR HOURNAME MINUTE
## 0 0 0 0 0 0
## MINUTENAME NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 0 0 0 0 0 0
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER
## 26997 0 0 0 0 0
## RD_OWNERNAME MILEPT MILEPTNAME LATITUDE LATITUDENAME LONGITUD
## 0 0 0 0 0 0
## LONGITUDNAME SP_JUR SP_JURNAME HARM_EV HARM_EVNAME MAN_COLL
## 0 0 0 0 0 0
## MAN_COLLNAME RELJCT1 RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT
## 0 0 0 0 0 0
## TYP_INTNAME WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 0 0 0 0 0 0
## LGT_CONDNAME WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 0 0 0 0 0 0
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME ARR_HOUR
## 0 0 0 0 0 0
## ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR HOSP_HRNAME HOSP_MN
## 0 0 0 0 0 0
## HOSP_MNNAME FATALS DRUNK_DR
## 0 0 0Vemos que no hay valores nulos en los datos. También verificamos si existen campos llenos de espacios en blanco. En este caso sí encontramos el campo TWAY_ID2 con 26997 valores en blanco. Valoramos no hacer ninguna acción de eliminar registros puesto que este campo no lo utilizaremos.
Vamos a crear histogramas y describir los valores para ver los datos en general de estos atributos para hacer una primera aproximación a los datos:
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')
if(!require('Rmisc')) install.packages('Rmisc'); library('Rmisc')
if(!require('dplyr')) install.packages('dplyr'); library('dplyr')
if(!require('xfun')) install.packages('xfun'); library('xfun')
summary(accidentData[c("FATALS","DRUNK_DR")])## FATALS DRUNK_DR
## Min. :1.000 Min. :0.0000
## 1st Qu.:1.000 1st Qu.:0.0000
## Median :1.000 Median :0.0000
## Mean :1.085 Mean :0.2664
## 3rd Qu.:1.000 3rd Qu.:1.0000
## Max. :8.000 Max. :4.0000histList<- list()
n = c("FATALS","DRUNK_DR")
accidentDataAux= accidentData %>% select(all_of(n))
for(y in 1:ncol(accidentDataAux)){
col <- names(accidentDataAux)[y]
ggp <- ggplot(accidentDataAux, aes_string(x = col)) +
geom_histogram(bins = 30, fill = "cornflowerblue", color = "black",ggtittle = "Contador de ocurrencias por variable")
histList[[y]] <- ggp # añadimos cada plot a la lista vacía
}
multiplot(plotlist = histList, coles = 1)
## [1] 1Observaciones:
Número de muertes: Todos los accidentes recogidos en estos datos reportan una muerte como mínimo. Siendo el accidente más grave con ocho víctimas y vemos que la distribución se acumula de forma muy evidente en una muerte por accidente.
Conductores bebidos involucrados en el accidente: Analizaremos con más detalle este dato más adelante para derivar un nuevo dato: Accidentes donde el alcohol está presente o no. De entrada, la media es de 0.26% de accidentes donde interviene un conductor bebido. La franja de conductores bebidos por accidente va de un conductor como mínimo a cuatro como máximo.
summary(accidentData[c("VE_TOTAL","VE_FORMS","PVH_INVL")])## VE_TOTAL VE_FORMS PVH_INVL
## Min. : 1.00 Min. : 1.000 Min. : 0.00000
## 1st Qu.: 1.00 1st Qu.: 1.000 1st Qu.: 0.00000
## Median : 1.00 Median : 1.000 Median : 0.00000
## Mean : 1.56 Mean : 1.517 Mean : 0.04269
## 3rd Qu.: 2.00 3rd Qu.: 2.000 3rd Qu.: 0.00000
## Max. :15.00 Max. :15.000 Max. :10.00000#Crearemos una lista para mostrar los atributos que interesan.
histList<- list()
n = c("VE_TOTAL","VE_FORMS","PVH_INVL")
accidentDataAux= accidentData %>% select(all_of(n))
for(y in 1:ncol(accidentDataAux)){
col <- names(accidentDataAux)[y]
ggp <- ggplot(accidentDataAux, aes_string(x = col)) +
geom_histogram(bins = 30, fill = "cornflowerblue", color = "black")
histList[[y]] <- ggp # añadimos cada plot a la lista vacía
}
multiplot(plotlist = histList, coles = 1)
## [1] 1Observaciones en cuanto a los vehículos implicados.
Número de vehículos implicados (VE_TOTAL) en total está en la franja de 1 hasta 59 siendo este el valor máximo y una media de 1.5. Número de vehículos en movimiento implicados (VE_FORMS), el valor más habitual es 1 con un valor máximo también de 59. Prevemos que hay un valor extremo que habrá que tratar para poder sacar más información de esta variable. Número de vehículos estacionados implicados (PVH_INVL): Por lo que respecta a esta variable lo habitual es que no haya vehículos estacionados en los incidentes recogidos en estos datos. Con todo aparecen casos aislados donde incluso había 10 coches estacionados.
summary(accidentData[c("PEDS","PERSONS","PERMVIT","PERNOTMVIT")])## PEDS PERSONS PERMVIT PERNOTMVIT
## Min. :0.0000 Min. : 0.000 Min. : 0.000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.: 1.000 1st Qu.: 1.000 1st Qu.:0.0000
## Median :0.0000 Median : 2.000 Median : 2.000 Median :0.0000
## Mean :0.2285 Mean : 2.173 Mean : 2.163 Mean :0.2387
## 3rd Qu.:0.0000 3rd Qu.: 3.000 3rd Qu.: 3.000 3rd Qu.:0.0000
## Max. :8.0000 Max. :61.000 Max. :61.000 Max. :9.0000#Crearemos una lista para mostrar los atributos que interesan.
histList<- list()
n = c("PEDS","PERSONS","PERMVIT","PERNOTMVIT")
accidentDataAux= accidentData %>% select(all_of(n))
for(y in 1:ncol(accidentDataAux)){
col <- names(accidentDataAux)[y]
ggp <- ggplot(accidentDataAux, aes_string(x = col)) +
geom_histogram(bins = 30, fill = "cornflowerblue", color = "black")
histList[[y]] <- ggp # añadimos cada plot a la lista vacía
}
multiplot(plotlist = histList, coles = 1)
## [1] 1Observaciones en cuanto a las personas implicadas en un accidente.
El número de peatones implicados (PEDS) es muy bajo siendo coherente con el tipo de vía que se estudia y dónde no es habitual que haya gente andando. Con todo el valor como media de 0.22 y máximo de 8 obliga a investigar más este dato. (PERSONS) El número de ocupantes de vehículo implicados se sitúa como media en 2.1 (PERMVIT) El número conductores y ocupantes de vehículos en circulación implicados tiene un valor de media de 2.1. Estas dos variables recogen la misma información, pero la cuantifican de diferente manera. El accidente con el mayor número de ocupantes es de 61 personas. Por lo que respecta al número peatones, ciclistas, personas en vehículos aparcados y otros (PERNOTMVIT) vemos que aumenta un poco la media respecto a peatón puesto que entendemos que se incluyen más casos.
Vamos a profundizar un poco en el tema de la relación del alcohol en los conductores y el número de accidentes.
accidentData$alcohol <- ifelse(accidentData$DRUNK_DR %in% c(0), 0, 1)
counts <- table(accidentData$alcohol)
barplot(prop.table(counts),col=c("green","red"), main="Accidentes con conductor bebido", legend.texto=c("No Alcohol","Sí Alcohol"),xlab ="Presencia Alcohol", ylab = "Porcentaje",ylim=c(0,0.8) )## Warning in plot.window(xlim, ylim, log = log, ...): "legend.texto" es un
## parámetro gráfico inválido## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "legend.texto" es un parámetro gráfico inválido## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "legend.texto" es un parámetro gráfico inválido## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "legend.texto"
## es un parámetro gráfico inválido
Vemos que porcentualmente, en la gran mayoría de accidentes, alrededor del 75% no hay presencia de alcohol en el conductor. Los conductores que dan positivo están alrededor de un 22%. Hemos buscado contrastar el dato con otros países y estarían en un valor central donde los valores extremos máximo por país superan el 50% y los mínimos están sobre el 10%
Observamos ahora como se distribuyen las muertes por accidente.
df1 <- accidentData %>%
group_by(accidentData$FATALS) %>%
dplyr::summarise(counts = n())
df1## # A tibble: 8 × 2
## `accidentData$FATALS` counts
## <int> <int>
## 1 1 33226
## 2 2 2154
## 3 3 289
## 4 4 71
## 5 5 20
## 6 6 4
## 7 7 1
## 8 8 1counts <- table(accidentData$FATALS)
barplot(prop.table(counts),col=("red"),ylim=c(0,0.99),main="Distribución de la mortalidad a los accidentes",xlab ="Número de muertos", ylab = "Porcentaje")
Vemos que la mayoría de los accidentes tienen como mínimo un muerto. Vamos ahora a relacionar mortalidad y alcohol.
counts <- table(accidentData$alcohol, accidentData$FATALS)
colors <- c("green", "red")
barplot(prop.table(counts), beside = TRUE, col = colors,
ylim = c(0, 1), axes = TRUE,
xlab = "Número de muertos",
ylab = "Porcentaje",
main = "Accidentes por muertes y conductores positivos en alcohol",
legend = c("No Alcohol", "Alcohol"),
fill = colors)## Warning in plot.window(xlim, ylim, log = log, ...): "fill" es un parámetro
## gráfico inválido## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "fill" es un parámetro gráfico inválido## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...): "fill"
## es un parámetro gráfico inválido## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "fill" es un
## parámetro gráfico inválido
Probaremos ahora si hay relación entre el estado donde ha pasado el accidente y el número de conductores bebidos. Filtramos los cinco estados donde hay más accidentes.
accidentDataST5=subset(accidentData, accidentData$STATENAME == "California" | accidentData$STATENAME == "Texas" | accidentData$STATENAME == "Florida" | accidentData$STATENAME == "Georgia" | accidentData$STATENAME == "North Carolina")
files=dim(accidentDataST5)[1]
ggplot(data=accidentDataST5[1:files,],aes(x=DRUNK_DR,fill=STATENAME))+geom_bar()+ggtitle("Relación entre conductor bebido y Estado")+labs(x="Número de conductores bebidos implicados")
Como reflexión este gráfico tiene que pasar por el filtro de percentuar el número de accidentes por estado y la población del estado para no sacar conclusiones apresuradas.
Veamos ahora como en un mismo gráfico de frecuencias podemos trabajar con 3 variables: FATALS, STATENAME y WEATHERNAME.
ggplot(data = accidentDataST5[1:files,],aes(x=FATALS,fill=STATENAME))+geom_bar(position="fill")+facet_wrap(~WEATHERNAME)+ggtitle("Número de muertes en accidente por Estado y clima")+labs(x="Número de muertes")
Esta gráfica está bien como mecánica de construcción, pero los resultados los ponemos en entredicho. Está bien como paso inicial, pero hay que profundizar mucho más.
Vamos a buscar las correlaciones en función de las muertes y unas variables elegidas que creemos que pueden ayudar a explicar el aumento de muertes por accidente:
DRUNK_DR conductores bebidos VE_TOTAL número de vehículos implicados en total VE_FORMS número de vehículos en movimiento implicados PVH_INVL número de vehículos estacionados implicados PEDS número de peatón implicados PERSONS número de ocupante de vehículo implicados PERMVIT número conductores y ocupantes implicados PERNOTMVIT número peatones, ciclistas… Cualquier cosa menos vehículo motorizado
# Utilizamos esta librería para usar la funcio multiplot()
if(!require('Rmisc')) install.packages('Rmisc'); library('Rmisc')
n = c("DRUNK_DR","VE_TOTAL","VE_FORMS","PVH_INVL","PEDS","PERSONS","PERMVIT","PERNOTMVIT")
accidentDataAux= accidentData %>% select(all_of(n))
histList2<- vector('list', ncol(accidentDataAux))
for(i in seq_along(accidentDataAux)){
message(i)
histList2[[i]]<-local({
i<-i
col <-log(accidentDataAux[[i]])
ggp<- ggplot(data = accidentDataAux, aes(x = accidentData$FATALS, y=col)) +
geom_point(color = "gray30") + geom_smooth(method = lm,color = "firebrick") +
theme_bw() + xlab("Muertes") + ylab(names(accidentDataAux)[i])
})
}
multiplot(plotlist = histList2, cols = 3)
Podemos ver que:
De forma general cualquier aumento en las variables elegidas implica un aumento de las muertes en el accidente.
El factor que hace aumentar más el número de víctimas son las variables relacionadas con los peatones y pasajeros de los coches involucrados en el accidente.
Utilizamos las columnas que nos interesa para hacer la matriz y la visualizaremos utilizando la función corrplot.
# https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
if(!require("corrplot")) install.packages("corrplot"); library("corrplot")
n = c("FATALS","DRUNK_DR","VE_TOTAL","VE_FORMS","PVH_INVL","PEDS","PERSONS","PERMVIT","PERNOTMVIT")
factores= accidentData %>% select(all_of(n))
res<-cor(factores)
corrplot(res,method="color",tl.col="black", tl.srt=30, order = "AOE",
number.cex=0.75,sig.level = 0.01, addCoef.col = "black")
No vemos que haya una correlación negativa significativa entre dos variables y sí una muy buena correlación ya previsible entre los peatones implicados y personas involucradas en el accidente que no van en coche (PEDS y PERNOTMVIT) Lo mismo podemos observar en cuanto al número de conductores y ocupantes implicados (PERMVIT) y el número de vehículos implicados en movimiento (VE_FORMS) o el total de vehículos (VE_TOTAL).
Vamos a probar si hay una correlación entre personas implicadas en el accidente y el número de muertes.
if (!require('tidyverse')) install.packages('tidyverse'); library('tidyverse')## Cargando paquete requerido: tidyverse## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'tidyr' was built under R version 4.4.3## Warning: package 'purrr' was built under R version 4.4.3## Warning: package 'forcats' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.4 ✔ tidyr 1.3.1
## ✔ readr 2.1.5## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::arrange() masks plyr::arrange()
## ✖ purrr::compact() masks plyr::compact()
## ✖ dplyr::count() masks plyr::count()
## ✖ dplyr::desc() masks plyr::desc()
## ✖ dplyr::failwith() masks plyr::failwith()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::id() masks plyr::id()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::mutate() masks plyr::mutate()
## ✖ dplyr::rename() masks plyr::rename()
## ✖ dplyr::summarise() masks plyr::summarise()
## ✖ dplyr::summarize() masks plyr::summarize()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorscor.test(x = accidentData$PERSONS, y = accidentData$FATALS, method = "kendall")##
## Kendall's rank correlation tau
##
## data: accidentData$PERSONS and accidentData$FATALS
## z = 53.008, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2558174ggplot(data = accidentData, aes(x = PERSONS, y = log(FATALS))) + geom_point(color = "gray30") + geom_smooth(color = "firebrick") + theme_bw() +ggtitle("Correlación entre personas implicadas en el accidente y número de muertes")## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
De la observación de este gráfico podemos concluir que efectivamente el número de muertes aumenta en función de las personas implicadas en un accidente pero que la correlación no es tan elevada ni continúa como se podía prever.
Si dos variables están altamente correlacionadas obviamente darán casi exactamente la misma información en un modelo de regresión, por ejemplo. Pero, al incluir las dos variables, en realidad estamos debilitando el modelo. No estamos añadiendo información incremental. En lugar de esto, estamos haciendo un modelo ruidoso. No es una buena idea.
Cómo hemos visto antes tenemos una correlación muy grande entre PEDS y PERNOTMVIT, por lo tanto, podríamos eliminar la columna de peatones (PEDS) y dejar el total de peatones y otros reflejado a PERNOTMVIT.
# accidentData$PEDS<- NULL
str(accidentData)## 'data.frame': 35766 obs. of 82 variables:
## $ STATE : int 1 1 1 1 1 1 1 1 1 1 ...
## $ STATENAME : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ ST_CASE : int 10001 10002 10003 10004 10005 10006 10007 10008 10009 10010 ...
## $ VE_TOTAL : int 1 4 2 1 1 2 1 2 2 2 ...
## $ VE_FORMS : int 1 4 2 1 1 2 1 2 2 2 ...
## $ PVH_INVL : int 0 0 0 0 0 0 0 0 0 0 ...
## $ PEDS : int 0 0 0 0 0 0 1 0 0 0 ...
## $ PERSONS : int 4 6 2 5 1 3 1 2 4 3 ...
## $ PERMVIT : int 4 6 2 5 1 3 1 2 4 3 ...
## $ PERNOTMVIT : int 0 0 0 0 0 0 1 0 0 0 ...
## $ COUNTY : int 51 73 117 15 37 103 73 25 45 95 ...
## $ COUNTYNAME : chr "ELMORE (51)" "JEFFERSON (73)" "SHELBY (117)" "CALHOUN (15)" ...
## $ CITY : int 0 350 0 0 0 0 330 0 0 1500 ...
## $ CITYNAME : chr "NOT APPLICABLE" "BIRMINGHAM" "NOT APPLICABLE" "NOT APPLICABLE" ...
## $ DAY : int 1 2 2 3 4 4 7 8 9 10 ...
## $ DAYNAME : int 1 2 2 3 4 4 7 8 9 10 ...
## $ MONTH : int 1 1 1 1 1 1 1 1 1 1 ...
## $ MONTHNAME : chr "January" "January" "January" "January" ...
## $ YEAR : int 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ DAY_WEEK : int 4 5 5 6 7 7 3 4 5 6 ...
## $ DAY_WEEKNAME: chr "Wednesday" "Thursday" "Thursday" "Friday" ...
## $ HOUR : int 2 17 14 15 0 16 19 7 20 10 ...
## $ HOURNAME : chr "2:00am-2:59am" "5:00pm-5:59pm" "2:00pm-2:59pm" "3:00pm-3:59pm" ...
## $ MINUTE : int 58 18 55 20 45 55 23 15 0 2 ...
## $ MINUTENAME : chr "58" "18" "55" "20" ...
## $ NHS : int 0 0 0 0 0 0 0 0 0 1 ...
## $ NHSNAME : chr "This section IS NOT on the NHS" "This section IS NOT on the NHS" "This section IS NOT on the NHS" "This section IS NOT on the NHS" ...
## $ ROUTE : int 4 6 3 4 4 3 4 4 4 2 ...
## $ ROUTENAME : chr "County Road" "Local Street - Municipality" "State Highway" "County Road" ...
## $ TWAY_ID : chr "cr-4" "martin luther king jr dr" "sr-76" "CR-ALEXANDRIA WELLINGTON RD" ...
## $ TWAY_ID2 : chr "" "" "us-280" "" ...
## $ RUR_URB : int 1 2 1 1 1 1 2 1 1 1 ...
## $ RUR_URBNAME : chr "Rural" "Urban" "Rural" "Rural" ...
## $ FUNC_SYS : int 5 4 4 7 5 4 4 5 5 3 ...
## $ FUNC_SYSNAME: chr "Major Collector" "Minor Arterial" "Minor Arterial" "Local" ...
## $ RD_OWNER : int 2 4 1 2 2 1 4 2 2 1 ...
## $ RD_OWNERNAME: chr "County Highway Agency" "City or Municipal Highway Agency" "State Highway Agency" "County Highway Agency" ...
## $ MILEPT : int 0 0 49 0 0 390 0 0 0 3019 ...
## $ MILEPTNAME : chr "None" "None" "49" "None" ...
## $ LATITUDE : num 32.4 33.5 33.3 33.8 32.8 ...
## $ LATITUDENAME: chr "32.43313333" "33.48465833" "33.29994167" "33.79507222" ...
## $ LONGITUD : num -86.1 -86.8 -86.4 -85.9 -86.1 ...
## $ LONGITUDNAME: chr "-86.09485" "-86.83954444" "-86.36964167" "-85.88348611" ...
## $ SP_JUR : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SP_JURNAME : chr "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" ...
## $ HARM_EV : int 42 12 34 42 42 12 8 12 12 12 ...
## $ HARM_EVNAME : chr "Tree (Standing Only)" "Motor Vehicle In-Transport" "Ditch" "Tree (Standing Only)" ...
## $ MAN_COLL : int 0 6 0 0 0 2 0 1 1 2 ...
## $ MAN_COLLNAME: chr "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" "Angle" "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" ...
## $ RELJCT1 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ RELJCT1NAME : chr "No" "No" "No" "No" ...
## $ RELJCT2 : int 1 1 3 1 1 1 3 1 8 1 ...
## $ RELJCT2NAME : chr "Non-Junction" "Non-Junction" "Intersection-Related" "Non-Junction" ...
## $ TYP_INT : int 1 1 3 1 1 1 2 1 1 1 ...
## $ TYP_INTNAME : chr "Not an Intersection" "Not an Intersection" "T-Intersection" "Not an Intersection" ...
## $ WRK_ZONE : int 0 0 0 0 0 0 0 0 0 0 ...
## $ WRK_ZONENAME: chr "None" "None" "None" "None" ...
## $ REL_ROAD : int 4 1 4 4 4 1 1 1 1 1 ...
## $ REL_ROADNAME: chr "On Roadside" "On Roadway" "On Roadside" "On Roadside" ...
## $ LGT_COND : int 2 3 1 1 2 2 3 1 2 1 ...
## $ LGT_CONDNAME: chr "Dark - Not Lighted" "Dark - Lighted" "Daylight" "Daylight" ...
## $ WEATHER : int 1 2 2 10 2 1 1 1 10 10 ...
## $ WEATHERNAME : chr "Clear" "Rain" "Rain" "Cloudy" ...
## $ SCH_BUS : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SCH_BUSNAME : chr "No" "No" "No" "No" ...
## $ RAIL : chr "0000000" "0000000" "0000000" "0000000" ...
## $ RAILNAME : chr "Not Applicable" "Not Applicable" "Not Applicable" "Not Applicable" ...
## $ NOT_HOUR : int 99 17 14 99 0 17 19 7 20 10 ...
## $ NOT_HOURNAME: chr "Unknown" "5:00pm-5:59pm" "2:00pm-2:59pm" "Unknown" ...
## $ NOT_MIN : int 99 18 58 99 45 0 23 21 0 3 ...
## $ NOT_MINNAME : chr "Unknown" "18" "58" "Unknown" ...
## $ ARR_HOUR : int 3 17 15 99 0 17 19 7 20 10 ...
## $ ARR_HOURNAME: chr "3:00am-3:59am" "5:00pm-5:59pm" "3:00pm-3:59pm" "Unknown EMS Scene Arrival Hour" ...
## $ ARR_MIN : int 10 26 15 99 55 19 29 28 10 7 ...
## $ ARR_MINNAME : chr "10" "26" "15" "Unknown EMS Scene Arrival Minutes" ...
## $ HOSP_HR : int 99 99 99 99 88 18 88 88 99 10 ...
## $ HOSP_HRNAME : chr "Unknown" "Unknown" "Unknown" "Unknown" ...
## $ HOSP_MN : int 99 99 99 99 88 51 88 88 99 29 ...
## $ HOSP_MNNAME : chr "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" "Unknown EMS Hospital Arrival Time" ...
## $ FATALS : int 3 1 1 1 1 1 1 1 1 1 ...
## $ DRUNK_DR : int 1 0 0 0 0 0 0 0 0 0 ...
## $ alcohol : num 1 0 0 0 0 0 0 0 0 0 ...Seguidamente vayamos a asignar un 1 por accidentes que se producen de madrugada (01h a 06h en invierno) y un 0 para el resto de franja horaria, es decir, vamos a categorizar la variable HOUR y así tendremos una variable numérica que nos permitirá trabajar mejor en el futuro. La denominaremos madrugada. Después la utilizaremos para ver cómo se distribuyen los accidentes en las dos franjas horarias. Debemos tener en cuenta que la hora incluye el código 99 qué quiere decir que la hora no está informada. Miraremos de filtrar los registros con este valor para excluirlos.
accidentDataAux=subset(accidentData, accidentData$HOUR <= 24)
accidentData$madrugada <- NA
accidentData$madrugada[accidentDataAux$HOUR >=1 & accidentDataAux$HOUR <= 6] <- 1
accidentData$madrugada[accidentDataAux$HOUR ==0 | accidentDataAux$HOUR >6 ] <- 0
counts <- table(accidentData$madrugada)
barplot(prop.table(counts),col=c("green","red"),legend.texto=c("Resto del día","Madrugada"),ylim=c(0,1), main="Distribución de accidentes la madrugada y resto del día",xlab="0 Resto del día 1 Madrugada",ylab="Porcentaje" )## Warning in plot.window(xlim, ylim, log = log, ...): "legend.texto" es un
## parámetro gráfico inválido## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "legend.texto" es un parámetro gráfico inválido## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "legend.texto" es un parámetro gráfico inválido## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "legend.texto"
## es un parámetro gráfico inválido
Ahora añadiremos un campo nuevo a los datos. Este campo contendrá el valor de la hora del accidente discretizada con un método simple de intervalos de igual amplitud.
summary(accidentDataAux[,"HOUR"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.00 15.00 13.19 19.00 23.00Discretizamos con intervalos. Los criterios de corte están cogidos de la Web del Parlamento de Cataluña.
accidentDataAux["segmento_horario"] <- cut(accidentDataAux$HOUR, breaks = c(0,4,11,14,18,22), labels = c("Madrugada", "Mañana", "Mediodía", "Anochecer","Noche"))Observamos los datos discretizados y construimos un gráfico para analizar cómo se agrupan los accidentes.
head(accidentDataAux$segmento_horario)## [1] Madrugada Anochecer Mediodía Anochecer <NA> Anochecer
## Levels: Madrugada Mañana Mediodía Anochecer Nocheplot(accidentDataAux$segmento_horario,main="Número de accidentes por segmento horario",xlab="Segmento horario", ylab="Cantidad",col = "ivory")
Ahora vamos a discretizar la variable que contiene el número de vehículos implicados en un accidente (VE_TOTALS) puesto que era una de las variables en que las distancias entre sus valores eran muy grandes:
# Utilizaremos la función discretize de arules: This function implements several basic unsupervised methods to convert a continuous variable into a categorical variable (factor) using different binning strategies.
# https://cran.r-project.org/web/packages/arules/index.html
if (!require('arules')) install.packages('arules'); library('arules')## Cargando paquete requerido: arules## Warning: package 'arules' was built under R version 4.4.3## Cargando paquete requerido: Matrix##
## Adjuntando el paquete: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Adjuntando el paquete: 'arules'## The following object is masked from 'package:dplyr':
##
## recode## The following objects are masked from 'package:base':
##
## abbreviate, writeset.seed(2)
table(discretize(accidentData$VE_TOTAL, "cluster" ))##
## [1,1.57) [1.57,3.38) [3.38,15]
## 19972 14972 822hist(accidentData$VE_TOTAL, main="Número de accidentes por vehículos implicados con kmeans",xlab="Vehículos implicados", ylab="Cantidad",col = "ivory")
abline(v=discretize(accidentData$VE_TOTAL, method="cluster", onlycuts=TRUE),col="red")
Podemos observar que sin pasar ningún argumento y permitiendo que el algoritmo elija el conjunto de particiones se muestran tres clústeres que agrupan los vehículos implicados en las franjas mencionadas. Podemos asignar el propio clúster como una variable más al dataset para trabajar después.
accidentData$VE_TOTAL_KM<- (discretize(accidentData$VE_TOTAL, "cluster" ))
head(accidentData$VE_TOTAL_KM)## [1] [1,1.5) [2.73,15] [1.5,2.73) [1,1.5) [1,1.5) [1.5,2.73)
## Levels: [1,1.5) [1.5,2.73) [2.73,15]Ahora normalizaremos el número de muertes por el máximo añadiendo un nuevo campo a los datos que contendrá el valor.
accidentData$FATALS_NM<- (accidentData$FATALS/max(accidentData[,"FATALS"]))
head(accidentData$FATALS_NM)## [1] 0.375 0.125 0.125 0.125 0.125 0.125Supongamos que queremos normalizar por la diferencia para ubicar entre 0 y 1 la variable del número de muertes del accidente dado que el algoritmo de minería que utilizaremos así lo requiere. Observamos la distribución de la variable original y las generadas.
accidentData$FATALS_ND = (accidentData$FATALS-min(accidentData$FATALS))/(max(accidentData$FATALS)-min(accidentData$FATALS))
max(accidentData$FATALS)## [1] 8min(accidentData$FATALS)## [1] 1hist(accidentData$FATALS,xlab="Muertos", col="ivory",ylab="Cantidad", main="Número de muertes en accidente")
hist(accidentData$FATALS_NM,xlab="Muertos",ylab="Cantidad",col="ivory", main="Muertos normalizado por el máximo")
hist(accidentData$FATALS_ND,xlab="Muertos",ylab="Cantidad", col="ivory", main="Muertos normalizado por la diferencia")
A continuación, vamos a normalizar las otras columnas para asegurarnos que cada variable contribuye por igual en nuestro análisis.
# Definimos la función de normalización
nor <-function(x) { (x -min(x))/(max(x)-min(x))}
# Guardamos un nuevo dataset normalizado
accidentData$type<- NULL
n = c("FATALS","DRUNK_DR","VE_TOTAL","VE_FORMS","PVH_INVL","PEDS","PERSONS","PERMVIT","PERNOTMVIT")
accidentData<- accidentData %>% select(all_of(n))
accidentData_nor <- as.data.frame(lapply(accidentData, nor))
head(accidentData_nor)## FATALS DRUNK_DR VE_TOTAL VE_FORMS PVH_INVL PEDS PERSONS PERMVIT
## 1 0.2857143 0.25 0.00000000 0.00000000 0 0 0.06557377 0.06557377
## 2 0.0000000 0.00 0.21428571 0.21428571 0 0 0.09836066 0.09836066
## 3 0.0000000 0.00 0.07142857 0.07142857 0 0 0.03278689 0.03278689
## 4 0.0000000 0.00 0.00000000 0.00000000 0 0 0.08196721 0.08196721
## 5 0.0000000 0.00 0.00000000 0.00000000 0 0 0.01639344 0.01639344
## 6 0.0000000 0.00 0.07142857 0.07142857 0 0 0.04918033 0.04918033
## PERNOTMVIT
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0Tanto el análisis de componentes principales, principal componente analysis (PCA) en inglés, como la descomposición de valores singulares, singular value decomposition (SVD) en inglés, son técnicas que nos permitan trabajar con nuevas características llamadas componentes, que ciertamente son independientes entre sí. En realidad, estas dos técnicas nos permiten representar el juego de datos en un nuevo sistema de coordenadas que denominamos componentes principales. Este sistema está mejor adaptado a la distribución del juego de datos, de forma que recoge mejor su variabilidad.
Aplicamos el análisis de componentes principales al dataset. Empezamos ejecutando la función prcomp().
pca.acc <- prcomp(accidentData_nor)
summary(pca.acc)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 0.1168 0.08755 0.0690 0.0484 0.03269 0.02437 0.01072
## Proportion of Variance 0.4520 0.25389 0.1577 0.0776 0.03539 0.01967 0.00381
## Cumulative Proportion 0.4520 0.70584 0.8635 0.9411 0.97652 0.99619 1.00000
## PC8 PC9
## Standard deviation 1.143e-14 2.491e-15
## Proportion of Variance 0.000e+00 0.000e+00
## Cumulative Proportion 1.000e+00 1.000e+00Como se puede observar la función summary, nos devuelve la proporción de varianza aplicada al conjunto total de cada atributo. Gracias a esto, el atributo 1 explica el 0.452 de variabilidad del total de datos; en cambio, el atributo 7 explica solo el 0.000381.
A continuación, se muestra un histograma para ver el peso de cada atributo sobre el conjunto total de datos:
if (!require('factoextra')) install.packages('factoextra'); library('factoextra')## Cargando paquete requerido: factoextra## Warning: package 'factoextra' was built under R version 4.4.3## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa#Los valores propios corresponden a la cantidad de variación explicada por cada componente principal (PC).
ev= get_eig(pca.acc)
ev## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 1.364366e-02 4.519469e+01 45.19469
## Dim.2 7.664656e-03 2.538921e+01 70.58391
## Dim.3 4.760727e-03 1.576993e+01 86.35384
## Dim.4 2.342530e-03 7.759642e+00 94.11348
## Dim.5 1.068327e-03 3.538839e+00 97.65232
## Dim.6 5.937915e-04 1.966938e+00 99.61926
## Dim.7 1.149407e-04 3.807416e-01 100.00000
## Dim.8 1.305543e-28 4.324617e-25 100.00000
## Dim.9 6.203717e-30 2.054984e-26 100.00000fviz_eig(pca.acc)
En este ejercicio se decidió utilizar el método de Káiser para decidir cuales de las variables obtenidas serán escogidas. Este criterio mantendrá todas aquellas variables cuya varianza sea superior a 1.
# Calculamos la varianza de los componentes principales a partir de la desviación estándar
var_acc <- pca.acc$sdev^2
var_acc## [1] 1.364366e-02 7.664656e-03 4.760727e-03 2.342530e-03 1.068327e-03
## [6] 5.937915e-04 1.149407e-04 1.305543e-28 6.203717e-30Con los resultados obtenidos es muy complicado decidir cuáles son los componentes principales componentes a escoger. Este hecho podría estar causado por no haber escalado los datos previamente. Por lo tanto, el siguiente paso es escalar los datos y volver a calcular la varianza para ver qué datos selecciona.
# Escalamos los datos
acc_scale <- scale(accidentData_nor)
# Calculamos las componentes principales
pca.acc_scale <- prcomp(acc_scale)
# Mostramos la varianza de dichas variables:
var_acc_scale <- pca.acc_scale$sdev^2
head(var_acc_scale)## [1] 3.5632826 1.8581532 1.1186798 1.0185157 0.8423082 0.5603141Después de analizar la varianza y aplicando el criterio de Káiser nos quedaremos con los componentes principales 1,2,3 y 4 que son los superiores a 1. Este criterio tiene el problema de sobreestimar el número de factores, pero a pesar de ello es el que aplicaremos para analizar los resultados.
Mostramos el histograma de porcentaje de varianza explicado con los datos escalados:
fviz_eig(pca.acc_scale)
ev = get_eig(pca.acc_scale)
ev## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.563283e+00 3.959203e+01 39.59203
## Dim.2 1.858153e+00 2.064615e+01 60.23817
## Dim.3 1.118680e+00 1.242978e+01 72.66795
## Dim.4 1.018516e+00 1.131684e+01 83.98479
## Dim.5 8.423082e-01 9.358980e+00 93.34377
## Dim.6 5.603141e-01 6.225712e+00 99.56948
## Dim.7 3.874645e-02 4.305161e-01 100.00000
## Dim.8 3.781577e-26 4.201752e-25 100.00000
## Dim.9 2.016132e-26 2.240146e-25 100.00000Los valores propios se pueden utilizar para determinar el número de componentes principales a retener después de la PCA (Kaiser 1961):
Un valor propio > 1 indica que los PCs representan más varianza de la que representa una de las variables originales de los datos estandarizados. Esto se utiliza habitualmente como punto de corte para el cual se conservan los PCs. Esto solo es cierto cuando los datos están estandarizados.
También podemos limitar el número de componentes a este número que representa una determinada fracción de la varianza total. Por ejemplo, si estamos satisfecho con el 80% de la varianza total explicada, usamos el número de componentes para conseguirlo que son los 4 componentes principales vistos antes.
Continuamos con el análisis de los componentes principales. Después de aplicar el método Káiser se han seleccionado los 4 componentes principales.
var <- get_pca_var(pca.acc_scale)
var## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"Los componentes de get_pca_var() se pueden utilizar en el diagrama de variables de la siguiente manera:
#Utilizamos los 4 componentes principales encontrados antes
head(var$coord[,1:4],11)## Dim.1 Dim.2 Dim.3 Dim.4
## FATALS -0.30725642 0.07606945 0.28817514 -0.73967833
## DRUNK_DR -0.04921565 -0.38170799 -0.26454706 -0.58497489
## VE_TOTAL -0.83708697 0.28706095 -0.34390477 0.13727724
## VE_FORMS -0.86483356 0.18416941 -0.02398033 0.21866655
## PVH_INVL -0.06062570 0.30358100 -0.85844575 -0.18336810
## PEDS 0.48249027 0.82661079 0.13673775 -0.08353380
## PERSONS -0.88514448 0.21755791 0.19732665 -0.07582156
## PERMVIT -0.88736748 0.19749105 0.22352113 -0.06716534
## PERNOTMVIT 0.45872440 0.85355861 0.04773257 -0.10804664La calidad de representación de las variables en el mapa de factores se denomina cos2 (coseno cuadrado, coordenadas cuadradas). Podemos acceder al cos2 de la siguiente manera:
head(var$cos2[,1:4],11)## Dim.1 Dim.2 Dim.3 Dim.4
## FATALS 0.094406509 0.005786562 0.083044911 0.547124037
## DRUNK_DR 0.002422180 0.145700988 0.069985146 0.342195618
## VE_TOTAL 0.700714591 0.082403990 0.118270488 0.018845040
## VE_FORMS 0.747937080 0.033918370 0.000575056 0.047815060
## PVH_INVL 0.003675475 0.092161422 0.736929104 0.033623861
## PEDS 0.232796860 0.683285394 0.018697211 0.006977896
## PERSONS 0.783480750 0.047331446 0.038937806 0.005748910
## PERMVIT 0.787421038 0.039002717 0.049961697 0.004511183
## PERNOTMVIT 0.210428078 0.728562294 0.002278398 0.011674076corrplot(var$cos2[,1:4], is.corre=FALSE)## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "is.corre" es un parámetro gráfico inválido## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "is.corre" es un parámetro gráfico inválido## Warning in title(title, ...): "is.corre" es un parámetro gráfico inválido
También es posible crear un diagrama de barras de variables cos2 mediante la función fviz_cos2():
fviz_cos2(pca.acc_scale, choice = "var", axes = 1:2)
Un cos2 elevado indica una buena representación de la variable en el componente principal. En este caso, la variable se coloca cerca de la circunferencia del círculo de correlación.
Un cos2 bajo indica que la variable no está perfectamente representada por los PC. En este caso, la variable está cerca del centro del círculo.
Para una variable dada, la suma del cos2 de todos los componentes principales es igual a uno.
Si una variable está perfectamente representada por solo dos componentes principales (Dim.1 y Dim.2), la suma del cos2 en estos dos PCs es igual a uno. En este caso las variables se colocarán en el círculo de correlaciones.
Para algunas de las variables, pueden ser necesarios más de 2 componentes para representar perfectamente los datos. En este caso las variables se sitúan dentro del círculo de correlaciones.
En resumen:
fviz_pca_var(pca.acc_scale,
col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
Las contribuciones de las variables en la contabilización de la variabilidad de un determinado componente principal se expresan en porcentaje.
Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importantes para explicar la variabilidad en el conjunto de datos.
Las variables que no están correlacionadas con ningún PC o con las últimas dimensiones son variables con una contribución baja y se pueden eliminar para simplificar el análisis global.
La contribución de las variables se puede extraer de la siguiente manera:
head(var$contrib[,1:4],11)## Dim.1 Dim.2 Dim.3 Dim.4
## FATALS 2.64942527 0.3114147 7.42347448 53.7177824
## DRUNK_DR 0.06797608 7.8411720 6.25604798 33.5974816
## VE_TOTAL 19.66486180 4.4347254 10.57232696 1.8502454
## VE_FORMS 20.99011424 1.8253807 0.05140488 4.6945826
## PVH_INVL 0.10314857 4.9598398 65.87489045 3.3012610
## PEDS 6.53321358 36.7722855 1.67136397 0.6851044
## PERSONS 21.98761218 2.5472306 3.48069265 0.5644400
## PERMVIT 22.09819245 2.0990044 4.46613018 0.4429174
## PERNOTMVIT 5.90545583 39.2089469 0.20366845 1.1461852Cuando más grande sea el valor de la contribución, más contribución habrá al componente.
corrplot(var$contrib[,1:4], is.cor=FALSE)
Las variables más importantes (que más contribuyen) se pueden resaltar a la gráfica de correlación de la siguiente manera:
fviz_pca_var(pca.acc_scale, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
Las variables correlacionadas positivas apuntan al mismo lado de la trama. Las variables correlacionadas negativas apuntan a lados opuestos del gráfico. Por ejemplo, vemos que las personas involucradas en un accidente (PVH_INVL) y conductor bebido (DRUNK_DR) apuntan direcciones opuestas por tanto no están nada correlacionas, además lo hemos visto antes puesto que tienen un coeficiente de correlación de -0.01.
Se observa que las variables que más aportan a las componentes principales son PEDS y PERNOTMVIT por un lado y VE_TOTAL, VE_FORMS, PERSONS y PERMVIT por otro. Esto es debido al hecho que están correlacionadas. En concreto por el diagrama de correlación de antes de que PEDS está muy bien correlacionada con PERNOTMVIT. De otra banda VE_TOTAL, VE_FORMS, PERSONS y PERMVIT están también bastante correlacionadas. La correlación de FATALS con este grupo de variables no es elevada, pero apunta en la misma dirección.
Podrían ahora rehacer los componentes excluyendo las variables que no aportan información. Una vez rehechas estas nuevas variables sustituyen a las originales que las forman y se podrían utilizar por ejemplo como un indicador de gravedad de accidente puesto que incluye vehículo en movimiento, parado, peatones, conductores y otros implicados en una sola variable.
Los datos estudiados contemplan accidentes de tráfico con víctimas en las redes de autopistas en los EUUU a largo del 2020. Todos los registros tienen un identificador único de accidente y una serie de hechos principales como número de muertos, número de conductores bebidos, vehículos y personas implicadas. Tenemos que añadir otras variables que los caracterizan agrupadas por ubicación geográfica, temporal, condiciones específicas del accidente, meteorológicas, la intervención del servicio de emergencias y otros factores.
Revisados los datos parecen bien informados. Los datos están bastante limpios y bien documentados. No plantean graves problemas de campos con valores nulos o vacíos y tienen bastante potencial para generar nuevos indicadores a partir de los datos.
Podemos afirmar que a lo largo del 2020 en las autopistas de EE. UU. sucedieron 35766 accidentes en los que perdieron la vida 38824 personas. Pretendiamos extraer relaciones entre la presencia de alcohol en los conductores y el número de accidentes, pero las conclusiones no fueron claras. Las relaciones más obvias comprobadas son el incremento de muertes en función del incremento del número de vehículos, pasajeros y peatones implicados.
Habría que profundizar mucho más. Sí que podemos perfilar cómo son los accidentes típicos en cuanto al número de vehículos y personas, conductores o peatones implicados.
El más habitual es un muerto por accidente incrementándose este valor en función de las variables relacionadas. Los conductores bebidos aparecen en uno de cada cuatro accidentes mortales aproximadamente. Los vehículos implicados en los accidentes son típicamente uno pudiéndose incrementar a dos en los casos más típicos. El número de peatones implicados es relativamente bajo dado el tipo de vía que estamos estudiando.
En cuanto al consumo de alcohol, con el grado de profundidad estudiado, no se observa un estado donde las proporcionalidades del número de conductor con presencia de alcohol sean superiores a otros estados.
Estudiando el número de muertes en accidente en relación con el estado donde ha sucedido y la condición climática, vemos necesario profundizar con técnicas por ejemplo de agregación para ver cómo se agrupan y poder obtener un perfil.
Se ha estudiado la franja horaria de la madrugada para ver si acumula un mayor número de accidentes, no siendo así. Parece que el número de accidentes mantiene también proporcionalidad respecto las franjas horarias. Se ha estudiado el número de accidentes por segmento horario con una discretización fijada en intervalo arbitrarios. La mayor presencia de accidentes en horario por la mañana y anochecer (ir y volver al trabajo) hace pensar en que queda pendiente estudiar dado el tipo de vía la distribución horaria de los accidentes lunes a viernes respecto a los fines de semana y festivos para ver si hay horas donde se acumulan más accidentes mortales.
Finalmente, con la técnica de los componentes principales hemos generado una nueva variable que combina otras variables con una correlación inicial que se podría considerar como índice de gravedad del accidente.
A partir de un juego de datos que te parezca interesante, propón un
proyecto completo de minería de datos. La organización de la respuesta
tiene que coincidir en las fases típicas del ciclo de vida de un
proyecto de minería de datos.
No hay que realizar las tareas de la fase, este es un ejercicio
teórico.
Se espera por lo tanto que se responda de forma argumentada a las siguientes preguntas (Metodología CRISP-DM):
Para cada fase indica cuál es el objetivo de la fase y el producto que se obtendrá. Utiliza ejemplos de qué y cómo podrían ser las tareas. Si hay alguna característica que hace diferente el ciclo de vida de un proyecto de minería respecto a otros proyectos indícalo.
Extensión mínima: 400 palabras
Para la realización de esta actividad se ha utilizado una base de datos de Kaggle, en la que se recopilan las ventas de chocolate de una determinada empresa. A partir de ahora la denominaré “ChocoDelight”. El dataset puede encontrarse en el siguiente link: https://www.kaggle.com/datasets/atharvasoundankar/chocolate-sales
A continuación, presento la propuesta de proyecto de minería de datos siguiendo la metodología CRISP-DM y sus fases.
- Comprensión del negocio
ChocoDelight comercializa una variedad de tipos de chocolate distinto: (negro, blanco y con leche). Estas ventas se realizan a través de distintos canales (como el retail y mayoristas). El objetivo de la compañía es el de ptimizar las ventas al mismo tiempo se aumenta el márgen sobre ventas. El negocio necesita responder las siguientes preguntas:
¿Qué clase de chocolate (Dark, Milk o White) es el que más factura en cada canal de clientes? ¿Es la demanda de chocolate estacional? Saber esto permitiría programar mejor las campañas de marketing. ¿Hay regiones geográficas o localizaciones más propensas al consumo de chocolate? Con estos interrogantes se pretende identificar patrones de compra que pueden servir de ayuda para orientar la estrategia comercial. Como diseñar campañas específicas para cada segmento de clientes, o poder planificar de antemano el stock en fechas y localizaciones clave.
- Comprensión de los datos
El dataset se conforma de 7 columnas:
Date: Fecha de la transacción de venta. Product Name: Nombre del producto de chocolate vendido (podrían existir varios nombres como “Dark Hazelnut Bar”, “Milk Truffle Bag”, etc.). Category: Tipo de chocolate (Dark, Milk, White). Quantity Sold: Cantidad de chocolates vendidos por transacción. Revenue: Ingreso total generado en la transacción (en moneda local). Customer Segment: Tipo de cliente (Retail, Wholesale). Location: Región o tienda donde se efectuó la venta.
Primeras hipotesis y verificaciones:
Se revisará si las fechas tienen un formato uniforme y coherente. Se revisará que los revenues están correctamente calculados, no existen valores negativos. Se verificará si existen o no valores ausentes (NA) en campos que son clave, como Quantity Sold o Revenue, pues estas podrían reflejar ventas no registradas o datos incompletos. El “product name” podría contener alguna variación lingüística o de escritura, lo que generaría clases adicionales y exigiría un proceso de normalización.
- Preparación de los datos
Integración y selección:
A pesar de que únicamente se dispone de un único archivo .csv, se podría cruzar esta información con datos externos, (festividades nacionales, o promociones), de este modo el análisis podría verse enriquecido.
Se seleccionarán las columnas relevantes para nuestros objetivos: puesto que se trata de un dataset muy acotado. Todas las columnas son relevantes y cada una de las variables pueden información no redundante.
Limpieza y selección:
Se deberán corregir los posibles pedidos duplicados debido a errores en la introducción Se deberán revisar posibles valores extremos en “Quantity sold” (por ejemplo 20.000 unidades en un único pedido) Se crearán variables derivadas a partir de “Date” como “Month” o “Year”, de este modo se modrá estudiar tendencias estacionales o mensuales. Posteriormente se normalizarán los datos si necesitamos sensibles a las escalas.
Codificación y discretización:
La variable location puede agruparse por regiones (Norte, Sur, etc)
- Modelado
Podrán llevarse a cabo distintos tipos de modelado:
Modelos de predicción de Ventas: mediante regresión lineal, o series temporales. Puede tratar de predecirse la facturación futura en función de la fecha.
Modelos de segmentación: algoritmos de clustering (k-means o DBSCAN) con tal de encontrar patrones de ventasimilares en diferentes localizaciones, o también “Category + customer segment”. Con ello se pretende encontrar nichos subexplotados.
Reglas de asociación: si de la propia empresa se extranjeras datos más detallados, como qué otros productos incluyen las transacciones, aplicando un apriori o FP-Growth podría estudiarse con qué tipos de artículos se suele vender el chocolate.
- Evaluación
Validación de los modelos: Se dividirá el conjunto de datos en uno de entreno y otro de test para medir la capacidad de generalización que tiene el modelo. Se realizará una prueba de validación cruzada. Se compararán diferentes técnicas (por ejemplo, la regresión lineal vs un random forest) eligiendo la que tenga mejor equilibrio entre su precisión y su interpretabilidad.
- Implementación
Despliegue e integración: El modelo final podría integrarse en una aplicación interna. De ese modo, el equipo de ventas podrá consultar las predicciones de la demanda y planificar los stocks en tiempo real. Se construirán dashboards o informes periódicos para observar la evolución de los KPIs. (Revenue total, unidades vendidas, márgenes, etc)
Mantenimiento y actualización Se definirán ciclos de refresco de datos, (diario o semanal) de tal modo que el modelo siempre se entrene con información nueva y las proyecciones de ventas se mantengan actualizadas. A través de la monitorización se estudiará que el modelo no sufra degradación a lo largo del tiempo.
Perspectivas a futuro Una de las posibles líneas de mejora es la incorporación de datos externos (como promociones y tendencias en redes sociales) para así afinar el modelo. Otra posibilidad es la introducción de datos de la competencias, estableciendo así un análisis de precios dinámicos.
A partir del juego de datos utilizado en el ejemplo de la PEC, realiza las tareas previas a la generación de un modelo de minería de datos explicadas en los módulos “El proceso de minería de datos” y “Preprocesado de los datos y gestión de características”.
Puedes utilizar de referencia el ejemplo de la PEC, pero procura cambiar el enfoque y analizar los datos en función de las diferentes dimensiones que presentan los datos. Así, no se puede utilizar la combinación de variables utilizada en el ejemplo: “FATALS”,“DRUNK_DR”,“VE_TOTAL”,“VE_FORMS”,“PVH_INVL”,“PEDS”,“PERSONS”,“PERMVIT”,“PERNOTMVIT”. Se debe analizar cualquier otra combinación que puede incluir (o no) algunas de estas variables con otras nuevas.
Opcionalmente y valorable se pueden añadir al estudio datos de otros años para realizar comparaciones temporales (https://www.nhtsa.gov/file-downloads?p=nhtsa/downloads/FARS/) o añadir otros hechos a estudiar relacionados, por ejemplo, el consumo de drogas en los accidentes (https://static.nhtsa.gov/nhtsa/downloads/FARS/2020/National/FARS2020NationalCSV.zip)
Para la resolución de este ejercicio, de la página oficial del FARS, descargaremos los registros de accidentes y de drogas de los años 2019, 2020 y 2021. De esta forma, obtenemos un componente temporal interesante, (comparaciones de antes y durante pandemia), al mismo tiempo que mantenemos un análisis acotado y más manejable y evitando un volumen de datos para procesar demasiado elevado.
Se cambia el nombre de accident.csv por accident2020.csv y se descarga el resto de archivos.
# Paso 1: Cargamos librerías necesarias
if(!require(dplyr)) install.packages("dplyr"); library(dplyr)
# Paso 2: Cargamos cada CSV de accident
accident_2019 <- read.csv("accident_2019.csv", row.names=NULL)
accident_2020 <- read.csv("accident_2020.csv", row.names=NULL)
accident_2021 <- read.csv("accident_2021.csv", row.names=NULL)
# Paso 3: Unificamos el tipo de las columnas conflictivas para poder ejecutar el bind_rows sin problemas
accident_2019 <- mutate(accident_2019,
LATITUDENAME = as.character(LATITUDENAME),
LONGITUDNAME = as.character(LONGITUDNAME))
accident_2020 <- mutate(accident_2020,
LATITUDENAME = as.character(LATITUDENAME),
LONGITUDNAME = as.character(LONGITUDNAME))
accident_2021 <- mutate(accident_2021,
LATITUDENAME = as.character(LATITUDENAME),
LONGITUDNAME = as.character(LONGITUDNAME))
# UNIFICAR LOS DATOS DE ACCIDENT
# Paso 4: Se unifican todos los datos de los 3 años
accident_all <- bind_rows(accident_2019, accident_2020, accident_2021)# Paso 4: Antes de Combinar accident_all y drugs_all, se comprueba la estructura y que las columnas clave State y Statename no tengan valores perdidos, ya que son clave foranea.
str(accident_all)## 'data.frame': 82338 obs. of 91 variables:
## $ STATE : int 1 1 1 1 1 1 1 1 1 1 ...
## $ STATENAME : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ ST_CASE : int 10001 10002 10003 10004 10005 10006 10007 10008 10009 10010 ...
## $ VE_TOTAL : int 2 2 3 1 1 2 1 1 1 1 ...
## $ VE_FORMS : int 2 2 3 1 1 2 1 1 1 1 ...
## $ PVH_INVL : int 0 0 0 0 0 0 0 0 0 0 ...
## $ PEDS : int 0 0 0 1 0 0 0 0 0 1 ...
## $ PERSONS : int 3 2 4 1 1 2 5 1 1 1 ...
## $ PERMVIT : int 3 2 4 1 1 2 5 1 1 1 ...
## $ PERNOTMVIT : int 0 0 0 1 0 0 0 0 0 1 ...
## $ COUNTY : int 81 55 29 55 3 85 55 69 3 101 ...
## $ COUNTYNAME : chr "LEE (81)" "ETOWAH (55)" "CLEBURNE (29)" "ETOWAH (55)" ...
## $ CITY : int 2340 1280 0 2562 0 0 0 790 1493 2130 ...
## $ CITYNAME : chr "OPELIKA" "GADSDEN" "NOT APPLICABLE" "RAINBOW CITY" ...
## $ DAY : int 7 23 22 22 18 7 2 9 8 7 ...
## $ DAYNAME : int 7 23 22 22 18 7 2 9 8 7 ...
## $ MONTH : int 2 1 1 1 1 1 1 2 2 2 ...
## $ MONTHNAME : chr "February" "January" "January" "January" ...
## $ YEAR : int 2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
## $ DAY_WEEK : int 5 4 3 3 6 2 4 7 6 5 ...
## $ DAY_WEEKNAME: chr "Thursday" "Wednesday" "Tuesday" "Tuesday" ...
## $ HOUR : int 12 18 19 3 5 12 9 21 7 5 ...
## $ HOURNAME : chr "12:00pm-12:59pm" "6:00pm-6:59pm" "7:00pm-7:59pm" "3:00am-3:59am" ...
## $ MINUTE : int 54 3 0 15 50 30 25 2 53 56 ...
## $ MINUTENAME : chr "54" "3" "0" "15" ...
## $ NHS : int 1 1 1 1 1 1 1 0 0 0 ...
## $ NHSNAME : chr "This section IS ON the NHS" "This section IS ON the NHS" "This section IS ON the NHS" "This section IS ON the NHS" ...
## $ ROUTE : int 1 1 1 1 1 1 1 4 6 6 ...
## $ ROUTENAME : chr "Interstate" "Interstate" "Interstate" "Interstate" ...
## $ TWAY_ID : chr "I-85" "I-759" "I-20" "I-59" ...
## $ TWAY_ID2 : chr "" "" "" "" ...
## $ RUR_URB : int 2 2 1 1 2 1 1 1 2 2 ...
## $ RUR_URBNAME : chr "Urban" "Urban" "Rural" "Rural" ...
## $ FUNC_SYS : int 1 1 1 1 1 1 1 7 4 5 ...
## $ FUNC_SYSNAME: chr "Interstate" "Interstate" "Interstate" "Interstate" ...
## $ RD_OWNER : int 1 1 1 1 1 1 1 2 4 4 ...
## $ RD_OWNERNAME: chr "State Highway Agency" "State Highway Agency" "State Highway Agency" "State Highway Agency" ...
## $ MILEPT : int 641 15 2118 1775 413 1567 1826 0 0 0 ...
## $ MILEPTNAME : chr "641" "15" "2118" "1775" ...
## $ LATITUDE : num 32.7 34 33.7 34 30.7 ...
## $ LATITUDENAME: chr "32.66622222" "33.99782778" "33.66084167" "33.95647222" ...
## $ LONGITUD : num -85.3 -86.1 -85.4 -86.1 -87.8 ...
## $ LONGITUDNAME: chr "-85.33665833" "-86.05399722" "-85.39101111" "-86.14052222" ...
## $ SP_JUR : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SP_JURNAME : chr "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" "No Special Jurisdiction" ...
## $ HARM_EV : int 12 12 12 8 1 12 32 1 3 8 ...
## $ HARM_EVNAME : chr "Motor Vehicle In-Transport" "Motor Vehicle In-Transport" "Motor Vehicle In-Transport" "Pedestrian" ...
## $ MAN_COLL : int 1 1 1 0 0 7 0 0 0 0 ...
## $ MAN_COLLNAME: chr "Front-to-Rear" "Front-to-Rear" "Front-to-Rear" "The First Harmful Event was Not a Collision with a Motor Vehicle in Transport" ...
## $ RELJCT1 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ RELJCT1NAME : chr "No" "No" "No" "No" ...
## $ RELJCT2 : int 1 1 1 1 1 1 1 1 3 1 ...
## $ RELJCT2NAME : chr "Non-Junction" "Non-Junction" "Non-Junction" "Non-Junction" ...
## $ TYP_INT : int 1 1 1 1 1 1 1 1 3 1 ...
## $ TYP_INTNAME : chr "Not an Intersection" "Not an Intersection" "Not an Intersection" "Not an Intersection" ...
## $ WRK_ZONE : int 0 0 0 0 0 0 0 0 0 0 ...
## $ WRK_ZONENAME: chr "None" "None" "None" "None" ...
## $ REL_ROAD : int 1 1 1 1 4 1 3 1 4 1 ...
## $ REL_ROADNAME: chr "On Roadway" "On Roadway" "On Roadway" "On Roadway" ...
## $ LGT_COND : int 1 2 2 2 2 1 1 2 1 3 ...
## $ LGT_CONDNAME: chr "Daylight" "Dark - Not Lighted" "Dark - Not Lighted" "Dark - Not Lighted" ...
## $ WEATHER1 : int 1 2 10 1 5 1 10 1 1 1 ...
## $ WEATHER1NAME: chr "Clear" "Rain" "Cloudy" "Clear" ...
## $ WEATHER2 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ WEATHER2NAME: chr "No Additional Atmospheric Conditions" "No Additional Atmospheric Conditions" "No Additional Atmospheric Conditions" "No Additional Atmospheric Conditions" ...
## $ WEATHER : int 1 2 10 1 5 1 10 1 1 1 ...
## $ WEATHERNAME : chr "Clear" "Rain" "Cloudy" "Clear" ...
## $ SCH_BUS : int 0 0 0 0 0 0 0 0 0 0 ...
## $ SCH_BUSNAME : chr "No" "No" "No" "No" ...
## $ RAIL : chr "0000000" "0000000" "0000000" "0000000" ...
## $ RAILNAME : chr "Not Applicable" "Not Applicable" "Not Applicable" "Not Applicable" ...
## $ NOT_HOUR : int 12 18 19 3 99 12 9 88 7 88 ...
## $ NOT_HOURNAME: chr "12:00pm-12:59pm" "6:00pm-6:59pm" "7:00pm-7:59pm" "3:00am-3:59am" ...
## $ NOT_MIN : int 59 3 2 4 99 30 21 88 53 88 ...
## $ NOT_MINNAME : chr "59" "3" "2" "4" ...
## $ ARR_HOUR : int 13 18 19 3 6 12 9 88 7 88 ...
## $ ARR_HOURNAME: chr "1:00pm-1:59pm" "6:00pm-6:59pm" "7:00pm-7:59pm" "3:00am-3:59am" ...
## $ ARR_MIN : int 9 7 12 11 0 54 24 88 59 88 ...
## $ ARR_MINNAME : chr "9" "7" "12" "11" ...
## $ HOSP_HR : int 13 99 20 88 88 88 99 88 88 88 ...
## $ HOSP_HRNAME : chr "1:00pm-1:59pm" "Unknown" "8:00pm-8:59pm" "Not Applicable (Not Transported)" ...
## $ HOSP_MN : int 27 99 5 88 88 88 99 88 88 88 ...
## $ HOSP_MNNAME : chr "27" "Unknown EMS Hospital Arrival Time" "5" "Not Applicable (Not Transported)" ...
## $ CF1 : int 0 0 14 0 0 0 0 0 0 0 ...
## $ CF1NAME : chr "None" "None" "Motor Vehicle struck by falling cargo,or something that came loose from or something that was set in motion by a vehicle" "None" ...
## $ CF2 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ CF2NAME : chr "None" "None" "None" "None" ...
## $ CF3 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ CF3NAME : chr "None" "None" "None" "None" ...
## $ FATALS : int 1 1 1 1 1 1 1 1 1 1 ...
## $ DRUNK_DR : int 1 0 0 0 1 0 0 1 0 0 ...print('NA')## [1] "NA"colSums(is.na(accident_all))## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS PVH_INVL
## 0 0 0 0 0 0
## PEDS PERSONS PERMVIT PERNOTMVIT COUNTY COUNTYNAME
## 0 0 0 0 0 0
## CITY CITYNAME DAY DAYNAME MONTH MONTHNAME
## 0 0 1 1 1 0
## YEAR DAY_WEEK DAY_WEEKNAME HOUR HOURNAME MINUTE
## 1 1 0 1 0 1
## MINUTENAME NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 0 1 0 1 0 0
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER
## 0 1 0 1 0 1
## RD_OWNERNAME MILEPT MILEPTNAME LATITUDE LATITUDENAME LONGITUD
## 0 1 0 1 1 1
## LONGITUDNAME SP_JUR SP_JURNAME HARM_EV HARM_EVNAME MAN_COLL
## 1 1 0 1 0 1
## MAN_COLLNAME RELJCT1 RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT
## 0 1 0 1 0 1
## TYP_INTNAME WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 0 1 0 1 0 1
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME WEATHER
## 0 48851 48851 48851 48851 1
## WEATHERNAME SCH_BUS SCH_BUSNAME RAIL RAILNAME NOT_HOUR
## 0 1 0 0 0 1
## NOT_HOURNAME NOT_MIN NOT_MINNAME ARR_HOUR ARR_HOURNAME ARR_MIN
## 0 1 0 1 0 1
## ARR_MINNAME HOSP_HR HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1
## 0 1 0 1 0 48851
## CF1NAME CF2 CF2NAME CF3 CF3NAME FATALS
## 48851 48851 48851 48851 48851 1
## DRUNK_DR
## 13085print('Blancos')## [1] "Blancos"colSums(accident_all=="")## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS PVH_INVL
## 0 0 0 0 0 0
## PEDS PERSONS PERMVIT PERNOTMVIT COUNTY COUNTYNAME
## 0 0 0 0 0 0
## CITY CITYNAME DAY DAYNAME MONTH MONTHNAME
## 0 0 NA NA NA 1
## YEAR DAY_WEEK DAY_WEEKNAME HOUR HOURNAME MINUTE
## NA NA 1 NA 1 NA
## MINUTENAME NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 1 NA 1 NA 1 1
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER
## 61541 NA 1 NA 1 NA
## RD_OWNERNAME MILEPT MILEPTNAME LATITUDE LATITUDENAME LONGITUD
## 1 NA 1 NA NA NA
## LONGITUDNAME SP_JUR SP_JURNAME HARM_EV HARM_EVNAME MAN_COLL
## NA NA 1 NA 1 NA
## MAN_COLLNAME RELJCT1 RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT
## 1 NA 1 NA 1 NA
## TYP_INTNAME WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 1 NA 1 NA 1 NA
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME WEATHER
## 1 NA NA NA NA NA
## WEATHERNAME SCH_BUS SCH_BUSNAME RAIL RAILNAME NOT_HOUR
## 1 NA 1 1 1 NA
## NOT_HOURNAME NOT_MIN NOT_MINNAME ARR_HOUR ARR_HOURNAME ARR_MIN
## 1 NA 1 NA 1 NA
## ARR_MINNAME HOSP_HR HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1
## 1 NA 1 NA 1 NA
## CF1NAME CF2 CF2NAME CF3 CF3NAME FATALS
## NA NA NA NA NA NA
## DRUNK_DR
## NA# UNIFICAR LOS DATOS DE DROGAS
# Paso 5: Se cargan los ficheros drugs de los años correspondientes
drugs_2019 <- read.csv("drugs_2019.csv")
drugs_2020 <- read.csv("drugs_2020.csv")
drugs_2021 <- read.csv("drugs_2021.csv")
# Paso 6: Se unifican todos los datos de drugs
drugs_all <- bind_rows(drugs_2019, drugs_2020, drugs_2021)# Paso 7: Se comprueba la estructura, y valores NA o "".
str(drugs_all)## 'data.frame': 325211 obs. of 9 variables:
## $ STATE : int 1 1 1 1 1 1 1 1 1 1 ...
## $ STATENAME : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ ST_CASE : int 10001 10001 10001 10002 10002 10002 10003 10003 10003 10003 ...
## $ VEH_NO : int 1 1 2 1 1 2 1 1 1 2 ...
## $ PER_NO : int 1 2 1 1 1 1 1 2 2 1 ...
## $ DRUGSPEC : int 0 0 0 1 1 0 1 1 1 1 ...
## $ DRUGSPECNAME: chr "Test Not Given" "Test Not Given" "Test Not Given" "Whole Blood" ...
## $ DRUGRES : int 0 0 0 401 417 0 1 600 605 1 ...
## $ DRUGRESNAME : chr "Test Not Given" "Test Not Given" "Test Not Given" "AMPHETAMINE" ...print('NA')## [1] "NA"colSums(is.na(drugs_all))## STATE STATENAME ST_CASE VEH_NO PER_NO DRUGSPEC
## 0 0 0 0 0 0
## DRUGSPECNAME DRUGRES DRUGRESNAME
## 0 0 0print('Blancos')## [1] "Blancos"colSums(drugs_all=="")## STATE STATENAME ST_CASE VEH_NO PER_NO DRUGSPEC
## 0 0 0 0 0 0
## DRUGSPECNAME DRUGRES DRUGRESNAME
## 0 0 0En este caso no existen valores Na ni ““.
# Paso 8: Creamos una variable 'drugs_involved' = 1 si hay al menos un resultado DRUGRES que indique positivo
# Según la codificación FARS, se indica que 0 es "test no administrado", "1 sería "negativo", y el resto son sustancias específicas. 996 sería "otras drogas".
# drugs_info sera un subsetting de drugs_all donde seleccionemos solo aquells casos donde hubo drogas involucradas.
drugs_info <- drugs_all %>%
group_by(STATE, ST_CASE) %>%
summarise(
drugs_involved = ifelse(any(!DRUGRES %in% c(0, 1)), 1, 0),
.groups = "drop" # para evitar el mensaje de agrupación
)
# COMBINAR TODA LA INFORMACIÓN DE ACCIDENT Y DRUGS
# Paso 9: Unimos a la tabla 'accident_all y drugs info'
accident_merged <- accident_all %>%
left_join(drugs_info, by = c("STATE", "ST_CASE"))
# Paso 10: Revisamos resultados
table(accident_merged$drugs_involved, useNA="ifany")##
## 0 1
## 12049 70289Se ha creado la variable drugs_involved como una variable booleana, en el que de 0 si no se ha realizado test o este ha dado negativo, y 1 en todos los que se haya detectado algún tipo de sustancia.
Se combinan las tablas accident_all con drugs_info “y no con drugs_all” porque si lo hiciéramos de ese modo, tendríamos múltiples filas para cada caso, ST_CASe, puesto que puede haber múltiples personas implicadas en el accidente.
Se comprobará qué drogas estaban más presentes directamente desde el dataset drugs_all.
Interesante comprobar que de todos los 82338 accidentes de accident_merged, en 70289 había presencia de drogas. El análisis preliminar con el alcohol podría haber llevado a conclusiones precipitadas hacia una falta de relación sobre el uso de drogas en la conducción.
# LIMPIEZA Y VERIFICACIÓN INICIAL
# Paso 11: Eliminación de duplicados (por si se ha repetido algún accidente)
accident_merged <- accident_merged %>% distinct()
# Paso 12: Valores ausentes de accident_merged
colSums(is.na(accident_merged))## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS
## 0 0 0 0 0
## PVH_INVL PEDS PERSONS PERMVIT PERNOTMVIT
## 0 0 0 0 0
## COUNTY COUNTYNAME CITY CITYNAME DAY
## 0 0 0 0 1
## DAYNAME MONTH MONTHNAME YEAR DAY_WEEK
## 1 1 0 1 1
## DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 0 1 0 1 0
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 1 0 1 0 0
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME
## 0 1 0 1 0
## RD_OWNER RD_OWNERNAME MILEPT MILEPTNAME LATITUDE
## 1 0 1 0 1
## LATITUDENAME LONGITUD LONGITUDNAME SP_JUR SP_JURNAME
## 1 1 1 1 0
## HARM_EV HARM_EVNAME MAN_COLL MAN_COLLNAME RELJCT1
## 1 0 1 0 1
## RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME
## 0 1 0 1 0
## WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 1 0 1 0 1
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME
## 0 48851 48851 48851 48851
## WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 1 0 1 0 0
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME
## 0 1 0 1 0
## ARR_HOUR ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR
## 1 0 1 0 1
## HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1 CF1NAME
## 0 1 0 48851 48851
## CF2 CF2NAME CF3 CF3NAME FATALS
## 48851 48851 48851 48851 1
## DRUNK_DR drugs_involved
## 13085 0# Paso 13: Variables vacías ("")
colSums(accident_merged =="")## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS
## 0 0 0 0 0
## PVH_INVL PEDS PERSONS PERMVIT PERNOTMVIT
## 0 0 0 0 0
## COUNTY COUNTYNAME CITY CITYNAME DAY
## 0 0 0 0 NA
## DAYNAME MONTH MONTHNAME YEAR DAY_WEEK
## NA NA 1 NA NA
## DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 1 NA 1 NA 1
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## NA 1 NA 1 1
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME
## 61541 NA 1 NA 1
## RD_OWNER RD_OWNERNAME MILEPT MILEPTNAME LATITUDE
## NA 1 NA 1 NA
## LATITUDENAME LONGITUD LONGITUDNAME SP_JUR SP_JURNAME
## NA NA NA NA 1
## HARM_EV HARM_EVNAME MAN_COLL MAN_COLLNAME RELJCT1
## NA 1 NA 1 NA
## RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME
## 1 NA 1 NA 1
## WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## NA 1 NA 1 NA
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME
## 1 NA NA NA NA
## WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## NA 1 NA 1 1
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME
## 1 NA 1 NA 1
## ARR_HOUR ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR
## NA 1 NA 1 NA
## HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1 CF1NAME
## 1 NA 1 NA NA
## CF2 CF2NAME CF3 CF3NAME FATALS
## NA NA NA NA NA
## DRUNK_DR drugs_involved
## NA 0# Paso 14: Se unifican los valores vacíos "Así solo tendremos Na y no Na y "".
accident_merged[accident_merged == ""] <- NA
# Se vuelve a verificar
colSums(is.na(accident_merged))## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS
## 0 0 0 0 0
## PVH_INVL PEDS PERSONS PERMVIT PERNOTMVIT
## 0 0 0 0 0
## COUNTY COUNTYNAME CITY CITYNAME DAY
## 0 0 0 0 1
## DAYNAME MONTH MONTHNAME YEAR DAY_WEEK
## 1 1 1 1 1
## DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 1 1 1 1 1
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 1 1 1 1 1
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME
## 61541 1 1 1 1
## RD_OWNER RD_OWNERNAME MILEPT MILEPTNAME LATITUDE
## 1 1 1 1 1
## LATITUDENAME LONGITUD LONGITUDNAME SP_JUR SP_JURNAME
## 1 1 1 1 1
## HARM_EV HARM_EVNAME MAN_COLL MAN_COLLNAME RELJCT1
## 1 1 1 1 1
## RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME
## 1 1 1 1 1
## WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 1 1 1 1 1
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME
## 1 48851 48851 48851 48851
## WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 1 1 1 1 1
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME
## 1 1 1 1 1
## ARR_HOUR ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR
## 1 1 1 1 1
## HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1 CF1NAME
## 1 1 1 48851 48851
## CF2 CF2NAME CF3 CF3NAME FATALS
## 48851 48851 48851 48851 1
## DRUNK_DR drugs_involved
## 13085 0Parece ser que en muchas variables falta un único dato. Intuitivamente podemos suponer que debe haber un registro o fila a la que le faltan muchos datos. Vamos a localizarla y eliminarla.
accident_merged %>%
# Paso 15: Localización de la fila incompleta
# Para cada fila calculamos el número de columnas que son NA
mutate(count_na = rowSums(is.na(.))) %>%
# Ordenamos por número de NA en orden descendente
arrange(desc(count_na)) %>%
# Vemos las 5 filas con más NA
head(5)## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS PVH_INVL PEDS PERSONS PERMVIT
## 1 15 Hawaii 150004 1 1 0 0 1 1
## 2 1 Alabama 10001 2 2 0 0 3 3
## 3 1 Alabama 10002 1 1 0 0 2 2
## 4 1 Alabama 10003 1 1 0 1 1 1
## 5 1 Alabama 10004 1 1 0 0 1 1
## PERNOTMVIT COUNTY COUNTYNAME CITY CITYNAME DAY DAYNAME MONTH
## 1 0 1 HAWAII (1) 4150 KEA NA NA NA
## 2 0 115 ST. CLAIR (115) 0 NOT APPLICABLE 12 12 2
## 3 0 73 JEFFERSON (73) 0 NOT APPLICABLE 11 11 2
## 4 1 73 JEFFERSON (73) 0 NOT APPLICABLE 7 7 2
## 5 0 117 SHELBY (117) 0 NOT APPLICABLE 3 3 2
## MONTHNAME YEAR DAY_WEEK DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 1 <NA> NA NA <NA> NA <NA> NA <NA>
## 2 February 2021 6 Friday 22 10:00pm-10:59pm 10 10
## 3 February 2021 5 Thursday 18 6:00pm-6:59pm 0 0
## 4 February 2021 1 Sunday 0 0:00am-0:59am 20 20
## 5 February 2021 4 Wednesday 16 4:00pm-4:59pm 20 20
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID TWAY_ID2 RUR_URB
## 1 NA <NA> NA <NA> <NA> <NA> NA
## 2 1 This section IS ON the NHS 1 Interstate I-20 <NA> 2
## 3 1 This section IS ON the NHS 1 Interstate I-459 <NA> 2
## 4 1 This section IS ON the NHS 1 Interstate I-459 <NA> 2
## 5 1 This section IS ON the NHS 1 Interstate I-65 <NA> 2
## RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER RD_OWNERNAME MILEPT
## 1 <NA> NA <NA> NA <NA> NA
## 2 Urban 1 Interstate 1 State Highway Agency 1569
## 3 Urban 1 Interstate 1 State Highway Agency 293
## 4 Urban 1 Interstate 1 State Highway Agency 183
## 5 Urban 1 Interstate 1 State Highway Agency 2300
## MILEPTNAME LATITUDE LATITUDENAME LONGITUD LONGITUDNAME SP_JUR
## 1 <NA> NA <NA> NA <NA> NA
## 2 1569 33.60164 33.60164167 -86.31238 -86.31238333 0
## 3 293 33.54136 33.54136111 -86.64374 -86.64374444 0
## 4 183 33.41980 33.41979722 -86.75257 -86.75257222 0
## 5 2300 33.36089 33.36089444 -86.77714 -86.77713889 0
## SP_JURNAME HARM_EV HARM_EVNAME MAN_COLL
## 1 <NA> NA <NA> NA
## 2 No Special Jurisdiction 12 Motor Vehicle In-Transport 2
## 3 No Special Jurisdiction 25 Concrete Traffic Barrier 0
## 4 No Special Jurisdiction 8 Pedestrian 0
## 5 No Special Jurisdiction 34 Ditch 0
## MAN_COLLNAME
## 1 <NA>
## 2 Front-to-Front
## 3 The First Harmful Event was Not a Collision with a Motor Vehicle in Transport
## 4 The First Harmful Event was Not a Collision with a Motor Vehicle in Transport
## 5 The First Harmful Event was Not a Collision with a Motor Vehicle in Transport
## RELJCT1 RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME WRK_ZONE
## 1 NA <NA> NA <NA> NA <NA> NA
## 2 0 No 1 Non-Junction 1 Not an Intersection 0
## 3 0 No 1 Non-Junction 1 Not an Intersection 0
## 4 0 No 1 Non-Junction 1 Not an Intersection 0
## 5 0 No 1 Non-Junction 1 Not an Intersection 0
## WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND LGT_CONDNAME WEATHER1
## 1 <NA> NA <NA> NA <NA> NA
## 2 None 1 On Roadway 2 Dark - Not Lighted NA
## 3 None 3 On Median 2 Dark - Not Lighted NA
## 4 None 2 On Shoulder 2 Dark - Not Lighted NA
## 5 None 4 On Roadside 1 Daylight NA
## WEATHER1NAME WEATHER2 WEATHER2NAME WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME
## 1 <NA> NA <NA> NA <NA> NA <NA>
## 2 <NA> NA <NA> 2 Rain 0 No
## 3 <NA> NA <NA> 2 Rain 0 No
## 4 <NA> NA <NA> 2 Rain 0 No
## 5 <NA> NA <NA> 1 Clear 0 No
## RAIL RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME ARR_HOUR
## 1 <NA> <NA> NA <NA> NA <NA> NA
## 2 0000000 Not Applicable 22 10:00pm-10:59pm 13 13 22
## 3 0000000 Not Applicable 99 Unknown 99 Unknown 19
## 4 0000000 Not Applicable 9 9:00am-9:59am 29 29 9
## 5 0000000 Not Applicable 16 4:00pm-4:59pm 20 20 16
## ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR HOSP_HRNAME
## 1 <NA> NA <NA> NA <NA>
## 2 10:00pm-10:59pm 25 25 23 11:00pm-11:59pm
## 3 7:00pm-7:59pm 9 9 88 Not Applicable (Not Transported)
## 4 9:00am-9:59am 40 40 88 Not Applicable (Not Transported)
## 5 4:00pm-4:59pm 28 28 99 Unknown
## HOSP_MN HOSP_MNNAME CF1 CF1NAME CF2 CF2NAME CF3 CF3NAME
## 1 NA <NA> NA <NA> NA <NA> NA <NA>
## 2 2 2 NA <NA> NA <NA> NA <NA>
## 3 88 Not Applicable (Not Transported) NA <NA> NA <NA> NA <NA>
## 4 88 Not Applicable (Not Transported) NA <NA> NA <NA> NA <NA>
## 5 99 Unknown EMS Hospital Arrival Time NA <NA> NA <NA> NA <NA>
## FATALS DRUNK_DR drugs_involved count_na
## 1 NA NA 1 77
## 2 2 NA 1 12
## 3 2 NA 1 12
## 4 1 NA 1 12
## 5 1 NA 1 12# Paso 16: Eliminación del registro ST_CASE 150004
accident_merged <- accident_merged %>%
filter(!(ST_CASE == 150004))
# Verificamos otra vez
colSums(is.na(accident_merged))## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS
## 0 0 0 0 0
## PVH_INVL PEDS PERSONS PERMVIT PERNOTMVIT
## 0 0 0 0 0
## COUNTY COUNTYNAME CITY CITYNAME DAY
## 0 0 0 0 0
## DAYNAME MONTH MONTHNAME YEAR DAY_WEEK
## 0 0 0 0 0
## DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 0 0 0 0 0
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 0 0 0 0 0
## TWAY_ID2 RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME
## 61539 0 0 0 0
## RD_OWNER RD_OWNERNAME MILEPT MILEPTNAME LATITUDE
## 0 0 0 0 0
## LATITUDENAME LONGITUD LONGITUDNAME SP_JUR SP_JURNAME
## 0 0 0 0 0
## HARM_EV HARM_EVNAME MAN_COLL MAN_COLLNAME RELJCT1
## 0 0 0 0 0
## RELJCT1NAME RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME
## 0 0 0 0 0
## WRK_ZONE WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND
## 0 0 0 0 0
## LGT_CONDNAME WEATHER1 WEATHER1NAME WEATHER2 WEATHER2NAME
## 0 48849 48849 48849 48849
## WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 0 0 0 0 0
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME
## 0 0 0 0 0
## ARR_HOUR ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR
## 0 0 0 0 0
## HOSP_HRNAME HOSP_MN HOSP_MNNAME CF1 CF1NAME
## 0 0 0 48849 48849
## CF2 CF2NAME CF3 CF3NAME FATALS
## 48849 48849 48849 48849 0
## DRUNK_DR drugs_involved
## 13084 0OBSERVACIONES IMPORTANTES
Seguimos analizando valores Na, TWAY_ID2, es la variable que ya tenía valores faltantes en el dataset de ejemplo. Al tratarse de vías de tránsito TWAY_ID de IDs de vías de tránsito de 1982 y TWAI_ID IDs de vías de tránsito de 2004, es posible que muchas vías tengan su identificador original del 1982, pero no el de 2004. De todos modos, estas variables no se van a usar. Como todas las vías tienen el identificador de 1982, se prescindirá del de 2004, eliminándolo.
Respecto a las variables WEATHER1, WEATHER1NAME y WEATHER2, WEATHER2NAME2. De entre los 3 datasets, estas columnas solo aparecen en el dataset de 2019, además, la variable WEATHER está completa igualmente en todos los datasets. Estas variables se eliminarán, solo interesaría WEATHER y WEATHERNAME, que es común a todos los datasets.
Los 18084 valores faltantes de la variable DRUNK_DR coinciden con el número de registros del dataset accident_2021. Tras realizar una comprobación, la variable DRUNK_DR (conductor bebido) se dejó de usar a partir del 2021. Quizás por la poca significancia que tiene discriminar el uso de alcohol en la conducción del resto de drogas, que sí que parecen tener una relación mucho más directa en la cause de accidentes. Se entiende pues, que los casos de alcohol están incluidos como drogas en general. Como el ejemplo original ya analiza el alcohol, se eliminará DRUNK_DR.
Respecto a los valores NA faltantes de las variables CF1, CF1NAME, CF2, CF2NAME y CF3, CF3NAME. Estas variables se definen en la lista de variables del ejemplo de la PEC, pero si uno se fija, no están incluidas en el dataset original proporcionado de 2020. El número de valores faltantes, 48849 coincide perfectamente con la suma de las observaciones de las observaciones de accident_2020 + accident_2021. Esto refleja que estas variables fueron usadas por última vez en el dataset de accident_2019. Luego se eliminarán.
# Paso 17: Tras haber hecho un análisis pormenorizado de cada columna conflictiva, se procede a su eliminación.
accident_merged <- accident_merged %>%
select(where(~ all(!is.na(.))))
# Verificamos una vez más
colSums(is.na(accident_merged)) ## STATE STATENAME ST_CASE VE_TOTAL VE_FORMS
## 0 0 0 0 0
## PVH_INVL PEDS PERSONS PERMVIT PERNOTMVIT
## 0 0 0 0 0
## COUNTY COUNTYNAME CITY CITYNAME DAY
## 0 0 0 0 0
## DAYNAME MONTH MONTHNAME YEAR DAY_WEEK
## 0 0 0 0 0
## DAY_WEEKNAME HOUR HOURNAME MINUTE MINUTENAME
## 0 0 0 0 0
## NHS NHSNAME ROUTE ROUTENAME TWAY_ID
## 0 0 0 0 0
## RUR_URB RUR_URBNAME FUNC_SYS FUNC_SYSNAME RD_OWNER
## 0 0 0 0 0
## RD_OWNERNAME MILEPT MILEPTNAME LATITUDE LATITUDENAME
## 0 0 0 0 0
## LONGITUD LONGITUDNAME SP_JUR SP_JURNAME HARM_EV
## 0 0 0 0 0
## HARM_EVNAME MAN_COLL MAN_COLLNAME RELJCT1 RELJCT1NAME
## 0 0 0 0 0
## RELJCT2 RELJCT2NAME TYP_INT TYP_INTNAME WRK_ZONE
## 0 0 0 0 0
## WRK_ZONENAME REL_ROAD REL_ROADNAME LGT_COND LGT_CONDNAME
## 0 0 0 0 0
## WEATHER WEATHERNAME SCH_BUS SCH_BUSNAME RAIL
## 0 0 0 0 0
## RAILNAME NOT_HOUR NOT_HOURNAME NOT_MIN NOT_MINNAME
## 0 0 0 0 0
## ARR_HOUR ARR_HOURNAME ARR_MIN ARR_MINNAME HOSP_HR
## 0 0 0 0 0
## HOSP_HRNAME HOSP_MN HOSP_MNNAME FATALS drugs_involved
## 0 0 0 0 0Finalmente, tenemos el dataset sin ningún valor vacío, pues hemos eliminado la fila que estaba incompleta y hemos comprobado que podemos prescindir de las variables que tenían valores faltantes. Estamos perdiendo información, sí, sobretodo con los 13.084 registros de conductores bebidos, pero no va aplicarse en el resto del ejercicio.
# TRANSFORMACIONES Y REDUCCIÓN DE LA DIMENSIONALIDAD
# Paso 18: Discretizar la hora en 4 franjas (Madrugada, mañana, tarde y noche)
# se ha optado por un criterio distinto al del Parlamento de Cataluña, porque quedaba una franja horaria desde 22-24 sin cubrir.
# Se iscretizaa la hora en 4 franjas (Madrugada, Mañana, Tarde, Noche)
accident_merged$HOUR_BIN <- cut(
accident_merged$HOUR,
breaks = c(-1, 6, 12, 18, 24), # Se usa -1 para que 0 quede en primer intervalo
labels = c("Madrugada", "Mañana", "Tarde", "Noche")
)
# Paso 19: Crear variable binaria WEEKEND_FLAG (DAY_WEEK = Sábado(7) o Domingo(1) segun FARS)
accident_merged$WEEKEND_FLAG <- ifelse(accident_merged$DAY_WEEK %in% c(1,7), 1, 0)
# Paso 20: Análisis de componentes principales (PCA)
# Seleccionamos columnas numéricas relevantes VE_TOTAL, PERSONS, PEDS, FATALS
numeric_vars <- c("VE_TOTAL", "PERSONS", "PEDS", "FATALS", "drugs_involved")
accident_numeric <- accident_merged[ , numeric_vars]
accident_numeric <- na.omit(accident_numeric)
# Se escala (centra y escala) antes de PCA
accident_numeric_scaled <- scale(accident_numeric)
# Se realiza el PCA
pca_acc <- prcomp(accident_numeric_scaled, center = TRUE, scale. = TRUE)
summary(pca_acc)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.3718 1.0079 0.9652 0.9193 0.5705
## Proportion of Variance 0.3764 0.2032 0.1863 0.1690 0.0651
## Cumulative Proportion 0.3764 0.5795 0.7659 0.9349 1.0000# Visualización de la varianza explicada (se requiere factoextra)
if(!require("factoextra")) install.packages("factoextra"); library("factoextra")
fviz_eig(pca_acc, addlabels=TRUE, ylim=c(0,50))
# Se la contribución de cada variable a los componentes
var <- get_pca_var(pca_acc)
var$contrib[,1:2] # contribuciones a PC1 y PC2## Dim.1 Dim.2
## VE_TOTAL 3.682083e+01 0.003503558
## PERSONS 4.126314e+01 0.192736947
## PEDS 1.119952e+01 10.265511345
## FATALS 1.071650e+01 6.389814536
## drugs_involved 2.231785e-06 83.148433613En total, los dos primeros componentes PC1 y PC2 acumulan cercca del 58% de la variabilidad de los datos.
¿Cómo contribuye cada variable a los dos primeros componentes? En Dim1: PERSONS es la variable que más contribuye a PC1. ve_total también tiene un peso grande. drugs_involved(0%) no influye nada, a la fatalidad. “Porque partimos de una muestra sesgada donde todos los accidentes son fatales”.
En Dim 2: drugs_involved (83.15%) domina por completo este segundo eje. Este eje captura la presencia vs la ausencia de drogas.
if(!require(ggplot2)) install.packages("ggplot2"); library(ggplot2)
# VISUALIZACIONES
# Distribución de FATALS
ggplot(accident_merged, aes(x=FATALS)) +
geom_histogram(binwidth=1) +
labs(title="Distribución de Número de Muertes",
x="Número de Muertes", y="Frecuencia")
# Accidentes con presencia de drogas vs sin
ggplot(accident_merged, aes(x=factor(drugs_involved))) +
geom_bar() +
labs(title="Accidentes con y sin presencia de drogas",
x="drugs_involved", y="Cantidad")
Gráficamente es muy llamativo cómo en la inmensa mayoría de accidentes con víctimas mortales en carretera hay presencia de drogas.
library(dplyr)
# Accidentes con o sin drogas a lo largo de los días de la semana
accident_by_day <- accident_merged %>%
group_by(DAY_WEEKNAME, drugs_involved) %>%
summarise(count=n(), .groups="drop")
# Visualización
ggplot(accident_by_day, aes(x=DAY_WEEKNAME, y=count, fill=factor(drugs_involved))) +
geom_bar(stat="identity", position="dodge") +
labs(
title="Accidentes por día de la semana y presencia de drogas",
x="Día de la semana",
y="Cantidad de accidentes",
fill="Drogas involucradas"
) +
theme(axis.text.x = element_text(angle=45, hjust=1))
Los viernes sábados y domingos son los días que más accidentes hay, tanto con drogas como sin drogas. Sin embargo, el uso de drogas aumenta el contraste entre el fin de semaana y los días laborales.
# Evolución temporal a lo largo de los años
accident_by_year <- accident_merged %>%
group_by(YEAR) %>%
summarise(total_muertes = sum(FATALS, na.rm=TRUE), .groups="drop")
ggplot(accident_by_year, aes(x=YEAR, y=total_muertes)) +
geom_line() +
geom_point() +
labs(
title="Evolución de muertes por accidente a lo largo de los años",
x="Año",
y="Total de muertes"
)
En este gráfico puede observarse como se redujo drásticamente el número de accidentes debido al confinamiento durante pandemia.
drugs_all %>%
group_by(DRUGRESNAME) %>%
summarise(count=n()) %>%
ggplot(aes(x=reorder(DRUGRESNAME, -count), y=count)) +
geom_bar(stat="identity") +
labs(
title="Distribución de los resultados del test de drogas",
x="Resultado (DRUGRESNAME)",
y="Cantidad de personas"
) +
theme(axis.text.x = element_text(angle=45, hjust=1))
El número de sustancias es abrumador, el gráfico no nos deja discernir qué sustancias son las más utilizadas.
# Haremos un top 20 de las más detectadas
top_drugs <- drugs_all %>%
filter(!DRUGRESNAME %in% c("Test Not Given", "Tested, No Drugs Found/Negative")) %>%
group_by(DRUGRESNAME) %>%
summarise(count=n()) %>%
arrange(desc(count)) %>%
slice_head(n=20) # top 20
ggplot(top_drugs, aes(x=reorder(DRUGRESNAME, count), y=count)) +
geom_bar(stat="identity") +
coord_flip() + # Para que las etiquetas se lean mejor
labs(
title="Top 20 sustancias detectadas en accidentes",
x="Sustancia",
y="Conteo de personas"
)
Otras drogas encabeza la lista, quizá por una mala detección. Se encuentran varios derivados del Cannabis en la lista, que sumados serían la primera sustancia.
library(dplyr)
library(ggplot2)
# 1) Contar cuántos registros (personas) hay por cada estado, en 'drugs_all'
# Esto incluye "Test not given", "No Drugs Found", etc.
drugs_by_state <- drugs_all %>%
group_by(STATENAME) %>%
summarise(count = n(), .groups = "drop") %>%
arrange(desc(count))
# 2) Nos quedamos con los top 20
top_20_states <- drugs_by_state %>%
slice_head(n = 20)
# 3) Gráfica de barras en orden descendente
ggplot(top_20_states, aes(x = reorder(STATENAME, count), y = count)) +
geom_bar(stat = "identity", fill = "steelblue") +
coord_flip() + # giramos el eje para lectura más cómoda
labs(
title = "Top 20 estados con más registros de test de drogas",
x = "Estado",
y = "Cantidad de registros"
)
Estos son los estados donde más registros de presencia de drogas en los accidentes se reportan.
A partir del juego de datos utilizado en el ejemplo de la PEC, realiza un análisis exploratorio de datos con el paquete explore() de R y comenta las ventajas e inconvenientes que presenta respecto al análisis realizado en el ejercicio 2.
Puedes utilizar la documentación publicada del paquete explore() en https://github.com/rolkra/explore
if(!require('explore')) install.packages('explore'); library('explore')
# Utilizaremos el juego de datos "pero el que contiene los 3 años que hemos hecho, para hacerlo más interesante"
# Redacta aquí el código R para el análisis exploratorio del juego de datos
accident_merged |> explore()Se abre una ventana a parte, donde podemos, sin la necesidad de emplear código, hacer un análisis exploratorio de datos. No obstante, vamos a indagar un poco en la documentación para sacarle más partido.
library(explore)
accident_merged |> describe()## # A tibble: 82 × 8
## variable type na na_pct unique min mean max
## <chr> <chr> <int> <dbl> <int> <dbl> <dbl> <dbl>
## 1 STATE int 0 0 51 1 24.1 56
## 2 STATENAME chr 0 0 51 NA NA NA
## 3 ST_CASE int 0 0 38134 10001 242168. 560121
## 4 VE_TOTAL int 0 0 18 1 1.58 59
## 5 VE_FORMS int 0 0 18 1 1.54 59
## 6 PVH_INVL int 0 0 12 0 0.04 11
## 7 PEDS int 0 0 11 0 0.24 10
## 8 PERSONS int 0 0 35 0 2.21 74
## 9 PERMVIT int 0 0 35 0 2.2 74
## 10 PERNOTMVIT int 0 0 11 0 0.25 10
## # ℹ 72 more rows# Exploración de la variable fatals
accident_2020 |>
explore(FATALS)
Muy similar al gráfico proporcionado en el ejercicio 1.
# Explorar "FATALS" y "VE_TOTAL" (dos numéricas) con target "drugs_involved"
accident_merged |>
explore(FATALS, target = drugs_involved)
Es curioso observar que tanto en los casos más habituales como en los poco frecuentes, el número de víctimas es mayor cuando hay drogas involucradas.
accident_merged |> explore( VE_TOTAL, target = drugs_involved)
Lo mismo ocurre con el número de vehículos involucrados. Hay más cuando hay drogas involucradas.
accident_merged %>% explore(DAY_WEEKNAME, target = drugs_involved)
Sábado es el día donde se observan diferencias significativas.
# report of all variables
#accident_merged %>% report(output_file = "report.html", output_dir = tempdir())Se genera un documento html, con un resumen de las variables, que se adjuntará junto con la PEC.
Entre las ventajas del uso de explore():
Gran rapidez y sencillez de uso Con pocas líneas de R puede visualizarse la estructura global del dataset, ver correlaciones, histogramas, estadísticos descriptivos y distribuciones de todas las variables.
Automatización No es necesario pedir que genere los gráficos como barras e histogramas, ni estadísticos como mín, max, mediana, etc. Este lo hace de forma automática, sin necesidad de programar.
Esta herramienta es ideal para una primera inspección o para un usuario principiante. Es muy didáctico a nivel visual.
Cuando se busca un overview, para presentárselo rápidamente a gente no técnica, permite ahorrar mucho tiempo.
Posibles inconvenientes y limitaciones de explore():
El nivel de personalización y ajuste es menor No permite tanta flexibilidad como ggplot2. explore() provee una solución más genérica, ya que si se dea cambiar colores de los ejes, superponer capas o filtrar los datos de forma dinámica, deberá hacerse antes o después.
Falta de profundidad en el análisis explore() facilita un primer análisis, no obstante, en su nivel actual, no parece quepueda sustituir análisis estadísticos más avanzados, como el test de hipótesis, reducción de dimensionalidad, y algunas comparaciones específicas. Lo más probable es que explore() sea insuficiente para un proyecto de minería de datos profesional.
Rendimiento en grandes datasets Si se utilizaran conjuntos de datos más grandes, explore() podría tener problemas de rendimiento debido a la carga de todos los gráficos. No sería de extrañar que se quedara corto de memoria. Con el análisis manual podemos muestrear para evitar este tipo de problemas.
Respecto con el Ejercicio 2
El tiempo dedicado al Ejercicio 2 ha sido ampliamente superior. En él realicé un proceso completo de limpieza y preprocesamiento de datos que requiere de mucha investigación de documentación y esfuerzo intelectual para un estudiante. Todo ello me ha otorgado un control más fino y la posibilidad de aprender de cada paso y poder justificarlo. Con explore con un primer momento habría ido más rápido, pero no habría llegado al nivel de profundidad que he llegado sin utilizarlo.
Mi conclusión es que, definitivamente, se sigue precisando del enfoque manual.