---
title: "Take-Home Midterm"
date: "2022-10-02"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
    toc: yes
  html_document:
    df_print: paged
  word_document:
    toc: yes
  header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \usepackage{hyperref}
  - \usepackage{xcolor}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
  - \DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, include = TRUE, message = FALSE)
knitr::opts_chunk$set(fig.width=6, fig.height=5, fig.align="center") 
# library(knitr)
```


\def\red#1{\textcolor{red}{#1}}
\def\blue#1{\textcolor{blue}{#1}}
\let\bar\overline
\let\ul\underline

\frenchspacing


\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
Solutions
\end{center}


```{r}

library(arm)  ## includes lme4, MASS, Matrix
library(ggplot2); theme_set(theme_bw())
library(gridExtra) ## to arrange ggplots...
library(GGally) ## for ggpairs...

library(leaps) ## regsubsets(), summary(), coef()
library(car)   ## subsets(), mmps(), vif(), etc.

library(glmnet) ## for glmnet, cv.glmnet, etc.

set.seed(100000)     ## cross-validation for lasso chooses folds randomly, and 
                     ## settng the pseudo-random number generator seed manually 
                     ## like this ensure that the folds will be the same each
                     ## time I re-run this rmd file...
                    
```

# Problem 1.

Return again to the \texttt{beauty} data.  We will use the  variables in the reduced data set \texttt{beauty.red}, with the \texttt{profevaluaion} variable also removed, and whatever transformations you decided to use in HW04. 

## Problem 1(a).

Print a summary of the fitted model you obtained after completing problems 3(c) and (d) on HW04, and write a short paragraph interpreting the fitted model for a college dean who is trying to understand what factors, other than teaching quality, might affect course evaluations.

\em

First, I re-create the model that I obtained after completing 3(c) and 3(d) (your model may be different---that's ok!).

```{r}
## Setting up beauty.red as the problem requires...
beauty <- read.csv("ProfEvaltnsBeautyPublic.csv")
prof.loc <- grep("profnumber",names(beauty))
multiclass.loc <- grep("multipleclass",names(beauty))
class.locs <- grep("class",names(beauty))
profeval.loc <- grep("profevaluation",names(beauty))
beauty.red <- beauty[,-c(prof.loc,multiclass.loc,class.locs,profeval.loc)]
```


```{r}
## Removing individual beauty scores and keeping standardized average....
btystd.locs <- grep("btystd",names(beauty.red))[-c(1,8)] 
## we are keeping btystdave, which comes first in this list,
## and btystdvariance, which is eighth...
beauty.red <- beauty.red[,-btystd.locs]
```


```{r echo=FALSE,eval=FALSE}
names(beauty.red)
```


```{r}
trans.names <- names(beauty.red)[c(3:12,19:20,23)]
## indices slightly different since I already removed "provevaluation"...
```


```{r echo=FALSE,eval=FALSE}
bjpowers <- rep(NA,length(trans.names))
names(bjpowers) <- trans.names
bjpowers["age"] <- 1 ## already fairly symmetric
bjpowers["beautyf2upper"]   <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautyflowerdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautyfupperdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautym2upper"]   <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautymlowerdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautymupperdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["btystdave"] <- 1 ## some right skew but reasonably bounded between -2 and +2
bjpowers["courseevaluation"] <- 1   ## even though 2.2 is suggested, I'm leaving this 
                                    ## alone, to make the linear regression coefficients
                                    ## more interpretable
bjpowers["didevaluation"]   <- -0.5 ## I may try a log here too.
bjpowers["percentevaluating"] <- 2  ## close, perhaps more interpretable
bjpowers["profevaluation"]   <- 3   ## I may try 2 later..
bjpowers["students"]        <- -0.5 ## I may try a log here too.
bjpowers["btystdvariance"]  <- 0.5  ## 0.29 is called for but hard to talk about
```


```{r}
bjpowers <- c(1.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, -0.5, 2.0, -0.5, 0.5, 3.0)
names(bjpowers) <- trans.names
## See solutions to hw04 to see where these powers came from...
## They are not the only possible transformations; you may have come up with different ones!
```

```{r}
beauty.trans <- beauty.red
for (i in trans.names) {
  beauty.trans[,i] <- beauty.red[,i]^bjpowers[i]
  names(beauty.trans)[grep(i,names(beauty.trans))] <- paste("t",i,sep=".")
}
names(beauty.trans)[grep("t.btystdave",names(beauty.trans))] <- "btystdave" 
## didn't transform btystdave...
names(beauty.trans)[grep("t.courseevaluation",names(beauty.trans))] <- "courseevaluation" 
## didn't transform courseevaluation...
names(beauty.trans)[grep("t.age",names(beauty.trans))] <- "age" 
## didn't transform age...
```


```{r echo=FALSE,eval=FALSE}
names(beauty.trans)
```


```{r}
lm.1 <- lm(courseevaluation ~ ., data=beauty.trans)

lm.2 <- update(lm.1, . ~ . - t.beautyf2upper - t.beautyflowerdiv - t.beautyfupperdiv -
                 t.beautym2upper - t.beautymlowerdiv - t.beautymupperdiv - 
                 t.profevaluation - t.percentevaluating)
```

And finally, here is a summary of my model after 3(c) and 3(d) from hw04 [your model may be different!]:

```{r}
summary(lm.2)
```


```{r eval=FALSE,echo=FALSE}
tab <- cbind(summary(lm.2)$coef, vif=c(NA,vif(lm.2)))

