---
title: "Code for Session 4"
author: "Dr. Richard M. Crowley"
date: ""
output:
  html_notebook
---

Note that the directories used to store data are likely different on your computer, and such references will need to be changed before using any such code.

```{r helpers, warning=FALSE, message=F}
library(knitr)
library(kableExtra)
library(tidyverse)
library(broom)
html_df <- function(text, cols=NULL, col1=FALSE, full=F) {
  if(!length(cols)) {
    cols=colnames(text)
  }
  if(!col1) {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full)
  } else {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full) %>%
      column_spec(1,bold=T)
  }
}
```

```{r, message=F, warning=F}
# Import Walmart data from Session_Project
library(lubridate)
preprocess_data <- function(df) {
  # Merge the data together (Pulled from outside of function -- "scoping")
  df <- inner_join(df, weekly.stores)
  df <- inner_join(df, weekly.features[,1:11])
  
  # Compress the weird markdown information to 1 variable
  df$markdown <- 0
  df[!is.na(df$MarkDown1),]$markdown <- df[!is.na(df$MarkDown1),]$MarkDown1
  df[!is.na(df$MarkDown2),]$markdown <- df[!is.na(df$MarkDown2),]$MarkDown2
  df[!is.na(df$MarkDown3),]$markdown <- df[!is.na(df$MarkDown3),]$MarkDown3
  df[!is.na(df$MarkDown4),]$markdown <- df[!is.na(df$MarkDown4),]$MarkDown4
  df[!is.na(df$MarkDown5),]$markdown <- df[!is.na(df$MarkDown5),]$MarkDown5
  
  # Fix dates and add useful time variables
  df$date <- as.Date(df$Date)
  df$week <- week(df$date)
  df$year <- year(df$date)
  df$month <- month(df$date)
  
  df
}

weekly <- read.csv("../../Data/WMT_train.csv", stringsAsFactors=FALSE)
weekly.features <- read.csv("../../Data/WMT_features.csv", stringsAsFactors=FALSE)
weekly.stores <- read.csv("../../Data/WMT_stores.csv", stringsAsFactors=FALSE)
df <- preprocess_data(weekly)
#df$id <- df$Store *10000 + df$week * 100 + df$Dept
df <- df %>%
  group_by(Store, Dept) %>% 
  mutate(store_avg=mean(Weekly_Sales, rm.na=T)) %>%
  ungroup()
```

```{r}
# Create the binary variable from Walmart sales data
df$double <- ifelse(df$Weekly_Sales > df$store_avg*2,1,0)
fit <- glm(double ~ IsHoliday, data=df, family=binomial)
tidy(fit)
```

```{r}
# Automating the above:
exp(coef(fit))
```

```{r}
c(-3.44, -2.9)
```

```{r}
a <- exp(-3.44)
b <- exp(-2.9)
c(a/(1+a), b/(1+b))
```

```{r}
model2 <- glm(double ~ IsHoliday + Temperature + Fuel_Price, data=df, family=binomial)
summary(model2)
```

```{r}
# Odds
exp(coef(model2))
```

```{r}
# Typical September days
hday_sep    <- mean(predict(model2, filter(df, IsHoliday, month==9), type="response"))
no_hday_sep <-  mean(predict(model2, filter(df, !IsHoliday, month==9), type="response"))
# Typical December days
hday_dec    <- mean(predict(model2, filter(df, IsHoliday, month==12), type="response"))
no_hday_dec <-  mean(predict(model2, filter(df, !IsHoliday, month==12), type="response"))

html_df(data.frame(Month=c(9,9,12,12),
                   IsHoliday=c(FALSE,TRUE,FALSE,TRUE),
                   Probability=c(no_hday_sep, hday_sep, no_hday_dec, hday_dec)))
```

```{r, message=F}
percent <- function(x) {
  paste0(round(x*100,2),"%")
}
library(rlang)
# Another use for a quosure :)
# You could always just use quoted arguements here though
extract_margin <- function(m, x) {
  x <- enquo(x)
  percent(filter(summary(m), factor==quo_text(x))$AME)
}
```

```{r, message=F}
# Calculate AME marginal effects
library(margins)
m <- margins(model2)
m
```

```{r}
summary(m) %>%
  html_df()
```

```{r, fig.height=4, fig.width=6}
plot(m, which=summary(m)$factor)
```

```{r}
margins(model2, at = list(IsHoliday = c(TRUE, FALSE)),
        variables = c("Temperature", "Fuel_Price")) %>%
  summary() %>%
  html_df()
```

```{r}
margins(model2, at = list(Temperature = c(0, 25, 50, 75, 100)),
        variables = c("IsHoliday")) %>%
  summary() %>%
  html_df()
```

