---
title:
output:
  revealjs::revealjs_presentation:
    theme: white
    css: white-edits.css
    highlight: tango
    transition: slide
    self_contained: false
    reveal_plugins: ["chalkboard","search"]
    transitionSpead: fast
    mathjax: "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_SVG"
    reveal_options:
      width: 1024
      height: 768
      minScale: 1.0
      maxScale: 1.0
      slideNumber: true
      controlsTutorial: true
      chalkboard:
        theme: whiteboard
        readOnly: false
---

## ACCT 420: Logistic Regression {data-background="../Backgrounds/title.jpg" class="header no-toc-progress" data-state="no-toc-progress"}

<!-- "reveal_plugins:..." and "chalkboard:" are handled by prep.py -->
```{r helpers, echo=FALSE, 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)
  }
}
```


<h3 class="date"><br>Session 4</h3>
<h2 class="author"><br>Dr. Richard M. Crowley</h2>

# Front matter {data-background="../Backgrounds/section.jpg" class="header"}

## Learning objectives {data-background="../Backgrounds/default.jpg" class="default"}

<div class="clearfix">
<div class="left">
<img src="../Course_Diagrams/Map-4.png">
</div>
<div class="right objectives">
- <b>Theory:</b>
    - Understanding binary problems
- <b>Application:</b>
    - Detecting shipping delays caused by typhoons
- <b>Methodology:</b>
    - Logistic regression
    - Spatial visualization
</div>
<div class="clearfix">

## Datacamp {data-background="../Backgrounds/default.jpg" class="default"}

> - Explore on your own
> - No specific required class this week

## Assignment 2 {data-background="../Backgrounds/reference.jpg" class="default"}

- Looking at Singaporean retail firms
    - Mostly focused on time and cyclicality
    - Some visualization
    - A little of what we cover today
- Optional (but encouraged):
    - You can work in *pairs* on this assignment
        - If you choose to do this, please only make 1 submission and include both your names on the submission

# Binary outcomes {data-background="../Backgrounds/section.jpg" class="header"}

## What are binary outcomes? {data-background="../Backgrounds/default.jpg" class="default"}

- Thus far we have talked about events with continuous outcomes
    - Revenue: Some positive number
    - Earnings: Some number
    - ROA: Some percentage
- Binary outcomes only have two possible outcomes
    - Did something happen, *yes* or *no*?
    - Is a statement *true* or *false*?

## Accounting examples of binary outcomes {data-background="../Backgrounds/default.jpg" class="default"}

- Financial accounting:
    - Will the company's earnings meet analysts' expectations?
    - Will the company have positive earnings?
- Managerial accounting:
    - Will we have ___ problem with our supply chain?
    - Will our customer go bankrupt?
- Audit:
    - Is the company committing fraud?
- Taxation:
    - Is the company too aggressive in their tax positions?

> We can assign a probability to any of these

## Regression approach: Logistic regression {data-background="../Backgrounds/reference.jpg" class="default"}

- When modeling a binary outcome, we use logistic regression
    - A.k.a. logit model
- The *logit* function is $logit(x) = \text{log}\left(\frac{x}{1-x}\right)$
    - Also called *log odds*

$$
\text{log}\left(\frac{\text{Prob}(y=1|X)}{1-\text{Prob}(y=1|X)}\right)=\alpha + \beta_1 x_1 + \beta_2 x_2 + \ldots + \varepsilon
$$