round(tab[order(tab[,4]),],2)

bjpowers

```

\bf
I did not write a short paragraph here, in order to show in some detail the reasoning I might use.  A short paragraph that summarizes your overall conclusions would be fine here.

\rm

\em
The estimated model here is 
\rm
\begin{eqnarray*}
   \mbox{(courseevaluation)} & = & (4.18) +  (0.53)\mbox{(onecredit)} + (-0.18)\mbox{(female)} + (0.20)\mbox{(blkandwhite)} + \\&& (2.66)\frac{1}{\sqrt{\mbox{(students)}}} + (-0.31)\mbox{(nonenglish)} + (0.09)\mbox{(btystdave)} + \\&& (0.20)\mbox{(fulldept)} +  (-0.01)\mbox{(age)} + (0.14)\mbox{(formal)} + (-0.16)\mbox{(minority)} + \\&& (-1.54)\frac{1}{\sqrt{\mbox{(didevaluation)}}} + (-0.15)\mbox{(tenuretrack)} +  (0.07)\mbox{(tenured)} + \\&& (-0.04)\sqrt{\mbox{(t.btystdvariance)}} + (0.02)\mbox{(lower)}
\end{eqnarray*}
\em
where I have organized the predictors in decreasing order of statistical importance.  

This says, according to the fitted model, that before we know anything, each course "starts with" a baseline evaluation of 4.18 points, on the 5-point course evaluation scale; this is near the median course rating for the college (which is reassuring!), as shown in the left panel below.

```{r echo=FALSE,fig.height=3.5}
par(mfrow=c(1,2))
cat("Summary Statistics for Course Evaluation:\n")
summary(beauty.red$courseevaluation)
hist(beauty.red$courseevaluation,main="",xlab="Course Evaluation")

cat("Summary Statistics for Standardized Average Beauty Score:\n")
summary(beauty.red$btystdave)
hist(beauty.red$btystdave,main="",xlab="Standardized Average Beauty Score")
```
From this baseline score of 4.18, the course evaluation 
\begin{itemize}\itemsep=0pt\parskip=4pt

\item Increases by just over half a point (0.53) on average, for courses that are only one credit (the ``fun'' or ``easy'' courses?).

\item Decreases by roughly a fifth of a point ($-$0.18) on average, if the instructor is female.

\item Increases by a fifth of a point if \texttt{blkandwhite=1} (but we are not sure what \texttt{blkandwhite} is\ldots).

\item Increases as the number of students (course size) in the course decreases.  The relationship here is a little hard to give a simple interpretation to, but as an example, if the course size goes from 50 to 25, the average increase in course evaluation is about `r round(2.66*(1/sqrt(25) - 1/sqrt(50)),2)` points. (I would advise refitting the model with untransformed \texttt{students} (course size) variable, to get a better handle on this.).

\item Decreases by about 1/3 of a point on average ($-$0.31) if the instructor is a non-native English speaker.

\item Increases by about 1/10 of a point on average (0.09) for each 1-point increase in standardized average beauty rating of the instructor by students.  In the right panel above we see that 1 point would be quite a large jump in beauty rating, so even though this variable is statistically important, beauty rating does not seem to have a great practical effect on course rating.

\item Increases by 1/5 of a point on average if everyone in the instructor's department has their pictures on the web.

\item Decreases by 1/100 point on average for each additional year older the instructor is. Thus, even though it is statistically important, instructor age doesn't seem to have much effect on course rating.

\item Increases by about 0.14 if the instructor wears formal attire in their picture on the web.

\item Decreases by about 0.16 on average if the instructor is a member of a minority.

\end{itemize}

The remaining variables, \verb|didevaluation| (number of students in the class who evaluated the course) \verb|tenuretrack| (is the instructor in the tenure track, or just hired to teach), \verb|tenured| (is the instructor tenured), \verb|btystdvariance| (a measure of how disparate the student ratings of instructor beauty were) and \verb|lower| (whether this is a freshman/sophomore course) were not statistically important, and could be ignored.  

One exception to this advice might be \verb|tenured| and \verb|tenuretrack|.  If the instructor is tenured, then they are also in the tenure track. i.e.,  \verb|tenured=1| implies that \verb|tenuretrack=1|, as we can see from this table:
```{r}
with(beauty.red,table(tenured,tenuretrack))
```
Therefore it is likely that each of these variables are sapping some effect from the other one, due to this partial collinearity.  We could decide to remove one and keep the other, or we could recode the variables as follows:
```{r eval=FALSE}
tenured <- tenured  ## i.e., leave "tenured" alone...
ttrack.not.tenured <- tenuretrack - tenured
```
in order to explore whether instructors working toward tenure tend to get higher or lower course ratings than those that already have tenure.

Some overall conclusions that we could make include:
\begin{itemize}\itemsep=0pt\parskip=4pt

\item One-unit courses get a big boost in course rating.

\item Beauty rating and age don't seem to matter very much by themselves in the course rating.

\item Nevertheless, belonging to a department with all instructors' pictures on the web, and dressing well for one's own picture, has a positive effect on course ratings.

\item Smaller classes tend to get better course ratings.

\item Being female, a member of a minority, or a non-native speaker of English all have significant negative effects on course ratings, on average.

\end{itemize}
\rm

## Problem 1(b).

Using the version of \texttt{beauty.red} that you used to obtain the model in part (a) [including any transformations that you applied for HW04], apply the lasso to select variables for predicting \texttt{courseevaluation} (or a transformation of \texttt{courseevaluation} if that's what you used in HW04) for this data set. 
\begin{itemize}
\item Is it feasible to use shrinkage plots for the lasso, as in lectures 08 and 09?  If so, try it.  If not, explain why not.
\item The function \texttt{cv.glmnet} in \texttt{library(glmnet)} tries to find an optimal $\lambda$ by cross-validation using mean-squared prediction error.  Read the documentation and try variable selection using \texttt{cv.glmnet}.  Note that \texttt{cv.glmnet} produces both \texttt{lambda.min} (the best value found by cross-validation) and \texttt{lambda.1se} (the value of $\lambda$ that is one SE larger than \texttt{lambda.min}, which many people use to protect against capitalization on chance). 
\item You can compare the results of the two values of $\lambda$ with code like this:
\begin{verbatim}
result <- cv.glmnet(x,y)
plot(result)
c(lambda.1se=result$lambda.1se,lambda.min=result$lambda.min)
cbind(coef(result),coef(result,s=result$lambda.1se),coef(result,s=result$lambda.min))
\end{verbatim}
\end{itemize}

\em


First, we try glmnet's shrinkage plot, using the transformed variables (I have put them all in \verb|beauty.trans|:
```{r fig.height=4.5}
y <- beauty.trans$courseevaluation
X <- as.matrix(beauty.trans[,-grep("courseevaluation",names(beauty.trans))])
lasso.fits <- glmnet(X,y,alpha=1)
plot(lasso.fits,xvar="lambda",label=TRUE)
data.frame("curve number"=1:dim(X)[2],"variable name"=dimnames(X)[[2]])  
```
Since there are 22 variables, it is rather hard to figure out which curves go with which variables.  Thus the shrinkage plot is not all that useful.

Now we will try selecting $\lambda$ by cross-validation:
```{r}
result <- cv.glmnet(X,y)
plot(result)
## shows estimated cross-valudation MSE at various log-lambda values,
## along with estimated CI for each MSE 

