---
title: "ACCT 420: Logistic Regression"
author: 
  - name: Dr. Richard M. Crowley
    email: rcrowley@smu.edu.sg
    affiliation: https://rmc.link/
title-slide-attributes:
  data-background-video: "../Backgrounds/title.mp4"
  data-background-video-loop: "true"
  data-background-size: stretch
format:
  revealjs:
    theme: default
    css: quarto_white.css
    navigation-mode: vertical
    controls: true
    transition: slide
    background-transition: fade
    transition-speed: default
    width: 1536
    height: 864
    logo: ../Logos/logo2.png
    code-link: true
    code-line-numbers: false
    chalkboard:
      theme: whiteboard
    auto-stretch: false
revealjs-plugins:
  - attribution
execute:
  echo: true
  warning: false
  message: false
editor: source
---

# Front Matter {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

<!--# Note that the downlit and xml2 pacakges are required to compile -->
<!--# Installed extensions include `attribution`, `qrcode`, and `material-icons` -->

```{r helpers, echo=FALSE, warning=FALSE, message=F}
library(knitr)
library(kableExtra)

#' Pretty formatted HTML tables
#' 
#' @param df The dataframe to format as a table
#' @param names Names to use for columns if not using the df column names. Vector
#' @param highlight_cols Columns to highlight. Default is c(1) for highlighting the first row. Use c() to remove.
#' @param highlight_rows Rows to highlight. Default is c() to not highlight any rows.
#' @param color_rows Rows to highlight in color. Default is c() to not highlight any rows.
#' @param highlight_end_rows Rows to highlight. Uses a reverse indexing where 1 is the last row, 2 is second to last, etc.
#' @param full Use to span the full width of a slide. FALSE by default.
#' @param fixed_header Whether the first row should be treated as a header row. TRUE by default.
html_df <- function(df, names=NULL, highlight_cols=c(1), highlight_rows=c(), color_rows=c(), highlight_end_rows=c(), full=F, fixed_header=T) {
  if(!length(names)) {
    names=colnames(df)
  }
  kbl(df,"html", col.names = names, align = c("l",rep('c',length(cols)-1))) %>%
    kable_paper(c("striped","hover"), full_width=full, fixed_thead=F, html_font = "\"Source Sans Pro\", Helvetica, sans-serif") %>%
    {if(length(highlight_cols)) column_spec(., highlight_cols,bold=T, background = "#ffffff99") else .} %>%
    {if(length(highlight_rows)) row_spec(., highlight_rows,bold=T, background = "#ffffff99") else .} %>%
    {if(length(color_rows)) row_spec(., color_rows,bold=T, background = "#aaaaff99") else .} %>%
    {if(length(highlight_end_rows)) row_spec(., nrow(df) + 1 - highlight_end_rows,bold=T, background = "#ffffff99") else .}
}

library(tidyverse)
library(broom)
```

<!--# List packages here that are not being correctly highlighted -->
<!--# The usual cause is calling the library from a echo=FALSE cell -->

```{r}
#| include: false
options(downlit.attached = c("knitr", "kableExtra", "tidyverse"))
```

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

:::: columns

::: {.column width="50%"}

![](../Course_Diagrams/Map-4.png){fig-alt="Roadmap for the course"}

:::

::: {.column width="50%"}

- <b>Theory:</b>
    - Understanding binary problems
- <b>Application:</b>
    - Detecting shipping delays caused by typhoons
- <b>Methodology:</b>
    - Logistic regression
    - Spatial visualization

:::

::::

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

- 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-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

- 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

## Brainstorming... {data-background="../Backgrounds/rain.jpg" class="default present-only"}

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

## Brainstorming... {data-background="../Backgrounds/group.jpg" class="default present-not" visibility="hidden"}

> What types of business problems or outcomes are binary?

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

- 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"}

- 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"}

- 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"}

![](../Figures/odds.png){height="600px" fig-align="center"}

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

> 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 <- df %>%
    inner_join(weekly.stores) %>%
    inner_join(weekly.features[,1:11])
  
  # Compress the weird markdown information to 1 variable
  df <- df %>%
    mutate(markdown = 0,
           markdown = ifelse(!is.na(MarkDown1),MarkDown1,markdown),
           markdown = ifelse(!is.na(MarkDown1),MarkDown2,markdown),
           markdown = ifelse(!is.na(MarkDown1),MarkDown3,markdown),
           markdown = ifelse(!is.na(MarkDown1),MarkDown4,markdown),
           markdown = ifelse(!is.na(MarkDown1),MarkDown5,markdown))
  
  # Fix dates and add useful time variables
  df <- df %>%
    mutate(date = as.Date(Date),
           week = week(date),
           year = year(date),
           month = month(date))

  df
}

