\newif\ifsolutions
\solutionsfalse

``` {r include = FALSE}
solutions <- FALSE
set.seed(02138) #CGIS Knafel
```

```{r eval = !solutions, echo = FALSE}
lambda <- NA
lambda_y <- NA
lambda_y_auto <- NA
lambda_y <- NA
lambda_t <- NA
lambda_t_auto <- NA
auto_model <- NA
t_terms <- NA
y_terms <- NA
selected_vars <- list()
```


**Read the following instructions carefully:**

- All answers must be typed up. Your text, math, code, and
  figures should be included together for each problem. RMarkdown will
  do this automatically and is strongly encouraged.
 
- Late submission is until 2/25 midnight, your score will be discounted .
 
- You are encouraged to work in groups, but you should write up your
  answers alone and cite your peers.
 
- Upload   your compiled PDF to Gradescope.

\newpage

```{r include=FALSE, echo=FALSE, message=FALSE}
library(knitr)
opts_chunk$set(comment = NA)
opts_chunk$set(echo = solutions)
opts_chunk$set(eval = solutions)
opts_chunk$set(message = solutions)
opts_chunk$set(warning = FALSE)
```


```{r echo = FALSE, message = FALSE}
library(ggplot2)
library(RColorBrewer)
library(scales)
library(grid)
library(wesanderson)
theme_andy <- function() {
  
  # Generate the colors for the chart procedurally with RColorBrewer
  palette <- brewer.pal("Greys", n=9)
  color.background = palette[2]
  color.grid.major = palette[3]
  color.axis.text = palette[6]
  color.axis.title = palette[7]
  color.title = palette[9]
  
  # Begin construction of chart
  theme_bw(base_size=9)  + 
    
    # Set the entire chart region to a light gray color
    theme(panel.background=element_rect(fill=color.background, color=color.background))  + 
    theme(plot.background=element_rect(fill=color.background, color=color.background))  + 
    theme(panel.border=element_rect(color=color.background))  + 
    
    # Format the grid
    theme(panel.grid.major=element_line(color=color.grid.major,size=.25))  + 
    theme(panel.grid.minor=element_line(color=color.grid.major,size=.15))  + 
    theme(panel.grid.minor.y=element_blank())  + 
    theme(axis.ticks=element_blank())  + 
    
    # Format the legend
    theme(legend.background = element_rect(fill=color.background))  + 
    theme(legend.key = element_rect(fill=color.background))  + 
    theme(legend.text = element_text(size=7,color=color.axis.title))  + 
    
    # Set title and axis labels, and format these and tick marks
    theme(plot.title=element_text(color=color.title, size=10, vjust=1.25))  + 
    theme(axis.text.x=element_text(size=7,color=color.axis.text))  + 
    theme(axis.text.y=element_text(size=7,color=color.axis.text))  + 
    theme(axis.title.x=element_text(size=8,color=color.axis.title, vjust=0))  + 
    theme(axis.title.y=element_text(size=8,color=color.axis.title, vjust=1.25))  + 
    
    theme(strip.background = element_rect(colour = color.background, 
                                          fill = palette[3], size = 0.5)) +
    # Title fonts
    #theme(plot.title=element_text(family="OpenSans-CondensedBold", margin=margin(b=15)))  + 
    #theme(plot.subtitle=element_text(family="OpenSans-CondensedLightItalic"))  + 
    
    # Plot margins
    theme(plot.margin = unit(c(0.35, 0.2, 0.3, 0.35), "cm"))
}

#theme_set(theme_andy())
```


# Problem 1

The purpose of this problem is to place OLS, LASSO, and ridge
regression into a common framework---a family of estimators that
minimize an objective function based at least partly on the sum of
squared residuals (SSR).

