---
title: "Outliers and Robust Regression"
author: "Jeffrey B. Arnold"
date: "05/20/2015"
output:
  html_document:
    toc: true
    fig_height: 3
    fig_width: 5
---
```{r echo = FALSE, results = 'hide'}
knitr::opts_chunk$set(cache = TRUE, autodep = TRUE,
                  fig.height = 3, fig.width = 5)
```

This example works through diagnostics for outliers, as well as methods of robust
regression.

## Setup

This example will use the following

```{r message = FALSE}
library("MASS")
library("dplyr")
library("tidyr")
library("broom")
library("boot")
library("ggplot2")
```
This ensures that we always use the `select` function from **dplyr** rather than 
the one from **MASS**.
```{r}
select <- dplyr::select
```
For the **ggplot2** plots, we will the default theme settings here, so that we
can reuse them for all plots, and also, if we feel like changing them, we only
need to change them in one location.
```{r}
theme_local <- theme_minimal
```


## Iver and Soskice Data

This is an example of from Iversen and Soskice (2003).
That paper is interested in the relationship between party systems and redistributive efforts of the government. 

The party system is measured using the effective number of parties;
the redistributive efforts of the government is measured as the percent people lifted from poverty by taxes and transfers

First, let's load the data
```{r}
iver <- read.csv("http://pols503.github.io/pols_503_sp15/data/iver.csv")
glimpse(iver)
```
The variables of interest are `lnemp` (log effective number of parties),
and `povred` (poverty reduction).
Let's plot the relationship between them
```{r}
ggplot(iver, aes(x = lnenp, y = povred)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("log(Number of Effective parties)") +
  ylab("Poverty Reduction") +
  theme_local()
```

## Influential Observations

What are influential points in a regression? 
They are points that 
How much would the regression line change if we deleted a the point and reran the regression? 
```{r}
iver_mod1 <- lm(povred ~ lnenp, data = iver)

iver_loo_regs <-
  # Start with the iver data
  iver %>%
  # Group by country
  group_by(cty) %>%
  # For each country
  # Run the regression without that country and store the coefficient values
  do({
    tidy(lm(povred ~ lnenp, data = filter(iver, cty != .$cty))) %>%
      select(term, estimate)
  }) %>%
  # Reshape the dataset so that each coefficient is in a column
  spread(term, estimate) %>%
  # Calculate how much these slopes differ from the one with all the data
  mutate(diff_slope = lnenp - coef(iver_mod1)["lnenp"],
         abs_diff_slope = abs(diff_slope)) %>%
  # Sort by the difference in slopes
  arrange(- abs_diff_slope)
    
iver_loo_regs
```
Switzerland looks particularly problematic.
The effect of `lnenp` on `povred` is 7.

We could also plot these lines against the original data, to get a more intuitive sense of how much dropping one observation affects the regression slopes.
```{r}
ggplot() +
  geom_abline(data = iver_loo_regs, aes(intercept = `(Intercept)`,
                                        slope = lnenp)) + 
  geom_point(data = iver, aes(x = lnenp, y = povred)) +
  xlab("log(Number of Effective parties)") +
  ylab("Poverty Reduction") +
  theme_local()
```

Conveniently, in linear regression we can find which observations will have the largest influence on regression lines without rerunning the regression.
Three statistics are of interest:

- Cook's distance: a single number that summarizes how much dropping an observation changes **all** the regression coefficients.
- Studentized residual: The scaled residual of the observation. 
- Hat score: How far the observation is from the center of the data.

