---
title: "Regression III: Lab2 Answers"
output: html_document
---

```{r setup, echo=FALSE, include=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message = FALSE, fig.align="center")
library(DAMisc)
library(dplyr)
options(java.parameters = "-Xmx5g")
```

Both of the questions in this lab use state repression as the dependent variable.  In a general sense, state repression is the violation of human rights by the state.  In this case, the focus is on the set of "physical integrity rights" - the rights to be free from torture, political imprisonment, extrajudicial killing and forced disappearance. 

## Question 1

This question uses the `q1data.rda` file (which wil put an object in your workspace called `q1data`).  You can get this file by either downloading it from the course website https://quantoid.net/teachicpsr/regression3 or by doing the following in R: 

```{r}
q1 <- file("https://quantoid.net/files/reg3/q1data.rda")
load(q1)
close(q1)
```

This file has a number of variables.  You can find a short description of each with: 


```{r}
searchVarLabels(q1data, "")
```

The dependent variable is going to be `physint` the CIRI physical integrity rights index.  All of the other variables, except for `ccode` and `Year` will be the independent variables.  

  1. Use `alsosDV` in the `DAMisc` package to figure out whether the 9-category (0-8) variable could be treated as an interval-level variabled and modeled with OLS without other modifications.  

```{r}
library(DAMisc)
mod <- alsosDV(physint ~ polity2 + pop + all_phones_pc + gdppc_mp +         revols + riots + agdems + terror_incidents, data=q1data)
```

Then, we could make a plot of the measurement function: 

```{r, fig.align="center", out.width='65%'}
plot(mod$result)
```

This looks as though we could use this in as an interval level variable.  For those of you in the scaling class, this would not be surprisnig because Cingranelli and Richards motivate the creation of the physical integrity rights index from the point of a cumulative scale and that model fits the data relatively well. 
  
  2. Once you've done that, diagnose problems with non-linearity using the conventional methods we talked about early last week (e.g., CR Plots).  Implement simple fixes to the problems if they exist. 
  
```{r, fig.align="center", fig.height=12, fig.width=12, out.width='100%'}
library(car)
mod <- lm(physint ~ polity2 + pop + all_phones_pc + gdppc_mp + 
        revols + riots + agdems + terror_incidents, data=q1data)
crPlots(model=mod, layout=c(3,3))
```

I'll make the following changes to the model and the re-evaluate: 

```{r, fig.height=12, fig.width=12, out.width='100%'}
mod2 <- lm(physint ~ poly(polity2, 2) + log(pop) + 
          log(all_phones_pc) + log(gdppc_mp) + 
        revols + riots + agdems + terror_incidents, data=q1data)
summary(mod2)
crPlots(mod2, layout=c(3,3))
```

I think the only other modification worth making might be to dichotomize revolutions (or at least recode the two high values as threes).  Looking at the CR Plot, it appears that the real decrease is from zero to 1 and then other changes result in rouchly no change in the DV.  First, if we want to look at transformations, we could use the `boxTidwell` function.  You'll note if you do this will all three of the variables that look like they could be amenable to a non-linear transformation (`pop`, `all_phones_pc`, `gdppc_mp`) you'll notice that the model fails:  

```{r, message=TRUE, include=TRUE, error=TRUE}
boxTidwell(physint ~ gdppc_mp + pop + all_phones_pc, 
           ~  poly(polity2, 2) + 
             I(revols > 0) + riots + agdems + terror_incidents,
           data=q1data)
```

With one additional argument, you can see where things go wrong. 

```{r, message=TRUE, include=TRUE, error=TRUE}
boxTidwell(physint ~ all_phones_pc + gdppc_mp + pop , 
           ~  poly(polity2, 2) + 
             I(revols > 0) + riots + agdems + terror_incidents,
           data=q1data, verbose=T)
```

The algorithm wants to send the power on `all_phones_pc` to really high values which result in: `r 70^2827`.  Taking out `all_phones_pc` highlights a different problem.  

```{r, message=TRUE, include=TRUE, error=TRUE}
boxTidwell(physint ~ gdppc_mp + pop , 
           ~  poly(polity2, 2) + all_phones_pc, 
             I(revols > 0) + riots + agdems + terror_incidents,
           data=q1data, verbose=T)
```

