NYC buses: company level predictors with R
[This article was first published on R Programming – DataScience+, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Categories
Tags
We are following on from our previously loaded data of NYC bus delays, and from our simple Cubist regression fit. We know that we need to add in more predictors, and we have access to company level data. A warning to the faint-hearted: this post is essentially all data manipulation and joins. No pretty graphs, no conclusions. I am documenting the slog because the slog is required. We will move on to better things afterwards.
Preparation of company and route level data
We read in all of our data and take a look:
#---II_Breakdown_and_Delays ##NR: this actually isn't what I want here. Where is the other data? in_csv <- "../data/II_Bus_Breakdown_and_Delays.csv" ii_breakdowns <- read_csv(in_csv) ii_breakdowns %>% glimpse ## Observations: 259,637 ## Variables: 21 ## $ School_Year <chr> "2015-2016", "2015-2016", "201... ## $ Busbreakdown_ID <int> 1212699, 1212700, 1212701, 121... ## $ Run_Type <chr> "Special Ed AM Run", "Special ... ## $ Bus_No <chr> "48186", "2518", "235", "2102"... ## $ Route_Number <chr> "N758", "L530", "K168", "K216"... ## $ Reason <chr> "Other", "Mechanical Problem",... ## $ Schools_Serviced <chr> "75485", "21854", "18366", "21... ## $ Occurred_On <dttm> 2015-09-02 06:27:00, 2015-09-... ## $ Created_On <dttm> 2015-09-02 06:29:00, 2015-09-... ## $ Boro <chr> "Nassau County", "Brooklyn", "... ## $ Bus_Company_Name <chr> "BORO TRANSIT, INC.", "RELIANT... ## $ How_Long_Delayed <chr> "25 minutes", NA, "30MINS", "2... ## $ Number_Of_Students_On_The_Bus <int> 0, 0, 0, 1, 0, 0, 0, 9, 0, 2, ... ## $ Has_Contractor_Notified_Schools <chr> "Yes", "Yes", "Yes", "Yes", "Y... ## $ Has_Contractor_Notified_Parents <chr> "No", "Yes", "Yes", "Yes", "Ye... ## $ Have_You_Alerted_OPT <chr> "No", "Yes", "No", "No", "No",... ## $ Informed_On <dttm> 2015-09-02 06:29:00, 2015-09-... ## $ Incident_Number <chr> NA, NA, NA, NA, NA, NA, NA, NA... ## $ Last_Updated_On <dttm> 2015-09-02 06:29:16, 2015-09-... ## $ Breakdown_or_Running_Late <chr> "Running Late", "Breakdown", "... ## $ School_Age_or_PreK <chr> "School-Age", "School-Age", "S... ##-------data loading #---III_Drivers_and_Attendants in_csv <- "../data/III_Drivers_and_Attendants.csv" iii_drivers <- read_csv(in_csv) iii_drivers %>% glimpse ## Observations: 832 ## Variables: 5 ## $ School_Year <chr> "2018-2019", "2018-2019", "2018-2019", "2018-... ## $ Vendor_Name <chr> "ACADEMY EXPRESS LLC", "E-Z BUS LLC (B2321)",... ## $ Type_of_Service <chr> "Coach", "School Age", "School Age", "School ... ## $ Active_Employees <int> 65, 5, 3, 2, 1, 3, 6, 226, 22, 18, 13, 2, 20,... ## $ Job_Type <chr> "Driver", "Attendant", "Driver", "Driver&Atte... #---IV_Routes in_csv <- "../data/IV_Routes.csv" iv_routes <- read_csv(in_csv) iv_routes %>% glimpse ## Observations: 42,114 ## Variables: 14 ## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016"... ## $ Route_Number <chr> "C911", "J499", "J500", "J501", "J502... ## $ Service_Type <chr> "D2D", "D2D", "D2D", "D2D", "D2D", "D... ## $ Vehicle_TypeDescription <chr> "Non-Wheelchair Accessible Alternativ... ## $ Route_Start_Date <chr> "03/21/2016", "09/09/2015", "09/01/20... ## $ Vendor_Code <chr> "VN", "RV", "HT", "HT", "HT", "HT", "... ## $ Vendor_Name <chr> "VAN TRANS LLC (B2192)", "RELIANT TRA... ## $ Vendor_Affiliation <chr> "VAN TRANS LLC (B2192)", "RELIANT TRA... ## $ `Garage _Street_Address` <chr> "670 Hillside Road", "297 NORMAN AVEN... ## $ Garage_City <chr> "Pelham Manor", "Brooklyn", "BROOKLYN... ## $ Garage_State <chr> "NY", "NY", "NY", "NY", "NY", "NY", "... ## $ Garage_Zip <int> 10803, 11222, 11224, 11224, 11224, 11... ## $ XCoordinates <dbl> 1034574, 1000787, 983439, 983439, 983... ## $ YCoordinates <dbl> 265547.6, 204481.0, 148884.0, 148884.... #---V_Routes_by_Transportation_Sites in_csv <- "../data/V_Routes_by_Transportation_Sites.csv" v_rbt <- read_csv(in_csv) v_rbt %>% glimpse ## Observations: 98,506 ## Variables: 3 ## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016", "2015-2016... ## $ Route_Number <chr> "J698", "J699", "J700", "J700", "J701", "J701", "... ## $ OPT_Code <int> 75004, 75580, 75003, 75140, 75003, 75140, 75003, ... #---VI_Transportation_Sites in_csv <- "../data/VI_Transportation_Sites.csv" vi_transport <- read_csv(in_csv) vi_transport %>% glimpse ## Observations: 18,220 ## Variables: 18 ## $ School_Year <chr> "2016-2017", "2016-2017", "2016-... ## $ OPT_Code <chr> "01001", "01569", "01600", "0161... ## $ Name <chr> "Comprehensive Kids Devel. Scho"... ## $ Affiliation <chr> "Other Religion", "Public", "Pub... ## $ Site_Type <chr> "School", "School", "School", "S... ## $ Street_Address <chr> "101 Norfolk Street", "525 EAST ... ## $ City <chr> "Manhattan", "Manhattan", "Manha... ## $ State <chr> "NY", "NY", "NY", "NY", "NY", "N... ## $ Zip <int> 10002, 10002, 10002, 10002, 1000... ## $ Longitude <dbl> -73.98727, -73.97544, -73.99063,... ## $ Latitude <dbl> 40.71879, 40.71876, 40.72249, 40... ## $ Door_To_Door_Service <chr> "Yes", "Yes", "No", "Yes", "Yes"... ## $ Stop_To_School_Service <chr> "No", "No", "No", "No", "Yes", "... ## $ Common_Carrier_Svc_Metrocards <chr> "No", "Yes", "Yes", "Yes", "Yes"... ## $ Site_Reimbursement <chr> "No", "No", "No", "No", "No", "N... ## $ Mid_Day_Service <chr> "No", "No", "No", "Yes", "Yes", ... ## $ D2D_Late_Day_Programs <chr> "No", "No", "No", "No", "No", "N... ## $ S2S_Late_Day_Programs <chr> "No", "No", "No", "No", "No", "N... #---VII_Vehicles in_csv <- "../data/VII_Vehicles.csv" vii_vehicles <- read_csv(in_csv) vii_vehicles %>% glimpse ## Observations: 764 ## Variables: 4 ## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016",... ## $ Vendor_Name <chr> "ACADEMY EXPRESS LLC", "ACME BUS CORP.... ## $ Vehicle_TypeDescription <chr> "Motor Coach Bus", "Mini-Wagon", "Mini... ## $ Active_Vehicles <int> 22, 71, 6, 62, 1, 9, 71, 6, 63, 47, 67... #---VIII_PreK_Riders_by_Transportation_Site in_csv <- "../data/VIII_PreK_Riders_by_Transportation_Site.csv" viii_prt <- read_csv(in_csv) viii_prt %>% glimpse ## Observations: 1,268 ## Variables: 5 ## $ School_Year <chr> "2017-2018", "2017-2018", "2017-2018", "2017-... ## $ OPT_Code <chr> "C225", "E011", "E019", "E082", "E148", "E187... ## $ Site_Name <chr> "Highbridge Advisory Council", "Crossroads Sc... ## $ School_Name <chr> "Marshall England", "Crossroads School for Ch... ## $ Number_of_Riders <int> 50, 2, 1, 3, 3, 2, 90, 67, 22, 6, 7, 8, 4, 3,... #---IX_PreK_Vendors_by_Transportation_Site in_csv <- "../data/IX_PreK_Vendors_by_Transportation_Site.csv" ix_pvt <- read_csv(in_csv) ix_pvt %>% glimpse ## Observations: 1,833 ## Variables: 5 ## $ School_Year <chr> "2017-2018", "2017-2018", "2017-2018", "2017-2018"... ## $ OPT_Code <chr> "C053", "C104", "C149", "C153", "C157", "C158", "C... ## $ Site_Name <chr> "Bank Street Family Center", "PAL WORLD OF CREATIV... ## $ School_Name <chr> "Bank Street Family Center", "KIDS CENTRIC, INC", ... ## $ Vendor_Name <chr> "PHILLIPS BUS SERVICE", "L & M BUS CORP.", "L & M ...
And we begin to engage in some 'exploratory' data cleaning. That is, we aren't quite yet sure what we are going to use, and what we aren't going to use. But we know dirt when we see it.
First we try to clean up iv_routes a little:
#iv_routes
check_similar <- function(df){
df %>%
mutate(
diff= stringdist(Vendor_Affiliation, str_trunc(Vendor_Name,str_length(Vendor_Affiliation))),
flag_company_same = ifelse(diff <= str_length(Vendor_Affiliation)/6, 1L, 0L),
diff = NULL
)
}
iv_routes %>% check_similar() %>% filter(flag_company_same != 1) %>%
distinct(Vendor_Name, Vendor_Affiliation)
## # A tibble: 12 x 2
## Vendor_Name Vendor_Affiliation
## <chr> <chr>
## 1 EMPIRE STATE BUS CORP. ALLIED TRANSIT CORP.
## 2 LOGAN BUS COMPANY INC. LITTLE LINDA BUS CO.,INC.
## 3 LORINDA ENT. LTD. LITTLE LINDA BUS CO.,INC.
## 4 CONSOLIDATED BUS TRANS. INC. BORO TRANSIT, INC.
## 5 LORISSA BUS SERVICE INC. BOBBY`S BUS CO. INC.
## 6 GRANDPA`S BUS CO., INC. BOBBY`S BUS CO. INC.
## 7 LITTLE RICHIE BUS SERVICE LITTLE LINDA BUS CO.,INC.
## 8 THIRD AVENUE TRANSIT JOFAZ TRANSPORTATION INC.
## 9 LOGAN TRANSPORTATION SYSTEMS BOBBY`S BUS CO. INC.
## 10 CONSOLIDATED BUS TRANSIT, INC. BORO TRANSIT, INC.
## 11 LORINDA ENTERPRISES, LTD. LITTLE LINDA BUS CO.,INC.
## 12 THIRD AVENUE TRANSIT, INC JOFAZ TRANSPORTATION INC.
#cleaning that needs to be done prior to welding
iv_cleaned <- iv_routes %>%
mutate(Route_Start_Date = mdy(Route_Start_Date),
Garage_City = tolower(Garage_City)
) %>%
dplyr::rename(
"Garage_Street_Address" = 'Garage _Street_Address',
"Garage_XCoord" = "XCoordinates",
"Garage_YCoord" = "YCoordinates"
) %>%
check_similar()
iv_cleaned %>% glimpse
## Observations: 42,114
## Variables: 15
## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016",...
## $ Route_Number <chr> "C911", "J499", "J500", "J501", "J502"...
## $ Service_Type <chr> "D2D", "D2D", "D2D", "D2D", "D2D", "D2...
## $ Vehicle_TypeDescription <chr> "Non-Wheelchair Accessible Alternative...
## $ Route_Start_Date <date> 2016-03-21, 2015-09-09, 2015-09-01, 2...
## $ Vendor_Code <chr> "VN", "RV", "HT", "HT", "HT", "HT", "H...
## $ Vendor_Name <chr> "VAN TRANS LLC (B2192)", "RELIANT TRAN...
## $ Vendor_Affiliation <chr> "VAN TRANS LLC (B2192)", "RELIANT TRAN...
## $ Garage_Street_Address <chr> "670 Hillside Road", "297 NORMAN AVENU...
## $ Garage_City <chr> "pelham manor", "brooklyn", "brooklyn"...
## $ Garage_State <chr> "NY", "NY", "NY", "NY", "NY", "NY", "N...
## $ Garage_Zip <int> 10803, 11222, 11224, 11224, 11224, 112...
## $ Garage_XCoord <dbl> 1034574, 1000787, 983439, 983439, 9834...
## $ Garage_YCoord <dbl> 265547.6, 204481.0, 148884.0, 148884.0...
## $ flag_company_same <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
Next we join this to v_rbt, trying as best we can to figure out how they join, and to account for any observations which can't be joined as cleanly. 'School_Year' and 'Route_Number' appear to be our main joining headers, but this join leaves out 'OPT_Code' for some observations. These observations are not represented in the correct school year. We take a guess, and assign these an OPT_Code based only on the route, hoping that things haven't changed from year to year.
v_rbt %>% group_by(School_Year,Route_Number) %>% count %>% filter(n>1)
## # A tibble: 27,712 x 3
## # Groups: School_Year, Route_Number [27,712]
## School_Year Route_Number n
## <chr> <chr> <int>
## 1 2015-2016 J510 2
## 2 2015-2016 J690 2
## 3 2015-2016 J700 2
## 4 2015-2016 J701 2
## 5 2015-2016 J702 2
## 6 2015-2016 J703 2
## 7 2015-2016 J704 2
## 8 2015-2016 J705 2
## 9 2015-2016 J706 2
## 10 2015-2016 J707 2
## # ... with 27,702 more rows
# There are 27712 School_Year, Route_Number combinations in v_rbt which appear more than once. Why is that?
#Some year/route combinations have more than one OPT_Code. Why?
#OPT_Code represents an organization in a specific location. Each Route/School_Year combination thus can serve multiple organizations in multiple locations.
#This is going to cause us some problems later on, since we will have bus level data with multiple organizations attached. It will probably require some summary operations.
v_joined <- iv_cleaned %>%
left_join(v_rbt, by = c("School_Year", "Route_Number"))
#62 observations have no OPT_code after this join. Check if the Route_Numbers can be matched alone
iv_cleaned %>% anti_join(v_rbt, by = c("School_Year", "Route_Number")) %>%
distinct(Route_Number) #59 routes
## # A tibble: 59 x 1
## Route_Number
## <chr>
## 1 Z003
## 2 T003
## 3 T004
## 4 Q8227
## 5 T008
## 6 T009
## 7 Z001
## 8 Z002
## 9 Z004
## 10 Z005
## # ... with 49 more rows
#question: are all of these routes represented in the v_rbt, even if the year is wrong?
iv_cleaned %>% anti_join(v_rbt, by = c("School_Year", "Route_Number")) %>%
inner_join(v_rbt, by = "Route_Number") %>% distinct(Route_Number) #39 routes can be filled by ignoring year. Check.
## # A tibble: 39 x 1
## Route_Number
## <chr>
## 1 Z003
## 2 T004
## 3 T008
## 4 T009
## 5 Z001
## 6 Z002
## 7 Z004
## 8 Z005
## 9 X2353
## 10 X9353
## # ... with 29 more rows
iv_cleaned %>% anti_join(v_rbt, by = "Route_Number") %>% distinct(Route_Number) #Yes, 20 routes left for which we have no OPT information.
## # A tibble: 20 x 1
## Route_Number
## <chr>
## 1 T003
## 2 Q8227
## 3 X8304
## 4 X8338
## 5 X8339
## 6 X8406
## 7 X9874
## 8 X9875
## 9 X9881
## 10 M8068
## 11 M8246
## 12 M9068
## 13 M9246
## 14 X8621
## 15 X9621
## 16 D002
## 17 Q2455
## 18 Q9455
## 19 A002
## 20 I001
#OK, so I have 39 routes which need to join the OPT_Code on. But I only want to perform this join on the members of v_joined which are missing OPT data
v_sub <- v_joined %>%
filter(is.na(OPT_Code)) %>%
mutate(OPT_Code = NULL) %>%
left_join(v_rbt %>% select(Route_Number,OPT_Code), by = "Route_Number")
v_corrected <- v_joined %>%
filter(!is.na(OPT_Code)) %>%
rbind(v_sub) #and 20 routes remain without an OPT_Code
v_corrected %>% glimpse
## Observations: 98,632
## Variables: 16
## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016",...
## $ Route_Number <chr> "C911", "J499", "J500", "J501", "J502"...
## $ Service_Type <chr> "D2D", "D2D", "D2D", "D2D", "D2D", "D2...
## $ Vehicle_TypeDescription <chr> "Non-Wheelchair Accessible Alternative...
## $ Route_Start_Date <date> 2016-03-21, 2015-09-09, 2015-09-01, 2...
## $ Vendor_Code <chr> "VN", "RV", "HT", "HT", "HT", "HT", "H...
## $ Vendor_Name <chr> "VAN TRANS LLC (B2192)", "RELIANT TRAN...
## $ Vendor_Affiliation <chr> "VAN TRANS LLC (B2192)", "RELIANT TRAN...
## $ Garage_Street_Address <chr> "670 Hillside Road", "297 NORMAN AVENU...
## $ Garage_City <chr> "pelham manor", "brooklyn", "brooklyn"...
## $ Garage_State <chr> "NY", "NY", "NY", "NY", "NY", "NY", "N...
## $ Garage_Zip <int> 10803, 11222, 11224, 11224, 11224, 112...
## $ Garage_XCoord <dbl> 1034574, 1000787, 983439, 983439, 9834...
## $ Garage_YCoord <dbl> 265547.6, 204481.0, 148884.0, 148884.0...
## $ flag_company_same <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ OPT_Code <int> 75666, 75760, 75006, 75437, 75437, 759...
Next we work on vi_transport, including looking into how unique OPT_Codes really are. We transform the vehicle descriptions into counts of the various kinds of passengers and attributes of each vehicle. This approach is a kind of 'type agnostic' approach, so that a new vehicle type can be predicted based on its properties, rather than just being a new predictor which the model cannot account for.
vi_cleaned <- vi_transport %>%
mutate(City = tolower(City),
OPT_Code = as.integer(OPT_Code)
) %>%
dplyr::rename(
"location_city" = "City",
"location_state" = "State",
"location_zip" = "Zip",
"location_longitude" = "Longitude",
"location_latitude" = "Latitude"
) %>%
filter(!is.na(OPT_Code))
#OPT_Codes with non-duplicate lat/long
vi_code_and_first_location <- vi_cleaned %>%
group_by(OPT_Code, School_Year) %>%
distinct(location_longitude,location_latitude) %>%
count() %>%
filter(n>1) %>%
inner_join(vi_cleaned, by = c("OPT_Code","School_Year")) %>%
filter(row_number() == 1) %>%
select(OPT_Code, School_Year, location_longitude, location_latitude)
mismatches_location <- vi_cleaned %>%
inner_join(vi_code_and_first_location, by = c("OPT_Code", "School_Year")) %>%
select(OPT_Code, School_Year, location_longitude.x, location_longitude.y, location_latitude.x, location_latitude.y) %>%
mutate(
diff_longitude = location_longitude.x - location_longitude.y,
diff_latitude = location_latitude.x - location_latitude.y
) %>%
filter(diff_longitude > 0 | diff_latitude > 0)
#So there are 65 non unique OPT_Codes, it seems, with 85 that produce duplicates.
#When performed with OPT_Code and School_Year, it seems there are no matches.
#So I need make my match by OPT_Code / School_Year. OPT_Codes are changing from year to year.
vi_joined <- v_corrected %>%
left_join(vi_cleaned, by = c("OPT_Code", "School_Year"))
vi_joined %>%
distinct(Vehicle_TypeDescription)
## # A tibble: 10 x 1
## Vehicle_TypeDescription
## <chr>
## 1 Non-Wheelchair Accessible Alternative (NWC)
## 2 Mini-Wagon
## 3 Type A or B Flex Vehicle
## 4 Standard Bus (SE)
## 5 Standard Bus (GE)
## 6 Ambulance Transportation Service
## 7 Hydraulic Lift
## 8 Type C or D Flex
## 9 Ramp-Wagon
## 10 Wheelchair Accessible Alternative (WC)
#Vehicle_TypeDescription: Number_of_Riders, vehicle_reg_seats, vehicle_disabled_seats, vehicle_amblatory_seats, vehicle_lift, vehicle_attendant
vehicle_typedescription <- data.frame(Vehicle_TypeDescription =
c("Non-Wheelchair Accessible Alternative (NWC)", "Mini-Wagon", "Type A or B Flex Vehicle", "Standard Bus (SE)", "Standard Bus (GE)","Ambulance Transportation Service", "Hydraulic Lift", "Type C or D Flex","Ramp-Wagon","Wheelchair Accessible Alternative (WC)", "Motor Coach Bus"),
vehicle_max_riders = c(3,15,6,36,60,1,16,35,8,4, 55),
vehicle_reg_seats = c(3,15,6,36,60,0,0,35,0,0, 55),
vehicle_disabled_seats = c(0,0,6,0,0,0,8,10,4,2, 0),
vehicle_ambulatory_seats = c(0,0,0,0,0,1,8,0,4,2, 0),
vehicle_lift = as.integer(c(0,0,1,0,0,0,1,1,0,1, 0)),
vehicle_attendant = as.integer(c(1,1,1,0,0,0,1,1,1,1,0)),
vehicle_emt = c(0,0,0,0,0,2,0,0,0,0,0),
vehicle_bathroom = c(0,0,0,0,0,0,0,0,0,0,1)
) %>%
mutate(
Vehicle_TypeDescription = as.character(Vehicle_TypeDescription)
) %>%
arrange(vehicle_max_riders)
# c("Non-Wheelchair Accessible Alternative (NWC)", 3, 3, 0, 0, 0, 1),
# c("Mini-Wagon", 15, 15, 0, 0, 0, 1),
# c("Type A or B Flex Vehicle", 6, 6, 6, 0, 1, 1),
# c("Standard Bus (SE)", 36, 36, 0, 0, 0, 0),
# c("Standard Bus (GE)", 60, 60, 0, 0, 0, 0),
# c("Ambulance Transportation Service", 1, 0, 0, 1, 0, 1),
# c("Hydraulic Lift", 16, 0, 8, 8, 1, 1),
# c("Type C or D Flex", 35, 35, 10, 0, 1, 1),
# c("Ramp-Wagon", 8, 0, 4, 4, 0, 1),
# c("Wheelchair Accessible Alternative (WC)", 4, 0, 2, 2, 1, 1))
vi_vehicle_details <- vi_joined %>%
left_join(vehicle_typedescription, by = "Vehicle_TypeDescription")
#Worth noting that Number_of_Riders from ix_joined has no analogue in this data set. I have that data only for the PreK data, it seems
#Put the garage count by Vendor_Code, School_Year
garage_count_abs <- v_corrected %>%
select(School_Year, Vendor_Code, Garage_XCoord, Garage_YCoord) %>%
group_by(School_Year, Vendor_Code) %>%
distinct(Garage_XCoord, Garage_YCoord) %>%
count() %>%
mutate(
garage_max = n,
garage_min = n,
n = NULL
)
vi_garage_count <- vi_vehicle_details %>%
left_join(garage_count_abs, by = c("School_Year", "Vendor_Code")) %>%
mutate(
OPT_Code = as.character(OPT_Code)
)
vi_garage_count %>% glimpse
## Observations: 98,632
## Variables: 42
## $ School_Year <chr> "2015-2016", "2015-2016", "2015-...
## $ Route_Number <chr> "C911", "J499", "J500", "J501", ...
## $ Service_Type <chr> "D2D", "D2D", "D2D", "D2D", "D2D...
## $ Vehicle_TypeDescription <chr> "Non-Wheelchair Accessible Alter...
## $ Route_Start_Date <date> 2016-03-21, 2015-09-09, 2015-09...
## $ Vendor_Code <chr> "VN", "RV", "HT", "HT", "HT", "H...
## $ Vendor_Name <chr> "VAN TRANS LLC (B2192)", "RELIAN...
## $ Vendor_Affiliation <chr> "VAN TRANS LLC (B2192)", "RELIAN...
## $ Garage_Street_Address <chr> "670 Hillside Road", "297 NORMAN...
## $ Garage_City <chr> "pelham manor", "brooklyn", "bro...
## $ Garage_State <chr> "NY", "NY", "NY", "NY", "NY", "N...
## $ Garage_Zip <int> 10803, 11222, 11224, 11224, 1122...
## $ Garage_XCoord <dbl> 1034574, 1000787, 983439, 983439...
## $ Garage_YCoord <dbl> 265547.6, 204481.0, 148884.0, 14...
## $ flag_company_same <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ OPT_Code <chr> "75666", "75760", "75006", "7543...
## $ Name <chr> "CARMEL ACADEMY", "THE SINAI SCH...
## $ Affiliation <chr> "Yeshiva", "Yeshiva", "Nonsectar...
## $ Site_Type <chr> "School", "School", "School", "S...
## $ Street_Address <chr> "270 LAKE AVENUE", "110 SOUTH OR...
## $ location_city <chr> "connecticut", "new jersey", "ne...
## $ location_state <chr> "CT", "NJ", "NJ", "NJ", "NJ", "N...
## $ location_zip <int> 6830, 7039, 7083, 7666, 7666, 87...
## $ location_longitude <dbl> -73.63757, -74.36126, -74.27827,...
## $ location_latitude <dbl> 41.04900, 40.77635, 40.70391, 40...
## $ Door_To_Door_Service <chr> "Yes", "Yes", "Yes", "Yes", "Yes...
## $ Stop_To_School_Service <chr> "No", "No", "No", "No", "No", "N...
## $ Common_Carrier_Svc_Metrocards <chr> "Yes", "No", "No", "No", "No", "...
## $ Site_Reimbursement <chr> "No", "No", "No", "No", "No", "N...
## $ Mid_Day_Service <chr> "No", "No", "No", "No", "No", "N...
## $ D2D_Late_Day_Programs <chr> "No", "No", "No", "No", "No", "N...
## $ S2S_Late_Day_Programs <chr> "No", "No", "No", "No", "No", "N...
## $ vehicle_max_riders <dbl> 3, 15, 3, 3, 3, 3, 3, 15, 3, 3, ...
## $ vehicle_reg_seats <dbl> 3, 15, 3, 3, 3, 3, 3, 15, 3, 3, ...
## $ vehicle_disabled_seats <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ vehicle_ambulatory_seats <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ vehicle_lift <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ vehicle_attendant <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ vehicle_emt <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ vehicle_bathroom <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ garage_max <int> 3, 4, 2, 2, 2, 2, 2, 4, 2, 2, 2,...
## $ garage_min <int> 3, 4, 2, 2, 2, 2, 2, 4, 2, 2, 2,...
When I initially worked through this data, I tinkered with the pre-K data, ix_pvt, quite a lot. But in the end, it turned out that this data set only really existed as a way to join vehicle and staff metadata to the main data set via OPT_Code. The issue is that the school- and pre-K-level data don't share headers. So pulling, say, the garage information from the school-level data leaves NA values in the pre-K.
The driver and vehicle level data is, in the end, what I am after. The school and pre-K data is really just a way of attaching that metadata to the main data. Here is my cleaning work.
#---drivers and attendants, and vehicles, by vendor name---
#NR: next, I want to join on the driver numbers and vehicle numbers metadata.
#This must be joined by Vendor Name and School Year, and might require fuzzy joining.
iii_drivers %>% glimpse #832 obs
## Observations: 832
## Variables: 5
## $ School_Year <chr> "2018-2019", "2018-2019", "2018-2019", "2018-...
## $ Vendor_Name <chr> "ACADEMY EXPRESS LLC", "E-Z BUS LLC (B2321)",...
## $ Type_of_Service <chr> "Coach", "School Age", "School Age", "School ...
## $ Active_Employees <int> 65, 5, 3, 2, 1, 3, 6, 226, 22, 18, 13, 2, 20,...
## $ Job_Type <chr> "Driver", "Attendant", "Driver", "Driver&Atte...
vendors_unique <- vi_garage_count %>%
select(Vendor_Name) %>%
rbind(ix_pvt %>% select(Vendor_Name)) %>%
distinct(Vendor_Name) %>%
arrange(Vendor_Name) #87 unique vendors
vendors_unique %>%
anti_join(iii_drivers %>% filter(Type_of_Service != "Coach"), by = "Vendor_Name") %>%
glimpse
## Observations: 1
## Variables: 1
## $ Vendor_Name <chr> "ROBIN TRANSPORTATION"
#only "Robin Transportation" isn't covered by either iii or vii. Check check.
#All vendors are covered by both iii_ and vii_ *except* ROBIN. Goddamn it, Robin.
#iii_drivers
#needs turning into tidy data. What do I mean by this?
#We are after company level information.
#Per company and school year, we want the following counts
#number of coach services
#number of school services
#number of pre-k services
#number of ambulance services
#number of attendant services
#number of drivers (count doubles)
#number of attendants (count doubles)
#number of EMTs (count doubles)
#number of paramedics (count doubles)
#service counts
iii_service_count <- iii_drivers %>%
group_by(Vendor_Name, School_Year) %>%
summarise(
drivers_numServ_coach = sum(Type_of_Service == "Coach"),
drivers_numServ_school = sum(Type_of_Service == "School Age"),
drivers_numServ_prek = sum(Type_of_Service == "Pre-K"),
drivers_numServ_amb = sum(Type_of_Service == "Ambulance"),
drivers_numServ_att = sum(Type_of_Service == "Attendant")
)
#staff counts
iii_staff_count <- iii_drivers %>%
spread(Job_Type, Active_Employees, fill = 0) %>%
group_by(Vendor_Name, School_Year) %>%
summarise(
drivers_total_driver = sum(Driver + `Driver&Attendant` + `Driver&EMT`+ `Driver&Paramedic`),
drivers_total_attendant = sum(`Driver&Attendant` + Attendant),
drivers_total_EMT = sum(`Driver&EMT` + EMT),
drivers_total_paramedic = sum(`Driver&Paramedic` + Paramedic)
)
#counts left: staff per service, and each category of staff per service.
#I'm only going to count drivers per service
#drivers per service
iii_drivers_per_service <- iii_drivers %>%
spread(Job_Type, Active_Employees, fill = 0) %>%
mutate(
drivers_total_driver = Driver + `Driver&Attendant` + `Driver&EMT`+ `Driver&Paramedic`
) %>%
group_by(Vendor_Name, School_Year, Type_of_Service) %>%
summarise(
drivers_dummy = sum(drivers_total_driver)
) %>%
ungroup() %>%
spread(Type_of_Service, drivers_dummy, fill = 0) %>%
dplyr::rename(
"drivers_num_servAmb" = "Ambulance",
"drivers_num_servAttend" = "Attendant",
"drivers_num_servCoach" = "Coach",
"drivers_num_servPreK" = "Pre-K",
"drivers_num_servSchool" = "School Age"
)
#total staff per service
iii_total_staff_per_service <- iii_drivers %>%
group_by(Vendor_Name, School_Year, Type_of_Service) %>%
summarise(
drivers_total_dummy = sum(Active_Employees)
) %>%
spread(Type_of_Service, drivers_total_dummy, fill = 0) %>%
dplyr::rename(
"drivers_staff_servAmb" = "Ambulance",
"drivers_staff_servAttend" = "Attendant",
"drivers_staff_servCoach" = "Coach",
"drivers_staff_servPreK" = "Pre-K",
"drivers_staff_servSchool" = "School Age"
)
iii_joined <- iii_service_count %>%
left_join(iii_staff_count, by = c("Vendor_Name", "School_Year")) %>%
left_join(iii_drivers_per_service, by = c("Vendor_Name", "School_Year")) %>%
left_join(iii_total_staff_per_service, by = c("Vendor_Name", "School_Year"))
#vii_vehicles
vii_vehicles %>% glimpse #764 obs
## Observations: 764
## Variables: 4
## $ School_Year <chr> "2015-2016", "2015-2016", "2015-2016",...
## $ Vendor_Name <chr> "ACADEMY EXPRESS LLC", "ACME BUS CORP....
## $ Vehicle_TypeDescription <chr> "Motor Coach Bus", "Mini-Wagon", "Mini...
## $ Active_Vehicles <int> 22, 71, 6, 62, 1, 9, 71, 6, 63, 47, 67...
#total vehicles of each kind is actually not relevant.
#total seats of each kind
vii_vehicle_counts <- vii_vehicles %>%
left_join(vehicle_typedescription, by = "Vehicle_TypeDescription") %>%
mutate(
vehicle_total_max_riders = Active_Vehicles * vehicle_max_riders,
vehicle_total_reg_seats = Active_Vehicles * vehicle_reg_seats,
vehicle_total_disabled_seats = Active_Vehicles * vehicle_disabled_seats,
vehicle_total_ambulatory_seats = Active_Vehicles * vehicle_ambulatory_seats,
vehicle_total_with_lifts = Active_Vehicles * vehicle_lift,
vehicle_total_with_attendants = Active_Vehicles * vehicle_attendant,
vehicle_total_with_emts = Active_Vehicles * vehicle_emt,
vehicle_total_with_bathrooms = Active_Vehicles * vehicle_bathroom
) %>%
mutate(
vehicle_max_riders = NULL,
vehicle_reg_seats = NULL,
vehicle_disabled_seats = NULL,
vehicle_ambulatory_seats = NULL,
vehicle_lift = NULL,
vehicle_attendant = NULL,
vehicle_emt = NULL,
vehicle_bathroom = NULL
) %>%
#spread(Vehicle_TypeDescription, Active_Vehicles, fill = 0) %>%
group_by(Vendor_Name, School_Year) %>%
summarise_if(is.numeric,
sum
)
In my actual workflow, I wrote out all of my intermediate outputs as files that I could keep, and joined them in a separate workbook. In lieu of that process, the code below the fold simply renames some of the data sets that we have developed above. Apologies for any confusion: I promise that my project organization is a little less chaotic than this reporting markdown presentation.
vi_school <- vi_garage_count iii_staff <- iii_joined vii_vehicles <- vii_vehicle_counts
Joining the data
After a lot of messing around, we are at the stage where we can join the company level staff and vehicle data onto the main data set. We also use the times to make variables that tell us how close to either rush hour we are.
in_csv <- "../output/intermediate/ii_spread.csv"
ii_spread <- read_csv(in_csv)
#Vendor_Name, School_Year for iii_staff and vii_vehicles
ix_module <- ix_pvt %>%
mutate(Service_Type = "D2D") %>%
select(Vendor_Name, School_Year, OPT_Code, School_Name, Service_Type) %>%
inner_join(iii_staff, by = c("Vendor_Name", "School_Year")) %>%
inner_join(vii_vehicles, by = c("Vendor_Name", "School_Year"))
#only Robin transportation can't be joined from ix
#The MONTAUK data needs porting over to 2018-2019 from 2017-2018.
vii_vehicles <- vii_vehicles %>%
filter(str_detect(Vendor_Name,"MONTAUK STUDENT TRANS, INC.") & School_Year == "2017-2018") %>%
mutate(School_Year = "2018-2019") %>%
rbind(vii_vehicles)
iii_staff <- iii_staff %>%
filter(str_detect(Vendor_Name,"MONTAUK STUDENT TRANS, INC.") & School_Year == "2017-2018") %>%
mutate(School_Year = "2018-2019") %>%
rbind(iii_staff)
vi_module <- vi_school %>%
distinct(Route_Number, Vendor_Name, School_Year, Service_Type) %>%
inner_join(iii_staff, by = c("Vendor_Name", "School_Year")) %>%
inner_join(vii_vehicles, by = c("Vendor_Name", "School_Year"))
#Now I need to look at joining these onto the main data set.
#---main join for basic fit
ii_ix <- ii_spread %>%
inner_join(ix_module, by = c("Schools_Serviced" = "OPT_Code", "School_Year", "Bus_Company_Name" = "Vendor_Name")) #31 421
ii_vi <- ii_spread %>%
inner_join(vi_module, by = c("School_Year", "Route_Number")) #190 067
#I need to put time occured back in the data set.
#Anonymise company, school, route, and so on. This can be used in another, more detailed join some other time.
ii_joined <- ii_vi %>%
select(-Route_Number, -School_Year, -Schools_Serviced, -Bus_Company_Name, -Vendor_Name) %>%
rbind(
ii_ix %>%
select(-Route_Number, -School_Year, -Schools_Serviced, -Bus_Company_Name, -School_Name)
) %>%
mutate(
service_type_d2d = (Service_Type == "D2D")*1L,
Service_Type = NULL
)
tz(ii_joined$Occurred_On) <- "America/New_York" #it's set to UTC
#remove zero variance columns
zero_var <- function(dat) {
out <- lapply(dat, function(x) length(unique(x)))
keep <- rownames(data.frame(which(out > 1)))
dat %>%
select(keep)
}
ii_times <- ii_joined %>%
mutate(
time = strftime(Occurred_On, format="%H:%M:%S", tz = "America/New_York"),
time_am = lubridate::am(hms(time))*1L,
rush_within = ((hms(time) >= hms("07:00:00") & hms(time) <= hms("09:00:00")) | (hms(time) >= hms("16:00:00") & hms(time) <= hms("18:00:00")))*1L,
rush_between = (hms(time) >= hms("09:00:00") & hms(time) <= hms("16:00:00"))*1L,
rush_outside = (hms(time) <= hms("07:00:00") | hms(time) >= hms("18:00:00"))*1L,
#rush_time_from_peak = min(abs(as.numeric(hms(time) - hms("08:00:00"))), abs(as.numeric(hms(time) - hms("17:00:00"))))/60
rush_min_from_peak = pmin(abs(as.numeric(hms(time) - hms("08:00:00"))),abs(as.numeric(hms(time) - hms("17:00:00"))))/60,
time = NULL,
Occurred_On = NULL,
time_delayed = time_delayed - time_diff_report,
time_diff_report = NULL
) %>%
zero_var
ii_times %>% glimpse
## Observations: 221,375
## Variables: 48
## $ Busbreakdown_ID <int> 1212699, 1212701, 1212703, 121...
## $ Has_Contractor_Notified_Schools <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ Has_Contractor_Notified_Parents <int> 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, ...
## $ Have_You_Alerted_OPT <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
## $ Number_Of_Students_On_The_Bus <int> 0, 0, 1, 0, 9, 0, 2, 3, 3, 4, ...
## $ time_delayed <int> 23, 28, 13, 21, 80, 9, 9, 24, ...
## $ reported_before_resolved <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ School_Age <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ Reason_Accident <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Reason_DelayedbySchool <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Reason_FlatTire <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Reason_HeavyTraffic <int> 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, ...
## $ Reason_LatereturnfromFieldTrip <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Reason_MechanicalProblem <int> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, ...
## $ Reason_ProblemRun <int> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...
## $ Reason_WeatherConditions <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Reason_WontStart <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_Bronx <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_Brooklyn <int> 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, ...
## $ Boro_Connecticut <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_Manhattan <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
## $ Boro_NassauCounty <int> 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, ...
## $ Boro_NewJersey <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_Queens <int> 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, ...
## $ Boro_RocklandCounty <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_StatenIsland <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ Boro_Westchester <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ drivers_numServ_school <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, ...
## $ drivers_numServ_prek <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ drivers_total_driver <dbl> 567, 185, 22, 567, 105, 567, 2...
## $ drivers_total_attendant <dbl> 383, 196, 14, 383, 131, 383, 3...
## $ drivers_num_servPreK <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ drivers_num_servSchool <dbl> 567, 185, 22, 567, 105, 567, 2...
## $ drivers_staff_servPreK <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ drivers_staff_servSchool <dbl> 949, 378, 35, 949, 220, 949, 5...
## $ Active_Vehicles <int> 699, 204, 20, 699, 104, 699, 3...
## $ vehicle_total_max_riders <dbl> 23024, 6810, 681, 23024, 1577,...
## $ vehicle_total_reg_seats <dbl> 20472, 6810, 681, 20472, 1305,...
## $ vehicle_total_disabled_seats <dbl> 1276, 126, 0, 1276, 136, 1276,...
## $ vehicle_total_ambulatory_seats <dbl> 1276, 0, 0, 1276, 136, 1276, 2...
## $ vehicle_total_with_lifts <int> 81, 21, 0, 81, 17, 81, 13, 0, ...
## $ vehicle_total_with_attendants <int> 258, 21, 3, 258, 104, 258, 321...
## $ service_type_d2d <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ time_am <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ rush_within <int> 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, ...
## $ rush_between <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ rush_outside <int> 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, ...
## $ rush_min_from_peak <dbl> 93, 75, 65, 65, 41, 20, 15, 5,...
Related Post
- Visualizations for correlation matrices in R
- Interpretation of the AUC
- Simple Experiments with Smoothed Scatterplots
- Understanding the Covariance Matrix
- Six Sigma DMAIC Series in R – Part 3
To leave a comment for the author, please follow the link and comment on their blog: R Programming – DataScience+.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.