---
title: "Homework 01 Solutions"
date: "9/2/2021"
output: 
  bookdown::pdf_document2:
    keep_tex: true
    fig_caption: true
    toc: false
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
csl: "american-statistical-association.csl"
link-citations: true
---

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

# (Sheather 2.8.6)
## (a)
\begin{proof}
\begin{align*}
(y_i - \hat{y}_i) &= (y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i)) && 
                  \text{since } \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i \\
                  &= (y_i - (\bar{y} - \hat{\beta}_1\bar{x} + \hat{\beta}_1 x_i)) &&
                  \text{since } \hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x}\\
                  &= (y_i - \bar{y}) - \hat{\beta}_1(x_i - \bar{x})
\end{align*}
\end{proof}

## (b)
\begin{proof}
\begin{align*}
(\hat{y}_i - \bar{y}) &= ((\hat{\beta}_0 + \hat{\beta}_1 x_i) - \bar{y}) \\
                  &= ((\bar{y} - \hat{\beta}_1\bar{x} + \hat{\beta}_1 x_i) - \bar{y}) \\
                  &= \hat{\beta}_1(x_i - \bar{x})
\end{align*}
\end{proof}

## (c)
\begin{proof}
Substituting the equalities from parts (a) and (b) above, we have:
\begin{align*}
  \sum_{i=1}^n (y_i - \hat{y}_i)(\hat{y}_i - \bar{y}) &= 
  \sum_{i=1}^n [(y_i - \bar{y}) - \hat{\beta}_1(x_i - \bar{x})]\hat{\beta}_1(x_i - \bar{x}) \\
  &= \hat{\beta}_1 \left[\sum_{i=1}^n (y_i-\bar{y})(x_i-\bar{x}) - \hat{\beta}_1 \sum_{i=1}^n (x_i - \bar{x})^2 \right] \\
  &= \hat{\beta}_1 \left[SXY - \frac{SXY}{SXX}SXX\right] \\
  &= \hat{\beta}_1 [SXY - SXY] \\
  &= 0
\end{align*}
\end{proof}

# (Gelman & Hill 2.6)
## (a)

```{r 2a}
var1 <- rnorm(1000, 0, 1)
var2 <- rnorm(1000, 0, 1)
model <- lm(var2 ~ var1)
summary(model)
```

The slope is not significantly different from 0, as we can see from the p-value in the `var1` row, `Pr(>|t|)` column in the model summary. 

## (b)

```{r 2b, message = F, warning = F}
library(arm)  # Needed for the se.coef() function.
z.scores <- rep(NA, 100)
for (k in 1:100) {
  var1 <- rnorm(1000, 0, 1)
  var2 <- rnorm(1000, 0, 1)
  fit <- lm(var2 ~ var1)
  z.scores[k] <- coef(fit)[2]/se.coef(fit)[2]
}

type1_errors <- sum(abs(z.scores) > 2)
print(type1_errors)
```

There were `r type1_errors` z-scores that were significant.

## (c)
The rule of thumb to consider a z-score significant if the absolute value exceeds 2 is roughly equivalent to setting $\alpha = 0.05$, since 95% of the standard normal
distribution falls within 1.96 standard deviations of the mean of 0. Therefore, we should expect the coefficients to be significant 5 out of 100 times on average, corresponding to a Type 1 error rate of 0.05.

# (Sheather 2.8.7)
The regression line is defined as $\mu(x) = E(Y|X = x)$. In other words, it's the
mean value of $Y$ at a given value of the covariates $X$. The confidence interval is
for the regression line itself; that is, it's a confidence interval 
around an estimated mean. It's entirely possible for the data points themselves to fall outside the 95% confidence interval for a mean. 

To gain some intuition, imagine that you have data from a normal
distribution and you want to estimate the mean of that distribution. The more
data you gather, the narrower the confidence interval around your estimate will
be. The underlying distribution hasn't changed, however, so the more data you
gather, the more of the data will fall outside the confidence interval for the
mean.

# (Sheather 3.4.8)
## Part 1
### (a)
Given the problem constraints, the only flexibility we have is in specifying
the intercept. We could (1) fix the intercept to some value, or (2) estimate the
intercept via least squares.

Option (1) typically means setting the intercept to
0, which can be nice for interpretability in cases where the value (0, 0) is
meaningful and where a different intercept would not make sense. For example,
suppose we were asked to regress the number of words typed in a paper against
the number of hours of work put in. Zero hours of work would necessarily produce zero words, so in that case there is a good reason to set the intercept to 0.

Here, however, the size of a diamond is presumably bounded away from 0. The notion
of a zero-carat diamond doesn't make sense, and no one will manufacture a ring with
a diamond that's so small as to be invisible to the naked eye. Indeed, the sizes of diamonds
in the data range from 0.12 to 0.35 carats. If this is the range we're
interested in modeling, then we may as well give the model the additional
flexibility of a non-zero intercept, since the y-value at the intercept will
not be meaningful regardless.

Hence, we take option (2) and estimate a model of the form
\begin{align*}
  \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i
\end{align*}