round(c(lambda.1se=result$lambda.1se,lambda.min=result$lambda.min),4)
round(log(c(log.lambda.1se=result$lambda.1se,log.lambda.min=result$lambda.min)),4)

cbind(coef(result,s=result$lambda.1se),coef(result,s=result$lambda.min))
```

Notice that many more variables are selected for \verb|lambda.min| than for \verb|lambda.1se|.  This suggests overfitting (capitalization on chance) and we might prefer the simpler model suggested for \verb|lambda.1se|.  

Let's look at the \verb|lambda.1se| model.  The variables kept in the model are: \verb|minority|, \verb|t.beautyfupperdiv|, \verb|female|, \verb|fulldept|, \verb|nonenglish|, \verb|onecredit|, \verb|t.percentevaluating|, and \verb|blkandwhite|.

If we refit the model with these variables on the full data set we get

```{r}
lm.3 <- lm(courseevaluation ~ minority + t.beautyfupperdiv + female + 
             fulldept + nonenglish + onecredit + t.percentevaluating + blkandwhite,
           data=beauty.trans)
summary(lm.3)
```
(Note: the results of \verb|cv.glmnet| can change from run to run, because \verb|cv.glmnet| is choosing the folds for cross-validation randomly.)

\rm

## Problem 1(c).

Compare the model in part (a) with the model in part (b): Make a table showing which variables remain in the final model for each of two models, and then write a brief paragraph saying which model you would use to help the college dean understand factors other than teaching quality that affect course evaluations, based on this table and any other evidence that seems relevant, and explain your reasoning. 

\em

```{r}

part.a.names <- rownames(summary(lm.2)$coef)  ## lm.2 is from part (a)

part.b.names <- rownames(summary(lm.3)$coef)  ## lm.3 is from part (b)

all.vars <- union(part.a.names,part.b.names)

tab <- matrix(NA,ncol=4,nrow=length(all.vars))
dimnames(tab) <- list(all.vars,c("lm.2 Est", "lm.2 pval", "     lm.3 Est", "lm.3 pval"))
tab[part.a.names,1] <- summary(lm.2)$coef[,1]
tab[part.a.names,2] <- summary(lm.2)$coef[,4]
tab[part.b.names,3] <- summary(lm.3)$coef[,1]
tab[part.b.names,4] <- summary(lm.3)$coef[,4]

round(tab,2)