For example, recall that the OLS estimator is a procedure that
searches the space of all possible linear regression coefficients,
$\mathbf{\beta}^\ast \in \mathbb{R}^k$ (where $k$ is the number of coefficients to
estimate), and picks the option that minimizes SSR. Thus, SSR is the
objective function in OLS. LASSO and ridge regression have modified
objective functions that include SSR, but also add a "penalty"
function that pulls estimates away from certain choices of
coefficients even when those choices might lead to a minimal SSR. In
this problem, we will visualize the shape of these penalty functions
and their effect on coefficient estimates. This problem uses a  
simulated dataset drawn from the following data-generating process:
$$y_i = 2 x_{1i} + 3 x_{2i} + \varepsilon_i,$$
where $X_k \sim N(0, 1) \text{ for } k = 1,2$ $\varepsilon_i \sim N(0, 1)$.

We provide code to visualize how SSR varies with $\mathbf{\beta}^\ast$. In the
corresponding plot, darker areas represent portions of the coefficient
search space that lead to higher SSR. The lightest point, marked with
a $\beta$, represents the solution of the optimization problem. Circular
contour lines (gray) indicate "level curves," or sub-regions in the
coefficient space that would lead to the same SSR.

You will do the same for LASSO and ridge objective functions, then
examine how the results change with the tuning parameter $\lambda$.

```{r include = TRUE, echo = TRUE, message = FALSE}
library(ggplot2)
library(dplyr)

## simulation dataset
set.seed(02139)
n <- 1000
d <- data.frame(x_1 = rnorm(n, mean=0, sd=1),
               x_2 = rnorm(n, mean=0, sd=1))
d$y <- 2 * d$x_1 + 3 * d$x_2 + rnorm(n, mean=0, sd=1)

## grid search for optimal beta
## We usually find beta hat analytically: (X'X)^{-1}X'y
## Here, though, we're exhaustively searching through combinations of the 
##  different possible betas, finding the error, and empirically finding the
##  best beta hats.
beta_grid <- expand.grid(beta_1=seq(-1,5,.1),beta_2=seq(-1,5,.1))

calc_ssr <- function(beta_star, d){
  X <- as.matrix(d[,c("x_1","x_2")])
  e <- d$y - X%*% beta_star
  return(crossprod(e)) 
}

## calculate the error for all combinations of beta1 and beta2
beta_grid$ssr <- apply(beta_grid, 1, calc_ssr, d)

## beta_ols is the the vector of two betas that empirically minimizes the
##  sum of squared error.
beta_ols <- beta_grid[which.min(beta_grid$ssr),]


## ssr contours and ols estimate
ggplot(beta_grid, aes(x=beta_1, y=beta_2)) + 
  geom_tile(aes(fill=log(ssr))) +
  stat_contour(aes(z=ssr), color="gray30") +
  # add the ols estimate
  geom_text(data=beta_ols, label="beta[plain(SSR)]",
           parse=TRUE, hjust=0, vjust=0) +
  geom_vline(xintercept=0) + 
  geom_hline(yintercept=0) +
  scale_fill_gradient(low="#FEE8C8", high="#E34A33",
                      breaks=range(log(beta_grid$ssr)), 
                      labels=c("low","high"), name="SSR") +
  xlab(expression(beta[1])) + 
  ylab(expression(beta[2])) + 
  coord_fixed() + 
  labs(title = "OLS objective function (SSR)")
```

**(a)** Calculate the penalty that ridge and lasso add to the OLS loss function, 
ignoring for now the sum of squared residuals. Visualize the shape of the penalties with
contours and shading as above, and in your own words, compare the different
shapes of the constraints. Can you reframe the penalty function in terms of a
prior on coefficients?

**Hint:** Assume $\lambda=1$ for now---the exact value of the penalty
function is less important than its shape.
 
**(b)** Next, visualize the full loss function for ridge (i.e., sum of SSR
plus the ridge penalty). On each graph, mark both ridge and OLS
coefficient estimates. Do this for $\lambda$ values of 250, 500, 1,000,
and 2,000. What sort of path do your ridge estimates follow as the
tuning parameter $\lambda$ increases?

