---
title: "Lab 4"
author: "Sergio Garcia-Rios - Jeff Arnold"
date: "Friday, April 24, 2015"
---

For this lab we will use the replication [data]( https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/17976) from [Michael Ross](http://www.sscnet.ucla.edu/polisci/faculty/ross/)'s "The Oil Curse: How Petroleum Wealth Shapes the Development of Nations."

We will be exploring the relationship between oil dependency and democracy.

## Initial Setup


This lab will use some libraries you've seen before and we should load them now

```{r, message=FALSE, warning=FALSE}
library(dplyr)
library(ggplot2)
library(broom)
```

## Reading in the Data

Read in the Ross data

```{r, results='hide'}
rossdata <- read.csv("http://UW-POLS503.github.io/pols_503_sp15/data/ross_2012.csv",
                     stringsAsFactors = FALSE)

head (rossdata)
glimpse(rossdata)
```

This is a pretty big data-set, we do not need all the variables let's subset our data to include only the variables we will use.

### Challenge

Create a new data-set containing only:

-  `cty`
-  `year`
-  `polity` 
-  `logoil_gdpcap2000_sup_Q_1`
-  `logGDPcap`
-  `oecd`

Call this new data-frame `data`

```{r}
data<-rossdata %>% 
  tbl_df() %>% 
  select(cty, year, polity, logoil_gdpcap2000_sup_Q_1, logGDPcap, oecd)
```

Note that I am putting this new data-frame in a `tbl_df`. You don't have to do it but let's try to use as much `dplyr` as possible.


## Data Management


Some of those names are too long, we should change it to something meaningful and short. This can be done easily using `dplyr` and `rename`.

```{r}
data<-data %>%
  rename(oil=logoil_gdpcap2000_sup_Q_1, gdp=logGDPcap)
```

This data-frame is way easier to glimpse at:
```{r}
glimpse(data)

data
```

A lot of missing values here. Let's omit them and then look at a summary of the data set
```{r}
data<-na.omit(data)
data %>%
  summary()
```

Finally! we are ready to start data-analyzing..

## Sacatterplots


We are going to be exploring the relationship between Democracy level (polity) and other covariates.

Let's explore these relationships with plots.

### Challenge

Create a plot that explores the relationship between democracy level and at least another variable but try to include more than two covariates using different colors and shapes.

We begin simple... 
```{r}
ggplot(data, aes(x = gdp, y = polity)) +
  geom_point(position = position_jitter(height = .5),  size = 3) + 
  theme_bw() 

```

Unfortunately, a simple scatter plot makes it hard to detect any relationship. However, `ggplot2` makes it easy to add different colors and shapes which might help identify trends.

```{r}
ggplot(data, aes(x = gdp, y = polity, colour = oil, shape=factor(oecd))) +
  geom_point(position = position_jitter(height = .5),  size = 3) + 
  theme_bw() 
```

```{r}
ggplot(data, aes(x = gdp, y = polity, colour = factor(oecd))) +
  geom_point(aes(size=(oil)), position = position_jitter(height = .5)) + 
  theme_bw() 
```

It seems like the higher the GDP the more democratic countries are, except if you are a high oil producer or a non-OECD. Let's explore these relationships using a regression.

```{r}
model1<-lm(polity ~ oil, data=data)
summary(model1)
```
Let's now include controls for GDP per capita and OECD membership
```{r}
model2<-lm(polity ~ gdp + oil + oecd, data=data)
summary(model2)
```


-  Which has a larger impact on the level of democracy: oil dependence or OECD membership?
-  Which has a larger impact on the level of democracy: oil dependence or GDP per capital?

Recall the OECD membership clustering? Let's try an interaction

```{r}
model3<-lm(polity ~ gdp + oil*oecd, data=data)
summary(model3)
```

How would you interpret these results?



## Visualiazing Regression Results

`broom` has three main functions, all of which return data frames (not lists, numeric vectors, or other types of object). `glance` returns a data frame with a single row summary of the model:

```{r}
glance(model2)
```
`tidy` returns a data frame with a row for each coefficient estimate:
```{r}
tidy(model2)
```
`augment` returns the original data frame used in the model with additional columns for fitted values, the standard errors of those fitted values, residuals, etc.
```{r}
head(augment(model2))
````

How about a coefficient plot, roppeladder... etc.

```{r}
ggplot(tidy(model2) %>% filter(term != "(Intercept)"), 
       aes(x = term, y = estimate, 
           ymin = estimate - 2 * std.error, 
           ymax = estimate + 2 * std.error)) + 
  geom_pointrange() + 
  coord_flip()

```

We can also use `coefplot` 
```{r, warning=FALSE, eval = FALSE}
library("coefplot")
coefplot(model2, coefficients = c("oecd", "oil", "gdp"))
```


### Challenge
-  What is wrong with this plot?
-  Is it useful?
-  Why? Why not?

## Regression tables

Several packages (`stargazer`, `texreg`, `apsrtable`) are useful for creating publication type regression tables. `stargazer` and `texreg` are the most complete package. Both allow output to either LaTeX or HTML tables for many types of statistical models. We'll use *stargazer* here:

```{r results = 'asis', warning=FALSE, message=FALSE}
library(stargazer)
stargazer(model1, model2, model3, type = "html")
```



## Predicted Values


We are going to use `predict` to get predicted values. We first have to set up a `newdata` 
```{r}
xnew <- list(gdp=5.9, oil=0, oecd=1)

predict(model2, newdata=xnew, interval="confidence")

xnew2 <- list(gdp=5.9, oil=0, oecd=0)

predict(model2, newdata=xnew2, interval="confidence")

```

What is  this really doing?
```{r}
model2
names(model2)

pe2<-model2$coefficients
pe2
1*pe2[1] + 5.9*pe2[2] + 0*pe2[3] + 1*pe2[4]


1*pe2[1] + 5.9*pe2[2] + 0*pe2[3] + 0*pe2[4]

```

We can create a matrix of hypothetical data to obtain predictions for a range of values:

```{r}

# create a vector of hypothetical values of GDP per capita

gdp.hyp <-seq(4,11,by=1)

#create a matrix containing all hypothetical values, which are constant for the other covariates:


xnew <- list(gdp=gdp.hyp, oil=rep(0, length(gdp.hyp)), oecd=rep(1, length(gdp.hyp)))

xnew
```

Now we feed this new data into `predict`
```{r}
pred.res <- predict(model2, newdata=xnew, interval="confidence")

pred.res
```

### Ploting Predicted Values

To plot these predicted values we have to create a data frame containing both the predicted values generated by predict and the data used to generate those values 
```{r}
mod2_predicted <-as.data.frame(pred.res)
mod2_pred_df <- cbind(xnew, mod2_predicted)

mod2_pred_df
```

We have now have a data-frame that can easily be taken by `ggplot`
```{r}
ggplot() +
  geom_line(data = mod2_pred_df, mapping = aes(x = gdp, y = fit)) +
  geom_ribbon(data = mod2_pred_df, mapping = aes(x = gdp, ymin = lwr, ymax = upr),
              alpha = 0.2)+
  ylab("Democracy")+
  theme_bw()
```

Now the model with the interaction.

```{r}
summary(model3)

oil.hyp <- seq(0,11,by=1)

xnew1 <- list(oil=oil.hyp, gdp=rep(mean(data$gdp), length(oil.hyp)), oecd=rep(0, length(oil.hyp)))

xnew1

pred.res1 <- predict(model3, newdata=xnew1, interval="confidence")
pred.res1


xnew2 <- list(oil=oil.hyp, gdp=rep(mean(data$gdp), length(oil.hyp)), oecd=rep(1, length(oil.hyp)))

xnew2


pred.res2 <- predict(model3, newdata=xnew2, interval="confidence")
```

### Challenge

How would you construct this data-frame to be used with `ggplot`?

```{r}
mod3_predicted_1 <- as.data.frame(pred.res1)
mod3_pred_df_1  <- cbind(xnew1, mod3_predicted_1)

mod3_predicted_2 <- as.data.frame(pred.res2)
mod3_pred_df_2 <- cbind(xnew2, mod3_predicted_2)

mod3_pred_df <- bind_rows(mod3_pred_df_1,mod3_pred_df_2)

mod3_pred_df
```

Now we can ggplot it
```{r}
ggplot(mod3_pred_df, aes(x =oil , y = fit, 
                         ymin = lwr, ymax = upr)) +
  geom_line(mapping = aes(colour = factor(oecd))) +
  geom_ribbon(mapping = aes(fill = factor(oecd)), alpha = 0.7) +
  ylab("Democracy") + 
  ggtitle("Predicted values of Democracy by GDP and OECD membership")+
  scale_fill_discrete("OECD")+
  scale_colour_discrete("OECD")+
  theme_bw()
```