```{r}

# load data
df <- read.csv("../../Data/Shipping/shipping_dataset.csv", stringsAsFactors = F, na="NA")
df_ships <- read.csv("../../Data/Shipping/index.csv", stringsAsFactors = F, na="NA")
df_ports <- read.csv("../../Data/Shipping/port_dataset.csv", stringsAsFactors = F, na="NA")
typhoon <- read.csv("../../Data/Shipping/typhoon.csv", stringsAsFactors = F)

# rename lat and lon to avoid conflicts
df_ports <- df_ports %>% rename(port_lat=lat, port_lon=lon)

# Join in ships and ports
df <- left_join(df, df_ships)
df <- df %>% filter(!imo == 9528029 | reported_destination == "AU ADL")
df <- left_join(df,df_ports, by=c("left_port"="port"))
df <- df %>% rename(left_lat=port_lat, left_lon=port_lon)
df <- left_join(df,df_ports, by=c("arrive_port"="port"))
df <- df %>% rename(arrive_lat=port_lat, arrive_lon=port_lon)

# fix dates
df$last_update = as.POSIXct(df$last_update, tz="UTC",origin="1970-01-01")
df$left_time = as.POSIXct(df$left_time, tz="UTC",origin="1970-01-01")
df$arrive_time = as.POSIXct(df$arrive_time, tz="UTC",origin="1970-01-01")

# Fix typhoon dates
typhoon$date <- as.POSIXct(paste(typhoon$date,typhoon$time), format="%Y%b%d %H%M%S", tz="UTC", origin="1970-01-01")
# Fix typhoon longitude -- they have lon * -1 in the data
typhoon$lon <- typhoon$lon * -1

typhoon_all <- typhoon
# filter to dates in sample
typhoon <- typhoon %>% filter(date > as.POSIXct("2018-08-30", tz="UTC"))

df <- df %>% rename(frame=run)

df_all <- df[df$frame != 32,]
df <- df %>% filter(last_update > as.POSIXct("2018-08-30", tz="UTC"))

library(plotly)  # needed for the toRGB() function
geo <- list(
  showland = TRUE,
  showlakes = TRUE,
  showcountries = TRUE,
  showocean = TRUE,
  countrywidth = 0.5,
  landcolor = toRGB("grey90"),
  lakecolor = toRGB("aliceblue"),
  oceancolor = toRGB("aliceblue"),
  projection = list(
    type = 'orthographic',  # detailed at https://plot.ly/r/reference/#layout-geo-projection
    rotation = list(
      lon = 100,
      lat = 1,
      roll = 0
    )
  ),
  lonaxis = list(
    showgrid = TRUE,
    gridcolor = toRGB("gray40"),
    gridwidth = 0.5
  ),
  lataxis = list(
    showgrid = TRUE,
    gridcolor = toRGB("gray40"),
    gridwidth = 0.5
  )
)

typhoon_Aug31 <- typhoon[typhoon$date > as.POSIXct("2018-08-31 00:00:00", tz="UTC") & typhoon$date < as.POSIXct("2018-09-01 00:00:00", tz="UTC"),]
typhoon_Jebi <- typhoon_all[typhoon_all$typhoon_name == "JEBI_Y",]
df_Aug31 <- df[df$frame == 1,]
```

```{r, warning=F, message=F, fig.height=5, fig.width=8}
library(plotly)  # for plotting
library(RColorBrewer)  # for colors
# plot with boats, ports, and typhoons
# Note: geo is defined in the appendix -- it controls layout
palette = brewer.pal(8, "Dark2")[c(1,8,3,2)]
p <- plot_geo(colors=palette) %>%
  add_markers(data=df_ports, x = ~port_lon, y = ~port_lat, color = "Port") %>%
  add_markers(data=df_Aug31, x = ~lon, y = ~lat, color = ~ship_type,
              text=~paste('Ship name',shipname)) %>%
  add_markers(data=typhoon_Aug31, x = ~lon, y = ~lat, color="TYPHOON",
              text=~paste("Name", typhoon_name)) %>%
   layout(showlegend = TRUE, geo = geo,
    title = 'Singaporean owned container and tanker ships, August 31, 2018')
p
```

```{r, warning=F, message=F, fig.height=5, fig.width=8}
# plot with boats and typhoons
palette = brewer.pal(8, "Dark2")[c(1,3,2)]
p <- plot_geo(colors=palette) %>%
  add_markers(data=df_all[df_all$frame == 14,], x = ~lon, y = ~lat, color = ~ship_type,
              text=~paste('Ship name',shipname)) %>%
  add_markers(data=typhoon_Jebi, x = ~lon, y = ~lat, color="Typhoon Jebi", text=~paste("Name", typhoon_name, "</br>Time: ", date)) %>%
   layout(
    showlegend = TRUE, geo = geo,
    title = 'Singaporean container/tanker ships, September 4, 2018, evening')
p
```