Repeat for LASSO regression, with $\lambda$ set to 1,250, 2,500, 5,000,
and 10,000. What similarities and differences do you see?

 

# Problem 2

In this problem, we will compare the abilities of OLS and LASSO  to
recover causal effects in a simulation exercise. We will use each method to
the bias and variance of their respective $\hat{\tau}$s and generate a graphical
comparision of the two

**(a)**. Create a $1,000 \times 800$  matrix of "control" variables $\mathbf{X}$, with 1,000 samples for each $X_k$ drawn from $X_k \stackrel{iid}{\sim} N(0, 1) \text{ for } k = 1,...,800$. Draw a (confounded) treatment $T_i = X_{44i} + \alpha_i$, where $\alpha_i \sim N(0, 1)$  Assume a DGP of 

\begin{equation}
Y_i = 2 T_i - 12X_{2i}  - 15 X_{15i} + 9 X_{44i} + 13 X_{300i} + \epsilon_i,
\end{equation}

where $\epsilon_i \sim N(0, 1)$. Let $\tau = 2$ be the effect of $T$ on $Y$.

 
**(b)**. First, we will investigate the ability of OLS to recover the true
effect in the presence of these control variables. Write two functions to
generate OLS estimates of the treatment effect: one that includes all control
variables and one that assumes a correctly specificed model (i.e. that only
includes variables with an effect). 
Generate at least 200 simulated datasets using the DGP above and use your your
function to find $\hat{\tau}_{OLS}$ for both the true and mis-specified models.
Report the bias and variance of the estimates.

```{r message = FALSE}
# Parallelize the simulations below
library(foreach)
library(doParallel)
cl <- makeCluster(detectCores())
registerDoParallel(cl)

# Set the global number of replications for all simuations
REPS <- 200
```

 

**(c)** Next, we will use lasso to estimate $\tau$. Repeat the simulation above,
but this time using `glmnet` instead of `lm` to estimate $\tau$ as a function of all X. Unlike OLS,
lasso has a tunable parameter, $\lambda$, which must be set by the researcher. 
One way of picking a good value for $\lambda$ is through "grid search" and
cross-validation. Write a function that uses 5-fold cross-validation to pick the
value for $\lambda$ out of a range of possible $\lambda$s to minimize the
cross-validated mean squared prediction error. You should consider at least 100
possible valus for $\lambda$ between 0 and the minimum value of $\lambda$ for
which no features are selected. Code is provided below to find the maximum
$\lambda$ and to generate the "grid" of $\lambda$s to consider:

```{r, eval = FALSE, include = TRUE, echo = TRUE}
X_both <- cbind(treat, X)
max_lambda <- max(abs(colSums(X_both*y)))/ncol(X_both)
lambda_vec <- seq(0, max_lambda, length.out = GRID_N)
```

(Note: this method of finding a maximum $\lambda$ requires scaled and centered 
variables. Our variables are standard normals so no centering and scaling is
required in Problem 2.)

Compare your value of $\lambda$ to the value returned by `glmnet`'s built-in
function `cv.glmnet`. Report the variables that are selected using your value of
$\lambda$ and compare them to the variables in the true DGP. Then use your value
of $\lambda$ in your simulations. 


 
 
 

 
 

**(d)** Finally, graphically compare the simulated $\hat{\tau}$s of the 3 models (correct
OLS, realistic OLS and lasso ) in three histograms. What do you find?
Which model would you prefer in your own data analysis?

 
# Problem 3

  
In this question, we will compare several classification algorithms by applying them to United Nations Security Council (UNSC) resolutions. Resolution topics and text (parsed from PDF without post processing) are in  \texttt{unsc\_text.csv}

We would like to get an intuition for machine learning algorithms by visualizing them, but the document-term matrix (DTM) that summarizes our UNSC resolutions has over three thousand dimensions (columns, or unique words) even after stemming/stopping and trimming obscure words. To facilitate visualization, we will first pre-process the DTM using dimensional reduction. This step finds the best way to project, or summarize, the DTM into $\mathbb{R}^2$—a transformation that is easier to manipulate and understand, but still explains much of the variation in the original $\mathbb{R}^{3309}$ dataset.

