---
jupyter:
  jupytext:
    formats: Chapter11//ipynb,jupyterbook/markdown/Chapter11//md:myst,Chapter11//Rmd
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.16.1
  kernelspec:
    display_name: R
    language: R
    name: ir
---

<!-- #region -->


# Diagnostics II

## Download

- [Slides](http://www.stanford.edu/class/stats191/Chapter11/Diagnostics_II.html)

- RStudio: [RMarkdown](http://www.stanford.edu/class/stats191/Chapter11/Diagnostics_II.Rmd),[Quarto](http://www.stanford.edu/class/stats191/Chapter11/Diagnostics_II.qmd)    
    
- [Jupyter](https://nbviewer.org/urls/www.stanford.edu/class/stats191/Chapter11/Diagnostics_II.ipynb)    
    

<!-- #endregion -->

## Gesell Score Data Set {.scrollable}

The Gesell score is used to measure child development. 

How does it correlate with the age (in months) at which a child speaks their 
first word?

```{r editor-font-scale=0.5, webr=TRUE}
df.child <- read.csv("https://dlsun.github.io/data/gesell.csv")
head(df.child)
```

## Influential Points

Which point do you think has the most "influence" on the fitted regression line?

```{r}
#| echo: false
df.child <- read.csv("https://dlsun.github.io/data/gesell.csv")
with(df.child, {
  plot(AgeAtFirstWord, GesellScore);
  abline(lm(GesellScore ~ AgeAtFirstWord));
  points(AgeAtFirstWord[18:19], GesellScore[18:19], col="red", pch=15);
  text(AgeAtFirstWord[18:19] - 0.7, GesellScore[18:19], labels=18:19, col="red")
})
```

## Influential Points {.scrollable}

Let's try deleting each point and seeing which one changes the regression
line more.

```{r editor-font-scale=0.5, webr=TRUE}
with(df.child, plot(AgeAtFirstWord, GesellScore))
abline(lm(GesellScore ~ AgeAtFirstWord, data=df.child))
abline(lm(GesellScore ~ AgeAtFirstWord, data=df.child[-18, ]),
       lty=2, col="blue")
```

<!-- #region -->
# Measuring Influence

How do we quantify how much "influence" a point has on the regression line?


## Mathematical Detour: The Hat Matrix

In thinking about influence, it is helpful to note that 
each fitted value $\hat Y_i$ can be expressed as a weighted sum of the 
$Y$ values:

$$ \hat Y_i = \sum_{j=1}^n h_{ij} Y_j. $$

. . .

Collectively, the weights $h_{ij}$ make up a matrix $H$.

. . .

$H$ is called the "hat matrix" because it puts the "hat" on $Y$.


## Diagnostic 1: Leverage

The **leverage** of an observation is the weight $h_{ii}$ that observation $i$ 
assigns to its _own_ $Y$ value.

::: {.incremental}
- For simple linear regression, there is an easy formula $h_{ii}$:
$$ h_{ii} = \frac{(X_i - \bar X)^2}{\sum_i (X_i - \bar X)^2} + \frac{1}{n}. $$
- For multiple regression, the formula requires matrix algebra.
:::

. . .

Either way, leverage measures how far observation $i$ is from the 
average, in the _explanatory_ variables.


## Diagnostic 1: Leverage {.smaller .scrollable}

In R, the `hatvalues` function calculates the leverage of each observation.
<!-- #endregion -->

```{r editor-font-scale=0.7, webr=TRUE}
model <- lm(GesellScore ~ AgeAtFirstWord, data=df.child)
hatvalues(model)
```

<!-- #region -->
. . .

Observation 18 has the highest leverage by far. This is no surprise because 
it had a very large $X$ value.

. . .

Some Properties of Leverage

::: {.incremental}
- $\frac{1}{n} \leq h_{ii} < 1$
- $\sum_i h_{ii} = p + 1$ (number of coefficients including the intercept), 
so the average leverage is $\bar h = \frac{p + 1}{n}$.
    - Some books recommend flagging observations with $h_{ii} > 2\bar h = 2\frac{p + 1}{n}$.
:::


## Diagnostic 2: DFFIT

Maybe we need a metric that looks more directly at how much the prediction 
changes when we delete observation $i$.

. . .

We compare the following predictions:

- $\hat Y$: predicted values from a model fit to all observations.
- $\hat Y^{(-i)}$: predicted values from a model fit to all observations, 
_except_ observation $i$.

. . .

A reasonable diagnostic is the "difference in fits":

$$ DFFIT_i = \hat Y_i - \hat Y^{(-i)}_i$$


## Diagnostic 2: DFFIT {.scrollable}

There is no built-in function to calculate DFFIT. (The function `dffits`
calculates a studentized version called DFFITS.)

. . .

However, it is not hard to implement DFFIT from scratch.
<!-- #endregion -->

```{r editor-font-scale=0.5, webr=TRUE}
i <- 18
model1 <- lm(GesellScore ~ AgeAtFirstWord, data=df.child)
model2 <- lm(GesellScore ~ AgeAtFirstWord, data=df.child[-i, ])
model1$fitted[i] - model2$fitted[i]
```

<!-- #region -->
## Diagnostic 3: Cook's Distance {.smaller}

:::: {.columns}
::: {.column width="30%"}
![](http://www.stanford.edu/class/stats191/Chapter11/figs/cook.jpeg)

Dennis Cook (1944 - )

Professor of Statistics, University of Minnesota
:::

::: {.column width="70%"}
**Cook's distance** measures how much the predictions for _all_ $n$ observations 
change, when observation $i$ is deleted:

$$ D_i = \frac{\sum_{k=1}^n (\hat Y_k - \hat Y^{(-i)}_k)^2}{p\hat\sigma^2}. $$

::: {.fragment}
Fortunately, this can be calculated without refitting the model
each time a different observation is deleted. An equivalent formula is:

$$ D_i = \frac{\text{residual}_i^2}{p\hat\sigma^2} \frac{h_{ii}}{(1 - h_{ii})^2}. $$
:::

:::
::::


## Diagnostic 3: Cook's Distance {.scrollable}

In R, the `cooks.distance` function calculates the Cook's distance of each observation.
<!-- #endregion -->

```{r editor-font-scale=0.5, webr=TRUE}
cooks.distance(model)
```

## Visual Comparison of Diagnostics

```{r}
#| echo: false
with(df.child, plot(AgeAtFirstWord, GesellScore))

model1 <- lm(GesellScore ~ AgeAtFirstWord, data=df.child)
model2 <- lm(GesellScore ~ AgeAtFirstWord, data=df.child[-18, ])

abline(model1)
abline(model2, lty=2)

arrows(df.child$AgeAtFirstWord[18], df.child$GesellScore[18],
       mean(df.child$AgeAtFirstWord), df.child$GesellScore[18],
       col="red", code=0, lwd=2)
arrows(df.child$AgeAtFirstWord[18] - 0.15, model1$fitted[18],
       df.child$AgeAtFirstWord[18] - 0.15, predict(model2, df.child)[18],
       col="blue", code=0, lwd=2)

for(i in 1:nrow(df.child)) {
  arrows(df.child$AgeAtFirstWord[i], model1$fitted[i],
         df.child$AgeAtFirstWord[i], predict(model2, df.child)[i],
         col="orange", code=0, lwd=2)
}

legend(x=35, y=120, 
       legend=c("Leverage", "DFFIT", "Cook's distance"), col=c("red", "blue", "orange"),
       lty=c(1, 1, 1), lwd=c(2, 2, 2))
```

## Putting It All Together {.scrollable}

R provides a function, `influence.measures`, that calculates these diagnostics
and more.

```{r editor-font-scale=0.5, webr=TRUE}
influence.measures(model)
```

<!-- #region -->
## How to Deal with Influential Observations

![](http://www.stanford.edu/class/stats191/Chapter11/figs/flowchart.png)


# Measuring Outliers

Even if observation 19 is not influential, it still has a large residual, 
so it may be an outlier.

How do we detect outliers?


## Studentized Residual {.smaller}

Recall the difference between $z$ and $t$ statistics. (See _Sleuth_, Chapter 2.2.3.)

:::: {.columns}
::: {.column width="50%"}
$$ z = \frac{\text{estimate} - \text{parameter}}{\text{SD(estimate)}} $$

Dividing by the SD is called _standardization_.

:::
::: {.column width="50%"}
$$ t = \frac{\text{estimate} - \text{parameter}}{\text{SE(estimate)}} $$

Dividing by the SE is called _studentization_.
:::
::::

. . .

What would a $t$-statistic for residuals look like?

::: {.incremental}
- Residuals have mean $0$.
- Residuals have SD $\sigma\sqrt{1 - h_{ii}}$.
:::

. . .

So the _studentized residuals_ are 
$\displaystyle t_i = \frac{\text{residual}_i - 0}{\hat\sigma\sqrt{1 - h_{ii}}}$.


## A Test for Outliers

Under the null hypothesis that observation $i$ comes from the model (i.e., 
is not an outlier), $t_i$ should follow a $t$-distribution.

::: {.incremental}
- For this reason, some books recommend flagging observation $i$ as an outlier 
if $|t_i| > 2$, since $t_{n-1}(.975) \approx 2$.
- But this ignores the fact that any one of the $n$ observations
could be an outlier! The familywise error rate is $\gg .05$.
- We can use the Bonferroni correction. Test each observation at 
$\alpha = .05 / n$ so that the familywise error rate is $\leq .05$.
:::


## Is Observation 19 an Outlier? {.smaller .scrollable}

First, let's calculate the studentized residual
$$ t_{19} = \frac{\text{residual}_{19} - 0}{\hat\sigma\sqrt{1 - h_{19,19}}}. $$
<!-- #endregion -->

```{r editor-font-scale=0.7, webr=TRUE}
stud.res <- (model$resid[19] /
             (summary(model)$sigma * sqrt(1 - hatvalues(model)[19])))
stud.res
```

. . .

Although this appears significant, we should account for the fact that there 
are 21 observations and any one of them could be an outlier. The actual cutoff
we should use is:

```{r editor-font-scale=0.7, webr=TRUE}
n <- 21
alpha <- .05 / n
qt(1 - alpha / 2, df=n - 1)
```

# Summary

- Influential Points
    - Leverage
    - DFFIT
    - Cook's distance
- Outliers
    - Studentized residual can be used as an outlier test, but need to be careful
    about multiple testing.
