Maps

1. EV charging station locations

d = readRDS('EVstations.rds')
colnames(d) = tolower(colnames(d))

d = d %>%
  select(-matches('cng|lng|lpg|hydrogen|e85|bd[.]blends|french|^ng[.]|plus4|level1')) %>%
  filter(fuel.type.code=='ELEC') %>%
  rename(lev2=ev.level2.evse.num, 
         lev3=ev.dc.fast.count, 
         network = ev.network, 
         lat = latitude, 
         lon = longitude) %>%
  mutate(status = case_when(status.code=='E' ~ 'avail', 
                            status.code=='P' ~ 'planned', 
                            status.code=='T' ~ 'temp.unavail', 
                            TRUE ~ ''), 
         lev2 = ifelse(is.na(lev2), 0, lev2), 
         lev3 = ifelse(is.na(lev3), 0, lev3)) %>%
  mutate(network = ifelse(network=='Tesla Destination', 'Tesla', network), 
         network = gsub('é', 'E', network), 
         network = gsub('É', 'E', network))

Plot the locations of Level 3 Charging Stations and this time make sure we can easily distinguish what network each charging station belongs to. Consider the top 5 networks (Tesla, Electrify America, ChargePoint, eVgo, and Non-Networked) separately and group all other networks together as “Other”.

library(tidyverse)
library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
# Read and prepare the data
d <- readRDS('EVstations.rds')
d <- d %>%
  filter(`Fuel.Type.Code` == 'ELEC') %>%
  mutate(
    network = case_when(
      `EV.Network` %in% c('Tesla', 'Tesla Destination') ~ 'Tesla',
      `EV.Network` == 'Electrify America' ~ 'Electrify America',
      `EV.Network` == 'ChargePoint Network' ~ 'ChargePoint',
      `EV.Network` == 'eVgo Network' ~ 'eVgo',
      `EV.Network` == 'Non-Networked' ~ 'Non-Networked',
      TRUE ~ 'Other'  # Handles all other networks not specified
    ),
    latitude = as.numeric(`Latitude`),
    longitude = as.numeric(`Longitude`)
  ) %>%
  filter(`EV.DC.Fast.Count` > 0)

# Load US map data
us_map <- map_data("state")