```

In the table above, the column of variables on the left lists all the variables in either model.  The next two columns show, if the variable is a predictor in lm.2 (from part a), its estimated coefficeint and pvalue for testing whether it is significantly different from zero.  An NA is listed if that variable is not in lm.2.  The last two columns provide similar information for lm.3 (from part b).

Except for the fact that lm.3 has an individual student's beauty scores, whereas lm.2 has the average beauty score, the variables in lm.3 are a subset of the variables in lm.2 (even if we exclude the variables with nonsignificant p-values).  There is a good argument for keeping the model with more variables, since that gives you more to talk about with the dean (and if you chose a model with more variables for that reason, that would be fine). However, in this case I like the smaller model since it focuses on important variables, but I would modify it in two ways:  (1) replace the individual beauty rating with the \verb|btystdave| and (2) put \verb|tenured| and \verb|tenuretrack| back in the model so that we could have the discussion about tenured, vs tenure-track but not tenured, vs instructor-only. 

(I am not sure what I would do with the formal and fulldept variables, which seem to function differently in the two models.)

\rm



# Problem 2.

The file \texttt{cdi.dat}, in the same Canvas folder as this midterm assignment sheet, is taken from Kutner et al. (2005)\footnote{Kutner, M.H., Nachsheim, C.J., Neter, J. \& Li, W. (2005) \textit{Applied Linear Statistical Models, Fifth Edition}.  NY: McGraw-Hill/Irwin.}: It provides selected county demographic information (CDI) for 440 of the most populous counties in the United States. Each line of the data set has an identification number with a county name and state abbreviation and provides information on 14 variables for a single county. Counties with missing data were deleted from the data set. The information generally pertains to the years 1990 and 1992.  The definitions of the variables are given in Table~\ref{t:variables}.

\begin{table}[hb]
\caption{Variable definitions for CDI data from Kutner et al. (2005).  \textit{Original source:}  Geospatial and Statistical Data Center, University of Virginia.}\label{t:variables}
\begin{center}
\sf\small
\let\mc\multicolumn
\begin{tabular}{rp{1.5in}p{4in}}
\hline
Variable\\ Number &  \mc{1}{c}{Variable Name} & \mc{1}{c}{Description} \\
\hline
1 & Identification number             &  1--440\\
2	&	County														&	County name\\
3	&	State															&	Two-letter state abbreviation\\
4	&	Land area													&	Land area (square miles)\\
5	&	Total population									&	Estimated 1990 population\\
6	&	Percent of population aged 18--34 & Percent of 1990 CDI population aged 18--34\\
7	&	Percent of population 65 or older	&	Percent of 1990 CDI population aged 65 or old\\
8	&	Number of active physicians				&	Number of professionally active nonfederal physicians during 1990\\
9	& Number of hospital beds						&	Total number of beds, cribs, and bassinets during 1990\\
10&	Total serious crimes							&	Total number of serious crimes in 1990, including murder, rape, robbery, aggravated assault, burglary, larceny-theft, and motor vehicle theft, as reported by law enforcement agencies\\
11&	Percent high school graduates			&	Percent of adult population (persons 25 years old or older) who completed 12 or more years of school\\
12&	Percent bachelor's degrees				&	Percent of adult population (persons 25 years old or older) with bachelor's degree\\
13&	Percent below poverty level				&	Percent of 1990 CDI population with income below poverty level\\
14&	Percent unemployment							&	Percent of 1990 CDI population that is unemployed\\
15&	Per capita income 								&	Per-capita income (i.e. average income per person) of 1990 CDI population (in dollars)\\
16&	Total personal income							&	Total personal income of 1990 CDI population (in millions of dollars)\\
17&	Geographic region									&	Geographic region classification used by the US Bureau of the Census, NE (northeast region of the US), NC (north-central region of the US), S (southern region of the US), and W (Western region of the US)\\
\hline
\end{tabular}
\end{center}
\end{table}

\clearpage

## Problem 2(a).

Data description.
\begin{itemize}
\item Make a table or tables showing appropriate summary statistics for each variable in the data set.  Note that summary statistics for continuous variables should be different from the summary statistics for categorical variables.
\item 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.   
\item Make some appropriate descriptive EDA plots to illustrate any important features of the variables or possible important relationships among them.
\end{itemize}

\em

```{r}
cdi <- read.table("cdi.dat",header=T)

apply(cdi,2,function(x) any(is.na(x)))
```
As we can see above, there is no missing data in the data set, as promised in the problem statement.

I will not summarize the variables \verb|id| and \verb|county|:
\begin{itemize}
\item \verb|id| is just the row number of each observation in the data set and is not interesting.
\item \verb|county| is nearly unique for each row (the combination "county, state" is unique for each row; some county names occur in two or three different states), and so no summary would be helpful.
\end{itemize}
The variables \verb|state| and \verb|region| are categorical and I will just make a table of counts for each:
```{r}
table(cdi$state)
length(unique(cdi$state))
table(cdi$region)
```
Whereas the data is fairly evenly spread across regions (except that West is a little light), there is a lot of variation in the number of counties included from each state, and two states appear to have zero counties (they are not included in the data at all).  The \verb|state| variable causes problems in some regression models, and is difficult to incorporate into \verb|regsubsets| and \verb|glmnet| variable selection and so I will leave it out of regression analyses below.

The other variables are continuous and I will provide 5-number summaries, means and SD's to sumamrize them:
```{r}
summary.with.sd <- function(x) {
  su <- summary(x)
  return(c(su[1:4],SD=sd(x),su[5:6]))
}
old.width <- options()$width
options(width=200)
round(t(apply(cdi[,-c(1:3,17)],2,summary.with.sd)),2)
options(width=old.width)
```

It seems valuable to make a scatter plot matrix with density plots for all the data except for \verb|id| and \verb|county| (I will also skip \verb|state| since it has too many levels for \verb|ggpairs| to deal with):

```{r}

ggpairs(cdi[,-c(1:3)], upper = list(continuous = wrap("cor", size = 2))) +
  theme(text=element_text(size=6))  ## theme affects the labels at the sides...