> There are other ways to model this though, such as [probit](https://en.wikipedia.org/wiki/Probit_model)


## Implementation: Logistic regression {data-background="../Backgrounds/reference.jpg" class="default"}

- The logistic model is related to our previous linear models as such:

> Both linear and logit models are under the class of General Linear Models (GLMs)

- To regress a GLM, we use the `glm()` command.
    - In fact, the `lm()` command we have been using is actually `glm()` when you specify the option `family=gaussian`
- To run a logit regression:

```{r, eval=F}
mod <- glm(y ~ x1 + x2 + x3 + ..., data=df, family=binomial)

summary(mod)
```

> `family=binomial` is what sets the model to be a logit

## Interpreting logit values {data-background="../Backgrounds/reference.jpg" class="default"}

- The **sign** of the coefficients means the same as before
    - **+**: *increases* the likelihood of $y$ occurring
    - **-**: *decreases* the likelihood of $y$ occurring
- The level of a coefficient is different
    - The relationship isn't linear between $x_i$ and $y$ now
    - Instead, coefficients are in log odds
        - Thus, $e^{\beta_i}$ gives you the *odds*, $o$
- You can interpret the odds for a coefficient
    - Increased by $[o-1]\%$
- You need to sum all relevant log odds before converting to a probability!

## Odds vs probability {data-background="../Backgrounds/reference.jpg" class="default"}

<img src="../Figures/odds.png" height="600px">

## Example logit regression {data-background="../Backgrounds/case.jpg" class="default"}

> Do holidays increase the likelihood that a department more than doubles its store's average weekly sales across departments?

```{r, include=F, echo=F, 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)
```

> Holidays increase the odds...  but by how much?

# Logistic regression interpretation {data-background="../Backgrounds/section.jpg" class="header"}

## A simple interpretation {data-background="../Backgrounds/default.jpg" class="default"}

- The model we just saw the following model:

$$
logodds(Double~sales) = -3.45 + 0.54 IsHoliday
$$

- There are two ways to interpret this:
    1. Coefficient by coefficient
    2. In total

## Interpretting specific coefficients {data-background="../Backgrounds/default.jpg" class="default"}

$$
logodds(Double~sales) = -3.45 + 0.54 IsHoliday
$$

- Interpreting specific coefficients is easiest done manually
- Odds for the $IsHoliday$ coefficient are `exp(0.54) = 1.72`
    - This means that having a holiday modifies the baseline (i.e., non-Holiday) odds by 1.72 to 1
        - Where 1 to 1 is considered no change
    - Baseline is 0.032 to 1

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

## Interpretting in total {data-background="../Backgrounds/default.jpg" class="default"}

- It is important to note that log odds are additive
    - So, calculate a new log odd by plugging in values for variables and adding it all up
        - Holiday: $-3.45 + 0.54 * 1 = -2.89$
        - No holiday: $-3.45 + 0.54 * 0 = -3.45$
- Then calculate odds and log odds like before
    - With holiday: `exp(-2.89) = 0.056`
    - Without holiday: `exp(-3.45) = 0.032`
    - Ratio of holiday to without: 1.72!
        - This is the individual log odds for holiday

> We need to specify values to calculate log odds in total

## Converting to probabilities {data-background="../Backgrounds/default.jpg" class="default"}

- We can calculate a probability at any given point using the log odds

$$
Probability = \frac{odds}{odds + 1}
$$

- Probability of double sales...
    - With a holiday: `0.056 / (0.056 + 1) = 0.052`
    - Without a holiday: `0.032 / (0.032 + 1) = 0.031`

> These are easier to interpret, but require specifying values for each model input to calculate

## Using predict() to simplify it {data-background="../Backgrounds/default.jpg" class="default"}

- `predict()` can calculate log odds and probabilities for us with minimal effort
    - Specify `type="response"` to get probabilities

```{r, eval=F}
test_data <- as.data.frame(IsHoliday = c(0,1))
predict(model, test_data)  # log odds
```
```{r, echo=F}
c(-3.44, -2.9)
```
```{r, eval=F}
predict(model, test_data, type="response")  #probabilities
```
```{r, echo=F}
a <- exp(-3.44)
b <- exp(-2.9)
c(a/(1+a), b/(1+b))
```

- Here, we see the baseline probability is 3.1%
- The probability of doubling sales on a holiday is higher, at 5.2%

## R practice: Logit {data-background="../Backgrounds/coding.jpg" class="default"}

- A continuation of last week's practices answering:
    - Is Walmart more likely to see a year over year decrease in quarterly revenue during a recession?
- Practice using `mutate()` and `glm()`
- Do exercises 1 and 2 in today's practice file
    - <a href="http://rmc.link/Slides/acct420v3/Session_4/Session_4_R.html">R Practice</a>
    - Shortlink: <a href="http://rmc.link/Slides/acct420v3/Session_4/Session_4_R.html">rmc.link/420r4</a>

# Logistic regression interpretation redux {data-background="../Backgrounds/section.jpg" class="header"}

## What about more complex models? {data-background="../Backgrounds/default.jpg" class="default"}

- Continuous inputs in the model
    - What values do we pick to determine probabilities?
- Multiple inputs?
    - We can scale up what we did, but things get messy
        - Mathematically, the inputs get interacted within the inner workings of logit...
            - So the impact of each input depends on the values of the others!

## Consider this model {data-background="../Backgrounds/default.jpg" class="default"}

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

## Odds and probabilities {data-background="../Backgrounds/default.jpg" class="default"}

```{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)))
```

## A bit easier: Marginal effects {data-background="../Backgrounds/default.jpg" class="default"}

> Marginal effects tell us the *average* change in our output for a change of 1 to an input

- The above definition is very similar to how we interpret linear regression coefficients
    - The only difference is the word *average* -- the effect changes a bit depending on the input data
- Using `package:margins`, we can calculate marginal effects
- There are a few types that we could calculate:
    - An *Average Marginal Effect* tells us what the average effect of an input is across all values in our data
        - This is the default method in the package
    - We can also specify a specific value to calculate marginal effects at (like with our probabilities last slides)

## Marginal effects in action {data-background="../Backgrounds/default.jpg" class="default"}

```{r, echo=F, include=F, 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
```

- A holiday increase the probability of doubling by a flat 1.33%
    - Not too bad when you consider that the probability of doubling is `r percent(mean(df$double))`
- If the temperature goes up by 1°F (0.55°C), the probability of doubling changes by `r extract_margin(m,Temperature)`
- If the fuel price increases by 1 USD for 1 gallon of gas, the probability of doubling changes by `r extract_margin(m,Fuel_Price)`

## `package:margins` niceties {data-background="../Backgrounds/default.jpg" class="default"}

- We can get some extra information about our marginal effects through `summary()`:

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

- Those p-values work just like with our linear models
- We also get a confidence interval
    - Which we can plot!

## Plotting marginal effects {data-background="../Backgrounds/default.jpg" class="default"}

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

> Note: The `which...` part is absolutely necessary at the moment due to a bug in the package

## Marginal effects at a specified value {data-background="../Backgrounds/default.jpg" class="default"}

```{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, 20, 40, 60, 80, 100)),
        variables = c("IsHoliday")) %>%
  summary() %>%
  html_df()
```

# Today's Application: Shipping delays {data-background="../Backgrounds/section.jpg" class="header"}

## The question {data-background="../Backgrounds/default.jpg" class="default"}

> Can we leverage global weather data to predict shipping delays?

<img src="../Figures/Singapore_port.jpg" height="400px">

## Formalization {data-background="../Backgrounds/default.jpg" class="default"}

1. Question
    - How can predict naval shipping delays?
2. Hypothesis (just the alternative ones)
    1. Global weather data helps to predict shipping delays
3. Prediction
    - Use Logistic regression and $z$-tests for coefficients
    - No hold out sample this week -- too little data

## A bit about shipping data {data-background="../Backgrounds/default.jpg" class="default"}

- WRDS doesn't have shipping data
- There are, however, vendors for shipping data, such as:

<a href="https://www.marinetraffic.com/en/ais/details/ports/290/Singapore_port:SINGAPORE"><img src="../Figures/Marine_traffic.png" height="120px"></a>

- They pretty much have any data you could need:
    - Over 650,000 ships tracked using ground and satellite based AIS 
        - AIS: Automatic Identification System
    - Live mapping
    - Weather data
    - Fleet tracking
    - Port congestion
    - Inmarsat support for ship operators

## What can we see from naval data? {data-background="../Backgrounds/case_digital.jpg" class="default"}

<div class="clearfix">
<div class="left">
Yachts in the Mediterranean

<img src="../Figures/Ships_med_yachts.png" width="400px">
</div>
<div class="right">
Oil tankers in the Persian gulf

<img src="../Figures/Ships_persiangulf_oil.png" width="400px">
</div>
</div>

## What can we see from naval data? {data-background="../Backgrounds/case_digital.jpg" class="default"}

<div class="clearfix">
<div class="left">
Shipping route via the Panama canal

<img src="../Figures/Ships_panama_canal.png" width="400px">
</div>
<div class="right">
River shipping on the Mississippi river, USA

<img src="../Figures/Ships_USA_river.png" width="250px">
</div>
</div>

## What can we see from naval data? {data-background="../Backgrounds/case_digital.jpg" class="default"}

<div class="clearfix">
<div class="left">
Busiest ports by containers and tons (Shanghai & Ningbo-Zhoushan, China)
<img src="../Figures/Ships_Ningbo-Zhoushan_busy-tonnage.png" width="400px">
</div>
<div class="right">
Busiest port for transshipment (Singapore)

<img src="../Figures/Ships_Singapore_busy-trans.png" width="400px">
</div>
</div>

## Examining Singaporean owned ships {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, include=F}

# 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, echo=F, 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
```

## Code for last slide's map {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, eval=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
```

- `plot_geo()` is from `package:plotly`
- `add_markers()` adds points to the map
- `layout()` adjusts the layout
- Within geo, a list, the following makes the map a globe
    - `projection=list(type="orthographic")`

## Singaporean ship movement {data-background="../Backgrounds/case_digital.jpg" class="default"}

<!-- Disabled due to slow down.  Output to a different file -->
```{r animated_ships, eval=F, echo=F, warning=F, message=F, fig.height=8, fig.width=8}
library(sf)  # Note: very difficult to install except on Windows
library(maps)
# Requires separately installing "maptools" and "rgeos" as well
# This graph requires ~7GB of RAM to render
world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))

df_all2 <- df_all[df_all$frame > 1 & df_all$frame < 30,]
df_all2 <- df_all2 %>% arrange(frame, imo)

p <- ggplot(data = world1) + 
  geom_sf() + 
  geom_point(data = df_all2, aes(x = lon, y = lat, frame=frame,
                                text=paste("name:",shipname)))
ggplotly(p) %>%
   animation_opts(
    1000, easing = "linear", redraw = FALSE)
```

<a href="http://rmc.link/Slides/acct420v3/Session_4/SG_ships_animation.html">Link to ship movement animation</a>

## Code for last slide's map {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, eval=F, fig.height=8, fig.width=8}
library(sf)  # Note: very difficult to install except on Windows
library(maps)
# Requires separately installing "maptools" and "rgeos" as well
# This graph requires ~7GB of RAM to render
world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))