The `gdppc_mp` variable also gets unusually large negative power which essentially makes all values of the transformed variable zero.  A little bit of trial and error shows that `all_phones_pc` could actually be back in the model, so long as `gdppc_mp` is not. 
```{r, message=TRUE, include=TRUE, error=TRUE}
boxTidwell(physint ~ all_phones_pc +  pop , 
           ~  gdppc_mp + poly(polity2, 2) + 
             I(revols > 0) + riots + agdems + terror_incidents,
           data=q1data)
```

We could include polynomials in `polity2` and `gdppc_mp` to capture their relative assumed forms. 

```{r}
mod3 <- lm(physint ~ poly(polity2, 2) + log(pop) + 
          I(all_phones_pc^3.5) + poly(gdppc_mp, 3) + 
          I(revols>0) + riots + agdems + terror_incidents,   
          data=q1data)
summary(mod3)
```

The Clarke test would show us that the model with the transformations and polynomials is better:

```{r}
library(games)
clarke(mod2, mod3)
```

For the polynomial terms, we might want to know whether we got the right degrees of freedom. 

```{r, fig.align="center", out.width='65%'}
NKnots(physint ~ polity2 + log(pop) + 
          I(all_phones_pc^3.5) + poly(gdppc_mp, 3) + 
          I(revols>0) + riots + agdems + terror_incidents,   max.knots=6, 
          data=q1data, var="polity2", plot=T, criterion="BIC", includePoly=T)
```


```{r, fig.align="center", out.width='65%'}
library(splines)
NKnots(physint ~ gdppc_mp + log(pop) + 
          I(all_phones_pc^3.5) + bs(polity2, df=4) + 
          I(revols>0) + riots + agdems + terror_incidents,   max.knots=6, 
          data=q1data, var="gdppc_mp", plot=T, criterion="BIC", includePoly=T)
```

Note that the we use the BIC to try to identify parsimonious models.  The 2-df solution is almost the same as the 6-df solution, so we use the more parsimonious specification. This would suggest the following model: 

```{r}
mod4 <- lm(physint ~ bs(polity2, df=4) + log(pop) + 
          I(all_phones_pc^3.5) + poly(gdppc_mp, 2) + 
          I(revols>0) + riots + agdems + terror_incidents,   
          data=q1data)
summary(mod4)
```


We could plot the variable effects out with: 

```{r, fig.align="center", fig.height=12, fig.width=12, out.width='100%', out.height='200%'}
library(effects)
plot(allEffects(mod4, xlevels=25))
```

Since the ALSOS result depends on the RHS specification of the model, we could run it again just to make sure

```{r, fig.align="center", out.width='65%'}
mod4a <- alsosDV(physint ~ bs(polity2, df=4) + log(pop) + 
          I(all_phones_pc^3.5) + poly(gdppc_mp, 2) + 
          I(revols>0) + riots + agdems + terror_incidents,   
          data=q1data)
plot(mod4a$result)
```

Now, if we wanted a relatively simple test, we could use a GAM with smoothers on all of the transformed variables and test against the alternative. 

```{r, fig.align="center", fig.height=9, fig.width=9, out.width='100%'}
library(mgcv)
mod4b <- gam(physint ~ s(polity2, bs="cs") + 
           s(log(pop), bs="cs") + s(all_phones_pc, bs="cs") +
           s(gdppc_mp, bs="cs") + I(revols > 0) + riots +
           agdems + terror_incidents, data=q1data)
anova(mod4, mod4b, test="F")
plot(mod4b, pages=1)
```

This suggests that the parametric model we estimated was not quite flexible enough to fit all of the nuances of the data. 

This isn't required to answer the question, but some of you have asked about time-series cross-sectional data.  If you don't want to estimate the model with a bunch of extra dummy variables in it, you could do the within transformation of all of the variables and then estimate models on the within-transformed data to get the fixed-effects estimator. 



