---
title: "Homework 06 Solutions"
header-includes: 
- \usepackage{float} 
- \floatplacement{figure}{H} 
output: 
  bookdown::pdf_document2:
    toc: false
    keep_tex: true
---


# Problem 1

```{r message = F, warning = F}
library(tidyverse)
library(kableExtra)
library(GGally)
library(grid)
library(gridExtra)
library(ggplotify)
library(reshape2)

cdidata <- read.table("cdi.dat",header=T)
```

__Note:__ There is more detail in the solutions for problem #1 than I would expect to see in a technical appendix.  I have included extra detail so that you can see some of my thinking, and some of my approaches to problem-solving in applied regression.  When you begin to convert your revised solutions to problem #1 into a technical appendix for the project 1 IDMRAD paper, you should include technical material (and explanations) that help support the claims that you will be making in the main paper.  You know you have done the technical appendix right, if you can tell the reader after each claim you make, which page(s) of the technical appendix has more evidence, calculations, etc., to back up that claim.  If you do not refer to part of your technical appendix in the main body of your paper, you may not need to include that part in your technical appendix at all.


```{r message = F, warning = F, include=F, echo=F}

## some things you might do even before what the question is asking for...

dim(cdidata)

unlist(lapply(cdidata,class))

str(cdidata)

```

## Part a

### Make a table or tables showing appropriate summary statistics for each variable in the data set. Note that summary statistics for continuous variables will be different from the summary statistics for categorical variables.

First I'll just look at the "head" of all the variables (see Tables 1 and 2); this is a bit like looking at `str()`, which is also worth looking at (I split the "head" into two parts, because there are too many variables to fit horizontally on the page).  
```{r}
head(cdidata[,1:10]) %>% kbl(booktabs=T,caption=" ") %>% kable_classic()
head(cdidata[,c(1,11:17)]) %>% kbl(booktabs=T,caption=" ") %>% kable_classic()
```

We can also check to see how many unique values each variable has (Table 3; this is especially relevant for the categorical variables)

```{r}
apply(cdidata,2,function(x) {length(unique(x))}) %>% 
  kbl(booktabs=T,col.names="unique values",caption=" ") %>% 
  kable_classic(full_width=F)
```

It looks like `id` is just the same as the row number for each row of the data frame -- not useful for data analysis!

`state` has 48 values -- that is a lot, and so I may set this variable aside until later in my analyses

`county` is a categorical variable with nearly as many unique values (373) as rows in the cdidata data frame (440).  A little more exploration (see R code just below) shows that if I combine `county` with `state`, I get 440 unique values: some counties in different states have the same name. The upshot is that we only have one observation per *unique* county (what we might have thought were multiple observations per county are really single observations from counties with the same name in different states), and so `county` is not a useful variable to include in models (though it might be useful to identify outliers and other interesting data points in our analyses)

```{r}
county.state <- with(cdidata,paste(county,state))
tmp <- as.data.frame(matrix(sort(county.state),ncol=4))
names(tmp) <- paste("Counties",c("1-110","111-220","221-330","331-440"))
tmp[1:30,] %>% kbl(booktabs=T,longtable=T,caption=" ") %>% kable_classic(full_width=F)
```

(I have only listed the first 30 rows, to give a sense of the list of county+state combinations; you can delete the `[1:30,]` in the R code to see the whole table...)

Scanning through this table of county & state names, we see that Allen is the first county name that appears in multiple states.  And a little further down, both "Baltimore MD" and "Baltimore_City MD" are listed; this makes me wonder whether these two data points are really independent (Baltimore the city is inside Baltimore the county!), but I will ignore possible dependence like this in my analyses below.

Next let's make a table with the usual one-dimensional summary statistics (Table 5)

```{r}
cdinumeric <- cdidata[,-c(1,2,3,17)] ## get rid of id, county, state and (for now) region
apply(cdinumeric,2,function(x) c(summary(x),SD=sd(x))) %>% as.data.frame %>% t() %>% 
  round(digits=2) %>% kbl(booktabs=T,caption=" ") %>% kable_classic()
```

and a separate table for the one categorical variable that I'll be using a lot (Table 6)

```{r}

tmp <- rbind(with(cdidata,table(region)))
row.names(tmp) <- "Freq"
tmp %>% kbl(booktabs=T,caption=" ") %>% kable_classic(full_width=F)

```

I don't see anything really exciting in the continuous variables (Table 5).  There are several variables with Mean substantially larger than Median (`land.area`, `pop`, `doctors`, `hosp.beds`, `crimes`, `per.cap.income`, and `total.income`), indicating possible right-skewing.  There are no variables with Mean substantially smaller than Median.  

For the `region` variable (Table 6), it might be of some interest that the most counties are in the South (region 'S') and the least are in the west (region 'W'). The low number of counties in the West could be indicative of a lack of sampling in the West, or it could be that counties are just larger (in land area) in the West, so there are fewer counties to sample from.  Similarly, the high number of counties in the South could be indicative of over-sampling, or perhap the South simply has a lot of counties that cover only small land areas.

\clearpage

### Indicate where (in which variables) there is missing data (NA’s), if any, how much there is (in each variable) and why it might be there.