df_all <- df_all %>% arrange(run, imo)

p <- ggplot(data = world1) + 
  geom_sf() + 
  geom_point(data = df_all, aes(x = lon, y = lat, frame=frame,
                                text=paste("name:",shipname)))
ggplotly(p) %>%
   animation_opts(
    1000, easing = "linear", redraw = FALSE)
```

- `world1` contains the map data
- `geom_sf()` plots map data passed to `ggplot()`
- `geom_point()` plots ship locations as longitude and latitude
- `ggplotly()` converts the graph to html and animates it
    - Animation follows the `frame` aesthetic

## What might matter for shipping? {data-background="../Backgrounds/group.jpg" class="default"}

> What observable events or data might provide insight as to whether a naval shipment will be delayed or not?

## {data-background="../Backgrounds/group.jpg" class="default"}

<iframe allowfullscreen frameborder="0" mozallowfullscreen src="https://app.wooclap.com/events/GBJOST/questions/5f530b02008031118f39b6d1" width="100%" height="600px"></iframe>

## Typhoon Jebi {data-background="../Backgrounds/answer.jpg" class="default"}

<iframe width="560" height="315" src="https://www.youtube.com/embed/kH9o1CSs5Rk?rel=0" frameborder="0" allow="autoplay; encrypted-media" allowfullscreen></iframe>

<br>

- <a href="https://www.youtube.com/embed/kH9o1CSs5Rk?rel=0">link</a>
- <a href="https://earth.nullschool.net/#2018/09/03/2100Z/wind/surface/level/orthographic=-221.34,35.55,806/loc=133.763,27.561">Nullschool plot</a>

## Typhoons in the data {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, echo=F, 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
```