Use the **broom** function augment to add residuals and other diagnostic data to the original regression data.
See `help(influence)` for functions to get these diagnostics using base R.
```{r}
iver_mod1_aug <- augment(iver_mod1) %>%
  mutate(cty = iver$cty)
glimpse(iver_mod1_aug)
```
Oddly, `augment` calculates the *standardized residual*,
$$
\mathtt{.std.resid} = E'_i = \frac{E_i}{S_E \sqrt{1 - h_i}}
$$
which divides by the regression residual standard error, which is itself a function of the residual of $i$, $S_E = \sqrt{\frac{\sum_j E_j}{n - k - 1}}$.
What we want is the *studentized residual* which divides by the standard error of the regression calculated omitting observation $i$:
$$
\mathtt{.resid / .sigma * sqrt(1 - .hat)} = E^*_i = \frac{E_i}{S_{E_{(i)}} \sqrt{1 - h_i}}
$$
where $S_{E_(i)}$ is the standard error of the regression run without observation $i$.
It is called the Studentized residual, because it is distributed Student's $t$; the standardized residual is not.
Add a new variable called `.student.resid`, which we can calculate using the residual (`.resid`), standard error of the regression that omits that observation (`.sigma`), and the hat value (`.hat`):
```{r}
iver_mod1_aug <-
  iver_mod1_aug %>% 
  mutate(.student.resid = .resid / .sigma * sqrt(1 - .hat))
```
In base R, the function `rstudent` calculates the Studentized residuals, and `rstandard` calculates the standardized residuals:
```{r}
setNames(rstudent(iver_mod1), iver$cty)
setNames(rstandard(iver_mod1), iver$cty)
```


This scatterplot weights observations by their hat score.
Points further from the mean of `lnenp` have higher hat scores.
```{r}
ggplot(data = iver_mod1_aug, aes(x = lnenp, y = povred)) +
  geom_point(mapping = aes(size = .hat)) +
  geom_smooth(method = "lm") +
  theme_local()
```

This scatterplot weights observations by their absolute Studentized residuals.
Those observations furthest from the regression line *and* high hat values, have the highest residuals.
```{r}
ggplot(data = iver_mod1_aug, aes(x = lnenp, y = povred)) +
  geom_point(mapping = aes(size = abs(.student.resid))) +
  geom_smooth(method = "lm") +
  theme_local()
```
Cook's distance is a measure of the overall influence of points on the regression; the point's effect on *all* the parameters. 
This plot weights points by their Cook's distance. 
We can see that the two points on the bottom (Switzerland and the US) have the highest Cook's distance.
```{r}
ggplot(data = iver_mod1_aug, aes(x = lnenp, y = povred)) +
  geom_point(mapping = aes(size = .cooksd)) +
  geom_smooth(method = "lm") +
  theme_local()
```

A standard plot to assess outliers is the Influence Plot.
The x-axis is hat scores, the y-axis is Studentized residuals.
The points are sized by Cook's Distance.
Rules of thumb lines are drawn at -2 and 2 for Studentized residuals, and $\bar{h} + 2 sd(h)$ and $\bar{h} + 3 sd(h)$ for hat scores.
```{r fig.height = 6, fig.width = 6}
ggplot() +
  geom_point(data = iver_mod1_aug,
             mapping = aes(x = .hat, y = .student.resid, size = .cooksd)) +
  # add labels to points, but only those points that are flagged as outliers
  # for at least one of the diagnostics considered here
  geom_text(data =
              filter(iver_mod1_aug,
                     .cooksd > 4 / iver_mod1$df.residual
                     | abs(.student.resid) > 2
                     | .hat > mean(.hat) + 2 * sd(.hat)),
            mapping = aes(x = .hat, y = .student.resid, label = cty),
            hjust = 0, size = 4, colour = "red") +
  geom_hline(data = data.frame(yintercept = c(-2, 0, 2)),
             mapping = aes(yintercept = yintercept),
             colour = "blue", alpha = 0.4) +
  geom_vline(data = data.frame(xintercept = mean(iver_mod1_aug$.hat) +
                                 sd(iver_mod1_aug$.hat) * c(2, 3)),
             mapping = aes(xintercept = xintercept),
             colour = "blue", alpha = 0.4) +
  xlab("hat") +
  ylab("Studentized residuals") + 
  scale_size_continuous("Cook's Distance") + 
  theme_local()
```


