---
title: |
  Spring 2023 \
  GE 461 Introduction to Data Science
# title: |
pagetitle: GE 461 Introduction to Data Science
papersize: a4paper
author: Statistical Models by Savaş Dayanık
# author: Statistical Models by Savaş Dayanık
always_allow_html: true
linkcolor: red
output: 
  bookdown::html_document2:
    theme: readable
    number_sections: false
    code_folding: "hide"
    toc: true
  bookdown::pdf_document2:
    number_sections: false
bibliography: GE461.bib
link-citations: yes
---

----

<center> <h4>  Week 5: Advertising and Promotion </h2> </center>

----

```{r setup, include=FALSE}
library(magrittr)
library(tidyverse)
library(car)
library(knitr)
library(kableExtra)
library(pander)
opts_chunk$set(echo = TRUE)

options(knitr.kable.NA =".") 
kable_format <- if (is_html_output()) "html" else "latex"
options(scipen = 999)
```


The Dodgers is a professional baseball team and plays in the Major Baseball League. The team owns a 56,000-seat stadium and is interested in increasing the attendance of their fans during home games. *At the moment the team management would like to know if bobblehead promotions increase the attendance of the team's fans?* This is a case study based on @miller2014modeling[ Chapter 2].

```{r, out.width=c("33%","33%","13%"), fig.align='center', fig.show="hold", fig.cap="56,000-seat Dodgers stadium (left),   shirts and caps (middle),  bobblehead (right)"}
include_graphics(c("los_angeles-dodgers-stadium.jpg",
                 "Los-Angeles-Dodgers-Promo.jpg",
                 "adrian_bobble.jpg"))
```

    
The 2012 season data in the `events` table of SQLite database `data/dodgers.sqlite` contain for each of  81 home play the 

* month, 
* day, 
* weekday, 
* part of the day (day or night),
* attendance, 
* opponent, 
* temperature, 
* whether cap or shirt or bobblehead promotions were run, and 
* whether fireworks were present.

## Prerequisites

We will use `R`, `RStudio`, `R Markdown` for the next three weeks to fit statistical models to various data and analyze them. Read @wickham2017r online