```
We can see (either from the plot or from the numerical summaries) that several variables have substantial right-skew: \verb|land.area|, \verb|pop|, \verb|doctors|, \verb|hosp.beds|, \verb|crimes|, \verb|per.cap.income|, and \verb|tot.income|.  Some other variables have some right-skew as well but it doesn't seem as severe.  

\bf

What I have done above is enough to answer the question.  Below I explore transformations just a little bit more, but this is not required to answer the question.

\rm\em

From the plots below it appears that applying log(x) to these variables cleans them up about as much as we could expect, though some variables still have some severe outliers.

```{r fig.height=4}
skew.candidates <- c("land.area", "pop", "doctors", "hosp.beds", "crimes", 
                     "per.cap.income", "tot.income")
boxplot(cdi[,skew.candidates],horizontal=TRUE,las=1,cex.axis=0.5)
boxplot(log(cdi[,skew.candidates]),horizontal=TRUE,las=1,cex.axis=0.5)
```



\rm

## Problem 2(b).

Build a regression model that predicts per-capita income from crimes and region of the country (using ony these three variables, not the full set of variables in the data set).  Should there be any interactions in the model?  What does your model say about the relationship between per-capita income and crimes?  Do your answers change, depending on whether you use number of crimes, or "per-capita crime" = (number of crimes)/(population) as a measure of crimes?  If so, which one best describes the relationship between per-capita income and crimes?  Why?  Show the fitted model results and explain your answer to these questions in terms of those results (as well as any economics knowledge you may have).

\em

Below are the additive and interactive models predicting per.cap.income from crime and region. We can see from the F test (anova table), from AIC, or from BIC, that the additive model, \verb|lm.1| is preferred.  According to \verb|summary(lm.1)|, the NE (Northeast) region of the country has higher income than the North Central (NC, the missing/baseline category) region and the other regions don't have significantly different per.cap.income, on average, from the NC region.  Raw crime counts seem to have a strong association with per.cap.income.  

However, the residual diagnostic plots show that lm.1 has some extreme outliers, at least one of which is strongly influential according to Cook's Distance, that may be affecting these results.

```{r}
lm.1 <- lm(per.cap.income ~ crimes + region, data=cdi)
lm.2 <- lm(per.cap.income ~ crimes * region, data=cdi)

anova(lm.1,lm.2)

cbind(AIC(lm.1,lm.2),BIC(lm.1,lm.2))

summary(lm.1)
par(mfrow=c(2,2))
plot(lm.1)
```

And here are the additive and interactive models predicting per.cap.income from per.cap.crime and region.  Once again, the additive model "wins". From the \verb|summary(lm.3)| table, we see that the story for regions is the same---NE seems to have substantially higher incomes than the rest---but the story for per.cap.crime is now different: per-capita crime does not seem to be significantly associated with per-capita income.  

The residual diagnostic plots still show some high outliers, but only one observation, \#6, that is close to being strongly influential on the fit.  The fit is not yet perfect (perhaps transformations would help) but the story is interestingly different from what we got when we used crime count alone.
```{r}
per.cap.crime <- with(cdi,crimes/pop)

lm.3 <- lm(per.cap.income ~ per.cap.crime + region, data=cdi)
lm.4 <- lm(per.cap.income ~ per.cap.crime * region, data=cdi)

anova(lm.3,lm.4)

cbind(AIC(lm.3,lm.4),BIC(lm.3,lm.4))

summary(lm.3)
par(mfrow=c(2,2))
plot(lm.3)

```
Finally, to answer the question about which model best describes the relationship between crime and per.cap.income, actually both models say there is very little relationship:
\begin{itemize}
\item In \verb|lm.1|, the variable \verb|crimes| is a significant predictor, but the coefficient is 0.0089, which means for 100 additional crimes, per.cap.income would only be expected to increase by \$0.89 (89 cents).  So although \verb|crimes| is statistically significant, it doesn't appear to be practically significant.

\item In \verb|lm.3|, the variable \verb|per.cap.crime| is not even close to being a significant predictor.  

\end{itemize}
It is also worth noting that the $R^2$'s in both \verb|lm.1| and \verb|lm.3| are pitifully small, so neither crime (however it is measured) nor region is doing a very good job of accounting for the variation in per-capita income.


\rm

## Problem 2(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 diagnostic plots, fit indices, added-variable or marginal model plots, etc.) to find the multiple regression model predicting per-capita income from the other variables, that makes the ``best'' tradeoff between the following criteria:
\begin{itemize}
\item Reflects the social science and the meaning of the variables
\item Satisfies modeling assumptions
\item Clearly indicated by the data
\item Can be explained to someone who is more interested in social, economic and health factors than in mathematics and
statistics.
\end{itemize}
Organize your work so that it tells an interesting data analysis story that a statistican (like me!) might like, and be sure to explain why and how you arrived at a final model.

(\textbf{Note:} 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.)

\em

Following an approach we have used before, I will 
\begin{itemize}
\item First consider transformations to pull in outliers and symmetrize distributions somewhat, in order to reduce the likelihood of high leverage or influentual points;

\item Next try some general variable selection, beginning with a "kitchen sink" model, and proceeding to other methods.

\item Look at diagnostic plots to see if improvements in distribution or functional form are needed;

\item Think about the practical meanings of the variables and the model(s) to see if any further adjustments are needed.

\end{itemize}

### Transformations:

In general I am mentally biased against a lot of transformations, and also biased against complicated searches for really good transformations.

Since I found a few log-transformations "by eye" that seemed to help somewhat when I was looking at summaries of the data, I will try those transformations here.

Recall that, above, I defined 
```{r eval=FALSE}
skew.candidates <- c("land.area", "pop", "doctors", "hosp.beds", "crimes", 
                     "per.cap.income", "tot.income")