## Code for last slide's map {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, eval=F, 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
```

- This map is made the same way as the first map

## Typhoons in the data using leaflet {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, echo=F, 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)
```

## Code for last slide's map {data-background="../Backgrounds/coding.jpg" class="default"}

```{r, eval=F, warning=F, message=F, fig.height=5, fig.width=8}
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, 5,
    ifelse(typhoon_Jebi$intensity_vmax < 118/1.852, 2, 3)))

# 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=3, 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=2, color="red", label=~date) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax <=118/1.852,],
    lng=~lon, lat=~lat, stroke = TRUE, radius=1, color="red", label=~date) %>%
  addMarkers(data=df_all[df_all$frame == 14,], lng=~lon, lat=~lat,
             label=~shipname, icon=shipicons)
```

- `pulseIcons()`: pulsing icons from `package:leaflet.extras`
- `iconList()`: pulls icons stored on your computer
- `leaflet()`: start the map; `addTiles()` pulls from <a href="https://www.openstreetmap.org/">OpenStreetMap</a>
- `setView()`: sets the frame for the map
- `addPulseMarkers()`: adds pulsing markers
- `addCircleMarkers()`: adds circular markers

## R Practice on mapping {data-background="../Backgrounds/coding.jpg" class="default"}

- Practice mapping typhoon data
    - 1 map using `package:plotly`
    - 1 map using `package:leaflet`
- Practice using `package:plotly` and `package:leaflet`
    - No practice using `package:ggplot2` as `package:sf` is missing on DataCamp light
        - And `package:sf` can be tough to install for anyone on a Mac
