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.

    Categories

    1. Basic Statistics

    Tags

    1. Data Management
    2. R Programming

    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

    1. Visualizations for correlation matrices in R
    2. Interpretation of the AUC
    3. Simple Experiments with Smoothed Scatterplots
    4. Understanding the Covariance Matrix
    5. 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.

    Never miss an update!
    Subscribe to R-bloggers to receive
    e-mails with the latest R posts.
    (You will not see this message again.)

    Click here to close (This popup will not appear again)