```{r, fig.align='center', out.width='100%'}
cc <- q1data[complete.cases(q1data[,c("physint", "polity2", "pop", "all_phones_pc", "gdppc_mp", 
            "revols", "riots", "agdems", "terror_incidents")]), ]
X <- model.matrix(~ physint + polity2 + log(pop) + all_phones_pc + gdppc_mp + 
        I(revols > 0) + agdems + terror_incidents + riots, data=cc)
Xw <- lm(X ~ as.factor(ccode), data=cc)$residuals
Xw <- as.data.frame(Xw[,-1])
Xw$year <- cc$Year  
names(Xw)[6] <- "revols"
names(Xw)[3] <- "pop"
mod4c <- gam(physint ~ s(polity2, bs="ts") + 
           s(pop, bs="ts") + s(all_phones_pc, bs="ts") +
           s(gdppc_mp, bs="ts") + revols + riots +
           agdems + terror_incidents + s(year, bs="ts"), data=Xw)
anova(mod4, mod4b, test="F")
plot(mod4c, pages=1)
```

## Question 2


We're going to continue to investigate repression in this question with a different dependent variable and potentially many other independent variables.  This question uses the `q2data.rda` file (which wil put an object in your workspace called `q2data`).  You can get this file by either downloading it from the course website https://quantoid.net/teachicpsr/regression3 or by doing the following in R: 

```{r}
q2 <- file("https://quantoid.net/files/reg3/q2data.rda")
load(q2)
close(q2)
```

There are too many variables to print out the list, but you can generate it with: 

```{r, eval=FALSE}
searchVarLabels(q2data, "")
```

I have organized the data so the country identifiers are first, the DV `fariss_repress` is next, then variables related to democracy and rights, next variables related to conflict and then finally variables related to other characteristics of the country. Before you start your investigation, I want you to take 33% of your data out and save it for later. You can do that as follows: 

```{r}
samps <- sample(1:nrow(q2data), floor(.33*nrow(q2data)), replace=FALSE)
keep <- setdiff(1:nrow(q2data), samps)
dat66 <- q2data %>% filter(1:nrow(q2data) %in% keep)
hold <- q2data %>% filter(1:nrow(q2data) %in% samps)
train33 <- sample(1:nrow(dat66), floor(.5*nrow(dat66)), replace=FALSE)
test33 <- setdiff(1:nrow(dat66), train33)
train <- dat66 %>% filter(1:nrow(dat66) %in% train33)
test <- dat66 %>% filter(1:nrow(dat66) %in% test33)
```

We'll do the preliminary investigation on half of the `dat66` object and the preliminary testing on the other half.  Only once we get to a final model we want to try will we do so on the `hold` data object. 

1. Much of the work on repression suggests that conflict increases demand for repression and democratic institutions and behaviors tend to reduce it.  Estimate a parametric model that captures this idea along with whatever controls you think might be important.

First, I'll juse the `checks` and `parcomp` variables for democracy.  I'll also use `terror_binary`, `guerrilla` and `riots` as the conflict variables.  Finally, I'll include, `gdppc_mp`,  `pop`, and `elec_prod_kwh_pc` as the other variables in the model. 

```{r}
mod <- lm(fariss_repress ~ checks + parcomp + terror_binary + guerrilla + riots + 
            gdppc_mp + pop + elec_prod_kwh_pc + Year, data=train)
summary(mod)
```


2. Try some of the machine learning algorithms (CART, random forest, XGBOOST, BART) and evaluate the PDPs and/or ice plots.  What to do those things say about the nature of the relationship between the independent variables and repression?  Is there any evidence of an interaction? 

```{r, fig.height=6, fig.width=4, out.width='50%'}
library(rpart)
cart1 <- rpart(formula(mod), data=train)
plot(cart1)
text(cart1)
cart.test <- predict(cart1, newdata=test)
cart.train <- predict(cart1, newdata=train)
cor(cart.test, test$fariss_repress)^2
```

```{r}
library(randomForest)
rf1 <- randomForest(formula(mod), data=train)
rf.test <- predict(rf1, newdata=test)
rf.train <- predict(rf1, newdata=train)
cor(rf.test, test$fariss_repress)^2
```

