---
title: "3.9: Logarithmic Models - R Practice"
author: "Ryan Safner"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
    #theme: 
    toc: true 
    toc_depth: 3
    toc_float: true
    code_folding: show
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Set Up
To minimize confusion, I suggest creating a new `R Project` (e.g. `regression_practice`) and storing any data in that folder on your computer.

Alternatively, I have made a project in R Studio Cloud that you can use (and not worry about trading room computer limitations), with the data already inside (you will still need to assign it to an object).

- [<i class="fas fa-cloud"></i> View Project on R Studio Cloud](https://rstudio.cloud/project/739808)
- [<i class="fab fa-markdown"></i> `Download Markdown`](/r/21-r-practice-rmd.Rmd)

## Question 1

We are returning to the speeding tickets data that we began to explore in [R Practice 15 on Multivariate Regression](http://metricsf19.classes.ryansafner.com/practice/15-practice) and [R Practice 19 on Dummy Variables & Interaction Effects](http://metricsf19.classes.ryansafner.com/practice/19-r-practice). Download and read in (`read_csv`) the data below.

- [<i class="fas fa-table"></i> `speeding_tickets.csv`](/data/speeding_tickets.csv)

This data again comes from a paper by Makowsky and Strattman (2009) that we will examine later. Even though state law sets a formula for tickets based on how fast a person was driving, police officers in practice often deviate from that formula. This dataset includes information on all traffic stops. An amount for the fine is given only for observations in which the police officer decided to assess a fine. There are a number of variables in this dataset, but the one's we'll look at are:

| Variable | Description |
|----------|-------------|
| `Amount` | Amount of fine (in dollars) assessed for speeding |
| `Age` | Age of speeding driver (in years) |
| `MPHover` | Miles per hour over the speed limit |
| `Black` | Dummy $=1$ if driver was black, $=0$ if not |
| `Hispanic` | Dummy $=1$ if driver was Hispanic, $=0$ if not |
| `Female` | Dummy $=1$ if driver was female, $=0$ if not |
| `OutTown` | Dummy $=1$ if driver was not from local town, $=0$ if not |
| `OutState` | Dummy $=1$ if driver was not from local state, $=0$ if not |
| `StatePol` | Dummy $=1$ if driver was stopped by State Police, $=0$ if stopped by other (local) |

> We again want to explore who gets fines, and how much.

---

<!--ANSWER BELOW HERE-->

```{r}
library(tidyverse)

speed<-read_csv("https://metricsf19.classes.ryansafner.com/data/speeding_tickets.csv")
```

---

## Question 2

Run a regression of `Amount` on `Age`. Write out the estimated regression equation, and interpret the coefficient on Age.

---

<!--ANSWER BELOW HERE-->

```{r}
reg_linear<-lm(Amount ~ Age, data = speed)
summary(reg_linear)
```

$\widehat{\text{Amount}_i}=131.71-0.29\text{Age}_i$

For every year of age, expected fines decrease by $0.29. 

---

## Question 3

Is the effect of `Age` on `Amount` nonlinear? Let's run a quadratic regression. 

### Part A

Create a new variable for $Age^2$. Then run a quadratic regression.

---

<!--ANSWER BELOW HERE-->

```{r}
# make Age_sq variable
speed<-speed %>%
  mutate(Age_sq=Age^2)

# run quadratic regression
reg_quad<-lm(Amount ~ Age + Age_sq, data = speed) 
summary(reg_quad) 
```

---

### Part B

Try running the same regression using the alternate notation: `lm(Y~X+I(X^2))`. This method allows you to not have to create a new variable first. Do you get the same results?

---

<!--ANSWER BELOW HERE-->

```{r}
reg_quad_alt<-lm(Amount ~ Age + I(Age^2), data = speed) 
summary(reg_quad_alt) 
```

Yes, this is the same output.

---

### Part C

Write out the estimated regression equation.

---

<!--ANSWER BELOW HERE-->

$$\widehat{\text{Amount}_i}=146.75-1.17\text{Age}_i+0.01\text{Age}_i^2$$

---

### Part D

Is this model an improvement from the linear model?^[Check $R^2$.]

---

<!--ANSWER BELOW HERE-->

Yes, a slight improvement. $R^2$ went from 0.004 on the linear model to 0.006 on the quadratic model. 

---

### Part E

Write an equation for the marginal effect of `Age` on `Amount`.

---

<!--ANSWER BELOW HERE-->

The marginal effect is measured by the first derivative of the regression equation with respect to Age. But you can just remember the resulting formula and plug in the parameters:

$$\begin{align*}
	\frac{d \, Y}{d \, X} &= \beta_1+2\beta_2 X\\
	\frac{d \, Amount}{d \, Age} &= -1.17+2(0.01) \, Age\\
	&=-1.17+0.02 \, Age\\
\end{align*}$$

```{r}
# PUT CODE HERE
```

---

### Part F

Predict the marginal effect on `Amount` of being one year older when you are 18. How about when you are 40?

---

<!--ANSWER BELOW HERE-->

For 18 year olds: 

$$\begin{align*}
	\frac{d \, Amount}{d \, Age} &=-1.17+0.02(18)\\
 		&=-1.17+0.36\\
 		&=-\$0.81\\
\end{align*}$$

For 40 year olds:

$$\begin{align*}
	\frac{d \, Amount}{d \, Age} &=-1.17+0.02(40)\\
 		&=-1.17+0.80\\
 		&=-\$0.37\\
 	\end{align*}$$

```{r}
# Let's do this in R:

# we need broom
library(broom)
tidy_reg_quad<-tidy(reg_quad)

tidy_reg_quad

# save beta 1 
quad_beta_1<-tidy_reg_quad %>%
  filter(term=="Age") %>%
  pull(estimate)

# save beta 2
quad_beta_2<-tidy_reg_quad %>%
  filter(term=="Age_sq") %>%
  pull(estimate)

# create function to estimate marginal effects
marginal_effect<-function(x){
  return(quad_beta_1+2*quad_beta_2*x)
}

# run the function on the 18-year-old and the 40-year-old
marginal_effect(c(18,40))

# close enough, we had some rounding error
```

---

### Part G

Our quadratic function is a $U$-shape. According to the model, at what age is the amount of the fine minimized?

---

<!--ANSWER BELOW HERE-->

We can set the derivative equal to 0, or you can just remember the formula and plug in the parameters: 

$$\begin{align*}
	\frac{d Y}{d X}  &= \beta_1+2\beta_2 X\\
	0 &=\beta_1+2\beta_2 X\\
	-\beta_1&=2\beta_2 X\\
	-\frac{1}{2} \frac{\beta_1}{\beta_2}&=Age^*\\
	-\frac{1}{2} \frac{-1.17}{0.01} &= Age^*\\
	-\frac{1}{2} 117 & \approx Age^*\\
	58.5 & \approx Age ^*\\
\end{align*}$$


```{r}
# Let's do this in R:

-0.5*(quad_beta_1/quad_beta_2)

# again, some rounding error
```

---

### Part H

Create a scatterplot between `Amount` and `Age` and add a a layer of a linear regression (as always), and an additional layer of your predicted quadratic regression curve. The regression curve, just like any regression *line*, is a `geom_smooth()` layer on top of the `geom_point()` layer. We will need to customize `geom_smooth()` to `geom_smooth(method="lm", formula="y~x+I(x^2)` (copy/paste this verbatim)! This is the same as a regression line (`method="lm"`), but we are modifying the formula to a polynomial of degree 2 (quadratic): $y=a+bx+cx^2$.

---

<!--ANSWER BELOW HERE-->

```{r, fig.retina=3}
ggplot(data=speed)+
  aes(x = Age,
      y = Amount)+
  geom_point()+
  geom_smooth(method="lm",
              formula="y~poly(x,2)",
              color="red")
```

---

### Part I

It's quite hard to see the quadratic curve with all those data points. Redo another plot and this time, only keep the quadratic `stat_smooth()` layer and leave out the `geom_point()` layer. This will only plot the regression curve.

---

<!--ANSWER BELOW HERE-->

```{r, fig.retina=3}
ggplot(data=speed)+
  aes(x = Age,
      y = Amount)+
  geom_smooth(method="lm",
              formula="y~poly(x,2)",
              color="red")
```

---

## Question 4

Should we use a higher-order polynomial equation? Run a cubic regression, and determine whether it is necessary.

---

<!--ANSWER BELOW HERE-->

```{r}
reg_cube<-lm(Amount ~ Age + I(Age^2) + I(Age^3), data = speed)
summary(reg_cube)
```

$$\widehat{\text{Amount}_i}=151.95-1.61\text{Age}_i+0.02\text{Age}_i^2-0.00008\text{Age}^3$$

The $t$-statistic on Age$^3$ is small (-1.03) and the $p$-value is 0.31, so the cubic term does not have a significant impact on `Amount`. We should *not* include it.

Just for fun, would the cubic model *look* any better?

```{r, fig.retina=3}
ggplot(data=speed)+
  aes(x = Age,
      y = Amount)+
  geom_point()+
  geom_smooth(method="lm", formula="y~x+I(x^2)", color="red")+
  geom_smooth(method="lm", formula="y~x+I(x^2)+I(x^3)", color="orange")
```

```{r, fig.retina=3}
ggplot(data=speed)+
  aes(x = Age,
      y = Amount)+
  geom_smooth(method="lm", formula="y~x+I(x^2)", color="red")+
  geom_smooth(method="lm", formula="y~x+I(x^2)+I(x^3)", color="orange")
```

---

## Question 5

Run an $F$-test to check if a nonlinear model is appropriate. Your null hypothesis is $H_0: \beta_2=\beta_3=0$ from the regression in pert (h). The command is `linearHypothesis(reg_name, c("var1", "var2"))` where `reg_name` is the name of the `lm` object you saved your regression in, and `var1` and `var2` (or more) in quotes are the names of the variables you are testing. This function requires (installing and) loading the "`car`" package (additional regression tools).

---

<!--ANSWER BELOW HERE-->

```{r}
# install.packages("car") # install once if you don't have
library("car") # load car package
linearHypothesis(reg_cube, c("I(Age^2)", "I(Age^3)")) # F-test 
```

We get a large $F$ of 26.43, with a very small $p$-value. Therefore, we can reject the null hyothesis that the model is linear $(\beta_2=0, \beta_3=0)$. We should in fact *not* use a linear model. Note it does *not* tell us if the model should be *quadratic* or *cubic* (or even *logarithmic* of some sort), *only* that it is not linear. Remember, this was a *joint* hypothesis of all *non-linear* terms $(beta_2$ and $\beta_3)$! 


---

## Question 6

Now let's take a look at speed (`MPHover` the speed limit).

### Part A
Creating new variables as necessary, run a **linear-log** model of `Amount` on `MPHover`. Write down the estimated regression equation, and interpret the coefficient on `MPHover` $(\hat{\beta_1})$. Make a scatterplot with the regression line.^[Hint: The simple `geom_smooth(method="lm")` is sufficient, so long as you use the right variables on the plot!]

---

<!--ANSWER BELOW HERE-->

```{r}
# create log of MPHover
speed<-speed %>%
  mutate(log_mph=log(MPHover))

# Run linear-log regression
linear_log_reg<-lm(Amount ~ log_mph, data = speed)
summary(linear_log_reg)

# note we could have done this without creating the variable
# just take log() inside the regression:
linear_log_reg_alt<-lm(Amount ~ log(MPHover), data = speed)
summary(linear_log_reg)
```

$$\widehat{\text{Amount}_i}=-200.10+115.75\text{ln(MPHover}_i)$$

A 1% increase in speed (over the speed limit) increases the fine by $\frac{115.75}{100}\approx \$1.16$

---

### Part B

Creating new variables as necessary, run a **log-linear** model of `Amount` on `MPHover`. Write down the estimated regression equation, and interpret the coefficient on `MPHover` $(\hat{\beta_1})$. Make a scatterplot with the regression line.^[Hint: The simple `geom_smooth(method="lm")` is sufficient, so long as you use the right variables on the plot!]

---

<!--ANSWER BELOW HERE-->

```{r}
# create log of Amount
speed<-speed %>%
  mutate(log_Amount=log(Amount))

# Run log-linear regression
log_linear_reg<-lm(log_Amount ~ MPHover, data = speed)
summary(log_linear_reg)

# again we could have done this without creating the variable
# just take log() inside the regression:
log_linear_reg_alt<-lm(log(Amount) ~ MPHover, data = speed)
summary(log_linear_reg_alt)
```


$$\widehat{\text{ln(Amount}_i)}=3.87+0.05\text{MPHover}_i$$

For every 1 MPH in speed (over the speed limit), expected fine increases by $0.05 \times 100\%=5\%$

---

### Part C

Creating new variables as necessary, run a **log-log** model of `Amount` on `MPHover`. Write down the estimated regression equation, and interpret the coefficient on `MPHover` $(\hat{\beta_1})$. Make a scatterplot with the regression line.^[Hint: The simple `geom_smooth(method="lm")` is sufficient, so long as you use the right variables on the plot!]

---

<!--ANSWER BELOW HERE-->

```{r}
# Run log-log regression
log_log_reg<-lm(log_Amount ~ log_mph, data = speed)
summary(log_log_reg)

# again we could have done this just taking log()s inside the regression:
log_log_reg_alt<-lm(log(Amount) ~ log(MPHover), data = speed)
summary(log_log_reg_alt)
```

$$\widehat{\text{ln(Amount}_i)}=2.31+0.86\text{ln(MPHover}_i)$$

For every 1% increase in speed (over the speed limit), expected fine increases by 0.86%.

---

### Part D

Which of the three log models has the best fit?^[Hint: Check $R^2$]

---

<!--ANSWER BELOW HERE-->

We can compare the $R^2$'s of the three models or compare scatterplots with the regression lines. I will make a table of the three regressions with `huxreg` for easy comparison of fit: 


```{r}
library(huxtable)
huxreg("Linear-Log" = linear_log_reg,
       "Log-Linear" = log_linear_reg,
       "Log-Log" = log_log_reg,
       coefs = c("Constant" = "(Intercept)",
                 "MPH Over" = "MPHover",
                 "Log(MPH Over)" = "log_mph"),
       statistics = c("N" = "nobs",
                      "R-Squared" = "r.squared",
                      "SER" = "sigma"),
       number_format = 2)
```

It appears the linear-log model has the best fit with the highest $R^2$ out of the three, but not by very much. 

If we wanted to compare scatterplots:

```{r, fig.height=2}
ggplot(data = speed)+
  aes(x = log_mph,
      y = Amount)+
  geom_point()+
  geom_smooth(method="lm", color="red")+
  labs(title="Linear-Log Model")
```

```{r, fig.height=2}
ggplot(data = speed)+
  aes(x = MPHover,
      y = log_Amount)+
  geom_point()+
  geom_smooth(method="lm", color="red")+
  labs(title="Log-Linear Model")
```

```{r, fig.height=2}
ggplot(data = speed)+
  aes(x = log_mph,
      y = log_Amount)+
  geom_point()+
  geom_smooth(method="lm", color="red")+
  labs(title="Log-Log Model")
```

---

## Question 7

Return to the quadratic model. Run a quadratic regression of `Amount` on `Age`, `Age`$^2$, `MPHover`, and all of the race dummy variables. Test the null hypothesis: *"the race of the driver has no effect on Amount"*

---

<!--ANSWER BELOW HERE-->

```{r}
full_reg<-lm(Amount~Age+Age_sq+MPHover+Black+Hispanic, data=speed)
summary(full_reg)

library(car)
linearHypothesis(full_reg, c("Black", "Hispanic"))
```

---

## Question 8

Now let's try standardizing variables. Let's try running a regression of `Amount` on `Age` and `MPHover`, but standardizing each variable. 

### Part A

Create new standardized variables for `Amount`, `Age`, and `MPHover`.^[Hint: use the `scale()` function inside of `mutate()`]

---

<!--ANSWER BELOW HERE-->

```{r}
# make standardized variables
speed<-speed %>%
  mutate(std_Amount = scale(Amount),
         std_Age = scale(Age),
         std_MPH = scale(MPHover))
```

---

### Part B

Run a regression of standardized `Amount` on standardized `Age` and `MPHover`. Interpret $\hat{\beta_1}$ and $\hat{\beta_2}$ Which variable has a bigger effect on `Amount`? 

---

<!--ANSWER BELOW HERE-->

```{r}
# make standardized variables
std_reg<-lm(std_Amount ~ std_Age + std_MPH, data = speed)
summary(std_reg)
```