```

and found that log-transformations seemed to provide some relief from extreme right-skewing in these variables.

So...
```{r}
cdi.trans <- cdi
cdi.trans[,skew.candidates] <-  log(cdi[,skew.candidates])
skew.cand.locs <- names(cdi.trans) %in% skew.candidates
names(cdi.trans)[skew.cand.locs] <- paste0("log.",skew.candidates)
cdi.trans <- cdi.trans[,-grep("id",names(cdi.trans))]
cdi.trans <- cdi.trans[,-grep("county",names(cdi.trans))]
cdi.trans <- cdi.trans[,-grep("state",names(cdi.trans))]

```

### Variable Selection: The Initial Kitchen Sink Model

Let's start with the "kitchen sink" model...
```{r}
cdi.0 <- lm(log.per.cap.income ~ ., data=cdi.trans)
summary(cdi.0)
```
This is weird; look at some of these numbers: $R^2$ and $R_{adj}^2$ both equal 1, residual standard error is $s=0.0001072$, and only two variables matter, \verb|log.pop| and \verb|log.tot.income|, and their $\hat\beta$'s are essentially $-1$ and 1 with tiny SE's.  All the others are not even close to being significant so we can guess that more formal variable selection would lead to a fitted model like this (with essentially no error):
\[
      \mbox{(log.per.cap.income)} = 1.382 - \mbox{(log.pop)} + \mbox{(log.tot.income)}
\]
What the heck is going on???

If we exponentiate both sides, we get
\[
\mbox{(per-capita income in each county)} = (3.98) \cdot\frac{\mbox{(total income in the county)}}{\mbox{(population of the county)}}
\]
The very definition of per-capita income is total income over population size, (and the factor of 4 ($\approx 3.98$) is due to the fact that not all the units are the same: per-capita income is in dollars, total income is in millions of dollars, etc.).

So this is not a very useful model for the social scientist.  It just says we rediscovered the formula for per-capita income!  To proceed we should take one or both of the variables \verb|log.pop| and \verb|log.tot.income| out of the model. 

I looked at summaries of the regressions taking out only \verb|log.pop|, only \verb|log.tot.income|, or both.  There isn't a purely mathematical reason for preferring one or the other of these three models, but I did try to guess what might make sense.  For example, when I took out \verb|log.tot.income| but left in \verb|log.pop|,  then \verb|log.crimes| had a significant positive coefficient.  Having income go up with crimes doesn't really make sense (and doesn't agree with part (b) of this problem! This seems to reflect a collinearity between total number of crimes and total population and not give us a really interpretable model. 

Finally I decided to take them both out, since they are "encoded" in the response variable \verb|log.per.cap.income| anyway.  Since I'm not going to use them anymore, I also took them out of the data frame, to make further variable selection, etc., easier.  
```{r}
cdi.trans <- cdi.trans[,!(names(cdi.trans) %in% c("log.pop","log.tot.income"))]

```
So my initial model will be:
```{r}
cdi.1 <- lm(log.per.cap.income ~ ., data=cdi.trans)
summary(cdi.1)
vif(cdi.1)
par(mfrow=c(2,2))
plot(cdi.1)
```
The residual diagnostic plots do not look terrible, and, overall, this is not an awful model to start with.  Some coefficients make sense, and those that don't (e.g. per-capita income goes down as percent high school grads goes up, or goes up as percent unemployment goes up) could be due to some remaining collinearity.  Also, all of the coefficients seem quite small, although the $R^2$ is still a decent 85\%.  All these are reasons to try some additional variable selection to see if we can improve the model (i.e., keep a high $R^2$ but have coefficients that are larger and still significant, in a model that satisfies the assumptions of linear regression and is interpretable for the social scientist).

### Variable Selection: Other Methods

Since there are not many predictors I can try all-subsets variable selection, and I will also try a couple others: 

\begin{itemize}

\item All subsets, with AIC and BIC

\item Stepwise, with AIC and BIC

\item lasso, using crossvalidation to choose $\lambda$, as in problem \#1.

\end{itemize}

#### All subsets:

I will do all-subsets variable selection without the categorical \verb|region| variable, since \verb|regsubsets()| doesn't treat categorical variables correctly, and then try to add \verb|region| back in to see if it helps.
```{r}
pmax <- dim(cdi.trans)[2] -  2 ## for log.per.cap.income and region 
all.subsets <- regsubsets(log.per.cap.income ~ . - region, data=cdi.trans, nvmax=pmax)
```


```{r}
tmp <- summary(all.subsets)
p <- 1:dim(tmp$which)[1]
n <- dim(cdi.trans)[1]

attach(tmp)

results <- data.frame(which, BIC=bic, AIC=n*log(rss) + 2*(p+2))

detach()

results
```


```{r}
minimize <- function(res,column) {
  obj <- res[,column]
  j <- (1:length(obj))[obj==min(obj)]
  fla <- names(res)[2:(pmax+1)][unlist(res[j,2:(pmax+1)])]
  fla <- paste("log.per.cap.income ~", paste(fla, collapse=" + "))
  return(fla)  
}

```


```{r}
best.bic.fla <- minimize(results,"BIC")


best.aic.fla <- minimize(results,"AIC")
```


```{r}
summary(best.bic.model <- lm(best.bic.fla,data=cdi.trans))

summary(best.aic.model <- lm(best.aic.fla,data=cdi.trans))

