---
title: "POLS 503: Lab 3"
author: "Carolina Johnson"
date: "April 17, 2015"
---
```{r echo = FALSE, results = 'hide'}
knitr::opts_chunk$set(warning = TRUE, message = TRUE)
```

## Outline

1. Check in: Questions from Homework, other lingering confusion
2. Review: Reading in data and initial data exploration
3. Linear regression
    a. Model specification in R
    b. `lm()` object and model results
    c. (Pre)Review model interpretation
4. An aside on R fundamentals: object types
5. Returning to regression:
    d. `predict()`
    e. Plotting regression summary/expected values

```{r, message = FALSE}
library("dplyr")
library("tidyr")
library("ggplot2")
library("readr")
```


## Reading in data review:

1. Download the data [iver.csv](http://UW-POLS503.github.io/pols_503_sp15/data/iver.csv)
2. Generate summary statistics for all variables in the dataset (number of observations, min, max, median, mean, standard deviation)
3. Plot the distributions of each numeric variable

This dataset is cross sectional data on industrial democracies. Containing:

--------- ----------------------------------------------------------------
`povred`  Percent of citizens lifted out of poverty by taxes and transfers
`enp`     Natural log of effective number of parties
`lnenp`   Natural log of effective number of parties
`maj`     Majoritarian election system dummy
`pr`      Proportional representation dummy
`unam`    Unanimity government dummy (Switzerland)
--------- ----------------------------------------------------------------

Source of data and model Torben Iversen and David Soskice, 2002, ``Why do some democracies redistribute more than others?'' Harvard University.

```{r echo = FALSE, results = 'hide'}
iver <- read_csv("../data/iver.csv")
```
```{r eval=FALSE}
iver <- read_csv("iver.csv")
```

```{r}
dim(iver)
iver
head(iver)

# Summarise iver variables
iver %>%
  summarise_each(funs(mean, sd), - cty, - elec_sys)

mean <- summarise_each(iver, funs(mean), - cty, - elec_sys)
sd <- summarise_each(iver, funs(sd), - cty, - elec_sys)
min <- summarise_each(iver, funs(min), - cty, - elec_sys)
max <- summarise_each(iver, funs(max), - cty, - elec_sys)
median <- summarise_each(iver, funs(median), - cty, - elec_sys)
rbind(mean, sd, min, max, median) %>%
  mutate(statistic = c("mean", "sd","min", "max", "median"))

# using tidyr
# Method 1: make long, then summarize over each variable
iver_long <- iver %>%
  gather(variable, value, - cty, - elec_sys)

iver_summary <- iver_long %>%
  group_by(variable) %>%
  summarise(mean = mean(value),
            median = median(value),
            min = min(value),
            max = max(value),
            sd = sd(value))

# Method 2
# See what each step is doing
iver_summary2 <- iver %>%
  summarise_each(funs(mean, sd, min, max, median), - cty, - elec_sys) %>%
  gather(variable_stat, value) %>%
  separate(variable_stat, c("variable", "stat"), sep = "_") %>%
  spread(stat, value)


povred.plot <- ggplot(iver, aes(povred)) + geom_histogram()
povred.plot

enp.plot <- ggplot(iver, aes(enp)) + geom_histogram()
enp.plot

lnenp.plot <- ggplot(iver, aes(lnenp)) + geom_histogram()
lnenp.plot
```

## Regression!

The basic command for linear regression in R is `lm()`. A call to this function takes the generic form illustrated below:
```{r, eval=FALSE}
res <- lm(y ~ x1 + x2 + x3, data = your.dataframe)
```

To see a summary of the regression results use `summary()`:
```{r, eval=FALSE}
summary(res)
```


Let's run a regression using the Iverson and Soskice data.  We save the regression as an object and then print the summary of the results:

**A simple bivariate model:**
```{r}
lm_bivar <- lm(povred ~ lnenp, data = iver)
summary(lm_bivar)
```

*Challenge:*

* How do we interpret this output?
* What happens if you just type `lm.bivar`?

**A multivariate model:**
```{r}
lm_multi <- lm(povred ~ lnenp + maj + pr, data = iver)
summary(lm_multi)
```

You could also specify the same model using R's default treatment of categorical variables in the formula:
```{r}
lm_cat <- lm(povred ~ lnenp + elec_sys, data = iver)
summary(lm_cat)
```


### Aside on the formula

Note: the first argument to the function is an R formula.  Formulas appear throughout many R functions, and have some special features of their syntax, some of which are illustrated below.

In `lm`, the formula is used to generate the exact $X$ matrix that will be used to estimate the model.  To see the matrix being generated internally by `lm`, add the argument `x = TRUE` to the `lm()` call:
```{r}
lm_cat <- lm(povred ~ lnenp + elec_sys, data = iver, x = TRUE)
lm_cat$x
```

We'll look at this again in one of the more complicated model specifications below.

**A new model with multiple regressors and no constant:**
```{r}
lm_multi_noc <- lm(povred ~ -1 + lnenp + maj + pr + unam, data = iver)
summary(lm_multi_noc)
```

**A new model with multiple regressors and an interaction:**
```{r}
lm_multi_interact <- lm(povred ~ lnenp * elec_sys, data = iver)
summary(lm_multi_interact)
```

If you want to add an interaction separately from the individual variables being interacted:
```{r}
lm_multi_interact <- lm(povred ~ lnenp + elec_sys + lnenp:elec_sys, data = iver)
summary(lm_multi_interact)
```

Take a look at the help page for `formula`
```{r eval=FALSE}
?formula
```

### Second aside on formulas

To better understand what the formula is doing, let's look at the model matrix one of the more complex formulas above generates:
```{r}
lm_multi_interact <- lm(povred ~ lnenp * elec_sys, data = iver, x = TRUE)
lm_multi_interact$x
```

**A new model with multiple regressors and a transformation:**

(This transformation is just illustrating that you have the option of taking the log of a variable *inside* the formula, rather than creating a new variable prior to fitting).
```{r}
lm_multi_log <- lm(povred ~ log(enp) + elec_sys, data = iver)
summary(lm_multi_log)
```

You can also do other transformations such as taking the square of a variable (we'll talk more about this substantively later in the course):
To apply a mathematical function to a variable within the formula object, enclose it in `I()`.
```{r}
lm_multi_sq <- lm(povred ~ enp + I(enp ^ 2) + elec_sys, data = iver)
summary(lm_multi_sq)
```

But what is this `lm()` object? Is it a data frame?
```{r}
is.data.frame(lm_multi)
```
No!

## An aside on R fundamentals: Object types

R has many different types of objects, that store information in different ways and are treated differently by different functions.

Some of the most common types of objects are:

1. Vectors
2. Matrices
3. Data frames
4. Lists

(Much of this should be review if you completed the R Data Camp)

### Vectors

Vectors are collections of single values, with a length equal to the number of items in the set.

For example, here is a vector of names:
```{r}
persons <- c("Sarah", "Melina", "Jefferson", "Brad", "Ashley")
print(persons)
```

And here is a vector of numbers:
```{r}
numbers <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
print(numbers)
```
Another way to create the same vector could have been:
```{r}
number_sequence <- seq(from = 1, to = 10, by = 1)
```

You can extract an item from a vector using square brackets:
```{r}
numbers[5]
persons[2:3]
```

Single variables in our data frame are *also* vectors. This is how we can do operations on them such as calculating the mean.

```{r}
countries <- iver$cty
countries
countries[3:5]
iver$cty[3:5]
```


### Matrices

Matrices are collections of equal-length vectors in row and column arrangement.  Matrices store information in a rectangular format, so look like data frames, but are less flexible as all data must be the same type (you can't mix character and numeric data, for example).  At first when you start working in R, matrices will be in use behind the scenes more than something you work with much.

As an example, see what happens when we convert our data frame into a matrix:
```{r}
iver_matrix <- as.matrix(iver)
iver_matrix
```
Why is everything in quotation marks now?

You can index a matrix using square brackets, indicating first the row you want, then the column. This is the same way that you could directly extract specific values from a data frame

```{r}
iver_matrix[2:4, ]
```
A blank before or after the comma indicates all rows or columns, respectively
```{r}
iver_matrix[8, 4]
```

If we leave out the character vectors and convert the data frame to a matrix, then the
matrix will be numeric,
```{r}
as.matrix(select(iver, - cty, - elec_sys))
```


### Data frames

You've  seen these a lot! Now there are all kinds of tools to exploit the features of data frames, many of which you're familiar with.

For the sake of completeness, here's the equivalent "base R" way of pulling out (or indexing) a data frame to select rows or columns that meet certain criteria (you will likely see this in other code or in help files etc as you explore resources on your ow):

First, let's select all countries with a majoritarian system of government:
```{r}
iver[iver$maj == 1, ] # this selects all columns
```

Next let's subset the data to all majoritarian countries but only keep the columns povred and lnenp (for example)
```{r}
iver[iver$maj == 1, c("povred", "lnenp")]

#let's keep the country names too!
iver[iver$maj == 1, c("cty", "povred", "lnenp")]
```

Here's how you create a new data frame from scratch:
```{r eval=FALSE}
mydata <- data.frame(variable1 = a.vector,
                     variable2 = another.vector,
                     variable3 = yet another vector)
```
An alternative function to create a data frame is `data_frame` from the **dplyr** package.
It differs from `data.frame` in a few of its defaults: by default it does not convert character vectors to factors, and it does not not rename columns.
```{r eval=FALSE}
mydata <- data_frame(variable1 = a_vector,
                     variable2 = another_vector,
                     variable3 = yet_another_vector)
```

An example:
```{r}
mydata <- data_frame(somenumbers = seq(2, 16, by = 2),
                     somewords = rep(c("pounce", "bounce", "IPA", "crescent"),
                                     2),
                     a_number = 5,
                     is_silly = rep(c(TRUE, FALSE), each = 4))
mydata
```

Note, vectors must be either the same length or multiples of each other's length (shorter vectors will be repeated)

```{r eval = FALSE}
mydata <- data_frame(somenumber = seq(1, 50, by = 3),)
```


### Lists

Finally, another object type is the list, which can store different types of R objects (kind like a vector, but instead of values of a variable, they're objects we might want to save and come back to).  Lists are everywhere. Many of the functions we use (such as `lm()`) return lists.  data frames are even a special kind of list!

List elements have names. The easiest way to access an element of a list is to use the `$` and the name (just like looking at a variable in a data frame).

To find out what in a list us `str()` and/or `names()`

*Challenge:*

1. Use `names()` and `str()` to explore the contents of one of the lm objects you've created. (Look at the help file for `lm` for further details)
2. Extract and save as separate objects:
    a. The coefficients
    b. The residuals (what are the residuals?)
    c. The fitted values (what are the fitted values?)

To extract the coefficients,
```{r}
coefficients_multi <- lm_multi$coefficients
coefficients_multi
# or
coef(lm_multi)
```
To extract the residuals,
```{r}
residuals_multi <- lm_multi$residuals
residuals_multi
#or
resid(lm_multi)
```

```{r}
fitted_multi <- lm_multi$fitted.values
fitted_multi
#or
fitted(lm_multi)
```

*Challenge:*

1. What important information is missing from the `lm()` list?

To extract standard errors of the estimates:
```{r}
se_multi <- lm_multi %>% vcov() %>% diag %>% sqrt
```
This calculates the standard errors by calculating the square root of the diagonal of the variance-covariance matrix of the parameters of the model object.  `vcov()` is an example of a function that has a specific "method" for different types of objects: it knows we have an  `lm` object and acts accordingly.

## Returning to regression

### Fitted values and predictions

hypothetical data
Another way to get the fitted values is with `predict()`:
```{r}
predict(lm_multi)
```

The nice thing about predict is that it will actually let us calculate the expected values from our model for any set of real or hypothetical data with the same X variables:

Here's the general form of a call to predict, giving 95% confidence intervals:
```{r eval=FALSE}
predict(object, #lm object
        newdata, # a data frame with same x vars as data, but new values
        interval = "confidence",
        level = 0.95 #the default
)
```

Let's try this with our model.

*Challenge:*

1. What would we expect the level of poverty reduction to be for a majoritarian country with 2 parties?
2. What would we expect the level of poverty reduction to be for a PR country as it goes from 1 to 5 parties?

*hint (refer to data frame info above for how to create a new dataframe for newdata argument)*
```{r}
predict(lm_cat, newdata = data_frame(lnenp = log(2), elec_sys = "maj"),
        interval = "confidence")
predict(lm_cat, newdata = data_frame(lnenp = log(seq(1:5)), elec_sys = "maj"),
        interval = "confidence")
```
or
```{r}
xnew <-  data_frame(lnenp = log(seq(1:5)), elec_sys = "maj")
predict(lm_cat, newdata = xnew, interval="confidence")
```
or
```{r}
data_frame(lnenp = log(seq(1:5)), elec_sys = "maj") %>%
  predict(lm_cat, newdata=., interval="confidence")
```

### Plotting regression results

Plotting regression results can be even more informative. Information dense, and more intuitive than regression tables!

To plot a regression line (not just using the `lm` smoother in ****ggplot2**2**), you can either fit a line to the observed values of X and the fitted values and CIs from `predict`, or fit a line to hypothetical data to illustrate the estimated relationship (the latter can help you to have smoother confidence intervals where you have fewer observations).

Generate a range of hypothetical values for a key variable of interest:
```{r}
lnenp_hyp <- seq(min(iver$lnenp), max(iver$lnenp), 0.1)
```

Calculate expected values of `povred` for each observed level of `lnenp`, setting the covariates to an fixed level, illustrating the effect of a change in `lnenp`, all else equal (for a "typical" respondent, use the mean of the covariates you are keeping fixed).  Remember to keep variable names identical to those in the model!

In order to set variable levels to their mean (to create a line that summarizes model for all countries, not just majoritarian or PR), I've gone back to using the original lm_multi model object, with the dummy variables.

```{r}
yhyp <- data.frame(lnenp = iver$lnenp, maj = mean(iver$maj),
                   pr = mean(iver$pr)) %>%
  predict(lm_multi, newdata = . , interval = "confidence")
```

We'll use these values of y and observed values of x to plot the regression line over a scatterplot of the observed data using **ggplot2**:
```{r}
plot <- ggplot(iver, aes(x = lnenp, y = povred)) +
  geom_line(aes(x = lnenp, y = yhyp[ , 1]))
plot
```
Now add confidence intervals,
```{r}
plot <- plot +
  geom_ribbon(aes(ymin = yhyp[ , 2], ymax = yhyp[ , 3]), alpha = (1 / 3))
plot
```
Finally, add the original data with text labels,
```{r}
plot + geom_text(mapping = aes(colour = elec_sys, label = cty), size = 3) +
  theme_minimal()
```


For interest, here's a comparison with the lines you would get if you did (as we did in lab) just set `elec_sys` to `"maj"` or `"pr"`:
```{r}
yhyp_maj <- data.frame(lnenp = iver$lnenp, elec_sys = "maj") %>%
  predict(lm_cat, newdata = . , interval = "confidence")
yhyp_pr <- data.frame(lnenp = iver$lnenp, elec_sys = "pr") %>%
  predict(lm_cat, newdata = . , interval = "confidence")

plot_line_compare <- plot +
  geom_line(aes(x = lnenp, y = yhyp_maj[ , 1]), colour = "red") +
  geom_line(aes(x = lnenp, y = yhyp_pr[ , 1]), colour = "green")
plot_line_compare
```

And compare how the different lines overlay the different countries:
```{r}
plot_line_compare +
  geom_text(mapping = aes(colour = elec_sys, label = cty), size = 3) +
  theme_minimal()
```

Another way of doing this is as follows,
Create a new data frame with different values of `lnenp` for each category of `elec_sys`.
This is easily done using the function `expand.grid`.
```{r}
xhyp_cat <- expand.grid(lnenp = iver$lnenp, elec_sys = unique(iver$elec_sys))
```
Then predict `y` values for each combination of `x` using the `predict` function,
and add those columns to the `x` values using `cbind`.
```{r}
yhyp_cat <- cbind(xhyp_cat,
									predict(lm_cat, newdata = xhyp_cat, interval = "confidence"))
```

In this format, the lines can be easily plotted in ggplot using the `elec_sys`
variable in `colour` and `fill` aesthetics to draw separate shapes for each "maj",
"pr", and "unem".
```{r}
ggplot() +
	geom_line(data = yhyp_cat,
						mapping = aes(x = lnenp, y = fit, ymin = lwr, ymax = upr,
													colour = elec_sys)) +
	geom_ribbon(data = yhyp_cat, 
							aes(x = lnenp, y = fit, ymin = lwr, ymax = upr,
								  fill = elec_sys),
							alpha = 0.2) +
	geom_text(data = iver, mapping = aes(x = lnenp, y = povred, label = cty,
                                       colour = elec_sys)) +
	scale_y_continuous("povred")

```
Because different datasets are used by different geoms,
each geom had to have a `data` and `mapping` (aesthetics) argument.

```{r results='asis'}
cat(texreg::htmlreg(list(lm_cat),
                html.tag = FALSE,
                head.tag = FALSE,
                body.tag = FALSE,
                doctype = FALSE))
```

```{r results='asis'}
stargazer::stargazer(lm_cat, type = "html")
```