Although the problem dicatates that we use a linear model, you should always visualize the data before specifying a model to help you decide whether that model is appropriate. Since we only have two variables, let's simply
look at a scatter plot of Price against Size (Figure \@ref(fig:5a-scatter)).

```{r 5a-scatter, fig.cap="Scatter plot of diamond Price against Size"}
diamonds <- read.table("../diamonds.txt", header = TRUE)
plot(Price ~ Size, data = diamonds)
```

It looks like there's a more or less linear relationship, so let's fit the model.

```{r}
model5a <- lm(Price ~ Size, data = diamonds)
```

The coefficients of the model are 

$\hat{\beta}_0$: `r round(coef(model5a)[1], 2)`

$\hat{\beta}_1$: `r round(coef(model5a)[2], 2)`

Let's add the best fit line to the scatter plot (Figure \@ref(fig:5a-scatter-model)).
```{r 5a-scatter-model, fig.cap="Scatter plot of diamond Price against Size with least squares line"}
plot(Price ~ Size, data = diamonds)
abline(a = coef(model5a)[1], b = coef(model5a)[2], col = "red", lty = "dashed")
```

The model looks reasonable. We'll examine diagnostic plots in Part 2 below.

### (b)
One weakness of the model is that by definition, it cannot capture nonlinearities in the relationship. Additionally, it has only one predictor; additional predictors could
provide additional information. It also cannot deal with nonconstant variance.

## Part 2
### (a)
In order to decide whether and how to modify the model above, let's examine some diagnostic plots.
First, the residuals against the predictor (Figure \@ref(fig:5a-residuals)).

```{r 5a-residuals, fig.cap="Model 5a residuals against predictor"}
plot(model5a$residuals ~ diamonds$Size, xlab = "Size", ylab = "Residuals")
```

There's no clear relationship (e.g. the residuals don't seem to be increasing or decreasing with Size),
so that's good. Now let's examine the diagnostic plots that R produces for us automatically (Figure \@ref(fig:5a-diagnostics)).

```{r 5a-diagnostics, fig.cap="Diagnostic plots for model 5a"}
par(mfrow = c(2, 2), mar = c(4, 4.5, 2.5, 4)) # reducing the margins slightly
plot(model5a)                                 # generating diagnostic plots
```

Let's discuss each of these in turn:

* The top left plot shows no relationship between the residuals and the fitted values (the $\hat{y}_i$s), so that's good.
* The top right plot shows that the residuals are approximately normally distributed, which is also good.
* The nonparametric red line in the bottom left plot shows a slight increasing trend between the fitted values and the standardized residuals, but there aren't many data points on the right side of the plot, so I wouldn't take this too seriously. (Recall that the variance of the residuals $\hat{e}_i$ is nonconstant, \emph{even if the variance of the underlying errors $e_i$ is constant}. The standardized residuals have approximately constant variance, so it is useful to examine these instead of the raw residuals, particularly when points of high leverage exist.)
* The bottom right plot shows that there are a couple outliers in the sense defined in Sheather: points whose standardized residual falls outside the range [-2, 2]. R has helpfully labeled these points with their row numbers in the data, so we can examine them. Additionally, recall the rule of thumb that a point is "high leverage," or equivalently a "leverage point," if its leverage value $h_{ii}$ is greater than $4/n$, where $n$ is the number of data points. In this case, $4/n = 4/49 \approx 0.082$, so we have a couple of leverage points in the data. One of these, point 42, has a relatively high Cook's distance; it lies close to the contour line that represents a Cook's distance of 0.5. 

Since R has helpfully labeled points 4, 19, and 42 for us in the bottom right diagnostic plot, let's examine these on the original graph (Figure \@ref(fig:5a-leverage-points)).

```{r 5a-leverage-points, fig.cap="Scatter plot of diamond Price against Size with least squares line and points 4, 19, 42 highlighted\\label{fig:scatter4}"}
plot(Price ~ Size, data = diamonds)
abline(a = coef(model5a)[1], b = coef(model5a)[2], col = "red", lty = "dashed")
points(diamonds[c(4, 19, 42), ], col = "red")
```

None of these points is obviously problematic. Let's check one more diagnostic, the inverse response plot:

```{r 5a-inverse-response, fig.cap="Inverse response plot for the diamonds model5a", message = F, warning = F}
# install.packages("alr3")  # if you don't already have it installed
library(alr3)
invisible(inverseResponsePlot(model5a))
```

Based on this plot, it looks like the original model, represented by the light blue line ($\lambda = 1$) fits the data well. It doesn't look like we would do better transforming the response variable, at least not with a scaled power transformation.

Given these results, it seems reasonable to keep the original model in 5(a)

### (b)
Since the model hasn't changed, see 5(a).

## Part 3
The model did not change.


# 5. IMRAD and IDMRAD analysis

## (a)
### Abstract
**Yes,** *this section has the right content.* The section nicely summarizes the later sections (introduction, method, results and discussion) and provides a clear overview.

### Introduction
**Yes,** *this section has the right content.* This section does well at motivating why we might be interested in examining the relationship between a critic's rating and the wine's price. It also clearly highlights the 3 questions that the paper will attempt to answer.