*  [Section 1.1](https://r4ds.had.co.nz/introduction.html#prerequisites) for how to download and install `R` and `RStudio`,
* [Chapter 27](https://r4ds.had.co.nz/r-markdown.html) for how to use `R Markdown` to interact with `R` and conduct various predictive analyses.

All materials for the next three weeks will be available on [Google drive](https://drive.google.com/drive/folders/1fdrM55XeDEkUQlb6ASJl4fdx1ozQyHEx?usp=share_link).

## March 13: Exploratory data analysis

@. Connect to  `data/dodgers.sqlite`. Read table `events` into a variable in `R`.

    * Read @baumer2017modern[, Chapters 1, 4, 5, 15] ([Second edition online](https://mdsr-book.github.io/mdsr2e/)) for getting data from and writing them to various SQL databases. 
    
    * Because we do not want to hassle with user permissions, we will use SQLite for practice. I recommend `PostgreSQL` for real projects.
    
    * Open `RStudio` terminal, connect to database `dodgers.sqlite` with `sqlite3`. Explore it (there is only one table, `events`, at this time) with commands 
      - `.help`
      - `.databases `
      - `.tables`
      - `.schema <table_name>`
      - `.headers on`
      - `.mode column`
      - `SELECT ...`
      - `.quit`
      
    * Databases are great to store and retrieve large data, especially, when they are indexed with respect to variables/columns along with we do search and match extensively. 
    
    * `R` (likewise, `Python`) allows one to seeminglessly read from and write to databases. For fast analysis, keep data in a database, index tables for fast retrieval, use `R` or `Python` to fit models to data.

```{r}
library(RSQLite)  ## if package is not on the computer, then install it only once using Tools > Install packages...
con <- dbConnect(SQLite(), "../data/dodgers.sqlite") # read Modern Data Science with R for different ways to connect a database.

## dbListTables(con)

tbl(con, "events") %>% 
  collect() %>% 
  mutate(day_of_week = factor(day_of_week, levels = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")),
         month = factor(month, levels = c("APR","MAY","JUN","JUL","AUG","SEP","OCT"))) %>% 
  mutate_if(is.character, factor) %>% 
  mutate(temp = round((temp- 32)*5/9)) -> events
```

```{r}
events %>% 
  count(bobblehead, fireworks)
```

@. What are the number of plays on each week day and in each month of a year?

Table \@ref(tab:monthweekday) and \@ref(fig:barweekdaymonth) summarize the number of games played on each weekday and month.

```{r}
events %>% 
  count(day_of_week, month) %>% 
  pivot_wider(names_from = day_of_week, values_from = n) %>% 
  pander(caption = "(\\#tab:monthweekday) Number of games played in each weekday and month")
```

```{r barweekdaymonth, fig.cap = "Barplot of counts of games for each weekday and month"}
events %>% 
  ggplot(aes(day_of_week)) +
  geom_bar(aes(fill=month))
  
```
Figure \@ref(fig:heatmap) shows heatmap of total attendance versus weekday and month. The colors change from bright yellow to dark red as attendance increases. Default heatmap shuffles rows and columns so as to bring together weekdays and months with similar attendance. Here we see May, Aug, and Jul together within the months and Saturday, Friday, Sunday within the weekdays. Learn more about xtabs (cross-table) heatmap by typing `?xtabs` and `?heatmap` in the R console.

```{r heatmap, fig.cap = "Heatmap of attendance versus weekday and month."}
xtabs(attend ~ day_of_week + month, events) %>% 
  heatmap()
```

In Figure \@ref(fig:scatterattend), I have added one more aes (colour) to capture day_night information. To avoid overplotting, I replaced `geom_point()` with `geom_jitter()`. These plots were also illuminating. So let us thank your friend who suggested this one, too.

```{r scatterattend, out.width="50%", fig.show="hold", fig.cap="Average and total attendances on each weekday and month in each part of day"}

sum_attend <- events %>% 
  group_by(day_of_week, month, day_night) %>% 
  summarize(mean_attend = mean(attend),
            total_attend = sum(attend), .groups = "drop")

sum_attend %>% 
  ggplot(aes(day_of_week, month)) +
  geom_jitter(aes(size = mean_attend, col = day_night), width = .1, height = .1, alpha=0.7) +
  scale_size(labels = scales::comma) +
  labs(title = "Average attendance", size = "attendance", col = "part of day",
       x = "Weekday", y = "Month")

sum_attend %>% 
  ggplot(aes(day_of_week, month)) +
  geom_jitter(aes(size = total_attend, col = day_night), width = .1, height = .1, alpha=0.7) +
  labs(title = "Total attendance", size = "attendance", col = "part of day",
       x = "Weekday", y = "Month") +
  scale_size(labels = scales::comma) +
  guides(col = guide_legend(order = 1), 
         shape = guide_legend(order = 2))
```


@. Check the orders of the levels of the `day_of_week` and `month` factors. If necessary, put them in the logical order.

```{r}
levels(events$day_of_week)
levels(events$month)
```

@. How many times were bobblehead promotions run on each week day?

```{r}
events %>% 
  count(day_of_week, bobblehead) %>% 
  pivot_wider(names_from = bobblehead, values_from = n) %>% 
  replace_na(list(YES = 0)) %>% 
  mutate(Total = YES + NO) %>% 
  select(-NO) %>% 
  rename(Bobblehead = YES)
```

@. How did the attendance vary across week days? Draw boxplots. On which day of week was the attendance the highest on average?

```{r}
events %>% 
  ggplot(aes(day_of_week, attend)) +
  geom_boxplot()
```

```{r}
events %>% 
  slice_max(order_by = attend, n=5)
```


@. Is there an association between attendance and 
    - whether the game is played in day light or night? 
    - whether skies are clear or cloudy?  

```{r}
events %>% 
  ggplot(aes(day_night, attend)) +
  geom_boxplot()
```

```{r}
t.test(x=events$attend[events$day_night=="Day"],
       y=events$attend[events$day_night=="Night"])
```

Since p-value (0.67) is large (greater than 0.05), we cannot reject null, which means there is no statistical difference between average attendance of games played in day and night.

```{r}
events %>% 
  ggplot(aes(skies, attend)) +
  geom_boxplot()
```

```{r}
t.test(x=events$attend[events$skies=="Clear"],
       y=events$attend[events$skies=="Cloudy"])
```

We do not see a statisticall significant difference between the average attendance of the games played under clear and cloudy skies.
  
@. Is there an association between attendance and temperature? 
    - If yes, is there a positive or negative association? 
    - Do the associations differ on clear and cloud days or day or night times? 

```{r, fig.width = 10}
events %>% 
  ggplot(aes(temp, attend)) +
  geom_jitter(aes(col = day_of_week)) +
  geom_text(aes(label = str_sub(day_of_week, 1,3), col = day_of_week)) +
  geom_smooth(se = FALSE)
```

```{r, fig.width = 10}
events %>% 
  ggplot(aes(temp, attend)) +
  geom_jitter(aes(col = month)) +
  geom_text(aes(label = str_sub(month, 1,3), col = month)) +
  geom_smooth(se = FALSE)
```

```{r}
events %>% 
  ggplot(aes(month, temp)) +
  geom_boxplot(aes(fill = day_night))
```

```{r, fig.width = 10}
events %>% 
  ggplot(aes(temp, attend)) +
  geom_jitter(aes(col = bobblehead)) +
  geom_smooth(se = FALSE)
```


```{r, fig.width = 10}
events %>% 
  ggplot(aes(temp, attend)) +
  geom_jitter(aes(col = opponent)) +
  geom_text(aes(label = str_sub(opponent, 1,4), col = opponent)) +
  geom_smooth(se = FALSE)
# +
#   geom_smooth(se=FALSE, method = "lm", formula = y ~ x, col = "red", lty = "dashed") +
#   geom_smooth(se=FALSE, method = "lm", formula = y ~ x + pmax(x-24,0), col = "red")
```

$$
attend = \beta_0 + \beta_1 temp + \beta_2 (temp - 23)^+ + \varepsilon_i
$$

```{r}
lmod0 <- lm(attend ~ temp , data = events) 
lmod0 %>% summary()
```
Null hypothesis tested with F-stat at the very bottom of the outplot of lm(): **"Null: There is no relation between attendance and temperature"** We say the data are consistent with the null hypothesis if p-value is large [always between 0 and 1]. In this case, pvlue is 0.38, which is large [small means < 0.05 in which case we reject the null]. Since p-value is large in this case, We conclude that the null is consistent with the data. We arrived at wrpng conclusion because the model was not checked through. Use diagnostics to make sure the model captures the features of data 

```{r}
plot(lmod0)
```



```{r}
lmod0 %>% summary()
lmod01 <- lm(attend ~ temp + pmax(0, temp - 24), data = events)
lmod01 %>% summary()
```
Here F-stat p-value is small (small if < 0.05). Therefore we reject the null (**null is equivalent to saying no relation between attendance and predictors in the right**). So find that temp is a useful predictor. 

On each row of table in the output, we test "null: corresponding variable has zero effect on attendance." The null is consistent with data is p-value on the last column is large (> 0.05). In this case both linear and broken line have nozero effects.


```{r}
plot(lmod0, which = 1)
plot(lmod01, which = 1)
```


**How does the break function help capture the increase/decrease in attendance with temperature?**

```{r, fig.asp=.8}
identity <- function(x) x
breakfun <-  function(x) pmax(x-24,0)
curve(identity(x), xlim = c(13,35), xlab = "temp", ylab = "attend", ylim = c(-10,40))
curve(breakfun(x), xlim = c(13,35), xlab = "temp", ylab = "attend", add = TRUE)
curve(identity(x) - 2* breakfun(x), xlim = c(13,35), xlab = "temp", ylab = "attend", add = TRUE, col ="red")
abline(v=24, lty = "dashed")
```

```{r}
x <- c(1, 10, 20, 23, 24, 25,30,35)
breakfun(x) %>% set_names(x)
```
```{r}
breakfun
```


```{r, eval =FALSE, echo=FALSE}
x <- sample(1:100, 10)
y <- sample(1:100, 10)
cbind(x,y)
pmax(x,y)
```

$$
attend = \beta_0 + \beta_1 temp + \beta_2 (temp-23)^+ + \varepsilon_i
$$

```{r}
events %>% 
  ggplot(aes(temp, attend)) +
  geom_jitter() +
  geom_smooth(se = FALSE) +
  geom_smooth(se = FALSE, method = "lm", 
              formula = y ~ x + pmax(x-24,0), col = "red") 
```

There is statistically significant relation between attendance and temperature.
    
\clearpage
## March 20: A linear regression model

Regress attendance on `month`, `day of the week`, and `bobblehead` promotion.

```{r}
lmod1 <- lm(attend ~ month + day_of_week + bobblehead, data = events)
```

Assuming that model is fine...

```{r}
summary(lmod1)
```
In the summary table, on the lines of bobbleheadYES, we test the null hypothesis: "Giving bobblehead has no effect on attendance". Null and data are consistent if p-value on the same line is large (> 0.05). If p-value is small, then data contradicts with the null, and therefore we reject the null. Since p-value = 0.0000359 is indeed small we reject the null. So we conclude that bobblehead can affect the attendance (in the positive direction because the estimate has positive value 10714.90). We expect the attendance to increase by 10715 on average.

But first check if the model is satisfactory

```{r}
plot(lmod1)
plot(lmod1, which = 5)
plot(lmod1, which = 4)
```
```{r}
plot(lmod1, which = 3)
car::ncvTest(lmod1) # null says the variance is constant. Since p is small (< 0.05), variance cannot be constant
plot(lmod1, which = 2)
shapiro.test(rstandard(lmod1)) # null says standardized residuals have normal distribution. Since p is small, the residuals do not seem to have normal distribution.
```

```{r}
car::residualPlots(lmod1)
```


```{r model2}
lmod2 <- lm(attend ~ month*day_of_week + bobblehead, events)
lmod2 %>% plot()
shapiro.test(rstandard(lmod2)) # residuals seem to have normal distr
car::ncvTest(lmod2) # variance seems to be constant
```
```{r model1}
plot(lmod1, which = 2)
shapiro.test(rstandard(lmod1))
car::ncvTest(lmod1) 
```

```{r}
summary(lmod2)
```

Compare simple and complex models using anova (useful for comparison of nested models with F-test)

```{r}
anova(lmod1, lmod2)
```

As far as we are concerned with the point predictions on the train data, these two models are not different. Null is that small model is correct. The null is consistent with data since p-value is large.

```{r}
AIC(lmod1, lmod2)
```

The lower the AIC, the better is the model. AIC picks complex model. F-test suggests simple model (model1). AIC (measures generalizibility to unseen data) and diagnostic plots recommend complex model. Since we want to predict attendance for future games, we take AIC more seriously. Therefore I will follow AIC + diagnostics recommendation and work with complex model.

**p-value is a measure of consistency of null hypothesis (=no relation between variable and attendance) with data.**

@. Is there any evidence for a relationship between attendance and other variables? Why or why not?

```{r}
lmod2 %>% summary()
```

Check F-statistic's p-value. If it is less than 0.05, then there is relation between attendance and predictors.
    
@. Does the bobblehead promotion have a statistically significant effect on the attendance?

```{r}
anova(update(lmod2, . ~ . - bobblehead), lmod2)
```
anova / F-test reject the small model (without bobblehead). Therefore we conclude that bobblehead is important to boost the attendance.

```{r}
AIC(update(lmod2, . ~ . - bobblehead), lmod2)
```

AIC confirms that bobblehead is important.

Check also with cross-validation.

@. Do month and day of week variables help to explain the number of attendants?

```{r}
anova(update(lmod2, . ~ . - month - month:day_of_week), lmod2)
```
Month is an important variable because small model's p-value is small -- therefore we reject the small model.

```{r}
update(lmod2, contrast = list(month = "contr.sum", bobblehead = "contr.sum")) %>% coef() 
```
```{r, fig.width=13}
effects::predictorEffect("month", lmod2) %>% plot()
```


@. How many fans are expected to be drawn alone by a bobblehead promotion to a home game? Give a 90\% confidence interval.

Expected number of additional fans drawn to a home game with a bobblehead promotion

```{r}
lmod2 %>% coef %>% .["bobbleheadYES"]
```

```{r}
confint(lmod2, parm = "bobbleheadYES", level = 0.95)
```

```{r}
confint(lmod2, parm = "bobbleheadYES", level = 0.80)
```

@. How good does the model fit to the data? Why? Comment on residual standard error and R$^2$. Plot observed attendance against predicted attendance. 

* Diagnostic plots to ensure the model assumptions
* Adj R2 `r round(100*summary(lmod2)$adj.r.squared)`% is mild -- try improving.
* RSE `r sigma(lmod2)` avg number of fans mistaken. Percentage error: `r round(100*sigma(lmod2)/56000)`%

@. Predict the number of attendees to a typical home game on a Wednesday in June if a bobblehead promotion is extended. Give a 90% prediction interval.

```{r, warning = FALSE}
predict(object = lmod2, newdata = data.frame(month = "JUN", day_of_week = "Wednesday", bobblehead = "YES"), interval = "prediction", level = 0.80)
```

## Project (will be graded)

Include **all variables** and conduct a full regression analysis of the problem. Submit your `R markdown` and `html` files to course homepage on moodle.

```{r}
events %>% names()
```

* Add other predictions
* Check if there are interacting variables
* Select useful variable (?step)
  - Identify important variables.
* Check if your final model fits well.
  - diagnostic plots
  - Adj R2, compare models with F-test, AIC
* After getting your final best model
  - does bobblehead still increase the attendance?
  - if so, predict attendance as in the last question with confidence intervals.
* Due date: Sunday 19:00 PM on April 23. Submit to Moodle. **Will not admit late projects.**
* Submit Rmd and html files // jupyter and html files. Write like a report to Dodgers management. Discussion of results in plain English supported with data analysis
* Groups up to three people are fine. Name your files with Bilkent ID's of group members. Write group member names and Bilkent ID's inside the report.
  
## Bibliography