# Base map layer
base_map <- ggplot() +
  geom_polygon(data = us_map, aes(x = long, y = lat, group = group), fill = "white", color = "black") +
  coord_fixed(1.3, xlim = c(-125, -66), ylim = c(25, 50)) +
  labs(title = "Locations of Level 3 EV Charging Stations",
       subtitle = "Color indicates network",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

# Overlay EV charging stations
# Make sure to specify data and aesthetics specifically for this layer
ev_map <- base_map +
  geom_point(data = d, aes(x = longitude, y = latitude, color = network), alpha = 0.6, size = 2) +
  scale_color_manual(values = c(
    'Tesla' = 'red',
    'Electrify America' = 'blue',
    'ChargePoint' = 'green',
    'eVgo' = 'yellow',
    'Non-Networked' = 'purple',
    'Other' = 'gray'
  ))

# Print the map with EV stations
print(ev_map)

2. Animation

Make a plot like above but facet by network, and animate it using open.date as the time variable. Since there are ~1458 unique open.dates, we probably don’t want one frame per date. So create a month column and use one frame per month. When a charging station appears, we want it to stay for the rest of the animation (not go away the next month). Upload your gif to Gradescope in the assignment denoted “GIF”, and at the end of your code below, include a screenshot of the animation taken sometime halfway through the animation. Tip for knitting your PDF: as discussed in class, comment out your animate code after creating a .gif so that your document knits more quickly.

Code

# library(tidyverse)
# library(lubridate)
# library(gganimate)
# library(maps)
# 
# # Read and prepare the data
# d <- readRDS('EVstations.rds')
# d <- d %>%
#   filter(`Fuel.Type.Code` == 'ELEC') %>%
#   mutate(
#     network = case_when(
#       `EV.Network` %in% c('Tesla', 'Tesla Destination') ~ 'Tesla',
#       `EV.Network` == 'Electrify America' ~ 'Electrify America',
#       `EV.Network` == 'ChargePoint Network' ~ 'ChargePoint',
#       `EV.Network` == 'eVgo Network' ~ 'eVgo',
#       `EV.Network` == 'Non-Networked' ~ 'Non-Networked',
#       TRUE ~ 'Other'
#     ),
#     latitude = as.numeric(`Latitude`),
#     longitude = as.numeric(`Longitude`),
#     open.date = as.Date(`Open.Date`, format = "%Y-%m-%d")
#   ) %>%
#   drop_na(open.date) %>%
#   mutate(month = floor_date(open.date, "month")) %>%
#   filter(`EV.DC.Fast.Count` > 0)
# 
# # Create a sequence of all months from the first to the last month data is available
# all_months <- seq(min(d$month, na.rm = TRUE), max(d$month, na.rm = TRUE), by = "1 month")
# 
# # Create an expanded grid of all months for each station, effectively duplicating each station for each month up to current
# d <- d %>%
#   group_by(id = row_number()) %>%
#   complete(month = all_months) %>%
#   fill(longitude, latitude, network, .direction = "down") %>%
#   filter(!is.na(longitude)) %>%
#   ungroup()
# 
# # Load US map data
# us_map <- map_data("state")
# 
# # Create base map layer
# base_map <- ggplot() +
#   geom_polygon(data = us_map, aes(x = long, y = lat, group = group), fill = "gray90", colour = "white") +
#   coord_fixed(1.3) +  # This keeps the aspect ratio of the map
#   theme_minimal()
# 
# # Combine base map with charging stations data
# p <- base_map +
#   geom_point(data = d, aes(x = longitude, y = latitude, color = network), size = 2) +
#   scale_color_manual(values = c(
#     'Tesla' = 'red',
#     'Electrify America' = 'blue',
#     'ChargePoint' = 'green',
#     'eVgo' = 'yellow',
#     'Non-Networked' = 'purple',
#     'Other' = 'gray'
#   )) +
#   facet_wrap(~network) +
#   transition_time(month) +
#   labs(title = "EV Charging Stations Opening Over Time: {closest_state}",
#        subtitle = "Data aggregated by month, cumulative visibility",
#        x = "Longitude", y = "Latitude") +
#   ease_aes('linear') +
#   shadow_mark(past = TRUE, future = FALSE)
# 
# # Render the animation
# anim <- animate(p, duration = 20, fps = 5, end_pause = 20, width = 800, height = 600, renderer = gifski_renderer())
# 
# # Save the animation as a GIF
# anim_save("EV_stations_animation.gif", animation = anim)
# 
# # Output animation object for preview in R environment
# anim

Screenshot here

gif image
gif image

3. Interactive Plot

Plot EV charging station locations on an interactive map with leaflet. Create labels that show the EV network, station name, street address, city, state, zip, and the number Level 2 and Level 3 charging stations. Size the dots by the number of chargers at that location. Since your computer might have issues showing all the points at once, plot just the locations in CT.

library(tidyverse)
library(leaflet)

# Read the data
d <- readRDS('EVstations.rds')

# Filter for stations in Connecticut (CT)
ct_stations <- d %>%
  filter(`State` == 'CT', `Fuel.Type.Code` == 'ELEC') %>%
  mutate(
    Total_Chargers = `EV.Level2.EVSE.Num` + `EV.DC.Fast.Count`,
    Label = sprintf(
      "<strong>Network:</strong> %s<br/>
       <strong>Station Name:</strong> %s<br/>
       <strong>Address:</strong> %s, %s, %s %s<br/>
       <strong>Level 2 Chargers:</strong> %s<br/>
       <strong>Level 3 Chargers:</strong> %s",
      `EV.Network`, `Station.Name`, `Street.Address`, `City`, `State`, `ZIP`,
      `EV.Level2.EVSE.Num`, `EV.DC.Fast.Count`
    )
  )

# Create the Leaflet map
map <- leaflet(ct_stations) %>%
  addTiles() %>%
  addCircleMarkers(
    ~Longitude, ~Latitude, radius = ~sqrt(Total_Chargers)*2,
    popup = ~Label, color = '#1E90FF', fillOpacity = 0.8, stroke = FALSE,
    labelOptions = labelOptions(style = list('font-weight' = 'normal', padding = '3px', 'border-radius' = '4px'))
  ) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  setView(lng = -72.7575, lat = 41.7759, zoom = 8)  # Center map on Connecticut

# Print the map
map

Shiny

4. Basic Shiny App

Do the following to create a basic shiny app and publish it to shinyapps.io

  • Configure R Studio to communicate with shinyapps.io by following the directions from this page https://shiny.rstudio.com/articles/shinyapps.html up to and including the “Method 1” section. You will show you how to install rsconnect, create a shinyapps.io account, and set up communication between R Studio and shinyapps.io. Note that you will only have to do this once.
  • Create a new shiny app like we did in class by clicking File, New File, Shiny Web App. Name the app hw03app1, choose Single File, choose an appropriate directory, and click Create. An app.R file will open in R Studio.
  • Click Run App to see the app. The app should appear in the Viewer tab in R Studio, or in its own window.
  • In the upper right of the app, click Publish. Make sure the correct files and account are selected. Click Publish in the lower right.

Paste the url for your app here. The url will be of the form https://yourusername.shinyapps.io/yourappname/

url here https://pantelicvv.shinyapps.io/hw03app1/

5. Shiny App for EV Charging Stations

Create a Shiny app that contains a leaflet map like in #3 and allows the user to select

  • a state (default: CT)
  • a date range (default: 2015-01-01 to present)
  • a network (default: all networks), and
  • Level 2 or Level 3 (default: Level 3)

and displays a leaflet plot of Level 2 or 3 EV charging stations for that state, for that network, and with an open.date in that date range. Note that you will have to copy the EVstations.rds file to the folder that contains the Shiny app.

Different types of user inputs can be found here https://shiny.rstudio.com/gallery/widget-gallery.html

Publish your shiny app to shinyapps.io, and paste the url for your app here:

Code

# library(shiny)
# library(leaflet)
# library(dplyr)
# library(lubridate)
# 
# # Load data outside the server function to improve performance and enable UI setup
# data <- readRDS("EVstations.rds")
# 
# # Get unique values for states and networks for the dropdown menus
# states <- unique(data$State)
# networks <- unique(data$`EV.Network`)
# networks <- c("All", networks[!is.na(networks)])  # Add "All" option and remove NA values
# 
# # Define UI for application
# ui <- fluidPage(
#   titlePanel("EV Charging Stations Interactive Map"),
#   
#   # Sidebar layout with input controls
#   sidebarLayout(
#     sidebarPanel(
#       selectInput("state", "Select a State:",
#                   choices = states, selected = "CT"),
#       dateRangeInput("dateRange", "Select Date Range:",
#                      start = "2015-01-01", end = Sys.Date()),
#       selectInput("network", "Select Network:",
#                   choices = networks, selected = "All"),
#       radioButtons("level", "Charger Level:",
#                    choices = c("Level 2", "Level 3"), selected = "Level 3")
#     ),
#     
#     # Main panel for displaying the leaflet map
#     mainPanel(
#       leafletOutput("map")
#     )
#   )
# )
# 
# # Define server logic
# server <- function(input, output) {
#   output$map <- renderLeaflet({
#     # Filter data based on user input
#     filtered_data <- data %>%
#       filter(`State` == input$state,
#              `Fuel.Type.Code` == 'ELEC',
#              as.Date(`Open.Date`) >= input$dateRange[1],
#              as.Date(`Open.Date`) <= input$dateRange[2]) %>%
#       filter(if(input$network != "All") `EV.Network` == input$network else TRUE) %>%
#       filter(if(input$level == "Level 2") `EV.Level2.EVSE.Num` > 0 else `EV.DC.Fast.Count` > 0)
#     
#     # Create the map
#     if (nrow(filtered_data) > 0) {
#       leaflet(filtered_data) %>%
#         addTiles() %>%
#         addCircleMarkers(
#           ~Longitude, ~Latitude, 
#           popup = ~paste(`Station.Name`, `Street.Address`, sep = "<br>"),
#           radius = ~sqrt(`EV.Level2.EVSE.Num` + `EV.DC.Fast.Count`) * 2,
#           color = ~ifelse(`EV.Network` == input$network, 'red', 'blue'),
#           fillOpacity = 0.8, stroke = FALSE
#         )
#     } else {
#       leaflet() %>% 
#         addTiles() %>%
#         addMarkers(lng = -72.675, lat = 41.765, popup = "No data available for this selection")
#     }
#   })
# }
# 
# # Run the application
# shinyApp(ui = ui, server = server)

url here

https://pantelicvv.shinyapps.io/sds361/

Bootstrapping s.e. for \(\hat{\lambda}\)

dc = readRDS('tracts.and.census.with.EV.stations.rds')
dc = dc@data ## keep just the data frame, not the polygons

dc = dc %>% 
  
  ## Change NAs to 0s, and 
  ## since charging stations seem to come in pairs, create new columns
  ## for pairs of charging stations
  mutate(lev2 = ifelse(is.na(lev2), 0, lev2), 
         lev3 = ifelse(is.na(lev3), 0, lev3), 
         lev2pairs = round(lev2/2),
         lev3pairs = round(lev3/2)) %>%
  
  ## Keep tracts with at least one charging station, 
  ## and get rid of a couple of outliers
  #filter((lev2!=0 | lev3!=0) & lev2<=50) %>%
  mutate(log.house.value = log(house.value))

6. s.e. for \(\hat{\lambda}\), all observations

In class we found the bootstrap s.e. for \(\hat{\lambda}\) when log.house.value = 12. Find the bootstrap s.e. for all values of log.house.value in the data. Hint: long format for a data frame is your friend. The s.e. for some values of log.house.value are small and some are large. Why?

# library(data.table)
# library(furrr)
# library(parallel)
# 
# # Set up parallel processing
# plan(multisession, workers = 4)
# 
# # Define the bootstrap function ensuring all necessary columns are visible
# bootstrap_analysis <- function(log_value, data, n_bootstrap = 1000) {
# 
#   results <- replicate(n_bootstrap, {
#     # Sample rows for the bootstrap sample
#     sample_rows <- sample(1:nrow(data), replace = TRUE)
#     # Create bootstrap sample
#     bootstrap_sample <- data[sample_rows, ]
#     
#     # Fit model
#     model <- glm(lev2 ~ log.house.value, data = bootstrap_sample, family = 'poisson')
#     
#     # Predict at the current log.house.value
#     predict_value <- predict(model, newdata = data.frame(log.house.value = log_value), type = "response")
#     
#     # Return results as a list
#     list(coef = coef(model)["log.house.value"], predict = predict_value)
#   }, simplify = FALSE)  # Use simplify = FALSE to handle different data types properly
# 
#   # Compute standard deviations
#   coef_se <- sd(sapply(results, function(x) x$coef))
#   predict_se <- sd(sapply(results, function(x) x$predict))
#   
#   list(coef_se = coef_se, predict_se = predict_se)
# }
# 
# # Example usage for a single log.house.value
# results <- bootstrap_analysis(log_value = 12, data = dc, n_bootstrap = 1000)
# 
# 
# results_list <- future_map_dfr(unique(dc$log.house.value), ~bootstrap_analysis(.x, dc, 1000))
# 
# print(results_list)

In terms of the standard error variability, multiple factors can cause it. For instance, high-end properties might include extra features and lead to greater variability in market pricing compared to more standardized, middle-range homes. In addition, if certain ranges of log.house.value have fewer data points, the bootstrap samples for these values will be less reliable and exhibit higher variability. Certain log.house.value values might act as outliers or leverage points and disproportionately influence the model’s fit. Bootstrap methods, by repeatedly sampling, could amplify the effects of these points and result in larger standard errors in their vicinity. Market dynamics in terms of economic factors, and demand-supply dynamics could also be affecting the variability in the standard errors.

7. Compare s.e. with bootstrapped s.e.

Plot \(\hat{\lambda} + 2 \, s.e.\) and \(\hat{\lambda} - 2 \, s.e.\) along with the curve for \(\hat{\lambda}\), and \(\pm 2 \,s.e.\) that you get from using predict(glm1, newdata = dc, type = 'response', se.fit = T) from the Poisson (glm1), Quasipoisson (glm2) and Negative binomial (glm3) models. Compare the results.

# Remove rows where log.house.value is NA or infinite
dc <- dc[!is.na(dc$log.house.value) & is.finite(dc$log.house.value), ]

# Now check and make sure that dc has no NA or non-finite values in log.house.value
summary(dc$log.house.value)  # This will give you a summary including min and max values
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.21   11.71   12.17   12.23   12.73   14.51
# Assuming the data is cleaned
newdata <- data.frame(log.house.value = seq(from = min(dc$log.house.value), to = max(dc$log.house.value), length.out = 100))

# Fit Poisson Model
glm1 <- glm(lev2 ~ log.house.value, data = dc, family = poisson())

# Fit Quasi-Poisson Model
glm2 <- glm(lev2 ~ log.house.value, data = dc, family = quasipoisson())

# Fit Negative Binomial Model
library(MASS)  # Load MASS for negative binomial
glm3 <- glm.nb(lev2 ~ log.house.value, data = dc)
# Assuming you want to make predictions across the same range as your original data
newdata <- data.frame(log.house.value = seq(min(dc$log.house.value), max(dc$log.house.value), length.out = 100))

# Predict with SE for Poisson
preds_poisson <- predict(glm1, newdata = newdata, type = "response", se.fit = TRUE)

# Predict with SE for Quasi-Poisson
# Note: Quasi-Poisson doesn't change the prediction; it adjusts the SE.
preds_quasipoisson <- predict(glm2, newdata = newdata, type = "response", se.fit = TRUE)

# Predict with SE for Negative Binomial
preds_negbin <- predict(glm3, newdata = newdata, type = "response", se.fit = TRUE)

# Calculate intervals
newdata$fit_poisson <- preds_poisson$fit
newdata$upper_poisson <- preds_poisson$fit + 2 * preds_poisson$se.fit
newdata$lower_poisson <- preds_poisson$fit - 2 * preds_poisson$se.fit

newdata$fit_quasipoisson <- preds_quasipoisson$fit
newdata$upper_quasipoisson <- preds_quasipoisson$fit + 2 * preds_quasipoisson$se.fit
newdata$lower_quasipoisson <- preds_quasipoisson$fit - 2 * preds_quasipoisson$se.fit

newdata$fit_negbin <- preds_negbin$fit
newdata$upper_negbin <- preds_negbin$fit + 2 * preds_negbin$se.fit
newdata$lower_negbin <- preds_negbin$fit - 2 * preds_negbin$se.fit
library(ggplot2)

g <- ggplot(newdata, aes(x = log.house.value)) +
  geom_line(aes(y = fit_poisson, color = "Poisson")) +
  geom_ribbon(aes(ymin = lower_poisson, ymax = upper_poisson), alpha = 0.2, fill = "blue") +
  geom_line(aes(y = fit_quasipoisson, color = "Quasi-Poisson")) +
  geom_ribbon(aes(ymin = lower_quasipoisson, ymax = upper_quasipoisson), alpha = 0.2, fill = "red") +
  geom_line(aes(y = fit_negbin, color = "Negative Binomial")) +
  geom_ribbon(aes(ymin = lower_negbin, ymax = upper_negbin), alpha = 0.2, fill = "green") +
  labs(title = "Model Predictions and Confidence Intervals",
       y = "Predicted Values (λ)",
       x = "Log of House Value") +
  scale_color_manual(values = c("Poisson" = "blue", "Quasi-Poisson" = "red", "Negative Binomial" = "green"))

print(g)

When assumptions of the model are closely met, the bootstrapped standard errors match the standard errors in the summary output. If the assumptions aren’t met, they can be different. For the Poisson model the assumption of mean=variance is not met, so the bootstrapped standard errors should be different from those given by the model summary output. The negative binomial model assumptions are more closely met, so they should match more closely in this case. The quasi-poisson model is a generalization of the Poisson model that allows for overdispersion, so the bootstrapped standard errors should closer to this model’s predictions than standard Poisson’s.

Bootstrapping

Suppose that we obtain a bootstrap sample from a set of \(n\) observations. We will now show that the probability that the \(j\)th observation is part of a bootstrap sample is 0.632 for large enough \(n\).

8. 1st and 2nd observation are in the sample

What is the probability that the 1st bootstrap observation is not the \(j\)th observation from the original sample? The 2nd bootstrap observation?

Each observation in the original dataset, including the jth observation, has an equal probability of being selected in each draw. Since there are n observations in total, the probability of selecting the jth observation in any single draw is \(1/n\). Therefore, the probability that the jth observation is not selected in the first bootstrap draw is \(1-(1/n)\). Since each draw in a bootstrap sample is independent (because of replacement), the probability that the jth observation is not selected in the second draw is the same as in the first draw, which is also \(1-(1/n)\).

9. jth observation is not in the sample

Argue that the probability that the \(j\)th observation is not in the bootstrap sample is \((1 - 1/n)^n\).

Since the bootstrap sample involves drawing with replacement, each draw is independent of the others. When calculating the probability that the jth observation is not selected across all n draws, we multiply the probabilities of not selecting it in each individual draw. This is because the outcome of each draw (whether or not we pick the jth observation) does not affect the outcomes of other draws. The probability that the jth observation is not selected in any of the n draws is then \((1 - 1/n)^n\). This arises from the principle of independence.

10. jth observation is in the sample

What is the probability that the \(j\)th observation is in the bootstrap sample? Find the limit of your answer as \(n\rightarrow \infty\). Argue that for large enough \(n\), the probability that the \(j\)th observation is in the bootstrap sample is around 0.632.

The probability that the jth observation is in the sample is the complement of it not being selected in any of the draws. Therefore: \(1-(1 - 1/n)^n\). To find the limit of this probability as n approaches infinity, we use a known limit from calculus: \(\lim_{n \to \infty} (1 - 1/n)^n = 1/e\). Therefore, the limit of the probability that the jth observation is in the sample as n approaches infinity is \(1-1/e \approx 0.632\). This is the probability that the jth observation is in the bootstrap sample for large enough n. As n becomes large, the expression closely approximates \(1/e\) as noted above. Therefore, the probability that the jth observation is in the sample asymptotically approaches this value, which is about 0.632. This is why, for large n, any specific observation has about a 63.2% chance of appearing in a bootstrap sample. This is a consequence of the properties of exponential decay and the uniform likelihood of selection in each independent draw.