Instead of a plot, we could find the id
Observations with high Cook's distance (greater than $4 / (n - k - 1)$):
```{r}
filter(iver_mod1_aug, .cooksd > (4 / iver_mod1$df.residual)) %>%
  select(cty, .cooksd, lnenp)
  
```
Observations with high hat scores (greater than 2 standard deviations than the mean hat score):
```{r}
filter(iver_mod1_aug, .hat > mean(.hat) + 2 * sd(.hat)) %>%
  select(cty, .hat, lnenp)
  
```
Observations with high Studentized residuals (+/- 2):
```{r}
filter(iver_mod1_aug, abs(.student.resid) > 2) %>%
  select(cty, .student.resid, lnenp)
```
Or combine these,
```{r}
filter(iver_mod1_aug,
       abs(.student.resid) > 2 |
         .hat > mean(.hat) + 2 * sd(.hat) |
         .cooksd > 4 / iver_mod1$df.residual) %>%
  select(cty, .cooksd, .hat, .student.resid, lnenp)
```

Also see `influencePlot` in **car**, and `influencePlot` in **simcf** for other implementations of this plot type.
One feature of those implementations is that they allow for the ability to identify the points on the plot.

Now that we've identified Switzerland as a problematic point, the question is what to do about it.
Checking the Switzerland data, it appears that it is correct and is not the result of data entry issues.
In general, we should avoid dropping points.
Perhaps the issue is that we have not accounted for different electoral systems. 
Let's try including 
```{r}
iver_mod2 <- lm(povred ~ lnenp + elec_sys, data = iver)
iver_mod2
iver_mod2_aug <- augment(iver_mod2) %>%
  mutate(.student.resid = .resid / (.sigma * sqrt(1 - .hat)),
         cty = iver$cty)
```

However, by including a categorical variable for electoral system in which Switzerland is the only
country with a unanamity government, we are effectively dropping Switzerland from the regression.
This means that we cannot calculate Cook's distance or studentized residuals, or hat scores
for Switzerland since a regression estimated *without* switzerland cannot estimate a coefficient
for the `unam` category, since Switzerland is the only member of that category.
```{r}
filter(iver_mod2_aug,
       abs(.student.resid) > 2 |
         .hat > mean(.hat) + 2 * sd(.hat) |
         .cooksd > 4 / iver_mod1$df.residual) %>%
  select(cty, .cooksd, .hat, .student.resid, lnenp)
```
But now that we've ignored Switzerland, both Italy and the United States seem to be influential.
This is because now that there are fewer observations per group, in some sense it is easier for observations to be influentia.
But, although the US and Italy have high studentized residuals, neither of them exceed the rule of thumb for
Cooks distance.
```{r}
filter(iver_mod2_aug,
         .cooksd > 4 / iver_mod1$df.residual) %>%
  select(cty, .cooksd, .hat, .student.resid, lnenp)
```


```{r fig.height = 6, fig.width = 6}
ggplot() +
  geom_point(data = filter(iver_mod2_aug, .cooksd < Inf),
             mapping = aes(x = .hat, y = .student.resid, size = .cooksd)) +
  # add labels to points, but only those points that are flagged as outliers
  # for at least one of the diagnostics considered here
  geom_text(data =
              filter(iver_mod2_aug,
                     .cooksd > 4 / iver_mod2$df.residual
                     | abs(.student.resid) > 2
                     | .hat > mean(.hat) + 2 * sd(.hat),
                     .cooksd < Inf),
            mapping = aes(x = .hat, y = .student.resid, label = cty),
            hjust = 0, size = 4, colour = "red") +
  geom_hline(data = data.frame(yintercept = c(-2, 0, 2)),
             mapping = aes(yintercept = yintercept),
             colour = "blue", alpha = 0.4) +
  geom_vline(data = data.frame(xintercept = mean(iver_mod2_aug$.hat) +
                                 sd(iver_mod2_aug$.hat) * c(2, 3)),
             mapping = aes(xintercept = xintercept),
             colour = "blue", alpha = 0.4) +
  xlab("hat") +
  ylab("Studentized residuals") + 
  scale_size_continuous("Cook's Distance") + 
  theme_local()
```