```{r, results='hide'}
library(xgboost)
Xm <- model.matrix(formula(mod), data=train)
Xmtest <- model.matrix(formula(mod), data=test)
y <- model.response(model.frame(formula(mod), data=train))
ytest <- model.response(model.frame(formula(mod), data=test))
dtrain <- xgb.DMatrix(data=Xm, info=list(label=y))
xgb1c <- xgb.cv(data=dtrain, nfold=5, ncross=25, nrounds=50)
xgb1 <- xgboost(data=dtrain, nrounds=50, params=list(objective="reg:linear"))
xgb.test <- predict(xgb1, newdata=Xmtest, nrounds=50)
xgb.train <- predict(xgb1, newdata=Xm, nrounds=50)
```

For the `bartMachine` analysis, is used `bartMachineCV` (which takes a reasonably long time), but it identified the following set of parameters as those that minimize the out of sample error: `k=2`, `nu=10`, `q=.75`, `num_trees=200`, so I implement that advice below. 


```{r}
library(bartMachine)
Xmdf <- as.data.frame(model.matrix(formula(mod), data=train))
Xmtestdf <- as.data.frame(model.matrix(formula(mod), data=test))
bm1 <- bartMachine(Xmdf, y, k=2, nu=10, q=.75, num_trees=200)
bm1.test <- predict(bm1, new_data=Xmtestdf)
bm1.train <- bm1$y_hat_train
```

Now that we've done all of those models, we can look at the in- and out-of-sample performance. 

```{r}
iscor <- cor(cbind(bm1.train, xgb.train, rf.train, cart.train), y)^2
oscor <- cor(cbind(bm1.test, xgb.test, rf.test, cart.test), ytest)^2
cors <- cbind( iscor, oscor)
colnames(cors) <- c("In", "Out")
rownames(cors) <- c("BART", "XGB", "RF", "CART")
cors
```

We could look to see what the important features are: 

```{r}
bm1.imp <- investigate_var_importance(bm1, type="splits", plot=F)
bm1.imp <- bm1.imp[[1]]
rf1.imp <- rf1$importance/sum(rf1$importance)
xgb.imp <- xgb.importance(model=xgb1)
cart.imp <- cart1$variable.importance/sum(cart1$variable.importance)
imports <-  xgb.imp[,c(1,4)]
names(imports)[2] <- "XGB"
imports <- cbind(imports, RF=rf1.imp[match(imports[[1]], rownames(rf1.imp))])
imports <- cbind(imports, CART=cart.imp[match(imports[[1]], names(cart.imp))])
imports <- cbind(imports, BART=bm1.imp[match(imports[[1]], names(bm1.imp))])
imports <- imports[,c(1,5,2,3,4)]
imports <- imports[order(rowMeans(imports[,-1]), decreasing = TRUE)]
print.data.frame(imports, digits=2)
```

In general, it looks like the four features at the top of the table are the most important.  We should make the various plots of interest for those variables.  First, let's look at the partial dependence plots. 

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
library(pdp)
cart.pdp <- partial(cart1, data=train, pred.var="parcomp", pred.grid=data.frame(parcomp=1:5))
rf.pdp <- partial(rf1, data=train, pred.var="parcomp", pred.grid=data.frame(parcomp=1:5))
## partial doesn't work with xgboost, so we're cooking one up on our own
# make new data that is just the original data
newdf <- lapply(1:5, function(x)Xm)
# for each one of the new datasets, replace parcomp with a different value from 1 to 5
for(i in 1:5){newdf[[i]][,"parcomp"] <- i}
# generate all of the predictions
xgb.preds1 <- lapply(newdf, function(x)predict(xgb1, data=dtrain, newdata=x))
xgb.partial <- sapply(xgb.preds1, mean)
bart.preds1 <- lapply(newdf, function(x)predict(bm1, new_data=as.data.frame(x)))
bart.partial <- sapply(bart.preds1, mean)
plot.df <- data.frame(
  preds = c(cart.pdp[,2], rf.pdp[,2], xgb.partial, bart.partial), 
  model = factor(rep(c("CART", "RF", "XGB", "BART"), each=5)), 
  parcomp = 1:5)
library(lattice)
xyplot(preds ~ parcomp, groups=model, data=plot.df, type="o", pch=16, 
       auto.key=list(space="top", columns=2, points=F, lines=T), 
       main = "PARCOMP PDP")