weekly <- read_csv("../../Data/WMT_train.csv")
weekly.features <- read_csv("../../Data/WMT_features.csv")
weekly.stores <- read_csv("../../Data/WMT_stores.csv")
df <- preprocess_data(weekly)

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 <- df %>% mutate(double = ifelse(Weekly_Sales > 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-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

$$
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"}

- 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"}

- 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"}

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

:::: {.columns}

::: {.column width="50%"}
**Outcomes**

```{r}
test_data <- data.frame(IsHoliday = c(F,T))
predict(fit, test_data)  # log odds
```
```{r}
# probabilities
predict(fit, test_data, type="response")
```
:::

::: {.column width="50%"}
**Converting from log odds to probabilities**

```{r, echo=T}
a <- exp(-3.446804)
b <- exp(-2.908164)
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"}

- 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/420r4">R Practice</a>
    - Shortlink: <a href="http://rmc.link/420r4">rmc.link/420r4</a>

# Logistic regression interpretation redux {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

- 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"}

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

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

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

```{r}
# Typical September days
hday_sep    <- mean(predict(fit2, filter(df, IsHoliday, month==9), type="response"))
no_hday_sep <-  mean(predict(fit2, filter(df, !IsHoliday, month==9), type="response"))
# Typical December days
hday_dec    <- mean(predict(fit2, filter(df, IsHoliday, month==12), type="response"))
no_hday_dec <-  mean(predict(fit2, 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"}

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

```{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 arguments 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(fit2)
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)`

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

- 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"}

```{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" .small}

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

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

# Today's Application: Shipping delays {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

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

![](../Figures/Singapore_port.jpg){height="400px" fig-align="center"}

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

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"}

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

[![](../Figures/Marine_traffic.png){height="120px" fig-align="center"}](https://www.marinetraffic.com/en/ais/details/ports/290/Singapore_port:SINGAPORE)

- 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"}

:::: {.columns}

::: {.column width="50%"}
Yachts in the Mediterranean

![](../Figures/Ships_med_yachts.png){width="600px" fig-align="center"}

:::

::: {.column width="50%"}
Oil tankers in the Persian gulf

![](../Figures/Ships_persiangulf_oil.png){width="600px" fig-align="center"}
:::

::::



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

:::: {.columns}

::: {.column width="50%"}
Shipping route via the Panama canal

![](../Figures/Ships_panama_canal.png){width="600px" fig-align="center"}

:::

::: {.column width="50%"}
River shipping on the Mississippi river, USA

![](../Figures/Ships_USA_river.png){height="500px" fig-align="center"}
:::

::::


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

:::: {.columns}

::: {.column width="50%"}
Busiest ports by containers and tons (Shanghai & Ningbo-Zhoushan, China)

![](../Figures/Ships_Ningbo-Zhoushan_busy-tonnage.png){width="500px" fig-align="center"}

:::

::: {.column width="50%"}
Busiest port for transshipment (Singapore)

![](../Figures/Ships_Singapore_busy-trans.png){width="600px" fig-align="center"}
:::

::::


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

```{r, include=F}

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

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

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

# fix dates
df <- df %>% 
  mutate(last_update = as.POSIXct(last_update, tz="UTC",origin="1970-01-01"),
         left_time = as.POSIXct(left_time, tz="UTC",origin="1970-01-01"),
         arrive_time = as.POSIXct(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=6, fig.width=12}
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"}

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

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

[Link to ship movement animation](http://rmc.link/Slides/acct420v5/Session_4/SG_ships_animation.html)

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

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

## Brainstorming... {data-background="../Backgrounds/rain.jpg" class="default present-only"}

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

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

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

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

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

<br>

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

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

```{r, echo=F, warning=F, message=F, fig.height=6, fig.width=10}
# 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"}

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

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

# typhoon icons
ty_icons <- pulseIcons(color='red',  #  supports hex colors like "#C0C0C0"
  heartbeat = 1.2,  #heartbeat is frequency of pulsing
  iconSize=3)

# ship icons
shipicons <- iconList(
  ship = makeIcon("../Figures/ship.png", NULL, 18, 18)
)

leaflet() %>%
  addTiles() %>% 
  setView(lng = 136, lat = 20, zoom=3) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax <= 150/1.852,], lng=~lon, lat=~lat,  # Not super typhoon
                   stroke = TRUE, radius=2, color="orange", label=~date) %>%
  addPulseMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax > 150/1.852,], lng=~lon, lat=~lat,  # Super typhoon
                   icon=ty_icons, 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"}

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

# typhoon icons
ty_icons <- pulseIcons(color='red',  #  supports hex colors like "#C0C0C0"
  heartbeat = 1.2,  #heartbeat is frequency of pulsing
  iconSize=3)

# ship icons
shipicons <- iconList(
  ship = makeIcon("../Figures/ship.png", NULL, 18, 18)
)

leaflet() %>%
  addTiles() %>% 
  setView(lng = 136, lat = 20, zoom=3) %>%
  addCircleMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax <= 150/1.852,], lng=~lon, lat=~lat,  # Not super typhoon
                   stroke = TRUE, radius=2, color="orange", label=~date) %>%
  addPulseMarkers(data=typhoon_Jebi[typhoon_Jebi$intensity_vmax > 150/1.852,], lng=~lon, lat=~lat,  # Super typhoon
                   icon=ty_icons, label=~date) %>%
  addMarkers(data=df_all[df_all$frame == 14,], lng=~lon, lat=~lat,
             label=~shipname, icon=shipicons)
```

- `pulseIcons()`: pulsing icons from `{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"}

- Practice mapping typhoon data
    - 1 map using `{plotly}`
    - 1 map using `{leaflet}`
- Practice using `{plotly}` and `{leaflet}`
- No practice using `{ggplot2}` with `{sf}`
    - The technical difficulties involved in setting up `{sf}`, along with the computational intensity, make it a bit tricky
- Do exercises 3 and 4 in today's practice file
    - <a href="http://rmc.link/420r4">R Practice</a>
    - Shortlink: <a href="http://rmc.link/420r4">rmc.link/420r4</a>

# Predicting delays due to typhoons {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

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

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

# Determine the main DV: delayed
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" .small}

- Use `st_distance()` from `{sf}` to get a quick and reasonably accurate calculation
    - Uses *great-circle distance*

```{r}
library(sf)

# sf can't handle NA values, so we need to filter them out
df3 <- df3 %>% mutate(n = row_number())
df_ty <- df3 %>% filter(!is.na(ty_lat), !is.na(ty_lon))

# Create lists of points on Earth
x <- st_as_sf(df_ty %>% select("lon", "lat"), coords = c("lon", "lat"), 
                 crs = 4326, agr = "constant")
y <- st_as_sf(df_ty %>% select("ty_lon", "ty_lat"), coords = c("ty_lon", "ty_lat"), 
                 crs = 4326, agr = "constant")

# Find the distance between the points
df_ty$dist_typhoon <- st_distance(x, y, by_element=T) / 1000  # distance in km

# Add the distance back to our main data
# Replaced NAs with a placeholder longer than any possible distance
df3 <- df3 %>%
  left_join(df_ty %>% select(n, dist_typhoon)) %>%
  mutate(dist_typhoon = ifelse(is.na(dist_typhoon), 20036, dist_typhoon)) %>%
  select(-n)
```

::: {.callout-caution}
## The `st_distance()` function is not 100% accurate

The *great-circle distance* function assumes the Earth is a sphere. As such, for some distances it is off by ~0.5%. A better approach is *Vincenty's Formula*, but the R package for calculating this is deprecated as of October 2023.

:::

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

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

![](../Figures/ifelse.png){width="600px" fig-align="center"}

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

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

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

- 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"}

```{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" .small}

```{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" .small}

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

> 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-video="../Backgrounds/section.mp4"}

## Wrap up {data-background="../Backgrounds/ending.jpg"}

- For next week:
    - Finish the first individual assignment
    - 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 this QR code:

{{< qrcode rmc.link/420survey4 qrsurvey width=200 height=200 >}}

## Packages used for these slides {data-background="../Backgrounds/code.jpg" visibility="uncounted"}

- `{broom}`
- `{DT}`
- `{downlit}`
- `{kableExtra}`
- `{knitr}`
- `{leaflet}`
- `{leaflet.extras}`
- `{margins}`
- `{maptools}`
- `{plotly}`
- `{quarto}`
- `{revealjs}`
- `{sf}`
- `{tidyverse}`

## Custom code {data-background="../Backgrounds/code.jpg" class="small" visibility="uncounted"}

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


## NOTES {class="never" visibility="hidden"}