We can check for NA's directly:

```{r}
apply(cdidata,2,function(x) any(is.na(x)) )
```

There do not appear to be any missing values in the data!  (In general we might also check for old-fashioned missing value codes like "9", "99", "98", etc., but there's no evidence of that in Table 5 (look at the Min and Max values - no "9's", "99's", etc.))


### Make some appropriate descriptive EDA plots to illustrate any important features of the variables or possible important relationships among them.

We want to capture (i) univariatoe distributions [primarily historgrams, boxplots, and the like] and (ii) two-variable relationships [primarily scatter plots].

We could try `ggpairs` from `library(GGally)` to capture both (i) and (ii), but there are too many variables to make this a legible plot.  Instead, we'll do (i) and (ii) separately.

For much of the rest of this problem, I want `region` to be in the data frame as well as the numerical variables, so I will create that now:
```{r}
cdigood <- data.frame(cdinumeric,region=cdidata$region)
```
First, univariate distributions:

```{r fig.width=7,fig.height=8.5, message = F, warning = F, fig.cap="Distributions of Variables"}

## You can get almost the same thing as below with 
##
## ggplot(gather(cdinumeric), aes(value)) + 
##    geom_histogram(bins=30) + 
##    facet_wrap(~key, scales = 'free_x')

hist.builder <- function(df) {  ## creates a list of graphs
    result <- NULL 
    for (var in names(df)) {
      d <- data.frame(dd=df[,var])
      if(mode(df[,var])=="numeric") {
        p <- ggplot(d,aes(x=dd)) + geom_histogram() +
               ggtitle(var) + xlab("") 
      } else {
          p <- ggplot(d,aes(x=dd)) + geom_bar() +
                 ggtitle(var) + xlab("") 
      }
    result <- c(result,list(p))
    }
    return(result)
}

grid.arrange(grobs=hist.builder(cdigood))  ## from library(grid) & library(gridExtra)
             
```

It looks from the histograms like the variables that will really need attention (because they are severely right-skewed) are `land.area`, `pop`, `doctors`, `hosp.beds`, `crimes`, and `tot.income`, and maybe `per.cap.income`.  These are the same variables that we identified from Table 5 above. 

And now for bivariate relationships...  We really only care about two things:

* checking for high correlations among the predictors

* checking to see if the predictors are linearly related to `per.cap.income`

When we have a large number of variables, it can be a more efficient use of the page to make a heatmap of the correlation matrix, rather than looking at the correlation matrix itself:

```{r fig.cap="Heat Map of Correlatons Among Variables"}

corgraph <- function(df) {
  cormat <- cor(df)
  melted_cormat <- melt(cormat)   ## need library(reshape2) for this...
  ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
    geom_tile() + 
    theme(axis.text.x = element_text(angle = 45,vjust=0.9,hjust=1)) +
    scale_fill_gradient2(low="red",mid="white",high="blue")
}

corgraph(cdinumeric)

```
We can make the following conclusions from the correlation matrix:

* `tot.income` and `pop` are highly correlated (no surprise there)

* both are reasonably highly correlated with `crimes`, `hosp.beds` and `doctors`

* the three variables `crimes`, `hosp.beds` and `doctors` seem strongly correlated with one another

* `per.cap.income` isn't really highly correlated with anything, but the best possibilities seem to be `pct.hs.grad`, `pct.bach.deg` (postively correlated with `per.cap.income`) and `pct.below.pov`, `pct.unemp` (negatively correlated with `per.cap.income`); all four of these variables are moderately highly correlated with one another

These observations suggest that we may run into multi-collinearity problems when we start fitting models, but there is also some hope that we can make a good model for `per.cap.income`.

Now we turn to scatter plots, but we are just going to concentrate on relationships with `per.cap.income`:


```{r fig.width=7,fig.height=8.5, echo=T, fig.cap="Scatter Plots with y = per.cap.income"}

## I'm not sure how I'd do this directly with ggplot...

scatter.builder <- function(df,yvar="per.cap.income") {
    result <- NULL 
    y.index <- grep(yvar,names(df))
    for (xvar in names(df)[-y.index]) {
      d <- data.frame(xx=df[,xvar],yy=df[,y.index])
      if(mode(df[,xvar])=="numeric") {
        p <- ggplot(d,aes(x=xx,y=yy)) + geom_point() +
               ggtitle("") + xlab(xvar) + ylab(yvar) 
      } else {
          p <- ggplot(d,aes(x=xx,y=yy)) + geom_boxplot(notch=F) +
                 ggtitle("") + xlab(xvar) + ylab(yvar) 
      }
    result <- c(result,list(p))
    }
    return(result)
}

grid.arrange(grobs=scatter.builder(cdigood))

```
The best possibilities for predicting `per.cap.income` are the same variables we identified from the correlation matrix: `pct.hs.grad`, `pct.bach.deg`,`pct.below.pov`, and `pct.unemp`.  The last plot shows how `per.cap.income` varies across the four regions of the country.  There is a lot of overlap in the boxplots, but the Northeast and the West seem to be doing a little better than the North Central and South regions.

### You were not asked to do any transofrmations at this point, but it is a good idea to stop and think about transformations for a minute anyway...

I am not going to work too hard on transformations for now.  I am going to take the logarithms of the variables I identified earlier, to address heavy skewing and potential leverage and influence issues.  But I will leave the other variables alone, even though there may be some skew in them.  I can explain the logarithms to the social scientist in terms of percent-change concepts, as we have done in class.  And the more untransformed variables there are, the easier it will be for the social scientist to think about the model(s) I present.

(Another plausible approach would be to observe that essentially all of the continuous variables are skewed to some extent, and so one could just log all of them for starters.)

```{r message = F, warning = F}
cdilogs <- cdigood

skewed.vars <- c("land.area", "pop", "doctors", "hosp.beds", "crimes", "tot.income", 
                 "per.cap.income")

for (tmp in skewed.vars) {
  loc <- grep(paste("^",tmp,"$",sep=""),names(cdilogs))
    ## note special characters "^" and "$" to anchor
    ## the search pattern at the begining and end of the 
    ## string
  cdilogs[,loc] <- log(cdilogs[,loc])
  names(cdilogs)[loc] <- paste("log.",names(cdilogs)[loc],sep="")
}

## Now I get to re-use those crazy functions I wrote...
```


```{r fig.width=7,fig.height=8.5, message = F, warning = F, fig.cap="Histograms after log transformations"}
grid.arrange(grobs=hist.builder(cdilogs))
```


```{r message = F, warning = F, fig.cap="Correlations after log transformations"}
corgraph(cdilogs[,-grep("region",names(cdilogs))])
```


```{r fig.width=7,fig.height=8.5, message = F, warning = F, fig.cap="Scatterplots after log transformations"}
grid.arrange(grobs=scatter.builder(cdilogs,"log.per.cap.income"))

```

Except for `log.pop` (which I probably won't use much anyway, since `per.cap.income` = `tot.income` / `pop` is a deterministic function of `pop`), the skewing seems to have largely been brought under control.  The correlations are similar to, but a little stronger than, the correlations for the untransformed variables, and we can see from the scatter plots that the linear relationships we thought were there seem stronger, and there may be some other relationships to pay attention to.

## Part B

__Build a regression model that predicts per-capita income from crime rate and region of the country. Should there be any interactions in the model? What does your model say about the relationship between per- capita income and crime rate? Do your answers change, depending on whether you use number of crimes, or “per-capita crime” = (number of crimes)/(population) as a crime rate measure? If so, which one best answers the question? Why? Show the fitted model results and explain your answer to these questions in terms of those results.__

For each version of the income and crime variables, there are essentially three models to think about.  Since logarithms cleaned up a lot of the skewing in the data, I will only use log-transformed variables.

For log.per.cap.income and log.crimes, the three models to consider are:

```{r}
ancova.01 <- lm(log.per.cap.income ~ log.crimes,data=cdilogs)

ancova.02 <- lm(log.per.cap.income ~ log.crimes + region,data=cdilogs)

ancova.03 <- lm(log.per.cap.income ~ log.crimes * region,data=cdilogs)
```

Residual plots appear in Figure 5 below.  None of them are awful, so we could trust F-tests to compare these models:
```{r}
anova(ancova.01,ancova.02, ancova.03)
```

It looks like ancova.02 (the additive, or "parallel lines" ANCOVA model with no interactions) is doing the best.

In order to compare this with a model involving “per-capita crime” = (number of crimes)/(population), we have to construct the new variable:

```{r}
attach(cdigood)
per.cap.crime <- crimes/pop
log.per.cap.crimes <-log(per.cap.crime)
detach()
```

Note that I could equally well have calculated

```{r}
attach(cdigood)
log.per.cap.crimes <-log(crimes) - log(pop)
detach()

```

Now let's look at the same three models with this new variable:
```{r}
ancova.04 <- lm(log.per.cap.income ~ log.per.cap.crimes,data=cdilogs)

ancova.05 <- lm(log.per.cap.income ~ log.per.cap.crimes + region,data=cdilogs)

ancova.06 <- lm(log.per.cap.income ~ log.per.cap.crimes * region,data=cdilogs)
```

Again, the residual plots in Figure 5 look OK, so we try F-tests to compare these nested models

```{r}
anova(ancova.04,ancova.05, ancova.06)
```
Now ancova.05 (once again, the additive or "parallel lines" model) is doing the best.

If we want to compare these two winners (ancova.02 vs. ancova.05), we need to use AIC or BIC, because the two winners are not nested models (you can't get one from the other by imposing one or more linear constraints).

```{r}
AIC(ancova.02,ancova.05)
BIC(ancova.02,ancova.05)
```

And so it seems ancova.02 is the best model, by either AIC or BIC (remember smaller is better!).  
```{r}

```

Recalling the model formula for ancova.02 and a brief summary of the model
```{r}
formula(ancova.02)
round(coef(summary(ancova.02)),2)
```
we can interpret the model as follows (using the handout "log xform and percent interpretation.pdf" from week03):

* All across the US, for every 1\% increase in crimes, we expect a 0.07\% increase in per-capita income, om average (this increase is statistically significant, but is it practically significant?).
```{r echo=F}
fmt <- function(x) {format(round(exp(x),2),scientific=F,big.mark=",")}
```

* Different regions of the country have different baseline per-capita incomes however: In the NC region, the baseline salary is `exp(9.19) =` $`r fmt(9.19)`.  In the NE it is `exp(9.19+0.010) =` $`r fmt(9.19+0.10)`, and so forth, so in the S it is $`r fmt(9.19-0.09)`, and in the W it is $`r fmt(9.19-0.06)`.  All of these region baselines are, according to the model, significantly different from the NC baseline.

* So, according to the model, the *level* of salary varies with region in the US, but the *way it is related to crime* does not.

There is an argument that per-capita crime is more comparable to, or at least on the same scale as, per-capita income, so we will briefly look at the second-best model, anova.05, to see how it compares to ancova.02:
```{r}
formula(ancova.05)
round(coef(summary(ancova.05)),2)
```

The interpretation of this model is similar to that of ancova.02:

* All across the US, for every 1\% increase in per-capita crime, there is an associated 0.04\% increase in per-capita income (so, slightly smaller effect, but still statistically significant; should that matter?).

* The regional baseline salaries are: NC: $`r fmt(9.94)`, NE: $`r fmt(9.94+0.11)`, S: is $`r fmt(9.94-0.07)`, and W: $`r fmt(9.94-0.02)`.  All but the W region have baselines that are, according to the model, significantly different from the NC baseline.

* Again, the level of salary varies with region, but not the way it varies with crime, according to the model.

One would have to choose which of these two models one would want to feature (and interpret!) in the main body of the IDMRAD paper...

\clearpage

ASIDE: 

We could also just have compared all 6 models with AIC or BIC.  In this case, we get the same result: the additive model ancova.02 has both the lowest AIC and the lowest BIC values.
```{r}
data.frame(AIC=AIC(ancova.01,ancova.02,ancova.03,ancova.04,ancova.05,ancova.06),
BIC=BIC(ancova.01,ancova.02,ancova.03,ancova.04,ancova.05,ancova.06))[,-3] %>% 
  kbl(booktabs=T,col.names=c("df","AIC","BIC")) %>% kable_classic(full_width=F)
```

RESIDUAL PLOTS: 

And finally, here are the residual disgnostic plots that justified our earlier use of Anova F-tests to compare nested models.

```{r fig.width=7,fig.height=8.5,fig.cap="Residual diagnostics for all 6 ANCOVA models, in order: ancova.01, ancova.02, ancova.03, ancova.04, ancova.05, ancova.06"}

oldmar <- par()$mar

par(mfrow=c(6,4))

par(mar=c(2,2,2,2))

invisible(lapply(list(ancova.01,ancova.02,ancova.03,ancova.04,ancova.05,ancova.06),
                function(x) plot(x,cex.main=0.5)))
  
par(mar=oldmar)
```


## Part C

__Use methods we have discussed in class and/or methods from Sheather Chapters 5, 6 & 7 (including, as needed: transformations, interactions, variable selection, residual analysis, fit indices, etc.) to find the multiple regression model predicting per-capita income from the other variables, that makes the “best” tradeoff between the following criteria:__

* __Reflects the social science and the meaning of the variables__

* __Satisfies modeling assumptions__

* __Clearly indicated by the data__

* __Can be explained to someone who is more interested in social, economic and health factors than in mathematics and statistics.__

__No matter what you do, you are likely to be unhappy with some or all of these criteria; the better you make one criterion, the worse another is likely to get. So you will have to find a compromise or tradeoff between these criteria. Explain how you decided to make the tradeoff(s) you made.__

I am reasonably happy with the set of variables in the `cdilogs` data frame, so for now I will do model selection using just these variables.  As shown in the first part of the answer to part (a), `id` and `county` are not useful predictors.  `state` might be, but I am going to see how far I can get without `state`, and then perhaps incorprate it at the end.

Before beginning, I need to take log.pop and log.tot.income out of consideration, since per.cap.income is a deterministic function of them (so if they are included, no other predictors can possibly matter, and so I won't learn anything about what is associated with per.cap.income).
```{r}
omit <- c(grep("log.pop",names(cdilogs)),grep("log.tot.income",names(cdilogs)))
cdilogred <- cdilogs[,-omit]
```

Because `region` is categorical, and almost none of our variable selection functions are smart about dealing with the group of indicator variables associated with a categorical variable, I am going to first work *without* the `region` variable, and then I will try to include the `region` variable in some intelligent way.
```{r}
cdilogred.cont <- cdilogred[,-grep("region",names(cdilogred))]
## a data set to work with that doesn't have 'region' in it -- only the continuous variables...
names(cdilogred.cont)
## there are 10 predictor variables and the response variable log.per.cap.income...
```

ALL-SUBSETS:

I'll start with all-subsets regression. 
```{r message=F,warning=F}
library(leaps)  ## for the regsubsets() function
library(car)    ## for the subsets() plotting function
```
A little preliminary exploration convinces me of two things:

* There are too many variables for the subsets() function to be an effective way to summarize the results of regsubsets()

* The documentation and a bit of googling shows that I have some other tools to work with:
  - plot() makes a *different* graphical summary of regsubsets results that *is useful*
  - summary() of regsubsets results can be stored in a variable, and there are several named components which can be used to locate the "best" (i.e. lowest-BIC) model
  - coef() works on regsubsets results, if you can identify which model in the sequence of models from regsubsets you want
  
So... onward!
```{r}
all.subsets.01 <- regsubsets(log.per.cap.income ~ ., data=cdilogred.cont,nvmax=10)
plot(all.subsets.01)
```
In the plot, the dark squares indicate which variables are in the model that has the BIC values on the left.  The darker the squares, the better the model (remember, small is good for BIC).  Rather than read them off the plot, though, we can use the components of summary() to get them directly.

We begin by locating the model with the lowest BIC...
```{r}
all.subsets.01.summary <- summary(all.subsets.01)
names(all.subsets.01.summary)
all.subsets.01.summary$bic
min(all.subsets.01.summary$bic)
print(best.model <- which.min(all.subsets.01.summary$bic))
coef(all.subsets.01,best.model)
```
Note that the variables listed here are the same as the ones at the top of the graph we made.  

Before we interpret the coefficients, let's refit the model so we can get standard errors, etc.  To make life easy, we'll set up a temporary data frame with just the variables from the lowest-BIC model.
```{r}
all.subsets.01.summary$which[best.model,]
tmp <- cdilogred.cont[,all.subsets.01.summary$which[best.model,][-1]]
all.subsets.01.final.model <- lm(log.per.cap.income ~ .,data=tmp)
summary(all.subsets.01.final.model)$coef
```

All the predictors have coefficients significantly different from zero.  However, most of the coefficients are small, and some seem to have the wrong sign (e.g. `pct.hs.grad` and `pct.unemp`).  Before diving into what might be going on, let's at least check VIFs and residual diagnostics...
```{r}
vif(all.subsets.01.final.model)
par(mfrow=c(2,2))
plot(all.subsets.01.final.model)
```
None of the VIFs seem excessively large, and the diagnostic plots don't show much except that the QQ plot suggests both the left and the right tails are a bit longer than expected for the normal distribution.

We can also look at marginal model plots, to see if we are missing a transformation, interaction, etc.
```{r}
mmps(all.subsets.01.final.model)
```
The marginal model plots look very good -- the blue data-based curves line up well with the red model-based curves.  We don't seem to be missing any important transformations, interactions, etc.

One last thing to try is to see if interaction with `region` helps in any way.
```{r}
tmp <- cbind(tmp,region=cdilogred$region)
all.subsets.01.final.with.region <- lm(log.per.cap.income ~ .*region,data=tmp)
summary(all.subsets.01.final.with.region)
```
A reasonable rule of thumb is: if any indicator for a categorical variable seems important (e.g. a statistically significant coefficient), then keep the whole categorical variable.  If none of them seem important, then drop the variable. The same thing works for interactions with categorical variables.

*(A more formally correct thing to do would be to check to see whether to include each group of variables involving `region` with partial F tests (the anova() function), AIC, and/or BIC, in a kind of stepwise regression procedure.)*

Using the rule of thumb, we would 

* Keep: `region`, `region:pct.hs.grad`, `region:pct.below.pov`, `region:pct.unemp`, and maybe `region:pct.bach.deg`  

* Drop: `region:log.land.area`, `region:pop.18_34`, `region:log.doctors`

Thus we arrive at
```{r}
all.subsets.01.final.with.some.region <- update(all.subsets.01.final.with.region,
                                                . ~ . - region:log.land.area - 
                                                  region:pop.18_34 - region:log.doctors)
summary(all.subsets.01.final.with.some.region)
                                                
```

This seems to be working... All of the main effects and interactions that involve `region` have at least one significant term, and the $R^2$ and residual SE hardly moved at all.

Let's quickly check VIFs and diagnostics, and then compare to the best BIC model that we obtained from stepwise without the `region` term.
```{r}
vif(all.subsets.01.final.with.some.region)
par(mfrow=c(2,2))
plot(all.subsets.01.final.with.some.region)
```

As usual, adding interactions usually creates collinearities (the interactions are necessarily somewhat collinear with their main effects).  We worry less about this, because the collinearities do not appear to be causing the t-statistics and p-values in the model summary to behave strangely.

The diagnostic plots look about as good or as bad as the last set we looked at. It seems to me that in this analysis, the residual diagnostic plots are always going the toughest to feel totally satisfied about.

Now let's compare our model with some region terms to the model we got from all-subsets, but with region not included:
```{r}
anova(all.subsets.01.final.model,all.subsets.01.final.with.some.region)
AIC(all.subsets.01.final.model,all.subsets.01.final.with.some.region)
BIC(all.subsets.01.final.model,all.subsets.01.final.with.some.region)

```
The anova (F test) and AIC really like the model with the region terms in it!  On the other hand, BIC prefers the simpler model.  This often happens: BIC chooses a simpler model than other methods, because it is going for a parsimonius explanatory model, rather than a highly predictive model.

 Finally, let's examine the table of estimated coefficients again for the model with some region terms in it, and offer some interpretations...
 
```{r}
formula(all.subsets.01.final.with.some.region)
round(summary(all.subsets.01.final.with.some.region)$coef,2)
```
Here are some things we can say (this is only a few examples, there are more things to say!):

* For every 1% increase in a county's land area, there is a 0.03% decrease in expected per-capita income.  (We might conjecture that this could be due to an urban-rural contrast: rural counties tend to be bigger than urban ones).

* For every 1% increase in the number of doctors in a county, the expected per-capita income increases by about 0.06%.  That makes sense; doctors are well-paid and could be big contributors to the per-capita average income.

* For every 1 percentage point increase in the percent of the population aged 18--34, there iz an expected 2% drop in per-capita income. (We might conjecture that this is because 18--34 year olds are not at peak earning capacity yet and so perhaps their lower incomes drags down the per-capita average).

* percent of the population that are high school graduates doesn't have much effect, except in the South, where a one percentage point increase in hs graduates induces an expected 2% decrease in per-capita income.  I don't hae a good explanation for this: It might depend on whether college graduates are counted as a subset of hs graduates rather than counting them separately, or it might have something to do with some unique feature of economics in the southern region of the US.

* In the main effect for region, and in several of the interactions for region, the West shows up as deviating significantly from the North Central part of the US
 
* etc.


STEPWISE REGRESSION:

Now let's try to see how far we get with stepwise regression.  I am going to take the same strategy: first try stepwise without `region`, and then see if adding `region` helps.

```{r message=F,warning=F}
library(MASS)
```


```{r}
stepwise.base <- lm(log.per.cap.income ~ .,data=cdilogred.cont)
## try to duplicate all-subsets with BIC
step.result.01.bic <- stepAIC(stepwise.base,
                       scope=list(lower = ~ 1, upper = ~ .), 
                       k=log(dim(cdilogred.cont)[1]),
                       trace=F)

anova(all.subsets.01.final.model,step.result.01.bic)
```

We can see that stepwise regression using the BIC criterion actually found exactly the all-subsets model, which is great.

```{r}
## now try with AIC
step.result.01.aic <- stepAIC(stepwise.base,
                       scope=list(lower = ~ 1, upper = ~ .), 
                       k=2,
                       trace=F)

anova(all.subsets.01.final.model,step.result.01.aic)

```
Stepwise regression using the AIC criterion added one more predictor, namely `pop.65_plus`.  Looking at the coefficient estimates below, we can say that, although the coefficient on `pop.65_plus` is significantly different from zero, its effect on expected per-capita income appears to be quite small.  Without knowing more about the economics, it's hard to make a decision about whether to keep `pop.65_plus` or not.

```{r}
round(summary(step.result.01.aic)$coef,2)
```

Well, now let's go for broke and see what happens if we just explore 2-way interactions briefly.
```{r}
step.result.02.bic <- stepAIC(stepwise.base,scope=list(lower = ~ 1, upper = ~ .^2), 
                      k=log(dim(cdilogred.cont)[1]),                         ## BIC penalty.
                      trace=F)

step.result.02.aic <- stepAIC(stepwise.base,scope=list(lower = ~ 1, upper = ~ .^2), 
                      k=2,                                                   ## AIC penalty.
                      trace=F)

comparison <- cbind(
AIC(all.subsets.01.final.model,step.result.01.aic,step.result.01.bic,
    step.result.02.aic,step.result.02.bic),
BIC(all.subsets.01.final.model,step.result.01.aic,step.result.01.bic,
    step.result.02.aic,step.result.02.bic))
comparison <- comparison[,-3]
names(comparison) <- c("df","AIC","BIC")
comparison %>% kbl(booktabs=T) %>% kable_classic()

```
As you can see in the table, both AIC and BIC like models with some interactions.  This is a quandary: now the model is getting complicated to explain to someone.  

\newpage

Let's look at the two models with interactons and see what we can say:
```{r}
round(summary(step.result.02.bic)$coef,2)
cat("\nR2 = ",summary(step.result.02.bic)$r.squared)
cat("\nR2adj = ",summary(step.result.02.bic)$adj.r.squared)
```

```{r}
round(summary(step.result.02.aic)$coef,2)
cat("\nR2 = ",summary(step.result.02.aic)$r.squared)
cat("\nR2adj = ",summary(step.result.02.aic)$adj.r.squared)
```
In the BIC-based model, there are two interactions that appear to be statistically significant, but their practical effect is almost zero on per-capita income.

The story is similar in the AIC model, which has many more interaction terms, but only two with an effect as large as 0.01. 

Although both interactions models produced big jumps in AIC and BIC (much bigger than 10!), the improvement in $R^2$ and $R^2_{adj}$ is pretty small, for all the terms that have been added to the models. 

For these reasons I might be willing to discuss these interactions with the social scientist, but I am disinclined to include them in a final model.

If I stick with the model found by all-subsets and stepAIC with a BIC penalty, then my conclusions about adding interactions with `region` will also be the same, and I will once again be led to `all.subsets.01.final.with.some.region`, which has some interesting and mostly-interpretable structure.

LASSO:

The last variable selection method I'll consider is the lasso.

```{r message=F, warning=F}
library(glmnet)
library(arm)

```

```{r}

loc <- grep("log.per.cap.income",names(cdilogred.cont))

y <- cdilogred.cont[,loc]

X <- apply(as.matrix(cdilogred.cont[,-loc]),2,function(x) rescale(x,"full"))

Xnames <- dimnames(X)[[2]]

lasso.result <- glmnet(X,y)

plot(lasso.result,xvar="lambda",xlim=c(-9,0))

abline(h=0,lty=2)

legend('bottomright',lty=1,col=1:length(Xnames),legend=Xnames,cex=0.75)

```
There is not such an obvious place to cut the shrinkage plot, to help determine what variables should be kept in the model.

Let's see if choosing lambda by cross-valiation works:
<!---------------------------------------------------------------------------------- -->

```{r}
cv.lasso.result <- cv.glmnet(X,y)

plot(cv.lasso.result)
c(lambda.1se=cv.lasso.result$lambda.1se,lambda.min=cv.lasso.result$lambda.min)
tmp <- cbind(coef(cv.lasso.result,s=cv.lasso.result$lambda.min),
             coef(cv.lasso.result,s=cv.lasso.result$lambda.1se)
            )
dimnames(tmp)[[2]] <- c("lambda(minMSE)","lambda(minMSE+1se)")

tmp
```

It's interesting to note that the model which minimizes 10-fold cross-validation error contains all 10 predictors (or sometimes 9---log.crimes can get left out---depending on what random folds we get), while the model that is 1 SE above that is again the same model that we first saw with all-subsets regression.  

Although the magnitudes of the coefficients are a little different from what we saw before (because they arise from mimimizing RSS penalized by the L1 penalty, rather than minimizing straight RSS), the signs are the same.

Once again, if we were to settle on this model, and try to bring in `region`, chances are we would end up with the same model as `all.subsets.01.final.with.some.region`.

```{r}
round(summary(all.subsets.01.final.with.some.region)$coef,2)
cat("\nR2 = ",summary(all.subsets.01.final.with.some.region)$r.squared)
cat("\nR2adj = ",summary(all.subsets.01.final.with.some.region)$adj.r.squared)

par(mfrow=c(2,2))
plot(all.subsets.01.final.with.some.region)

```

Here are  a few of the models plus's and minus's, and some tradeoffs I considered:

* __Plus's:__ 
  - The model is fairly parsimonious, and most of the estimated coefficients have the expected sign. 
  - The model is confirmed by stepwise and lasso procedures.  
  - Those procedures also found more complex models with somewhat better fit, but improvements in fit seemed small compared to the added complexity of the model. 
  - All of the variables are either in their original scale, or the have been replaced with their logarithm.  This facilitates explaining the models to anyone who is interested in and knowledgeable about the social science & economics but less knowledgeable about technical matters.

* __Minus's:__ 
   - The coefficient on `pct.unemp` seems to go the wrong way, and the coefficient on 'pct.hs.grad` is quite small, statistically and practically (it remains in the model because there is a noticable interaction that it participates in).
   - The residual diagnostic plots are just OK.
   - The fact that stepwise regression found some well-fitting models with interactions between continuous variables suggests exploring those more complex models in the future.
   - I did not have time to explore the state variable at all.  Some of the relationship between these demographic variables and per capita income might be explainable in terms of varying economic policy from one state to the next. (If one includes `state` in the model, one could take out `region` because the two are perfectly collinear (states are entirely nested within regions).
   
Finally, it would be very useful to have additional data to compare some of the models we found.  We are using reasonable methods for variable selection, but since it is all within-sample (our entire data set is our training sample), there is ample room for overfitting noise in the data.  Some of our inferences about which variables to leave in or take out may be based on overly optimistic standard error estimates, for example. If we were able to cross-validate on some new or hold-out data, we might be able to better distinguish what is the best model, at least in terms of prediction error.

# Problem 2

This will be much quicker and will mostly consist of code and comments written by a TA from an earlier version of this course.

```{r message = F, warning = F}
beauty <- read_csv("../ProfEvaltnsBeautyPublic.csv") %>%
  dplyr::select(-profevaluation, -profnumber, -multipleclass, 
                -starts_with("class"))
```

```{r warning =F , message = F , echo = F}
beauty2 <- read_csv("../ProfEvaltnsBeautyPublic.csv") %>%
  dplyr::select( -profnumber, -multipleclass, 
                -starts_with("class"))
```


## Part a

```{r}
library(leaps)
all_subsets <- regsubsets(courseevaluation ~.,
                          data = beauty,
                          nvmax = 10, intercept = T,
                          method = "exhaustive")

summ <- summary(all_subsets)

best_model_index <- which.min(summ$bic)

dropped <- summ$outmat %>% apply(2, function(col) {mean(col == " ") == 1})
all_subset_info <- cbind(summ$outmat[,!dropped], BIC = summ$bic) 
rownames(all_subset_info) <- 1:nrow(all_subset_info)
all_subset_info %>% data.frame %>%
  rownames_to_column(var = "num variables") %>%
  knitr::kable() %>% 
  kableExtra::kable_styling(latex_options = "scale_down") %>%
  row_spec(best_model_index, background = "gray")
```

```{r}
var_names <- colnames(all_subset_info)[all_subset_info[best_model_index,] == "*"]

all_subsets_model <- lm(as.formula(paste("courseevaluation ~",
                                         paste0(var_names, collapse = " + "))),
                        data = beauty)

summary(all_subsets_model)
```

## Part b

We'll show forward regression.

```{r}
forward_step <- step(lm(courseevaluation ~1, data = beauty),
  scope  = as.formula(paste("courseevaluation ~",
                      paste0(names(beauty)[names(beauty) != "courseevaluation"], 
                             collapse = " + "))),
                          direction = "forward", trace = 0)

summary(forward_step)
```

## Part c

```{r message = F, warning = F}
library(glmnet)
cv_lasso_selection <- cv.glmnet(x = as.matrix(beauty[,names(beauty) != "courseevaluation"]),
          y = beauty$courseevaluation, alpha = 1)

```
```{r}
lasso <- glmnet(x = as.matrix(beauty[,names(beauty) != "courseevaluation"]),
                y = beauty$courseevaluation,
                alpha = 1, lambda = cv_lasso_selection$lambda.1se)


coef_lasso <- as.matrix(coef(lasso))
coef_lasso_df <- data.frame(columns = names(coef_lasso[coef_lasso[,1] != 0,]),
         beta = coef_lasso[coef_lasso[,1] != 0,]) 
rownames(coef_lasso_df) <- NULL

coef_lasso_df
```

## Part d

```{r echo = T}
#basically coding up a document term matrix function for some reason

dtm_coef <- function(names_list){
  all_words <- unique(unlist(names_list))
  matrix_info <- matrix(nrow = length(names_list), ncol = length(all_words))
  
  r_idx <- 1
  for (col_names in names_list){
    matrix_info[r_idx,] <- all_words %in% col_names
    r_idx <- r_idx + 1
  }
  
  colnames(matrix_info) <- all_words
  rownames(matrix_info) <- names(names_list)
  
  return(matrix_info)
}

variables_in_model <- list(
  "best subset" = names(coef(all_subsets_model)),
  "lasso" = coef_lasso_df$columns,
  "forward stepwise" = names(forward_step$coefficients),
  "homework 05, 2c" = c(names(beauty[,names(beauty) != "courseevaluation"]), 
                        "profevaluation", "(Intercept)")
)
```

```{r tab-feature, echo = T}
t(dtm_coef(variables_in_model)) %>%
  data.frame(check.names = F) %>% 
  rownames_to_column(var = "features") %>%
  mutate(`best subset` = ifelse(`best subset`,  
                                cell_spec(`best subset`, "latex", color = "black"),
                                cell_spec(`best subset`, "latex", color = "red")),
         `lasso` = ifelse(`lasso`,  
                                cell_spec(`lasso`, "latex", color = "black"),
                                cell_spec(`lasso`, "latex", color = "red")),
         `forward stepwise` = ifelse(`forward stepwise`,  
                                cell_spec(`forward stepwise`, "latex", color = "black"),
                                cell_spec(`forward stepwise`, "latex", color = "red")),
         `homework 05, 2c` = ifelse(`homework 05, 2c`,  
                                cell_spec(`homework 05, 2c`, "latex", color = "black"),
                                cell_spec(`homework 05, 2c`, "latex", color = "red"))) %>%
  knitr::kable(format = "latex", caption = "Features include in each model.",
               escape = F) %>%
  kableExtra::kable_styling(latex_options = c("scale_down"),)
```

The number of features (excluding the intercept) in our models increase from 6 in the best subset model, lasso with 10, forward stepwise regression selecting 14 and our full model with 30 features in homework 05. If we take out the lasso model these models are nested inside each other. **It's hard to evaluate prediction models as we left nothing to evaluate them on.** We observe that the lasso had actually higher training mse than the best subset model - but that's not to crazy because the penalty means that we're get more training error allowed. 

Following the advice of Gelman and Hill from the end of lecture 11, it might also be useful to compare tables of estimated coefficients and their SE's, to see whether the direction and size of influence of each predictor on teacher ratings is (a) plausible and (b) consistent (or inconsistent) across models.  As usual, however, since we are selecting variables and fitting the model on the same data, we can be misled by overly optimistic standard errors, overfitting to noise in the data, etc.  So this approach ultimately requires some care.

```{r echo = F}
full_model <- lm(courseevaluation~., data = beauty2)
training_mse <- rep(NA, 4)
m_idx <- 1
for (model in list(all_subsets_model, lasso, forward_step, full_model)){
  if (inherits(lasso, "glmnet")){
    newx <- as.matrix(beauty[,names(beauty) != "courseevaluation"])
    yhat <- predict(object = model, newx = newx)
  } else {
    yhat <- predict(object = model)
  }
  y <- beauty$courseevaluation
  
  training_mse[m_idx] <- mean((y-yhat)^2)
  m_idx <- m_idx + 1
}

```


