```{r, message = FALSE}
library(mosaic)
#### For reproducibility
set.seed(4321)

#### Generate some artificial data, where Y depends
#### on a few Xs, and has no relationship to others

## N = number of data points to generate
## K = number of potential predictors to generate
N <- 15; K <- 10

## Generate random predictors (from a Normal distribution)
Predictors <- do(N) * rnorm(K, mean = 0, sd = 1)

## Name them X1, ..., XK (here K = 10)
names(Predictors) <- paste0("X", 1:K)
```

```{r, message = FALSE}
#### Define the "true" model

## set beta0 to 0
beta0 <- 0
## choose the other betas randomly between -1 and 1
beta <- runif(4, min = -1, max = 1)
beta
```

```{r, message =  FALSE}
## Create the "true" prediction function
f_true <- function(X, beta0, beta) 
{ 
  result = beta0
  for(k in 1:length(beta)) {
    result = result + beta[k] * X[,k] 
  }
  return(result)
}

## Create some "true" y-hat values, using the first 4 predictors
Yhat_true = f_true(Predictors[,1:4], beta0, beta)

## Create some sample residuals
epsilon = rnorm(N, mean = 0, sd = 1)

## Assemble the results into a data frame
FakeData <- cbind(
  Predictors, 
  Yhat_true, 
  epsilon, 
  Y = Yhat_true + epsilon)
```

```{r, message = FALSE}
### Fit some models, using some of the predictors

## The first few use too few predictors
m1 <- lm(Y ~ X1, data = FakeData)
m2 <- lm(Y ~ X1 + X2, data = FakeData)
m3 <- lm(Y ~ X1 + X2 + X3, data = FakeData)

## This one uses exactly the predictors that matter
m4 <- lm(Y ~ X1 + X2 + X3 + X4, data = FakeData)

## The rest include predictors that are unrelated to Y
m5 <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data = FakeData)
m6 <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = FakeData)
m7 <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data = FakeData)
m8 <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data = FakeData)
m9 <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9, data = FakeData)
m10 <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10, data = FakeData)
```

```{r}
### Defining a convenience function to extract adjusted R^2 
### from a fitted model
adj.rsquared <- function(model) {return(summary(model)$adj.r.squared)}
```

```{r, message = FALSE}
library(olsrr)

### Compute fit measures for our list of models
models <- list(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)

### Compute metrics for each model and assemble the 
### results into a data frame
results <- data.frame(
  model_id = 1:K,
  R2 = sapply(models, rsquared),
  adjR2 = sapply(models, adj.rsquared),
  Cp = sapply(models, ols_mallows_cp, m10),
  predR2 = sapply(models, ols_pred_rsq)
)
```

```{r}
### Plot the R^2 results
gf_point(
  R2 ~ factor(model_id), 
  data = results, 
  color = "red",
  xlab = "Number of Predictors Used",
  ylab = "R-squared (higher is better)") %>%
gf_point(
  adjR2 ~ factor(model_id), 
  data = results, 
  color = "blue") %>%
gf_point(
  predR2 ~ factor(model_id),
  data = results,
  color = "green")
```

```{r}
### Plot Mallow's Cp (which is in different units)
gf_point(
  Cp ~ factor(model_id), 
  data = results, 
  color = "orange",
  xlab = "Number of Predictors Used",
  ylab = "Mallows' Cp (lower is better)")
```