```
We can see that the best BIC model and the best AIC model differ only in that the best AIC model adds \verb|pop.65_plus| to the BIC model.

Let's see what happens when we add region back in:
```{r}

best.bic.plus.region <- update(best.bic.model, . ~ . + region)

best.aic.plus.region <- update(best.aic.model, . ~ . + region)

cbind(AIC(best.bic.model,
          best.bic.plus.region,
          best.aic.model,
          best.aic.plus.region),
      BIC(best.bic.model,
          best.bic.plus.region,
          best.aic.model,
          best.aic.plus.region)
)

```
So, we see that the minimum BIC model does not include region, but the minimum AIC model does.


#### Stepwise:

Fortunately, the \verb|stepAIC| procedure does know how to deal with categorical variables like \verb|region|, so no special treatment of region is required:
```{r}
summary(best.bic.stepwise <- stepAIC(cdi.1,k=log(n),trace=0))
summary(best.aic.stepwise <- stepAIC(cdi.1,k=2,trace=0))
```
The stepwise AIC model adds the variables \verb|pop.65_plus| and \verb|region| to the stepwise BIC model.  Otherwise they are the same.

#### Lasso:

Finally, we try the lasso, with $\lambda$ selected by cross-validation.  Because \verb|glmnet| cannot deal with categorical variables, we take \verb|region| out, and then check to see if it should be added in again, just as we did with the all-subsets procedure.

```{r}

y <- cdi.trans$log.per.cap.income
X <- as.matrix(cdi.trans[,!(names(cdi.trans) %in% c("log.per.cap.income","region"))])

result <- cv.glmnet(X,y)

plot(result)

c(lambda.1se=result$lambda.1se,lambda.min=result$lambda.min)

log(c(log.lambda.1se=result$lambda.1se,log.lambda.min=result$lambda.min))

cbind(coef(result,s=result$lambda.1se),coef(result,s=result$lambda.min))

```
We can see that the predictors for the lasso model that minimizes 10-fold cross-validation prediction error include all the available variables (possibly indicating overfit) while the predictors for the lasso model with prediction error one SE above the minimum includes only the variables \verb|log.land.area|, \verb|pop.18_34|, \verb|log.doctors|, \verb|pct.back.deg|, \verb|pct.below.pov|, and \verb|pct.unemp|.

We need to see if adding \verb|region| back in improves the model or not.  For this, we will refit the two lasso models without the lasso penalty, and compare (using AIC and BIC) with the same models but with \verb|region| added in.

```{r}
result.1se <- as.matrix(coef(result,s=result$lambda.1se))
names.1se <- dimnames(result.1se)[[1]][result.1se!=0][-1]

result.min <- as.matrix(coef(result,s=result$lambda.min))
names.min <- dimnames(result.min)[[1]][result.min!=0][-1]

lasso.1se.fla <- paste("log.per.cap.income ~",paste(names.1se,collapse=" + "))
lasso.min.fla <- paste("log.per.cap.income ~",paste(names.min,collapse=" + "))

lasso.1se.model <- lm(lasso.1se.fla,data=cdi.trans)
lasso.1se.plus.region <- update(lasso.1se.model, . ~ . + region)

lasso.min.model <- lm(lasso.min.fla,data=cdi.trans)
lasso.min.plus.region <- update(lasso.min.model, . ~ . + region)

cbind(AIC(lasso.1se.model,
          lasso.1se.plus.region,
          lasso.min.model,
          lasso.min.plus.region),
      BIC(lasso.1se.model,
          lasso.1se.plus.region,
          lasso.min.model,
          lasso.min.plus.region)
)


```
Interstingly, AIC likes the most complex model \verb|lasso.min.plus.region|, whereas BIC likes the simplest model \verb|lasso.1se.model|.

Let's summarize the models we have so far in a table:

```{r}
var.names <- c("(Intercept)",names(cdi.trans)[names(cdi.trans)!="log.per.cap.income"])
var.names <- c(var.names[!(var.names=="region")], "regionNE" ,"regionS", "regionW", "regionNC") 

tab <- matrix(NA,ncol=6,nrow=length(var.names))
dimnames(tab) <- list(var.names,c("AIC.all",
                                  "AIC.step",
                                  "lasso.AIC",
                                  "BIC.all",
                                  "BIC.step",
                                  "lasso.BIC"))

models <- list(best.aic.plus.region,
               best.aic.stepwise,
               lasso.min.plus.region,
               best.bic.model,
               best.bic.stepwise,
               lasso.1se.model
               )

for (m in 1:6) {
  coefs <- summary(models[[m]])$coef
  vars <- dimnames(coefs)[[1]]
  vals <- coefs[,1]
  vals <- ifelse(vals==0,NA,vals)
  tab[vars,m] <- vals
}
tab <- rbind(tab,AIC=sapply(models,AIC),BIC=sapply(models,BIC))

round(tab,4)