This introduction could be improved by more clearly highlighting the connection between the second and third questions to the first question (which is directly related to the motivation in the first paragraph)^[This comment could suggest that our answer should be actually "No, the introduction section does not contain the expected content".]. 

### Methods
**No,** *this section's content does not meet expectations.* This section doesn't highlight (in writing) a lot of the methods that were used to to answer the questions posed. Specifically, all 3 analytics tools are is only mentioned in passing in the final sentence of the section. No motivation is provided to suggest why such analysis was done, which should have at least occured for Figure 2's log-transformation of the price and scores. 

In terms of the organization advice given on slide 9 of lecture 02 (week01), the method section failed to use a "context, content, and conclusion" structure in presenting what methods were used. This is especially true for this section as the figures are presented first without any discussion of what was done.

As most of this section presents data information, it might be useful for the author to switch to an IDMRAD approach and put that information into a Data section. In the commentary on the **results** section of the paper in section \@ref(results), I will discuss parts of results that really should be in the methods section in more details.

### Results
**No,** *this section's content does not meet expectations.* As alluded to above, the results section contains several pieces that really should be in the **methods** section. Specifically, the descriptions of the methods / figures used should be in the method section. Additionally, the discussion of the motivation of using a log transformation should have introduced Figure 2 in the methods section.

Additionally, this section presents the idea that the slope of regression line (relating the critic's scores to the price point) is a valuable measure to compare the critic's impact even though the scales are different. This is hard to justify and muddies the water in terms of making sure the reader agrees with the conclusions (#7 of slide 9). Another, slightly nit-picky critique of this section, would be that the introduction failed to well motivate the later 2 questions and the author spent a lot of time discussion the results related to the questions without well motivating the results. This might lose a reader and relates to #9 of slide 9 - allocation the time of the reader wisely.

This section's commentary on the spread of the prices relative to the critic's score does attempt to well address the question relative to the methods used, and if I was an editor / advisor I would encourage that to be expanded on.

## Discussion
**Yes,** *this section has the right content.* Conditional on the comments above, the **discussion** section does highlight overarching conclusions, and this section spends a lot of time discussion limitations and future improvements. The author did convert potential criticisms of the work into limitations and future work ideas. 

If anything, this section spends the vast majority of presenting limitations and could have been better structured to capture that idea.

\pagebreak

## (b)
*A rewrite of the **Data** and **Methods** section. Stylistically, I tried to stay as close to the current write-up as possible.*

**Data**

The data for this study come from @parker2003  and @coates2004. The prices (in pounds sterling) include duty but exclude shipping and VAT in London in September 2003 [@sheather2009, pp. 8]. Parker’s rating system is 100-point based and the scale can be seen in Figure \@ref(fig:tab-parker). Whereas Coates’ rating system is 20-point based and can be seen in Figure \@ref(fig:tab-coates). The dataset contains prices for 72 wines from the 2000 vintage in Bordeaux, and all the variables in Figure \@ref(fig:tab-variables), and was obtained from supplementary material provided by @sheather2009.
 
```{r tab-parker, fig.cap="Parker's rating system for wine", fig.pos="h", fig.align="center", echo = F, out.width="50%"}
knitr::include_graphics("../Table01.jpg")
```

```{r tab-coates, fig.cap="Coates' rating system for wine", fig.pos="h", fig.align="center", echo = F, out.width="50%"}
knitr::include_graphics("../Table02.jpg")
```

```{r tab-variables, fig.cap="Description of dataset's variables", fig.pos="h", fig.align="center", echo = F, out.width="80%"}
knitr::include_graphics("../Table03.jpg")
```


**Methods**

In order to understand the relationship between the prices of wine and critics scores we creates a set of scatter plots relating these two variables together. Figure \@ref(fig:fig-matrixplot1) shows these relationships and the relationship between Parker and Coates' critic scores of the same wines. From \@ref(fig:fig-matrixplot1)'s plots relating both critic's scores to the price of wine, we saw that the data points demonstrate an exponential trend, as such, we created a log-log transformation of plots in Figure \@ref(fig:fig-matrixplot1) as can be seen in Figure \@ref(fig:fig-matrixplot2).

```{r fig-matrixplot1, fig.cap="Pairwise scatterplots of critic scores and prices of wines.", fig.pos="h", fig.align="center", echo = F, out.width="50%"}
knitr::include_graphics("../Fig01.jpg")
```

```{r fig-matrixplot2, fig.cap="Pairwise scatterplots of critic scores and prices of wines (log-log transformed).", fig.pos="h", fig.align="center", echo = F, out.width="50%"}
knitr::include_graphics("../Fig02.jpg")
```

Beyond critic's scores, we explored the relationship between the price of wine and some binary variables that might also effect the wine's price in Figure \@ref(fig:fig-boxplot) by comparing the distribution of prices conditional on these binary variables values.

```{r fig-boxplot, fig.cap="Boxplots comparing the price of wine conditional on additional features about the wine (see Figure \\ref{fig:tab-variables} for more details).", fig.pos="h", fig.align="center", echo = F, out.width="80%"}
knitr::include_graphics("../Fig03.jpg")
```