- Do exercises 3 and 4 in today's practice file
    - <a href="http://rmc.link/Slides/acct420v3/Session_4/Session_4_R.html">R Practice</a>
    - Shortlink: <a href="http://rmc.link/Slides/acct420v3/Session_4/Session_4_R.html">rmc.link/420r4</a>

# Predicting delays due to typhoons {data-background="../Backgrounds/section.jpg" class="header"}

## Data {data-background="../Backgrounds/case_digital.jpg" class="default"}

```{r, include=F}
# 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()
```

- If the ship will report a delay of at least 3 hours some time in the next 12-24 hours
- What we have:
    - Ship location
    - Typhoon location
    - Typhoon wind speed

> We need to calculate distance between ships and typhoons

## Distance for geo {data-background="../Backgrounds/case.jpg" class="default"}

- There are a number of formulas for this
    - *Haversine* for a simple calculation
    - *Vincenty's formulae* for a complex, incredibly accurate calculation
        - Accurate within **0.5mm**
- Use `distVincentyEllipsoid()` from `package:geosphere` to get a reasonably quick and accurate calculation
    - Calculates distance between two sets of points, x and y, structured as matrices
    - Matrices must have longitude in the first column and latitude in the second column
    - Provides distance in meters by default

```{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
```

## Clean up {data-background="../Backgrounds/case.jpg" class="default"}

- Some indicators to cleanly capture how far away the typhoon is

```{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)
```

<img src="../Figures/ifelse.png">

## Do typhoons delay shipments? {data-background="../Backgrounds/case.jpg" class="default"}

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

> It appears so!

## Interpretation of coefficients {data-background="../Backgrounds/case.jpg" class="default"}

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

- Ships 1,000 to 2,000 km from a typhoon have a 17% increased odds of having a delay

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

- Ships 1,000 to 2,000 km from a typhoon have an extra `r extract_margin(m1, typhoon_2000)` chance of having a delay (baseline of `r percent(mean(df3$delayed, na.rm=T))`)

## What about typhoon intensity? {data-background="../Backgrounds/default.jpg" class="default"}

- Hong Kong's typhoon classification: <a href="https://www.hko.gov.hk/informtc/class.htm">Official source</a>
    1. 41-62 km/h: Tropical depression
    2. 63-87 km/h: Tropical storm
    3. 88-117 km/h: Severe tropical storm
    4. 118-149 km/h: **Typhoon**
    5. 150-184 km/h: **Severe typhoon**
    6. 185+km/h: **Super typhoon**

```{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)
```

## Typhoon intensity and delays {data-background="../Backgrounds/case.jpg" class="default"}

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

> Moderate storms predict delays when within 500km

> Super typhoons predict delays when 1,000 to 2,000km away 

## Interpretation of coefficients {data-background="../Backgrounds/case.jpg" class="default"}

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

- Delays appear to be driven mostly by 2 factors:
    1. A typhoon 1,000 to 2,000 km away from the ship
    2. Weak typhoons

## Interpretating interactions {data-background="../Backgrounds/case.jpg" class="default"}

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

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

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

## What might matter for shipping? {data-background="../Backgrounds/group.jpg" class="default"}

> What other observable events or data might provide insight as to whether a naval shipment will be delayed or not?

- What is the reason that this event or data would be useful in predicting delays?
    - I.e., how does it fit into your mental model?

# End matter {data-background="../Backgrounds/ending.jpg" class="header"}

## For next week {data-background="../Backgrounds/ending.jpg" class="default"}

- For next week:
    - Second individual assignment
        - Finish by 2 classes from now 
        - Submit on eLearn
    - Think about who you want to work with for the project
- Survey on the class session at <a href="https://docs.google.com/forms/d/e/1FAIpQLSdMghdqvaZ9MiaFKsvV-iTE1WVwb4feghgZ7HjE2UDZKBDUVQ/viewform?usp=sf_link">rmc.link/420survey4</a>

## Packages used for these slides {data-background="../Backgrounds/coding.jpg" class="default"}

<div class="clearfix">
<div class="left">
- `package:broom`
- `package:geosphere`
- `package:kableExtra`
- `package:knitr`
- `package:leaflet`
- `package:leaflet.extras`
- `package:lubridate`
- `package:magrittr`
</div>
<div class="right">
- `package:margins`
- `package:maps`
- `package:maptools`
- `package:plotly`
- `package:revealjs`
- `package:rgeos`
- `package:sf`
- `package:tidyverse`
</div>
</div>

## Custom code {data-background="../Backgrounds/coding.jpg" class="default"}

<small style="width:100%">
```{r, eval=F}
# styling for plotly maps
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
  )
)
```

</small>

## NOTES {class="never"}

Add in some questions...