```

We could also do the same thing for the values of `gdppc_mp`.  Here, instead of the values `1:5`, we need to pick other values along the scale of `gdppc_mp`.  

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
s <- quantile(q2data$gdppc_mp, probs=seq(.05,.95, by=.1))
cart.pdp <- partial(cart1, data=train, pred.var="gdppc_mp", pred.grid=data.frame(gdppc_mp = s))
rf.pdp <- partial(rf1, data=train, pred.var="gdppc_mp", pred.grid=data.frame(gdppc_mp=s))
## partial doesn't work with xgboost, so we're cooking one up on our own
# make new data that is just the original data
newdf <- lapply(1:length(s), function(x)Xm)
# for each one of the new datasets, replace parcomp with a different value from 1 to 5
for(i in 1:length(s)){newdf[[i]][,"gdppc_mp"] <- s[i]}
# generate all of the predictions
xgb.preds2 <- lapply(newdf, function(x)predict(xgb1, data=dtrain, newdata=x))
xgb.partial <- sapply(xgb.preds2, mean)
bart.preds2 <- lapply(newdf, function(x)predict(bm1, new_data=as.data.frame(x)))
bart.partial <- sapply(bart.preds2, mean)
plot.df <- data.frame(
  preds = c(cart.pdp[,2], rf.pdp[,2], xgb.partial, bart.partial), 
  model = factor(rep(c("CART", "RF", "XGB", "BART"), each=length(s))), 
  gdppc_mp = s)
xyplot(preds ~ gdppc_mp, groups=model, data=plot.df, type="o", pch=16, 
       auto.key=list(space="top", columns=2, points=F, lines=T), 
       main = "GDP/capita PDP")
```

We could also make the ICE plots, which will allow us to evaluate whether or not there is much in the way of any interactions. 

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
library(RColorBrewer)
cols <- brewer.pal(5, "Set1")
library(ICEbox)
cart.ice <- ice(cart1, X=Xmdf, predictor="parcomp")
cice <- clusterICE(cart.ice, nClusters=1)
```

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
rf.ice <- ice(rf1, X=Xm, predictor="parcomp")
rice <- clusterICE(rf.ice, nClusters=5, plot_legend=T, colorvec=cols)
```


```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
xgb.ice <- ice(xgb1, X=Xm, predictor="parcomp")
xice <- clusterICE(xgb.ice, nClusters=5, plot_legend=T, colorvec=cols)
```

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
bart.ice <- ice(bm1, X=Xmdf, predictor="parcomp")
bice <- clusterICE(bart.ice, nClusters=5, plot_legend=T, colorvec=cols)
```


And, then we could do the same thing for GDP/capita. 

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
cart.ice <- ice(cart1, X=Xmdf, predictor="gdppc_mp")
cice <- clusterICE(cart.ice, nClusters=1)
```

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
rf.ice <- ice(rf1, X=Xm, predictor="gdppc_mp")
rice <- clusterICE(rf.ice, nClusters=5, plot_legend=T, colorvec=cols)
```


```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
xgb.ice <- ice(xgb1, X=Xm, predictor="gdppc_mp")
xice <- clusterICE(xgb.ice, nClusters=5, plot_legend=T, colorvec=cols)
```

```{r, fig.align="center", fig.height=6, fig.width=6, out.width='65%'}
bart.ice <- ice(bm1, X=Xmdf, predictor="gdppc_mp", num_grid_pts=10)
bice <- clusterICE(bart.ice, nClusters=5, plot_legend=T, colorvec=cols)
```


There is some sense from looking at the plots that non-linearity and interactions could be a problem.  

3. Use the `diagFun` function that we talked about in class to diagnose with your parametric model.  What, if any, problems exist. 

```{r}
source("https://quantoid.net/files/reg3/diagfun.r")
mod <- update(mod, data=dat66)
out <- diagFun(mod, data=dat66, splitSample=T, returnMods=T, partialFirst="X_res")