**(a)**  Load \texttt{unsc\_text.csv}. Using regular expressions, drop resolutions for which no text was found—i.e., that contain only whitespace in the `text` field. Clean whitespace, numbers, and punctuation, then stem the document. Reduce the size of the DTM by removing words that only occur in two or fewer documents.
 
Finally, normalize each row so that it sums to unity, to help account for different document lengths. (Why?)


**(b)**  Singular value decompositions (SVD) are a re-expression of a $ n × k$ matrix (with $k > n$) in terms of $n$ orthogonal unit vectors, usually ordered by their explanatory power. The SVD is a useful way to obtain a lower-dimensional summary of the DTM: the first $m$ left-singular vectors,  $dtm\_svd\$u\$[,1:k]$   , represent a $m$-dimensional projection of the data, or the $m$ most important “concepts” that are embodied in the DTM. They sometimes have a substantive interpretation, but often it is unclear exactly what they represent.

In this problem, we will use the second and third left-singular vectors of the DTM to visualize some classification algorithms in a two-dimensional space. (The first left-singular vector, which loads heavily onto words like ``united``, ``nation``, ``security``, ``council``, ``resolution``, and ``paragraph``, appears to represent ``bureaucratic-ness``; while it is important in explaining many of the words that appear in the resolutions, it contains little useful information for classification.)
Get the SVD of dtm and extract its left-singular vectors $svd(dtm)\$u$. This is a $n×n$ matrix, where each column $svd(dtm)\$u\$[,i]$  is a different weighted combination of our columns of old (pre-SVD) covariates. In your own words, interpret the rows of U and the columns of V (svd(dtm)\$v).

Visualize all resolutions in a plot with the second and third left-singular vectors (LSV2 and LSV3, or $dtm\_svd\$u\$[,2]$  and $dtm\_svd\$u\$[,3]$`, respectively) on the $x$- and $y$-axis, respectively. On this plot, highlight resolutions that are hand-coded (by the UN Bibliographic Information System) as pertaining to the “IRAQ-KUWAIT SITUATION” and “WESTERN SAHARA QUESTION” in different colors. What do you see? Note that the SVD knows nothing about the subject area of resolutions; only word counts were used in the decomposition. (Hint: Remember that there is more than one agenda topic per resolution. You may wish to use regular expressions.)

For the remainder of this problem, we will work only with resolutions coded as pertaining to “IRAQ-KUWAIT SITUATION” or “WESTERN SAHARA QUESTION”. Create a dataframe of these resolutions, with the second and third left-singular vectors as covariates and their topic as outcome.


**(c)**  Next, we will use OLS as a classifier. Run a linear regression (without interaction term) of topic on LSV2 and LSV3, and interpret its fitted values as predictions by cutting them at some suitable threshold. Create a plot with the following components:

 
- LSV2 and LSV3 on the $x$- and $y$-axis, respectively
- Points representing the training data, with plotting symbols that distinguish between
resolutions relating to `Iraq` and `Sahara`
- Shaded regions that indicate your model’s prediction for incoming test data (for test data projected to the subspace of RSV2 and RSV3).


Hint: You may find it helpful to use \texttt{expand.grid()}, \texttt{predict()}, and \texttt{geom.tile()}.

**(d)**   Repeat for the following models:


- OLS with interaction
- Linear SVM
- SVM with radical (Gaussian) Kernel and $\gamma=1$
- SVM with radical (Gaussian) Kernel and $\gamma=2$
- $1-$nearest neighbor
- $3-$nearest neighbor


While there are many good implementions of support vector machines and K-nearest neighbors, two options are \texttt{svm()} (package \texttt{e1071}) and \texttt{knn()} (package \texttt{class}).

