diff --git a/script/Assign5.Rmd b/script/Assign5.Rmd
new file mode 100644
index 0000000..4583fb6
--- /dev/null
+++ b/script/Assign5.Rmd
@@ -0,0 +1,678 @@
+---
+title: "Assign5"
+author: "Viva Wan, Yaohan Xu, Neve Zhang"
+date: "2024-04-25"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+```{r}
+library(tidyverse)
+library(sf)
+library(lubridate)
+library(tigris)
+library(tidycensus)
+library(viridis)
+library(riem)
+library(gridExtra)
+library(knitr)
+library(kableExtra)
+library(RSocrata)
+library(jsonlite)
+library(FNN)
+library(ggplot2)
+source("https://raw.githubusercontent.com/urbanSpatial/Public-Policy-Analytics-Landing/master/functions.r")
+
+palette6 <- c("#eff3ff","#bdd7e7","#6baed6","#3182bd","#08519c","#08306b")
+palette5 <- c("#eff3ff","#bdd7e7","#6baed6","#3182bd","#08519c")
+palette4 <- c("#D2FBD4","#92BCAB","#527D82","#123F5A")
+palette2 <- c("#6baed6","#08519c")
+```
+
+```{r}
+plotTheme <- theme(
+ plot.title =element_text(size=12),
+ plot.subtitle = element_text(size=8),
+ plot.caption = element_text(size = 6),
+ axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
+ axis.text.y = element_text(size = 10),
+ axis.title.y = element_text(size = 10),
+ # Set the entire chart region to blank
+ panel.background=element_blank(),
+ plot.background=element_blank(),
+ # panel.border=element_rect(colour="#F0F0F0"),
+ # Format the grid
+ panel.grid.major=element_line(colour="#D0D0D0",size=.2),
+ axis.ticks=element_blank())
+```
+
+```{r}
+mapTheme <- theme(plot.title =element_text(size=12),
+ plot.subtitle = element_text(size=8),
+ plot.caption = element_text(size = 6),
+ axis.line=element_blank(),
+ axis.text.x=element_blank(),
+ axis.text.y=element_blank(),
+ axis.ticks=element_blank(),
+ axis.title.x=element_blank(),
+ axis.title.y=element_blank(),
+ panel.background=element_rect(fill = "black"),
+ panel.border=element_blank(),
+ panel.grid.major=element_line(colour = 'transparent'),
+ panel.grid.minor=element_blank(),
+ legend.direction = "vertical",
+ legend.position = "right",
+ plot.margin = margin(1, 1, 1, 1, 'cm'),
+ legend.key.height = unit(1, "cm"), legend.key.width = unit(0.2, "cm"))
+```
+
+## Data Gathering
+
+### Bike Trips Data
+
+We used July 2022 Indego Data, Philadelphia
+
+Used `as.POSIXct` to convert start and end time record into proper formatting `%m/%d/%Y %H:%M`
+
+```{r read data}
+dat <- read.csv("indego-trips-2022-q3/indego-trips-2022-q3.csv")
+
+dat <- dat %>%
+ mutate(start_time = as.POSIXct(start_time, format = "%m/%d/%Y %H:%M"),
+ end_time = as.POSIXct(end_time, format = "%m/%d/%Y %H:%M"))
+
+```
+
+Let’s use some date parsing to bin” the data by 15 and 60 minute intervals by rounding.
+
+Notice we use the time format ymd_hms to denote year, month, day and hour, minute and seccond. We extract the week of the observation (ranging from 1-52 throughout the year) and the dotw for day of the week.
+
+```{r parsing data}
+#parsing data
+dat2 <- dat %>%
+ mutate(interval60 = floor_date(ymd_hms(start_time), unit = "hour"),
+ interval15 = floor_date(ymd_hms(start_time), unit = "15 mins"),
+ month = month(interval60),
+ week = week(interval60),
+ dotw = wday(interval60, label=TRUE)) %>%
+ filter(month == 7)
+```
+
+### Weather Data
+
+Import weather data from Philadelphia Airport (PHL) using `riem_measures`. We can `mutate` the data to get temperature, wind speed, precipitation on an hourly basis and plot the temperature and precipitation trends over our study period.
+
+```{r reading weather data}
+weather.Panel <-
+ riem_measures(station = "PHL", date_start = "2022-07-01", date_end = "2022-07-31")
+
+weather.Panel <- weather.Panel %>%
+ dplyr::select(valid, tmpf, p01i, sknt)%>%
+ replace(is.na(.), 0) %>%
+ mutate(interval60 = ymd_h(substr(valid,1,13))) %>%
+ mutate(week = week(interval60),
+ dotw = wday(interval60, label=TRUE)) %>%
+ group_by(interval60) %>%
+ summarize(Temperature = max(tmpf),
+ Precipitation = sum(p01i),
+ Wind_Speed = max(sknt)) %>%
+ mutate(Temperature = ifelse(Temperature == 0, 42, Temperature))
+
+glimpse(weather.Panel)
+```
+Visualized weather data in Philadelphia of July 2022
+
+```{r plot_weather, fig.height=8, catche = TRUE}
+grid.arrange(
+ ggplot(weather.Panel, aes(interval60,Precipitation)) + geom_line() +
+ labs(title="Percipitation", x="Hour", y="Perecipitation") + plotTheme,
+ ggplot(weather.Panel, aes(interval60,Wind_Speed)) + geom_line() +
+ labs(title="Wind Speed", x="Hour", y="Wind Speed") + plotTheme,
+ ggplot(weather.Panel, aes(interval60,Temperature)) + geom_line() +
+ labs(title="Temperature", x="Hour", y="Temperature") + plotTheme,
+ top="Weather Data - Philadelphia PHL - July, 2022")
+```
+
+### Census Data
+
+```{r}
+PHLCensus <-
+ get_acs(geography = "tract",
+ variables = c("B01003_001", "B19013_001",
+ "B02001_002", "B08013_001",
+ "B08012_001", "B08301_001",
+ "B08301_010", "B01002_001"),
+ year = 2022,
+ state = "PA",
+ geometry = TRUE,
+ county= "Philadelphia",
+ output = "wide") %>%
+ rename(Total_Pop = B01003_001E,
+ Med_Inc = B19013_001E,
+ Med_Age = B01002_001E,
+ White_Pop = B02001_002E,
+ Travel_Time = B08013_001E,
+ Num_Commuters = B08012_001E,
+ Means_of_Transport = B08301_001E,
+ Total_Public_Trans = B08301_010E) %>%
+ select(Total_Pop, Med_Inc, White_Pop, Travel_Time,
+ Means_of_Transport, Total_Public_Trans,
+ Med_Age,
+ GEOID, geometry) %>%
+ mutate(Percent_White = White_Pop / Total_Pop,
+ Mean_Commute_Time = Travel_Time / Total_Public_Trans,
+ Percent_Taking_Public_Trans = Total_Public_Trans / Means_of_Transport)
+
+PHLTracts <- PHLCensus %>%
+ as.data.frame() %>%
+ distinct(GEOID, .keep_all = TRUE) %>%
+ select(GEOID, geometry) %>%
+ st_sf() %>%
+ st_transform(crs = "EPSG:2272")
+```
+We add the spatial information to our rideshare data as origin and destination data, first joining the origin station, then the destination station to our census data. We don’t use the destination data in this exercise.
+
+*Census data only helps us with later testing, but not model building.* - unit of analysis is bike station.
+
+```{r}
+dat_census <- dat2 %>%
+ filter(!is.na(start_lat) & !is.na(start_lon)
+ & !is.na(end_lat) & !is.na(end_lon)) %>%
+ st_as_sf(coords = c("start_lon", "start_lat"), crs = 4326) %>%
+ st_transform(crs = "EPSG:2272") %>%
+ st_join(PHLTracts %>% st_transform(crs = "EPSG:2272"), join = st_intersects, left = TRUE) %>%
+ rename(Origin.Tract = GEOID) %>%
+ mutate(start_lon = unlist(map(geometry, 1)),
+ start_lat = unlist(map(geometry, 2))) %>%
+ select(-geometry) %>%
+ st_as_sf(coords = c("end_lon", "end_lat"), crs = 4326) %>%
+ st_transform(crs = "EPSG:2272") %>%
+ st_join(PHLTracts %>% st_transform(crs = "EPSG:2272"), join = st_intersects, left = TRUE) %>%
+ rename(Destination.Tract = GEOID) %>%
+ mutate(end_lon = unlist(map(geometry, 1)),
+ end_lat = unlist(map(geometry, 2)))#%>%
+ #as.data.frame() %>%
+ #select(-geometry)
+```
+
+### Amenities Data
+
+Retrieved additional amenities data from OpenDataPhilly per assignment request.
+
+links:
+
+- https://opendataphilly.org/datasets/bike-network/
+
+- https://opendataphilly.org/datasets/schools/
+
+- SEPTA Highspeed Station Shapefile: https://gis-septa.hub.arcgis.com/datasets/af52d74b872045d0abb4a6bbbb249453_0/explore?location=40.009318%2C-75.215731%2C11.68
+
+```{r reading amenities data}
+BikeNetwork <- st_read("Amenities/Bike_Network/Bike_Network.shp") %>%
+ st_transform(crs = "EPSG:2272")
+
+Schools <- st_read("Amenities/Schools/Schools.shp")%>%
+ st_transform(crs = "EPSG:2272")
+
+HighSpeed <- st_read("Amenities/Highspeed_Stations/Highspeed_Stations.shp")%>%
+ st_transform(crs = "EPSG:2272")
+```
+transformed bike network into points. Measured nearest neighbor distance from *station* to each type of amenities.
+
+```{r calculate nearest neighbor}
+BikeNetwork.pt <- st_coordinates(BikeNetwork) %>%
+ as.data.frame() %>%
+ st_as_sf(., coords = c("X", "Y"), crs = "EPSG:2272")
+
+dat_census <-
+ dat_census %>%
+ mutate(
+ bikenetwork_nn1 = nn_function(st_coordinates(dat_census),
+ st_coordinates(BikeNetwork.pt), k = 1),
+ highspeed_nn1 = nn_function(st_coordinates(dat_census),
+ st_coordinates(HighSpeed), k = 1),
+ schools_nn1 = nn_function(st_coordinates(dat_census),
+ st_coordinates(Schools), k = 1))
+
+dat_census <- dat_census %>%
+ as.data.frame() %>%
+ select(-geometry)
+```
+
+## Data Exploration
+
+See title of each plot for plotting outcome.
+
+```{r}
+ggplot(dat_census %>%
+ group_by(interval60) %>%
+ tally())+
+ geom_line(aes(x = interval60, y = n))+
+ labs(title="Bike share trips per hr. Philadelphia, July, 2022",
+ x="Date",
+ y="Number of trips")+
+ plotTheme
+```
+
+```{r}
+dat_census %>%
+ mutate(time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 18 ~ "Overnight",
+ hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
+ hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
+ hour(interval60) >= 15 & hour(interval60) <= 18 ~ "PM Rush"))%>%
+ group_by(interval60, start_station, time_of_day) %>%
+ tally()%>%
+ group_by(start_station, time_of_day)%>%
+ summarize(mean_trips = mean(n))%>%
+ ggplot()+
+ geom_histogram(aes(mean_trips), binwidth = 0.5)+
+ labs(title="Mean Number of Hourly Trips Per Station. Philadelphia, July, 2022",
+ x="Number of trips",
+ y="Frequency")+
+ facet_wrap(~time_of_day)+
+ plotTheme
+
+ggplot(dat_census %>%
+ group_by(interval60, start_station) %>%
+ tally())+
+ geom_histogram(aes(n), binwidth = 2)+
+ labs(title="Bike share trips per hr by station. Philadelphia, July, 2022",
+ x="Number of Stations",
+ y="Trip Counts")+
+ plotTheme
+
+ggplot(dat_census %>% mutate(hour = hour(start_time)))+
+ geom_freqpoly(aes(hour, color = dotw), binwidth = 1)+
+ labs(title="Bike share trips in Philadelphia, by day of the week, July, 2022",
+ x="Hour",
+ y="Trip Counts")+
+ plotTheme
+
+
+ggplot(dat_census %>%
+ mutate(hour = hour(start_time),
+ weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday")))+
+ geom_freqpoly(aes(hour, color = weekend), binwidth = 1)+
+ labs(title="Bike share trips in Philadelphia - weekend vs weekday, July, 2022",
+ x="Hour",
+ y="Trip Counts")+
+ plotTheme
+```
+```{r origin_map }
+ggplot()+
+ geom_sf(data = PHLTracts, fill = "black")+
+ geom_point(data = dat_census %>%
+ mutate(hour = hour(start_time),
+ weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday"),
+ time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 18 ~ "Overnight",
+ hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
+ hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
+ hour(interval60) >= 15 & hour(interval60) <= 18 ~ "PM Rush"))%>%
+ group_by(start_station, start_lat, start_lon,
+ weekend, time_of_day) %>%
+ tally(),
+ aes(x=start_lon, y = start_lat, color = n),
+ fill = "transparent", alpha = 0.4, size = 0.3)+
+ scale_colour_viridis(direction = -1,
+ discrete = FALSE, option = "D")+
+ ylim(min(dat_census$start_lat), max(dat_census$start_lat))+
+ xlim(min(dat_census$start_lon), max(dat_census$start_lon))+
+ facet_grid(weekend ~ time_of_day)+
+ labs(title="Bike share trips per hr by station. Philadelphia, July, 2022")+
+ mapTheme
+```
+Added one more section to test for correlation between amenity locations and start trips made at each station (monthly total).
+
+Schools variable does not show any correlation, hence eliminated from later model.
+
+```{r}
+dat_nn <- dat_census %>%
+ group_by(start_station) %>%
+ summarize(month_tot = n(),
+ bikenetwork_nn1 = mean(bikenetwork_nn1),
+ highspeed_nn1 = mean(highspeed_nn1),
+ schools_nn1 = mean(schools_nn1))
+
+correlation_matrix <- cor(dat_nn)
+
+ggcorrplot(
+ correlation_matrix,
+ p.mat = cor_pmat(dat_nn),
+ colors = c("blue", "white", "darkred"),
+ type = "lower",
+ insig = "blank"
+) +
+labs(title = "Correlation Test Amenities Variables")+
+theme(axis.text.x = element_text(size = 10))
+```
+
+## Create Space-Time Panel
+
+First we have to make sure each unique station and hour/day combo exists in our data set. This is done in order to *create a “panel” (e.g. a time-series) data set where each time period in the study is represented by a row - whether an observation took place then or not*. So if a station didn’t have any trips originating from it at a given hour, we still need a zero in that spot in the panel.
+
+```{r space-time panel}
+study.panel <-
+ expand.grid(interval60=unique(dat_census$interval60),
+ start_station = unique(dat_census$start_station)) %>%
+ left_join(., dat_census %>%
+ select(start_station, Origin.Tract,
+ start_lon, start_lat, bikenetwork_nn1, highspeed_nn1, schools_nn1)%>%
+ distinct() %>%
+ group_by(start_station) %>%
+ slice(1))
+
+```
+We create the full panel by *summarizing counts by station for each time interval*, keep census info and lat/lon information along for joining later to other data. We remove data for station IDs that are FALSE.
+
+Note amenities features are also added.
+
+```{r full panel}
+ride.panel <-
+ dat_census %>%
+ mutate(Trip_Counter = 1) %>%
+ right_join(study.panel) %>%
+ group_by(interval60, start_station, Origin.Tract, start_lon, start_lat,
+ bikenetwork_nn1, highspeed_nn1, schools_nn1) %>% #amenities feature added here
+ summarize(Trip_Count = sum(Trip_Counter, na.rm=T)) %>%
+ left_join(weather.Panel) %>%
+ ungroup() %>%
+ filter(is.na(start_station) == FALSE) %>%
+ mutate(week = week(interval60),
+ dotw = wday(interval60, label = TRUE)) %>%
+ filter(is.na(Origin.Tract) == FALSE)
+
+ride.panel <-
+ left_join(ride.panel, PHLCensus %>%
+ as.data.frame() %>%
+ select(-geometry), by = c("Origin.Tract" = "GEOID"))
+```
+
+## Create Time Lags
+
+Creating time lag variables will add additional nuance about the demand during a given time period - hours before and during that day.
+
+We can also try to control for the effects of holidays that disrupt the expected demand during a given weekend or weekday. (We have a holiday on May 28 - Memorial Day. For that three day weekend we could use some dummy variables indicating temporal proximity to the holiday.)
+
+**NEVE NOTE**: Replaced May 28 Memorial Day with July 4th Independence Day.
+
+Keep in mind, that unique fixed effects must be in your training set when you run your models.
+
+```{r}
+ride.panel <-
+ ride.panel %>%
+ arrange(start_station, interval60) %>%
+ mutate(lagHour = dplyr::lag(Trip_Count,1),
+ lag2Hours = dplyr::lag(Trip_Count,2),
+ lag3Hours = dplyr::lag(Trip_Count,3),
+ lag4Hours = dplyr::lag(Trip_Count,4),
+ lag12Hours = dplyr::lag(Trip_Count,12),
+ lag1day = dplyr::lag(Trip_Count,24),
+ holiday = ifelse(yday(interval60) == 185,1,0)) %>%
+ mutate(day = yday(interval60)) %>%
+ mutate(holidayLag = case_when(dplyr::lag(holiday, 1) == 1 ~ "PlusOneDay",
+ dplyr::lag(holiday, 2) == 1 ~ "PlustTwoDays",
+ dplyr::lag(holiday, 3) == 1 ~ "PlustThreeDays",
+ dplyr::lead(holiday, 1) == 1 ~ "MinusOneDay",
+ dplyr::lead(holiday, 2) == 1 ~ "MinusTwoDays",
+ dplyr::lead(holiday, 3) == 1 ~ "MinusThreeDays",
+ TRUE ~ NA_character_))
+
+ride.panel$holidayLag[is.na(ride.panel$holidayLag)] <- "0"
+```
+
+## Developing Models
+
+Splitted data based on whether the trip occurred after or before week 29, to keep temporal cohesiveness.
+
+Since holiday effect have to be accounted for by the training data (ride.Train), all data *before* week 29 is splitted into the training set.
+
+```{r}
+ride.Train <- filter(ride.panel, week < 29)
+ride.Test <- filter(ride.panel, week >= 29)
+```
+
+Testing out different models. reg1-reg5 follow class demo. reg6 incorporates added amenities feature per assignment request.
+
+```{r}
+reg1 <-
+ lm(Trip_Count ~ hour(interval60) + dotw + Temperature, data=ride.Train)
+
+reg2 <-
+ lm(Trip_Count ~ start_station + dotw + Temperature, data=ride.Train)
+
+reg3 <-
+ lm(Trip_Count ~ start_station + hour(interval60) + dotw + Temperature + Precipitation,
+ data=ride.Train)
+
+reg4 <-
+ lm(Trip_Count ~ start_station + hour(interval60) + dotw + Temperature + Precipitation +
+ lagHour + lag2Hours +lag3Hours + lag12Hours + lag1day,
+ data=ride.Train)
+
+reg5 <-
+ lm(Trip_Count ~ start_station + hour(interval60) + dotw + Temperature + Precipitation +
+ lagHour + lag2Hours +lag3Hours +lag12Hours + lag1day + holidayLag + holiday,
+ data=ride.Train)
+
+reg6 <-
+ lm(Trip_Count ~ start_station + hour(interval60) + dotw + Temperature + Precipitation +
+ lagHour + lag2Hours +lag3Hours +lag12Hours + lag1day + holidayLag + holiday +
+ bikenetwork_nn1 + highspeed_nn1, #did not include schools since there's no correlation
+ data=ride.Train)
+```
+
+### Effect of amenities?
+
+Turns out amenities are good predictors.
+
+```{r final model}
+summary(reg6)
+vif(reg6)
+```
+
+## Predict for Test Data
+
+create a nested data frame of test data by week. Nesting means that instead of merely having a “flat” file consisting of rows and columns, you have a matrix of other objects - imagine each cell in a matrix containing another matrix within it, or a list, or a list of lists.
+
+```{r create nested data frame}
+ride.Test.weekNest <-
+ ride.Test %>%
+ nest(data = -week)
+```
+
+think of `map` as visiting each dataframe in a nested data set and applies a function to it. We create a function called model_pred which we can then `map` onto each data frame in our nested structure.
+
+```{r function - map into nested data frame}
+model_pred <- function(dat, fit){
+ pred <- predict(fit, newdata = dat)}
+```
+
+```{r predictions}
+week_predictions <-
+ ride.Test.weekNest %>%
+ mutate(ATime_FE = map(.x = data, fit = reg1, .f = model_pred),
+ BSpace_FE = map(.x = data, fit = reg2, .f = model_pred),
+ CTime_Space_FE = map(.x = data, fit = reg3, .f = model_pred),
+ DTime_Space_FE_timeLags = map(.x = data, fit = reg4, .f = model_pred),
+ ETime_Space_FE_timeLags_holidayLags = map(.x = data, fit = reg5, .f = model_pred),
+ FTime_Space_FE_timeLags_holidayLags_amenities = map(.x = data, fit = reg6, .f = model_pred)) %>%
+ gather(Regression, Prediction, -data, -week) %>%
+ mutate(Observed = map(data, pull, Trip_Count),
+ Absolute_Error = map2(Observed, Prediction, ~ abs(.x - .y)),
+ MAE = map_dbl(Absolute_Error, mean, na.rm = TRUE),
+ sd_AE = map_dbl(Absolute_Error, sd, na.rm = TRUE))
+
+week_predictions
+```
+
+## Examining Model Performance
+
+The best models - the lag models, are accurate to less than an average of one ride per hour, at a glance, that’s pretty alright for overall accuracy. - In our project, the best model is reg6, which accounts for lags as well as amenities.
+
+**Which models perform best - and how would you describe their fit?**
+
+**Why don’t the holiday time lags seem to matter?**
+
+```{r MAPE by mode specification and week}
+
+week_predictions %>%
+ dplyr::select(week, Regression, MAE) %>%
+ gather(Variable, MAE, -Regression, -week) %>%
+ ggplot(aes(week, MAE)) +
+ geom_bar(aes(fill = Regression), position = "dodge", stat="identity") +
+ scale_fill_manual(values = palette6) +
+ labs(title = "Mean Absolute Errors by model specification and week") +
+ plotTheme
+
+week_predictions %>%
+ mutate(interval60 = map(data, pull, interval60),
+ start_station = map(data, pull, start_station)) %>%
+ dplyr::select(interval60, start_station, Observed, Prediction, Regression) %>%
+ unnest() %>%
+ gather(Variable, Value, -Regression, -interval60, -start_station) %>%
+ group_by(Regression, Variable, interval60) %>%
+ summarize(Value = sum(Value)) %>%
+ ggplot(aes(interval60, Value, colour=Variable)) +
+ geom_line(size = 1.1) +
+ facet_wrap(~Regression, ncol=1) +
+ labs(title = "Predicted/Observed bike share time series",
+ subtitle = "Philadelphia; A test set of 2 weeks",
+ x = "Hour", y= "Station Trips") +
+ plotTheme
+
+```
+
+We can look at our mean absolute errors by station - there seems to be a spatial pattern to our error (what is it?), but we need to go a bit further to get at the temporal element of the error.
+
+```{r}
+week_predictions %>%
+ mutate(interval60 = map(data, pull, interval60),
+ start_station = map(data, pull, start_station),
+ start_lat = map(data, pull, start_lat),
+ start_lon = map(data, pull, start_lon)) %>%
+ select(interval60, start_station, start_lon, start_lat, Observed, Prediction, Regression) %>%
+ unnest() %>%
+ filter(Regression == "FTime_Space_FE_timeLags_holidayLags_amenities") %>%
+ group_by(start_station, start_lon, start_lat) %>%
+ summarize(MAE = mean(abs(Observed-Prediction), na.rm = TRUE))%>%
+ggplot(.)+
+ geom_sf(data = PHLTracts, color = "grey", fill = "white")+
+ geom_point(aes(x = start_lon, y = start_lat, color = MAE),
+ fill = "transparent", alpha = 0.4)+
+ scale_colour_viridis(direction = -1,
+ discrete = FALSE, option = "D")+
+ ylim(min(dat_census$start_lat), max(dat_census$start_lat))+
+ xlim(min(dat_census$start_lon), max(dat_census$start_lon))+
+ labs(title="Mean Abs Error, Test Set, Best Model") #model 7
+```
+
+
+## Space-Time Error Evaluation - Pending Annotation
+
+If we plot observed vs. predicted for different times of day during the week and weekend, some patterns begin to emerge. We are certainly underpredicting in general, but *what do we begin to see about some of the outcomes that our model cannot explain?*
+
+```{r}
+week_predictions %>%
+ mutate(interval60 = map(data, pull, interval60),
+ start_station = map(data, pull, start_station),
+ start_lat = map(data, pull, start_lat),
+ start_lon = map(data, pull, start_lon),
+ dotw = map(data, pull, dotw)) %>%
+ select(interval60, start_station, start_lon,
+ start_lat, Observed, Prediction, Regression,
+ dotw) %>%
+ unnest() %>%
+ filter(Regression == "FTime_Space_FE_timeLags_holidayLags_amenities")%>%
+ mutate(weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday"),
+ time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 18 ~ "Overnight",
+ hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
+ hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
+ hour(interval60) >= 15 & hour(interval60) <= 18 ~ "PM Rush"))%>%
+ ggplot()+
+ geom_point(aes(x= Observed, y = Prediction))+
+ geom_smooth(aes(x= Observed, y= Prediction), method = "lm", se = FALSE, color = "red")+
+ geom_abline(slope = 1, intercept = 0)+
+ facet_grid(time_of_day~weekend)+
+ labs(title="Observed vs Predicted",
+ x="Observed trips",
+ y="Predicted trips")+
+ plotTheme
+```
+
+Is there a spatial pattern to these big errors? Let’s look at our errors on a map by weekend/weekday and time of day. **What is the implication for rebalancing with these errors - does it matter if they are higher volume locations?**
+
+```{r}
+week_predictions %>%
+ mutate(interval60 = map(data, pull, interval60),
+ start_station = map(data, pull, start_station),
+ start_lat = map(data, pull, start_lat),
+ start_lon = map(data, pull, start_lon),
+ dotw = map(data, pull, dotw) ) %>%
+ select(interval60, start_station, start_lon,
+ start_lat, Observed, Prediction, Regression,
+ dotw) %>%
+ unnest() %>%
+ filter(Regression == "FTime_Space_FE_timeLags_holidayLags_amenities")%>%
+ mutate(weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday"),
+ time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 18 ~ "Overnight",
+ hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
+ hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
+ hour(interval60) >= 15 & hour(interval60) <= 18 ~ "PM Rush")) %>%
+ group_by(start_station, weekend, time_of_day, start_lon, start_lat) %>%
+ summarize(MAE = mean(abs(Observed-Prediction), na.rm = TRUE))%>%
+ ggplot(.)+
+ geom_sf(data = PHLTracts, fill = "black")+
+ geom_point(aes(x = start_lon, y = start_lat, color = MAE),
+ fill = "transparent", size = 0.5, alpha = 0.4)+
+ scale_colour_viridis(direction = -1,
+ discrete = FALSE, option = "D")+
+ ylim(min(dat_census$start_lat), max(dat_census$start_lat))+
+ xlim(min(dat_census$start_lon), max(dat_census$start_lon))+
+ facet_grid(weekend~time_of_day)+
+ labs(title="Mean Absolute Errors, Test Set")+
+ mapTheme
+```
+
+(How the model error change with reference to census tracts' socio-economic conditions?) - this is where previously added census info come into place.
+
+```{r}
+week_predictions %>%
+ mutate(interval60 = map(data, pull, interval60),
+ start_station = map(data, pull, start_station),
+ start_lat = map(data, pull, start_lat),
+ start_lon = map(data, pull, start_lon),
+ dotw = map(data, pull, dotw),
+ Percent_Taking_Public_Trans = map(data, pull, Percent_Taking_Public_Trans),
+ Med_Inc = map(data, pull, Med_Inc),
+ Percent_White = map(data, pull, Percent_White)) %>%
+ select(interval60, start_station, start_lon,
+ start_lat, Observed, Prediction, Regression,
+ dotw, Percent_Taking_Public_Trans, Med_Inc, Percent_White) %>%
+ unnest() %>%
+ filter(Regression == "FTime_Space_FE_timeLags_holidayLags_amenities")%>%
+ mutate(weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday"),
+ time_of_day = case_when(hour(interval60) < 7 | hour(interval60) > 18 ~ "Overnight",
+ hour(interval60) >= 7 & hour(interval60) < 10 ~ "AM Rush",
+ hour(interval60) >= 10 & hour(interval60) < 15 ~ "Mid-Day",
+ hour(interval60) >= 15 & hour(interval60) <= 18 ~ "PM Rush")) %>%
+ filter(time_of_day == "AM Rush") %>%
+ group_by(start_station, Percent_Taking_Public_Trans, Med_Inc, Percent_White) %>%
+ summarize(MAE = mean(abs(Observed-Prediction), na.rm = TRUE))%>%
+ gather(-start_station, -MAE, key = "variable", value = "value")%>%
+ ggplot(.)+
+ geom_point(aes(x = value, y = MAE), alpha = 0.4)+
+ geom_smooth(aes(x = value, y = MAE), method = "lm", se= FALSE)+
+ facet_wrap(~variable, scales = "free")+
+ labs(title="Errors as a function of socio-economic variables",
+ y="Mean Absolute Error (Trips)")+
+ plotTheme
+```
+
+## Prediction Interpretation
+
+Assignment requirement: Conclude with how useful your algorithm is for the bike re-balancing plan.
+
+---
\ No newline at end of file
diff --git a/script/Assignment4.Rmd b/script/Assignment4.Rmd
index 6a80be6..aa665b9 100644
--- a/script/Assignment4.Rmd
+++ b/script/Assignment4.Rmd
@@ -576,10 +576,10 @@ filter(testProbs.thresholds, Threshold == .5) %>%
- Average cost of a crime committed by a recidivist: $50,000
- Average compensation for wrongful incarceration: $50,000
-**True negative Cost**: "An individual is predicted not to recidivate and indeed does not recidivate." - No cost
-**True positive Cost**: "An individual is predicted to recidivate and then does recidivate." - $40k
-**False negative Cost**: "An individual is predicted not to recidivate but then does recidivate." - $100k
-**False positive Cost**: "An individual is predicted to recidivate but does not actually recidivate." - $50k
+- **True negative Cost**: "An individual is predicted not to recidivate and indeed does not recidivate." - No cost
+- **True positive Cost**: "An individual is predicted to recidivate and then does recidivate." - $40k
+- **False negative Cost**: "An individual is predicted not to recidivate but then does recidivate." - $100k
+- **False positive Cost**: "An individual is predicted to recidivate but does not actually recidivate." - $50k
```{r cost}
whichThreshold.cost <- whichThreshold %>%
@@ -608,7 +608,7 @@ whichThreshold.cost %>%
threshold.cost <- whichThreshold.cost%>%
filter(Variable == 'Cost.all')
-#-> lowest cost when threshold = 0.49
+#-> lowest cost when threshold = 0.50
```
diff --git a/script/Assignment4.html b/script/Assignment4.html
index 1ec8811..37aac4c 100644
--- a/script/Assignment4.html
+++ b/script/Assignment4.html
@@ -2271,7 +2271,7 @@ ROC Curve
theme(legend.title = element_text(size = 9),
legend.text = element_text(size = 8)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
-
+
# choose the threshold:
# baseline (threshold =0.5): sensitivity = 0.8120290, specificity = 0.5394175
# a better one: sensitivity increase > specificity decrease
@@ -2302,7 +2302,7 @@ ROC Curve
theme(legend.title = element_text(size = 9),
legend.text = element_text(size = 8)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
-
+
#plot the changes
whichThreshold.nocost %>%
select(Sensitivity.change, Specificity.change, change.sum, Threshold)%>%
@@ -2318,7 +2318,7 @@ ROC Curve
theme(legend.title = element_text(size = 9),
legend.text = element_text(size = 8)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
-
+
# table the change
whichThreshold.nocost %>%
select(Threshold, Sensitivity.change, Specificity.change, change.sum)%>%
@@ -3838,20 +3838,23 @@ Error and Equity Discussion
Cost-Benifit Analysis
-- Average cost of keeping someone in prison for a year: $40,000
-
-- Average cost of arrest and trial for a recidivist: $10,000
-
-- Average cost of a crime committed by a recidivist: $50,000
-- Average compensation for wrongful incarceration: $50,000
+Average cost of keeping someone in prison for a year:
+$40,000
+
Average cost of arrest and trial for a recidivist: $10,000
+
Average cost of a crime committed by a recidivist:
+$50,000
Average compensation for wrongful incarceration: $50,000
True negative Cost: “An individual is predicted
+not to recidivate and indeed does not recidivate.” - No cost
True positive Cost: “An individual is predicted
+to recidivate and then does recidivate.” - $40k
False negative Cost: “An individual is predicted
+not to recidivate but then does recidivate.” - $100k
False positive Cost: “An individual is predicted
+to recidivate but does not actually recidivate.” - $50k
-
True negative Cost: “An individual is predicted not
-to recidivate and indeed does not recidivate.” - No cost True
-positive Cost: “An individual is predicted to recidivate and
-then does recidivate.” - $40k False negative Cost: “An
-individual is predicted not to recidivate but then does recidivate.” -
-$100k False positive Cost: “An individual is predicted
-to recidivate but does not actually recidivate.” - $50k
whichThreshold.cost <- whichThreshold %>%
dplyr::select(starts_with("Count"), Threshold) %>%
mutate(Cost.TN = Count_TN*0,
@@ -3874,11 +3877,11 @@ Cost-Benifit Analysis
theme(legend.title = element_text(size = 9),
legend.text = element_text(size = 8)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
-
+
threshold.cost <- whichThreshold.cost%>%
filter(Variable == 'Cost.all')
-#-> lowest cost when threshold = 0.49
+#-> lowest cost when threshold = 0.50
Conclusion
diff --git a/script/Assignment4.md b/script/Assignment4.md
index 0b3e5fb..f1472b6 100644
--- a/script/Assignment4.md
+++ b/script/Assignment4.md
@@ -1316,10 +1316,10 @@ filter(testProbs.thresholds, Threshold == .5) %>%
- Average cost of a crime committed by a recidivist: $50,000
- Average compensation for wrongful incarceration: $50,000
-**True negative Cost**: "An individual is predicted not to recidivate and indeed does not recidivate." - No cost
-**True positive Cost**: "An individual is predicted to recidivate and then does recidivate." - $40k
-**False negative Cost**: "An individual is predicted not to recidivate but then does recidivate." - $100k
-**False positive Cost**: "An individual is predicted to recidivate but does not actually recidivate." - $50k
+- **True negative Cost**: "An individual is predicted not to recidivate and indeed does not recidivate." - No cost
+- **True positive Cost**: "An individual is predicted to recidivate and then does recidivate." - $40k
+- **False negative Cost**: "An individual is predicted not to recidivate but then does recidivate." - $100k
+- **False positive Cost**: "An individual is predicted to recidivate but does not actually recidivate." - $50k
```r
@@ -1353,7 +1353,7 @@ whichThreshold.cost %>%
threshold.cost <- whichThreshold.cost%>%
filter(Variable == 'Cost.all')
-#-> lowest cost when threshold = 0.49
+#-> lowest cost when threshold = 0.50
```
diff --git a/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.RData b/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.RData
deleted file mode 100644
index 54e01f8..0000000
Binary files a/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.RData and /dev/null differ
diff --git a/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.RData b/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.RData
new file mode 100644
index 0000000..192c9be
Binary files /dev/null and b/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.RData differ
diff --git a/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.rdb b/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.rdb
similarity index 100%
rename from script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.rdb
rename to script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.rdb
diff --git a/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.rdx b/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.rdx
similarity index 91%
rename from script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.rdx
rename to script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.rdx
index 4258680..82d238a 100644
Binary files a/script/Assignment4_cache/html/cost_59f3f7e55f8a596554596643d0fb42b5.rdx and b/script/Assignment4_cache/html/cost_a91d77ba3ae2a37f0250edde7dc94389.rdx differ
diff --git a/script/Assignment4_files/figure-html/cost-1.png b/script/Assignment4_files/figure-html/cost-1.png
index a2d1e71..d460f0a 100644
Binary files a/script/Assignment4_files/figure-html/cost-1.png and b/script/Assignment4_files/figure-html/cost-1.png differ
diff --git a/script/Assignment4_files/figure-html/roc_curve-1.png b/script/Assignment4_files/figure-html/roc_curve-1.png
index 46302bb..84dc484 100644
Binary files a/script/Assignment4_files/figure-html/roc_curve-1.png and b/script/Assignment4_files/figure-html/roc_curve-1.png differ
diff --git a/script/Assignment4_files/figure-html/roc_curve-2.png b/script/Assignment4_files/figure-html/roc_curve-2.png
index 1d02b88..fe13b54 100644
Binary files a/script/Assignment4_files/figure-html/roc_curve-2.png and b/script/Assignment4_files/figure-html/roc_curve-2.png differ
diff --git a/script/Assignment4_files/figure-html/roc_curve-3.png b/script/Assignment4_files/figure-html/roc_curve-3.png
index 2db707d..dbb7597 100644
Binary files a/script/Assignment4_files/figure-html/roc_curve-3.png and b/script/Assignment4_files/figure-html/roc_curve-3.png differ