```
Note that the stepwise and all-subsets approaches led to the same model, using AIC as the measure; but the lasso gave the largest model (identical to cdi.1, and probably an overfitting model).  Similarly the all-subsets and stepwise approaches used with BIC led to just a single model, and the lasso gave a slightly smaller one (the results here may vary a bit, depending on the random selection of folds in the \verb|cv.glmnet| algorithm).

Since AIC is minimized for the AIC.step model (best.aic.stepwise) and BIC is minimized for the BIC.step model (best.bic.stepwise), I will just consider these two models going forward. 

I next considered two-way interactions, with code like this

```{r eval=FALSE}
summary(best.aic.stepwise.2 <- update(best.aic.stepwise, . ~ .^2))
summary(best.bic.stepwise.2 <- update(best.bic.stepwise, . ~ .^2))
```
(output not shown, to save space...)

While both updated models produced some significant interactions, they were largely uninterpretable (e.g. \verb|pop.65_plus:pct.bach.deg|), had very small coefficients, and sapped the main effects of statistical significance (this often happens, because interactions are collinear with their main effects).  Because of my lack of success with two-way interactions, I did not try three-way or higher interactions.

### Diagnostics

Let's look at diagnostic plots for the two models.

First, best.aic.stepwise:
```{r}
vif(best.aic.stepwise)
par(mfrow=c(2,2))
plot(best.aic.stepwise)
```


```{r fig.width=7.5,fig.height=9.5}
## I divided the added-variable and marginal model plots up into several pages each 
## to make them more legible...
avPlots(best.aic.stepwise, ~ log.land.area + pop.18_34 + pop.65_plus + log.doctors)
avPlots(best.aic.stepwise, ~ pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp)
mmps(best.aic.stepwise, ~ log.land.area + pop.18_34 + pop.65_plus + log.doctors)
mmps(best.aic.stepwise, ~ pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp)
```

Next, best.bic.stepwise:
```{r}
vif(best.bic.stepwise)
par(mfrow=c(2,2))
plot(best.bic.stepwise)
```


```{r fig.width=7.5,fig.height=9.5}
## I divided the added-variable and marginal model plots up into several pages each 
## to make them more legible...
avPlots(best.bic.stepwise,~ log.land.area + pop.18_34 + log.doctors)
avPlots(best.bic.stepwise,~ pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp)
mmps(best.bic.stepwise,~ log.land.area + pop.18_34 + log.doctors)
mmps(best.bic.stepwise,~ pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp)
```

I don't see much in these plots to differentiate the two models.  The casewise diagnostic plots asre not perfect, but they only show fairly mild violations of the modeing assumptions; the VIF'a are a little higher for the best.AIC.stepwise model, but that is not surprising since it has more predictors.  The avPlots and marginal model plots all look pretty good; there are no suggestions that further transformations are needed for either model.

For discussion with the social scientist I will take the bigger model, \verb|best.AIC.stepwise|, because it has more variables to have a conversation about. 

(You may have chosen a smaller model here.  That is fine as long as you gave a thoughtful reason for your choice.)


\rm

## Problem 2(d).  

Provide a careful and easy-to-follow interpretation of your final model for a client or collaborator who is more interested in social, economic and health factors than in mathematics and statistics.  This should be long enough that you hit all the important points, but not so long that your collaborator gets bored and stops reading.

\em

Here is a summary of the model.
```{r}

vif(best.aic.stepwise)
summary(best.aic.stepwise)

```

Since the model has a logarithm for $y$ and for some of the $x$'s, the interpretations in the handout "log xform and percent interpretation.pdf" from week03 will be important.  In particular
\begin{itemize}
\item When $\log(x)$ and $\log(y)$ are both in the model, a 1\% change in $x$ is associated with an approximately $\beta$\% change in $E[y]$.
\item When $x$ and $log(y)$ are in the model, a one-unit change in $x$ is associated with an approximately $\beta\cdot100$\% change in $E[y]$.
\end{itemize}
The model summarized above provides a pretty good account of variation in per-capita income across counties in the data set. Indeed $R^2\approx85\%$ of that variation is accounted for by the model, a remarkably high $R^2$ for social science data.

From the estimated coefficients we can infer that
\begin{itemize}

\item A one percent increase in land area is assocated with a 0.03\% decrease in avarage per-capita income. This seems counterintuitive until one realizes that major centers of industry and commerce are mostly located in states with smaller populations.

\item A 1\% increase in the number of doctors is associated with a 0.06\% increase in average per-capita income.  This makes sense, since doctors are usually on the high end of the earning scale.

\item A one-point increase in the percentage of bachelors' degree holders is associated with a roughly 0.02$\times$100\%$=$2\% increase on average in per capita increase.  This makes sense, since degree-holdrs tend to earl more mony than on-degree hopders.

\item A one-point increase in the percentage of persons graduating from high school is associated with a 0.5\% {\em drop} in average per-capita income.  This seems counterintuitive, but what might be going is this: since virtually all bachelors' degree holders are also hs graduates, quite often both predictors will change at the same time, and so it makes sense to look at the combined effect of $0.0159212-0.0055059=0.0104153$: a one-unit increase in both variables is associated with about a 1\% increase in expected per capita income.

\item A one point increase in the percent unemployed is associated with a 0.9\% increase in per-capita income.  This seems conterintuitive, but again this variable may be mildly collinear with pecent below poverty, which has a negative effect on log-income.  

\item The number of people aged 18--35 (and so not yet peak income earners) and the number of people 65 or odler (and hence past peak income years) are both negatively associated with per-capita income, and both of these make sense.  

\item Region matters, to the extent that per capita income in the South seems noticably lower than in the North-Central region (and there is no noticable difference in per-capita income between the other regions and the NE region). 

\end{itemize}
We can see from this list that changes in the predictors affect per capita income in the ways that one might expect.  The exceptions seem to be due to mathematical co-incidences in the relationships between the variables, rather than substantial counterintuitive effects.

\rm

