```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = 'G:\\Il mio Drive\\Teaching\\Data Science Lab\\Lecture 3')
```

# Lecture 3: Data Modeling

The goal of a model is to provide a simple, low-dimensional, interpretable summary of a dataset. Models are a really useful way to help you peel back layers of structure as you are exploring your dataset. Every statistical model can be "divided" in two parts:
\begin{enumerate}
\item a family of models that express a prece, but generic, pattern that you want to capture (i.e., the pattern can be a straight line or a quadratic curve);
\item a fitted model, that can be found by selecting the family of models that is the closest to your data.
\end{enumerate}

It is important to understand that a fitted model is just the closest model from a family of models. This implies that you have the "best" model according to some criteria and based on a set of assumptions. This does not imply that your model is a good model or that your model is "true". George Box, a famous british statistician, once said one of the most quoted statistical quotes: \textit{"all models are wrong, but some are useful"}.

It is worth reading the fuller context of the quote as it is quite illustrative of the philosophy behind any statistical model:
\textit{"Now it would be very remarkable if any system existing in the real world could be exactly represented by any simple model. However, cunningly chosen parsimonious models often do provide remarkably useful approximations. For example, the law PV = RT relating pressure P, volume V and temperature T of an "ideal" gas via a constant R is not exactly true for any real gas, but it frequently provides a useful approximation and furthermore its structure is informative since it springs from a physical view of the behavior of gas molecules. For such a model there is no need to ask the question "Is the model true?" If "truth" is to be the "whole truth" the answer must be "No." The only question of interest is \textbf{"Is the model illuminating and useful?"}}

This does not mean that all the models are wrong and, we should just go for the least wrong model. This quote should be interpeted as a call for careful laying down the assumptions on which the quality of the model is built on. As Berkeley statisticain Mark Van Der Laan stated in a recent article on "The statistical formulation and theory should define the algorithm" [source](https://magazine.amstat.org/blog/2015/02/01/statscience_feb2015/).


In this lecture we will go see how to perform in R two types of models:
\begin{enumerate}
\item linear regression models;
\item regularization and selection models.
\end{enumerate}


```{r, warning=FALSE}
library(tidyverse)
library(modelr)
library(hdm)
library(stabs)
library(AER)
library(sandwich)
library(lmtest)
library(broom)
library(lars)
library(glmnet)
```

```{r}
library(readxl)
data <- read_excel("G:\\Il mio Drive\\Econometrics Lab\\Data\\Compustat Data.xlsx")
data <- data[, !names(data) %in% c("Interest Expense - Total (Financial Services)",
                                   "Net Interest Income", "Nonperforming Assets - Total")]
data_clean <- na.omit(data)
```

```{r}
x <- data_clean$`Assets - Total`[which(data_clean$`Assets - Total`< 
                                 quantile(data_clean$`Assets - Total`, 0.95))]
y <- data_clean$`Sales/Turnover (Net)`[which(data_clean$`Assets - Total`<
                                quantile(data_clean$`Assets - Total`, 0.95))]

reg_data <- as.data.frame(cbind(x, y))
```

```{r}
ggplot(reg_data, aes(x, y)) +
  geom_point()
```

You can see a quite clear pattern in the data. Let's now use a model to capture the pattern and make it more explicit.

Let's first generate a set of random model an let's overlay them on the data.

```{r}
models <- tibble(
  beta1 = runif(length(x), 0, 200),
  beta2 = runif(length(x), -4, 4)
)
```

```{r}
ggplot(reg_data, aes(x, y)) +
  geom_abline(
    aes(intercept = beta1,
        slope = beta2),
    data = models, alpha = 1/15
  ) + 
  geom_point()
```

```{r}
model1 <- function(beta, data){
  beta[1] + data$x * beta[2]
}
```


```{r}
fitted.values <- model1(c(50, 1.5), reg_data)
```

```{r}
head(fitted.values)
```

Let's now get the residuals of our model.

```{r}
measure_distance <- function(mod, data) {
 diff <- data$y - model1(mod, data)
 sqrt(mean(diff ^ 2))
}
measure_distance(c(50, 1.5), reg_data)
```

We can use "purrr" to compute the distance for all the models defined previously. We will need a helper function because our distance expectes the model as a numeric vector of length 2.

```{r}
reg_data_dist <- function(beta1, beta2) {
 measure_distance(c(beta1, beta2), reg_data)
}
models <- models %>%
 mutate(dist = purrr::map2_dbl(beta1, beta2, reg_data_dist))
```

```{r}
models
```

We can now overlay the best 10 models on the data.

```{r}
ggplot(reg_data, aes(x, y)) +
 geom_point(size = 2, color = "grey30") +
 geom_abline(
 aes(intercept = beta1, slope = beta2, color = -dist),
 data = filter(models, rank(dist) <= 10)
 )

```

We can also think about these models as observations, and visualize
them with a scatterplot of beta1 versus beta2, again colored by -dist. We
can no longer directly see how the model compares to the data, but
we can see many models at once. Again, I've highlighted the 10 best
models, this time by drawing red circles underneath them:


```{r}
ggplot(models, aes(beta1, beta2)) +
 geom_point(
 data = filter(models, rank(dist) <= 1),
 size = 4, color = "red"
 ) +
 geom_point(aes(colour = -dist))
```

Instead of trying lots of random models, we could be more systematic and generate an evenly spaced grid of points (this is called a grid search). I picked the parameters of the grid roughly by looking at where the best models were in the preceding plot:

```{r}
grid <- expand.grid(
 beta1 = seq(0, 200, length = 50),
 beta2 = seq(-4, 4, length = 50)
 ) %>%
 mutate(dist = purrr::map2_dbl(beta1, beta2, reg_data_dist))
grid %>%
 ggplot(aes(beta1, beta2)) +
 geom_point(
 data = filter(grid, rank(dist) <= 1),
 size = 4, colour = "red"
 ) +
 geom_point(aes(color = -dist))

```

When you overlay the best 10 models back on the original data, they
all look pretty good:

```{r}
ggplot(reg_data, aes(x, y)) +
 geom_point(size = 2, color = "grey30") +
 geom_abline(
 aes(intercept = beta1, slope = beta2, color = -dist),
 data = filter(grid, rank(dist) <= 1)
 )
```

You could imagine iteratively making the grid finer and finer until you narrowed in on the best model. But there's a better way to tackle that problem: a numerical minimization tool called Newton-Raphson search. The intuition of Newton-Raphson is pretty simple: you pick a starting point and look around for the steepest slope. You then ski down that slope a little way, and then repeat again and again, until you can't go any lower. In R, we can do that with optim():

```{r}
best <- optim(c(0, 0), measure_distance, data = reg_data)
best$par
```

```{r}
ggplot(reg_data, aes(x, y)) +
 geom_point(size = 2, color = "grey30") +
 geom_abline(intercept = best$par[1], slope = best$par[2])
```

Don't worry too much about the details of how optim() works. It's the intuition that's important here. If you have a function that defines the distance between a model and a dataset, and an algorithm that can minimize that distance by modifying the parameters of the model, you can find the best model. The neat thing about this approach is that it will work for any family of models that you can
write an equation for. There's one more approach that we can use for this model, because it
is a special case of a broader family: linear models. A linear model has the general form
$y = a_1 + a_2 \cdot x_1 + a_3 \cdot x_2 + ... + a_n \cdot x_{(n - 1)}$. 
So this simple model is equivalent to a general linear model where n is 2 and $x_1$ is $x$. R has a tool specifically designed for fitting linear models called lm(). lm() has a special way to specify the model family: formulas. Formulas look like $y ~ x$, which lm() will translate to a function like $y = a_1 + a_2 * x$. We can fit the model and look at the output:

```{r}
model_1 <- lm(y ~ x, data = reg_data)
summary(model_1)
```

Now let's add an additional variable in the linear regression to compare the two different models.

```{r}
z <- data_clean$Employees[which(data_clean$`Assets - Total`<
                          quantile(data_clean$`Assets - Total`, 0.95))]
reg_data <- cbind(reg_data, z)
```

```{r}
model_2 <- lm(y ~ x + z, data = reg_data)
summary(model_2)
```

In R, you can either write down all the variables that you want to use as regressors in your model or you can just use $y \sim \: .$.

```{r}
model_3 <- lm(y ~ ., data = reg_data)
summary(model_3)
```

A very easy way to compare two different linear regressions is through the likelihood ratio test. In statistics, the likelihood-ratio test assesses the goodness of fit of two competing statistical models based on the ratio of their likelihoods.

```{r}
library(lmtest)
lrtest(model_1, model_2)
```

$p < 0.001$ indicates that the model with all predictors fits significantly better than the model with only one predictor. Another "goodness-of-fit" measure that can be used is the $R^2$:
\begin{equation}
R^2= 1 - \frac{ESS}{TSS}.
\end{equation}

```{r}
summary(model_1)$r.squared
summary(model_2)$r.squared
```

We can also get the fitted values of the model for any $x$ and $z$ by running the following chunck of code.

```{r}
coeffs = coefficients(model_2)
assets = 159 
employees = 2
y <- coeffs[1] +coeffs[2]*assets +coeffs[3]*employees
y
```

Or, equivalently:

```{r}
newdata <- data.frame(x = 159, z = 2)
predict(model_2, newdata) 

predict(model_2, newdata, interval="confidence") 
```

Once we fitted our favourite model, we can check the residuals from the model: $e_i = y_i - \hat{f}(x_i)$.

```{r}
model.res = resid(model_2)
plot(reg_data$y, model.res, ylab="Residuals", xlab="Sales", main="Residuals v. Sales") 
abline(0, 0)  
```

Moreover, we can standardize the residuals and plot them against normalized scores for the outcome variable. 

```{r}
model_2.stdres = rstandard(model_2)
qqnorm(model_2.stdres ,  ylab="Standardized Residuals",  xlab="Normal Scores",  main="Standardized Residuals v. Sales") 
qqline(model_2.stdres)
```

In R, you can introduce an interaction between the regressors by using $*$. Always remember to include also the single regressors in the formula.

```{r}
model_int<-lm(y ~ x  + z + x*z, data = reg_data)
summary(model_int)
```

You can't directly introduce a quadratic term in the regression formula. Hence, you need to create an additional variable with the square term and then you can include it in the regression.

```{r}
x2 <- x^2
model_squared<-lm(y ~ x  +  x2 + z + x*z, data = reg_data)
summary(model_squared)
```

## Variables Selection

Here, I am going to show an application based based on an article from Barro and Lee (1994). The hypothesis we want to test is if less developed countries, with lower GDP per capita, grow faster than developed countries. In other words, there is a catch up effect. The model equation is as follows:

\begin{equation}
 y_i=\alpha_0d_i+\sum_{i=1}^p\beta_jx_{i,j}+\varepsilon_i
\end{equation}

where $y_i$ is the GDP growth rate over a specific decade in country $i, d_i$ is the log of the GDP at the beginning of the decade, $x_{i,j}$ are controls that may affect the GDP. We want to know the effects of $d_i$ on $y_i$, which is measured by $\alpha_0$. If our catch up hypothesis is true, $\alpha_0$ must be positive and hopefully significant.

The dataset is available in the package. It has 62 variables and 90 observations. Each observation is a country, but the same country may have more than one observation if analysed in two different decades. The large number of variables will require some variable selection, and I will show what happens if we use a single LASSO selection and the Double Selection. The hdm package does all the DS steps in a single line of code, we do not need to estimate the two selection models and the Post-OLS individually. I will also run a naive OLS will all variables just for illustration.
This application can be found [here](https://www.r-bloggers.com/the-package-hdm-for-double-selection-inference-with-a-simple-example/).

```{r}
rm(list=ls())
data("GrowthData") # = use ?GrowthData for more information = #
dataset <- GrowthData[,-2] # = The second column is just a vector of ones = #
```

```{r}
# = Naive OLS with all variables = #
# = I will select only the summary line that contains the initial log GDP = #
summary(lm(Outcome ~., data = dataset))
```

```{r}
OLS <- summary(lm(Outcome ~., data = dataset))$coefficients[1, ]
OLS
```

```{r}
rlasso(Outcome~., data = dataset, post = FALSE)
# = Single step selection LASSO and Post-OLS = #
# = I will select only the summary line that contains the initial log GDP = #
lasso <- rlasso(Outcome~., data = dataset, post = FALSE) # = Run the Rigorous LASSO = #
selected <- which(coef(lasso)[-c(1:2)] !=0) # = Select relevant variables = #
selected
fm <- paste(c("Outcome ~ gdpsh465", names(selected)), collapse = "+")
SS <- summary(lm(fm, data = dataset))$coefficients[1, ]
SS
```

```{r}
# = Double Selection = #
X <- as.matrix(dataset[,-1])
y <- dataset$Outcome
DS <- rlassoEffects(X , y, I = ~ dataset$gdpsh465, data = dataset)
DS <- summary(DS)$coefficients[1,]
```

```{r}
results <- rbind(OLS,SS,DS)
results
```

The OLS estimate is positive, however the standard error is very big because we have only 90 observations for more than 60 variables. The Single Selection estimate is also positive and, in this case, significant. However, the Double Selection showed a negative and significant coefficient. If the DS is correct, our initial catch up hypothesis is wrong and poor countries grow less than rich countries. We can't say that the DS is correct for sure, but it is backed up by a strong theory and lots of simulations that show that the SS is problematic. It is very, very unlikely that the SS results are more accurate than the DS. It is very surprising how much the results can change from one case to the other. You should at least be skeptic when you see this type of modelling and the selection of controls is not clear.

The "hdm" package has several other implementations in this framework such as instrumental variables and logit models and there are also more examples in the package vignette.

## Stability Selection

```{r, warning=FALSE}
####################################################################
### using stability selection with Lasso methods:
stab.lasso <- stabsel(x = X, y = y,
                             fitfun = lars.lasso, cutoff = 0.75,
                             PFER = 1)

stab.stepwise <- stabsel(x = X, y = y,
                                fitfun = lars.stepwise, cutoff = 0.75,
                                PFER = 1)
plot(stab.lasso, main = "Lasso")
plot(stab.stepwise, main = "Stepwise Selection")
## --> stepwise selection seems to be quite unstable even in this low
##     dimensional example!
```

```{r}
## set seed (again to make results comparable)
set.seed(1234)
stab.glmnet <- stabsel(x = X, y = y,
                              fitfun = glmnet.lasso, cutoff = 0.75,
                              PFER = 1)
plot(stab.glmnet, main = "Lasso (glmnet)")
```

```{r}
## Select variables with maximum coefficients based on lasso estimate
set.seed(1234) 
 
## use cross-validated lambda 
lambda.min <- cv.glmnet(x = as.matrix(X), y = y)$lambda.min
stab.maxCoef <- stabsel(x = X, y = y,
                               fitfun = glmnet.lasso_maxCoef, 
                               # specify additional parameters to fitfun
                               args.fitfun = list(lambda = lambda.min),
                               cutoff = 0.75, PFER = 1)
                               
## WARNING: Using a fixed penalty (lambda) is usually not permitted and 
##          not sensible. See ?fitfun for details.
      
## now compare standard lasso with "maximal parameter estimates" from lasso
plot(stab.maxCoef, main = "Lasso (glmnet; Maximum Coefficients)")
plot(stab.glmnet, main = "Lasso (glmnet)")
## --> very different results.
```

## Instruments Selection

Reproduction of the analysis by Angrist and Krueger (1991).

```{r}
load("G:\\Il mio Drive\\Teaching\\Data Science Lab 2020\\angrist_krueger_1991.rda")
ak91 <- mutate(ak91,
               qob_fct = factor(qob),
               q4 = as.integer(qob == "4"),
               yob_fct = factor(yob))
```

Average years of schooling by quarter of birth for men born in 1930-39 in the 1980 US Census.

```{r}
ak91_age <- ak91 %>%
  group_by(qob, yob) %>%
  summarise(lnw = mean(lnw), s = mean(s)) %>%
  mutate(q4 = (qob == 4))
ggplot(ak91_age, aes(x = yob + (qob - 1) / 4, y = s)) +
  geom_line() +
  geom_label(mapping = aes(label = qob, color = q4)) +
  theme(legend.position = "none") +
  scale_x_continuous("Year of birth", breaks = 1930:1940) +
  scale_y_continuous("Years of Education", breaks = seq(12.2, 13.2, by = 0.2),
                     limits = c(12.2, 13.2))
```

Average log wages by quarter of birth for men born in 1930-39 in the 1980 US Census.

```{r}
ggplot(ak91_age, aes(x = yob + (qob - 1) / 4, y = lnw)) +
  geom_line() +
  geom_label(mapping = aes(label = qob, color = q4)) +
  scale_x_continuous("Year of birth", breaks = 1930:1940) +
  scale_y_continuous("Log weekly wages") +
  theme(legend.position = "none")
```


Regress log wages on 4th quarter.

```{r}
mod1 <- lm(lnw ~ q4, data = ak91)
coeftest(mod1, vcov = sandwich)
```

Regress years of schooling on 4th quarter.

```{r}
mod2 <- lm(s ~ q4, data = ak91)
coeftest(mod2, vcov = sandwich)
```

IV regression of log wages on years of schooling, with 4th quarter as an instrument for years of schooling.

```{r}
mod3 <- ivreg(lnw ~ s | q4, data = ak91)
coeftest(mod3, vcov = sandwich, diagnostics = TRUE)
```

IV reg using interaction between years and quarters as instruments. Controls for year of birth.

```{r}
mod4 <- ivreg(lnw ~ s | q4*yob_fct, data = ak91)
summary(mod4, vcov = sandwich, diagnostics = TRUE)
```

```{r}
lasso <- rlasso(s ~ q4*yob_fct, data = ak91, post = FALSE) # = Run the Rigorous LASSO = #
selected <- which(coef(lasso)[-c(1:2)] !=0) # = Select relevant variables = #
selected

ak91$yob_fct1936 <- ifelse(ak91$yob_fct==1936,1,0)
mod5 <- ivreg(lnw ~ s | q4*yob_fct1936,
                        data = ak91)
summary(mod5, vcov = sandwich, diagnostics = TRUE)
```