```{r, warning=F, message=F, fig.height=5, fig.width=9}
library(leaflet)
library(leaflet.extras)

# typhoon icons
icons <- pulseIcons(color='red',
  heartbeat = ifelse(typhoon_Jebi$intensity_vmax > 150/1.852, 0.8,
    ifelse(typhoon$intensity_vmax < 118/1.852, 1.6, 1.2)),
  iconSize=ifelse(typhoon_Jebi$intensity_vmax > 150/1.852, 12,
    ifelse(typhoon_Jebi$intensity_vmax < 118/1.852, 3, 8)))

# ship icons
shipicons <- iconList(
  ship = makeIcon("../Figures/ship.png", NULL, 18, 18)
)
leaflet() %>%
  addTiles() %>% 
  setView(lng = 136, lat = 34, zoom=4) %>%
  addPulseMarkers(data=typhoon_Jebi[seq(1,nrow(typhoon_Jebi),5),], lng=~lon,
                  lat=~lat, label=~date, icon=icons) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax > 150/1.852,],
    lng=~lon, lat=~lat,stroke = TRUE, radius=12, color="red", label=~date) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax <= 150/1.852 &
    typhoon_Jebi$intensity_vmax > 118/1.852,], lng=~lon, lat=~lat,
    stroke = TRUE, radius=12, color="red", label=~date) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax <=118/1.852,],
    lng=~lon, lat=~lat, stroke = TRUE, radius=12, color="red", label=~date) %>%
  addMarkers(data=df_all[df_all$frame == 14,], lng=~lon, lat=~lat,
             label=~shipname, icon=shipicons)
```

```{r}
# Load datas
df3 <- read.csv("../../Data/Shipping/combined.csv", stringsAsFactors = F)

library(geosphere)

# Calculate distance to and from ports
df3$dist_toport <- distVincentyEllipsoid(as.matrix(df3[,c("lon","lat")]), as.matrix(df3[,c("arrive_lon","arrive_lat")]))
df3$dist_fromport <- distVincentyEllipsoid(as.matrix(df3[,c("lon","lat")]), as.matrix(df3[,c("left_lon","left_lat")]))
df3 <- df3 %>%
  arrange(imo, last_update) %>%
  group_by(imo) %>%
  mutate(delayed=ifelse(difftime(arrive_time, lead(arrive_time), units="days") > 1/8 & arrive_time > last_update & left_port == lead(left_port),1,0)) %>%
  ungroup()
```

```{r}
library(geosphere)
x <- as.matrix(df3[,c("lon","lat")])  # ship location
y <- as.matrix(df3[,c("ty_lon","ty_lat")]) # typhoon location

df3$dist_typhoon <- distVincentyEllipsoid(x, y) / 1000
```

```{r}
df3$typhoon_500 = ifelse(df3$dist_typhoon < 500 & 
                         df3$dist_typhoon >= 0, 1, 0)
df3$typhoon_1000 = ifelse(df3$dist_typhoon < 1000 &
                          df3$dist_typhoon >= 500, 1, 0)
df3$typhoon_2000 = ifelse(df3$dist_typhoon < 2000 &
                          df3$dist_typhoon >= 1000, 1, 0)
```

```{r}
fit1 <- glm(delayed ~ typhoon_500 + typhoon_1000 + typhoon_2000, data=df3,
            family=binomial)
summary(fit1)
```

```{r}
odds1 <- exp(coef(fit1))
odds1
```

```{r}
m1 <- margins(fit1)
summary(m1)
```

```{r}
# Cut makes a categorical variable out of a numerical variable using specified bins
df3$Super <- ifelse(df3$intensity_vmax * 1.852 > 185, 1, 0)
df3$Moderate <- ifelse(df3$intensity_vmax * 1.852 >= 88 &
                         df3$intensity_vmax * 1.852 < 185, 1, 0)
df3$Weak <- ifelse(df3$intensity_vmax * 1.852 >= 41 & 
                   df3$intensity_vmax * 1.852 < 88, 1, 0)
df3$HK_intensity <- cut(df3$intensity_vmax * 1.852 ,c(-1,41, 62, 87, 117, 149, 999))
table(df3$HK_intensity)
```

```{r}
fit2 <- glm(delayed ~ (typhoon_500 + typhoon_1000 + typhoon_2000) : 
              (Weak + Moderate + Super), data=df3,
            family=binomial)
tidy(fit2)
```

```{r}
m2 <- margins(fit2)
summary(m2) %>%
  html_df()
```

```{r}
margins(fit2, at = list(Weak = c(1)),
        variables = c("typhoon_500", "typhoon_1000", "typhoon_2000")) %>%
  summary() %>%
  html_df()
```

```{r}
margins(fit2, at = list(Moderate = c(1)),
        variables = c("typhoon_500", "typhoon_1000", "typhoon_2000")) %>%
  summary() %>%
  html_df()
```

```{r}
margins(fit2, at = list(Super = c(1)),
        variables = c("typhoon_500", "typhoon_1000", "typhoon_2000")) %>%
  summary() %>%
  html_df()
```

