---
title: "Melbourne Pedestrian Traffic"
output:
  xaringan::moon_reader:
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
```

### What do we want to do?

* Investigate how weather impacts foot traffic around Melbourne.

### How are we doing this?

* Collect and combine data from pedestrian sensors and weather stations
* Explore the data a little bit
* Build and evalute a Poisson regression model

We'll need a few packages to help us out
```{r, warning=FALSE, message=FALSE}
library(tidyverse)
library(ggmap)
library(gridExtra)
library(readr)
library(knitr)
library(broom)
library(rwalkr)
```
---

class:center

### A Map of Melbourne
```{r echo=FALSE, warning=FALSE, message=FALSE, fig.width=10}
stations <- read_table("ghcnd-stations.txt", 
  col_names=c("ID", "lat", "lon", "elev", "state", "name", 
              "v1", "v2", "v3"), skip=353, n_max=17081)

oz <- map_data("world", xlim=range(stations$lon),
               ylim=range(stations$lat))
ped_loc <- read_csv("Pedestrian_Sensor_Locations.csv")
melb_stns <- stations %>% filter(lon > min(ped_loc$Longitude), 
                                 lon < max(ped_loc$Longitude), 
                                 lat > min(ped_loc$Latitude),
                                 lat < max(ped_loc$Latitude))

melb_stns <- melb_stns %>% mutate(type = "weather") %>%
  select(lat, lon, type)
ped_loc <- ped_loc %>% mutate(type = "ped") %>%
  rename(lat = Latitude, lon = Longitude) %>%
  select(lat, lon, type)

ped_loc_weather_stn <- bind_rows(ped_loc, melb_stns)
ped_loc_weather_stn$type <- factor(ped_loc_weather_stn$type)  

qmplot(lon, lat, data=ped_loc_weather_stn, 
       colour=type, shape=type, size=I(3)) +
  scale_colour_brewer(palette="Dark2")
```

---
### A subset of the weather data is contained in 'melb_ghcn.csv'.

```{r, message=FALSE, warning=FALSE}
melb_ghcn <- read_csv("melb_ghcn.csv")
melb_ghcn %>% 
  select(stn_id, date, variable, value) %>%
  head(4)
```
What do the variable and value columns mean?

The data is manipulated, here's a few important lines:
```{r, eval=FALSE}
mutate(value = value/10,
       high_prcp = ifelse(PRCP>5, "rain", "none"), 
       high_tmp = ifelse(TMAX>33, "hot", "not"), 
       low_tmp = ifelse(TMIN<6, "cold", "not"))
melb_ghcn_wide$PRCP[is.na(melb_ghcn_wide$PRCP)] <- 0
```
---

### Temperature over Time

```{r fig.width=12, echo=FALSE, warning=FALSE, message=FALSE}
melb_ghcn <- melb_ghcn %>% 
  select(stn_id, date, variable, value) %>%
  filter(variable %in% c("TMAX", "TMIN", "PRCP")) %>%
  mutate(year=substr(date, 1, 4), month_num=substr(date, 5, 6),
         mday=substr(date, 7, 8)) %>%
  mutate(date=as.Date(paste(mday, month_num, year, sep="-"),
                      format="%d-%m-%Y")) %>%
  mutate(mday=as.numeric(mday), month_num=as.numeric(month_num), 
         year=as.numeric(year), value = value/10)
melb_ghcn_wide <- melb_ghcn %>% spread(variable, value)
melb_ghcn_wide$PRCP[is.na(melb_ghcn_wide$PRCP)] <- 0
melb_ghcn_wide <- melb_ghcn_wide %>%
  mutate(high_prcp = ifelse(PRCP>5, "rain", "none"), 
         high_tmp = ifelse(TMAX>33, "hot", "not"), 
         low_tmp = ifelse(TMIN<6, "cold", "not"))
ggplot(melb_ghcn, aes(x=date, y=value)) + geom_point() + facet_wrap(~variable)
```

---

### Choose a sensor for your group.
```{r, warning=FALSE, message=FALSE}
# Read sensor counts
ped_sub <- read_csv("pedestrian_counts_sub.csv")
ped_sub <- ped_sub %>% 
  filter(year < 2015) %>%
  dplyr::arrange(sensor_id, date, time)
unique(ped_sub$sensor_name)
```

---
### We have to make a few changes to the code

```{r, eval=FALSE}
ped_run <- run_melb(year = 2013:2014,
                    sensor = "Bourke Street Mall (North)")
ped_run <- ped_run %>%
  rename(count = Count, time = Time, date = Date) %>%
  mutate(day = wday(date, label=TRUE),
         month = month(date, label=TRUE))

ped_weather <- left_join(ped_run, melb_ghcn_wide, by="date")
ped_weather <- ped_weather %>% 
  mutate(time = factor(ped_weather$time), 
         high_prcp = factor(ped_weather$high_prcp,
                            levels=c("none", "rain")), 
         high_tmp = factor(ped_weather$high_tmp,
                           levels=c("not", "hot")),
         low_tmp = factor(ped_weather$low_tmp,
                          levels=c("not", "cold")))
```

---


### Poisson Regression

$$\log(\mu_i) = \beta_0 + \sum_{j=1}^K \beta_j x_{i, j}$$
$$y_i = \mu_i + e_i$$

The dependent variable has a *poisson* distribution with mean $\mu_i$.

We want three way interactions between day of the week, month of the year, and hour of the day:
$$\log(\mu_i) = \mbox{Standard Variables} + \beta_{44}(\mbox{Day = Monday AND Time = 00:00}) + \dots$$ 
$$+ \beta_{??} (\mbox{Day = Thursday AND Time = 14:00 AND Month = June}) + \dots$$

In R:
```{r, eval=FALSE}
glm(count~day*time*month+high_tmp+low_tmp+high_prcp,
    data=ped_weather, family=poisson(link="log"))
```
There are 2019 coefficients to estimate, be patient.

---
### Predictions for the Poisson Regression

Many of the variables are factors, make sure to convert your prediction set to the right format.

```{r, eval = FALSE}
newdat <- data.frame() #Something needs to go here 
newdat$time <- factor(newdat$time, levels=0:23)
newdat$high_tmp <- factor(newdat$high_tmp, levels=c("not", "hot"))
newdat$low_tmp <- factor(newdat$low_tmp, levels=c("not", "cold"))
newdat$high_prcp <- factor(newdat$high_prcp,
                           levels=c("none", "rain"))
newdat$day <- factor(newdat$day,
                levels=c("Mon", "Tues", "Wed", "Thurs", 
                         "Fri", "Sat", "Sun"))
newdat$month <- factor(newdat$month,
                levels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
                         "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
```