out$r2
print(out$coefs)
```

It looks like the `marsNL` model is probably the right mix of flexibility and explanatory parsimony.  Because there are no interactions, the `plotmo` function should prove interesting to identify the nature of the non-linearities. 

```{r, fig.align='center', fig.height=8, fig.width=8, out.width='100%'}
plotmo(out$marsNL, pt.col=1)
```

There are a few big outliers on the `checks` variable, which seem to be driving the non-linearity in that variable.  The `parcomp` variable probably should be dummied out and represented with 4 regressors.  That's something we might have done from the beginning anyway. It appears that the main variation on `guerrilla` is between 0 and non-0, so that could be dichotomized.  It looks, from the plot, like `gdppc_mp` and `pop` could both be logged.  Let's look again at what would have happened if we had specified these ahead of time. 

```{r}
mod2 <- lm(fariss_repress ~ log(checks) + as.factor(parcomp) + terror_binary + I(guerrilla > 0) + riots + 
            poly(gdppc_mp,3) + log(pop) + elec_prod_kwh_pc + Year, data=dat66)
out <- diagFun(mod2, data=dat66, splitSample=T, returnMods=T, partialFirst="X_res")
out$r2
```
Specifying the model this way indicates that there is very little left over for the non-parametric model to do.  

4. Now, put all of the possible independent variables into the machine learning algorithms above.  How do the results of these models compare to the parameteric model and to the results from the machine learning model on a subset of the variables? 

```{r}
cart.all <- rpart(fariss_repress ~ ., data=train[,-c(1,3,4)])
rf.all <- randomForest(fariss_repress ~ ., data=train[,-c(1,3,4)])
X.all <- model.matrix(fariss_repress ~ ., data=train[,-c(1,3,4)])
X.all.test <- model.matrix(fariss_repress ~ ., data=test[,-c(1,3,4)])
y.all <- model.response(model.frame(fariss_repress ~ ., data=train[,-c(1,3,4)]))
#xgb.all <- xgb.cv(data=X.all, label=y.all, watchlist=list(data=X.all.test), nfold=10, nrounds=100)
xgb.all <- xgboost(data=X.all, label=y.all, nrounds=16)
bart.all <- bartMachine(X=as.data.frame(X.all), y.all, k=2, nu=10, q=.75, num_trees=200)

bart.imp <- investigate_var_importance(bart.all, type="splits", plot=F)
bart.imp <- bart.imp[[1]]
rf.imp <- rf.all$importance/sum(rf.all$importance)
xgb.imp <- xgb.importance(model=xgb.all)
cart.imp <- cart.all$variable.importance/sum(cart.all$variable.importance)
imports <-  xgb.imp[,c(1,4)]
names(imports)[2] <- "XGB"
imports <- cbind(imports, RF=rf.imp[match(imports[[1]], rownames(rf.imp))])
imports <- cbind(imports, CART=cart.imp[match(imports[[1]], names(cart.imp))])
imports <- cbind(imports, BART=bart.imp[match(imports[[1]], names(bart.imp))])
imports <- imports[,c(1,5,2,3,4)]
imports <- imports[order(rowMeans(imports[,-1]), decreasing = TRUE)]
head(imports, n=10)
```

It seems that some of the features I neglected are important in trying to predict the outcome. Here's what the in- and out-of-sample predictive capacity look like. 

```{r}
cart.testa <- predict(cart.all, newdata=test)
cart.traina <- predict(cart.all, newdata=train)
rf.testa <- predict(rf.all, newdata=test)
rf.traina <- predict(rf.all, newdata=train)
xgb.testa <- predict(xgb.all, newdata=X.all.test, nrounds=50)
xgb.traina <- predict(xgb.all, newdata=X.all, nrounds=50)
bm1.testa <- predict(bart.all, new_data=as.data.frame(X.all.test))
bm1.traina <- bm1$y_hat_train


iscora <- cor(cbind(bm1.traina, xgb.traina, rf.traina, cart.traina), y)^2
oscora <- cor(cbind(bm1.testa, xgb.testa, rf.testa, cart.testa), ytest)^2
corsa <- cbind( iscora, oscora)
colnames(corsa) <- c("In", "Out")
rownames(corsa) <- c("BART", "XGB", "RF", "CART")
corsa
```

Compare this to the in- and out-of-sample predictions from the models with just our selected few variables in them: 

```{r}
cors
```

What we get is some additional out-of-sample predictive capacity when we cast our net a bit wider.  