Although there are still a few observations with large residuals, and with a small dataset, it is almost inevitable that some observations will have outsized influence on the results, from an outlier perspective the new model seems less problematic.
However, we accomplished this at the cost of effectively ignoring Switzerland.
The model is able to estimate how different Switzerland is from what would be predicted, but by including a dummy variable
that is only 1 for Switzerland, we are treating Switzerland as *sui generis*.
Also note, that although the category is called `unam`, it would be inappropriate to interpret it as the effect of that
type of government since Switzerland is the only country in that category.
We cannot separate the effect of the government type from all the other things that make Switzerland unique.
It would be more appropriate to call it the "Switzerland" category in this instance.

## Robust and Resistant Regression Methods

Methods of dealing with outliers include robust and resistant regression methods.
Many forms of robust regression are available through the **MASS* library functions `lqs` and `rls`.
These include least median squares:
```{r}
library("MASS")
iver_lms <- lqs(povred ~ lnenp, data = iver, method = "lms")
iver_lms
```
least trimmed squares
```{r}
iver_lts <- lqs(povred ~ lnenp, data = iver, method = "lts")
iver_lts
```
M-method with Huber weighting,
```{r}
iver_huber <- rlm(povred ~ lnenp, data = iver, method = "M",
                  scale.est = "Huber")
iver_huber
```
MM-methods,
```{r}
iver_mm <- rlm(povred ~ lnenp, data = iver, method = "MM",
               scale.est = "Huber")
iver_mm
```

Now plot all of them together,
```{r fig.height=6,fig.width=6}
iver_line_compare <-
  bind_rows(data_frame(method = "OLS",
                       intercept = coef(iver_mod1)["(Intercept)"],
                       slope = coef(iver_mod1)["lnenp"]),
            data_frame(method = "LMS",
                       intercept = coef(iver_lms)["(Intercept)"],
                       slope = coef(iver_lms)["lnenp"]),
            data_frame(method = "LTS",
                       intercept = coef(iver_lts)["(Intercept)"],
                       slope = coef(iver_lts)["lnenp"]),
            data_frame(method = "Huber",
                       intercept = coef(iver_huber)["(Intercept)"],
                       slope = coef(iver_huber)["lnenp"]),
            data_frame(method = "MM",
                       intercept = coef(iver_mm)["(Intercept)"],
                       slope = coef(iver_mm)["lnenp"])
            )

print(iver_line_compare)

ggplot() +
  geom_abline(data = iver_line_compare, 
              mapping = aes(intercept = intercept, slope = slope,
                            colour = method)) +
  geom_point(data = iver, mapping = aes(x = lnenp, y = povred)) +
  scale_colour_discrete() +
  xlab("log(number of effective parties)") +
  ylab("Poverty reduction") +
  theme_local()

```

Note that these robust and resistant estimators do no include standard errors.
To get standard errors, we need to bootstrap these estimates.
The following code uses the `bootstrap` function combined `do` to generate bootstraps; see this [vignette](http://cran.r-project.org/web/packages/broom/vignettes/bootstrapping.html).
Then 
```{r}
iver %>%
  bootstrap(5000) %>%
  do({
    mod <- lqs(povred ~ lnenp, method = "lms",
               data = .)
    data.frame(term = names(coef(mod)),
               estimate = coef(mod))
  }) %>%
  group_by(term) %>%
  summarise(mean = mean(estimate),
            lb = quantile(estimate, 0.025),
            ub = quantile(estimate, 0.975))
  
  
```
We find that the standard errors are very large compared to those of the OLS.
```{r}
sqrt(diag(vcov(iver_mod1)))
```
This makes LTS not particularly useful in small datasets.

This is alternative code to calculate the bootstrap of LTS using the **boot** package.
`boot` is a more general and powerful method of bootstrapping, supporting many different sampling methods. 
However, it requires defining a function to calculate the statistic (in this case I write the function `leasttrimmed`), and does not return a data frame, but a special `boot` object.
```{r}
run_lqs <- function(d, i, ...) {
  coef(lqs(povred ~ lnenp, data = d, subset = i))
}

boot(iver[ , c("povred", "lnenp")], run_lqs, R = 1000)
```


Adapted from an example in Christopher Adolph (Spring 2014), "Outliers and Robust Regression Techniques" [lecture slides]. <http://faculty.washington.edu/cadolph/503/topic6.pw.pdf>.
