---
title: "Homework 03 Solutions"
date: "9/14/2020"
output:
  bookdown::pdf_document2:
    toc: no
    keep_tex: true
---

# Problem 1

## (a) EDA

```{r}
data <- read.table("../cars04.dat",header=T)  
num_of_NA_vals <- sum(is.na(data))
```

Our data contains `r nrow(data)` unique vehicles' information and 13 features about each vehicle, and the data is very clean, with `r num_of_NA_vals` missing data values across all vehicles. Table \@ref(tab:tab-variables) contains information about all the variables in the dataset.

Table: (\#tab:tab-variables) Details about all variables in `car04.dat` data frame

Variable name   | Description
----------------| :------------------------------------------------------
`Hybrid` | Is the vehicle a hybrid? (Yes = 1, No = 0)
`SuggestedRetailPrice` | Suggested retail price, what the manufacturer thinks the vehicle is worth, including adequate profit for the	automaker and the dealer (U.S. dollars)
`DealerCost` | Price dealer had to pay manufacturer (U.S. dollars)
`EngineSize` | Engine Size (liters)
`Cylinders` | Number of cylinders
`Horsepower` | Gross Hoursepower
`CityMPG` | City miles per gallon
`HighwayMPG` | Highway miles per gallon
`Weight` | Weight of vehicle (pounds)
`WheelBase` | Wheelbase - horizontal distance between front and rear wheels (inches)
`Length` | Length of vehicle (inches)
`Width` | Width of vehicle (inches)

```{r fig.width = 10, fig.height = 7, fig.cap = "Univariate Distribution of Variables in `car04.dat`"}
par(mfrow = c(3,4))
for (c_name in colnames(data)){
  if(length(unique(data[,c_name])) <= 6){
     c_vals <- names(table(data[,c_name]))
     c_counts <- table(data[,c_name])
     barplot(names.arg = c_vals,height= c_counts, main = c_name)
  }else {
    hist(data[,c_name], main = c_name)
  }
}

## or 
# summary(data)
# table(data$Hybrid)
# 
# par(mfrow=c(3,4))
# 
# for (i in 1:12) {
#     hist(data[,i],main=names(data)[i])
# }
# 
# for (i in 1:12) {
#     boxplot(data[,i],main=names(data)[i])
# }
```


```{r, figures-side, fig.show="hold", out.width="50%", fig.cap = "Distributions of Prices and log10(Prices) relative to each other"}
par(mar = c(4, 4, .1, .1))
library(psych)

pairs.panels(data[,2:3])
log_price <- log10(data[,2:3])
names(log_price) <- paste0("log10(",names(log_price), ")")
pairs.panels(log_price)
```


## (b) Model Fitting

```{r model1}
lm.1 <- lm(SuggestedRetailPrice ~ DealerCost, data=data)
```

```{r tab-summary}
summary(lm.1)
```

```{r fig-diagnostics}
par(mfrow = c(2,2))
plot(lm.1)
```

In examining the relationship in the Section \@ref(a-eda) there was a very linear relationship between the dealer's cost and the suggested retail price. After fitting, the very expensive cars (say with predicted sales cost above $60,000) seem to  underestimated. Additional, there appears to be natural clusters in the residual vs fitted plot and with future examination of some of these clusters (using `identify`) we find that certain clusters relate to vehicles made by the same brand - (e.g. lower left 5 points of the plot are all Suzuki cars). The residuals vs normal QQ plot suggests that the residuals have slightly heavier tails (i.e. more extreme values of residuals than might be expected). The scale-location plot shows us that as the estimated value increase we see an increase in the variability (suggesting heteroskedasticity of the errors)^[This observation can also be gathered for the residual-fitted visualization.] - though it should be acknowledged that there are less observations with high estimated fitted values. Our Cook's distance figure shows that we have the Mercedes-Benz CL600 2dr has a very high impact on the model (having a Cook's distance above 1 and a high leverage score - with the average expected score being 1/`r nrow(data)`.).

```{r echo = F, eval = F}
plot(lm.1$residuals~ lm.1$fitted.values)
identify(lm.1$residuals~ lm.1$fitted.values) #37 38 39 40 83
rownames(data)[c(37, 38, 39, 40, 83)]
grep(pattern = "Suzuki", rownames(data))
```

The fitted regression would suggest that the retail price is supposed to be around 8-9% larger than the dealer cost, and it is natural to see that as the price goes up the difference in predicting the truth also scales (think as a proportion of the predicted price).

# Problem 2 Box-Cox 

```{r message = F, warning = F}
library(car)
```


## (a) Transform Exploration

### i.

```{r fig-boxcox-x}
boxCox(lm(DealerCost~1,data=data))
```
```{r}
powerTransform(lm(DealerCost~1,data=data))
```

Based on the figure above we might transform `DealerCost` by raising it to a power of $-.25$.

### ii.

```{r fig-boxcox-y}
data$TransCost <- data$DealerCost^(-.25)

boxCox(lm(SuggestedRetailPrice~TransCost,data=data))

```

```{r}
powerTransform(lm(SuggestedRetailPrice~TransCost,data=data))
```


## (b) Transformed Model

```{r}
data$srp.20 <- data$SuggestedRetailPrice^(-0.20)
data$dc.25 <- data$DealerCost^(-0.25)

par(mfrow=c(1,1))
pairs.panels(data.frame(`SuggestedRetailPrice^(-.2)` = data$srp.20, 
                        `DealerCost^(-.25)` = data$dc.25,check.names = F))

lm.2 <- lm(srp.20 ~ dc.25, data = data)

summary(lm.2)

par(mfrow=c(2,2))
plot(lm.2)
```
Relative to the original data, the transformed `Suggested RetailPrice` and `DealerCost` have distributions that are much more symmetric.  In the residuals vs Fitted plot we still see some outliers^[Looks like 2 Suzuki branded cars] when the fitted value is higher but overall the distribution of residuals looks much better - it has a conditional average residual value around zero and the variability doesn't look too different as we look across the range of fitted values. The residuals now only have a heavy positive tail^[The heavy positive tail on the QQ plot is directly related to those large residuals in the residual vs fitted plot.]. The scale-Location plot confirms out observation from the residual vs Fitted that this new model with the transformations seems to have homoscedastic residuals. The Cook's distance plot also shows less impact from highly influential points, there is a point with leverage at .04 (compared with the average `r round(1/nrow(data),3)`, but this plot looks better than the original model.

I'm pretty pleased with this new model, most of the linear model assumptions look well met except for a few outliers (which where also identified in our examination of the first model).

# Problem 3: Log-Log

Note that we used log-based 10 (this is just a good think to be thinking about when you're transforming and are interested in interpretation).

```{r}
lm.3 <- lm(log10(SuggestedRetailPrice) ~ log10(DealerCost), data=data)

log_price <- log10(data[,2:3])
names(log_price) <- paste0("log10(",names(log_price), ")")
pairs.panels(log_price)

summary(lm.3)

par(mfrow=c(2,2))
plot(lm.3)
```

For our "Log-Log" model, there appears to be a slight trend in the average residual as we increase the fitted value, with a slightly positive average residual for the majority of the fitted value values with log(fitted) around 4.2 to 4.6 but negative residual for fitted values on the the extremes. There is also a set of cars with very low residuals with log(fitted) values between 4.2 and 4.6. The normal QQ plot looks decent, with a slightly heavy tail with negative residuals. For the extremely large fitted values we might also be seeming an indication for increased variability of residuals in the Scale-Location plot (although with so few points, this claim is probably premature). There is a group of points with larger leverage than the average, but their Cook's distances aren't that large.

This fit looks much better than the original regression, but there are still some odd groupings in the residual-fitted plot. Overall it looks like a decent model.

# Problem 4

## (a) Best Model (based on assumptions)

The Box-Cox transformed model seems to best meet the assumptions of the linear regression, specifically in the true linear relationship (examined with $E[\hat{e}|x] = 0$). Additionally, the second model seems well meet the assumption of constant noise variance. 

## (b) "Complete the Sentence"

- For the model in #1(b), a \$1 increase in `DealerCost` results in a \$`r round(coef(lm.1)[2],2)` increase in the expected value of `SuggestedRetailPrice`.
- For the model in #2(b), a 1 unit increase in $($`DealerCost`$)^{-.25}$ results in a `r round(coef(lm.2)[2],2)`% increase in the expected value of $($`SuggestedRetailPrice`$)^{-.2}$.
- For the model in #3, a 1% increase in `DealerCost` results in a `r round(coef(lm.3)[2],2)`% increase in the expected value of `SuggestedRetailPrice`.


## (c) Best Model (for consulting)

Giving the difficulting of translating the relationship between `DealerCost` to `SuggestedRetailPrice` for the Box-Cox transformed model, I would probably use the Log-Log model for a consulting project. The Log-Log model is still interpretable and more closely meets the OLS assumptions.

# Problem 5

```{r}
pyth_data <- read.table("../pyth.dat",header = T)
```

## (a)

\begin{align*}
M1: \texttt{y} = & \beta_0+\beta_1 \texttt{x1}+\varepsilon \\
M2: \texttt{y} = & \beta_0+\beta_1 \texttt{x2}+\varepsilon
\end{align*}

```{r}
lm.1 <- lm(y~x1, data=pyth_data)
lm.2 <- lm(y~x2, data=pyth_data)
summary(lm.1)
summary(lm.2)
```

```{r out.width = "50%", fig.show="hold"}
plot(y~x1, data = pyth_data)
plot(y~x2, data = pyth_data)

par(mfrow=c(2,2))

plot(lm.1)

plot(lm.2)
```

$M2$ is better, it reduces more of the variability of $y-\hat{y}$ by much more than $M1$.


## (b)

\begin{align*}
M3:  \texttt{y2} = & \beta_0+\beta_1 \texttt{x12}+\varepsilon \\
M4:  \texttt{y2} = & \beta_0+\beta_1 \texttt{x22}+\varepsilon
\end{align*}

```{r}

pyth_data$y2 <- pyth_data$y^2
pyth_data$x12 <- pyth_data$x1^2
pyth_data$x22 <- pyth_data$x2^2

## OR
# attach(pyth_data)
# 
# y2 <- y^2
# x12 <- x1^2
# x22 <- x2^2
# 
# detach()


lm.3 <- lm(y2~x12, data=pyth_data)

lm.4 <- lm(y2~x22, data=pyth_data)

summary(lm.3)
summary(lm.4)

```

```{r out.width = "50%", fig.show="hold"}
par(mfrow=c(2,2))

plot(lm.3)

plot(lm.4)
```


$M4$ is better, it reduces more of the variability of $y^2-\hat{\texttt{y2}}$ by much more than $M3$.

## (c)

\begin{align*}
M5: \texttt{y} = & \beta_0+\beta_1 \texttt{x1}+\beta_2\texttt{x2}+\varepsilon \\
M6: \texttt{y2} = & \beta_0+\beta_1 \texttt{x12}+\beta_2\texttt{x22}\varepsilon
\end{align*}

```{r}
lm.5 <- lm(y ~ x1 + x2, data=pyth_data)

lm.6 <- lm(y2 ~ x12 + x22, data=pyth_data)

summary(lm.5)
summary(lm.6)
```

```{r out.width = "50%", fig.show="hold"}

par(mfrow=c(2,2))

plot(lm.5)

plot(lm.6)
```

Both $M5$ and $M6$ represent a lot their respectively output, but $M6$ basically captures all variablity in $y^2$.

## (d)

\begin{align*}
M7: \texttt{y} = & \beta_0+\beta_1 \sqrt{\texttt{x12}+\texttt{x22}}+\varepsilon \\
\end{align*}

```{r}
pyth_data$x3 <- sqrt(pyth_data$x12 + pyth_data$x22)

lm.7 <- lm(y ~ x3,data=pyth_data)

summary(lm.7)

par(mfrow=c(2,2))

plot(lm.7)
```

Based on observations made on $M6$, $M7$ seems to capture basically all the variability of $y$. This relates to the fact that $M6$ was basically suggesting that $\texttt{y2} \approx 0 + 1 \texttt{x12} + 1 \texttt{x22}$, which could be re-expressed as 

\begin{align*}
y^2 &\approx 0 + 1 \texttt{x1}^2 + 1 \texttt{x2}^2 \\
\Leftrightarrow y &\approx \sqrt{\texttt{x1}^2 + \texttt{x2}^2} \;.\\
\end{align*}

