You can download the Rmd file [here](lab-10.Rmd).

```{r setup, include=FALSE}
# Install and load required packages
library(caret)
library(glmnet)
library(randomForest)
library(car)
library(tidyverse)
library(psych)

# Knit options
knitr::opts_chunk$set(echo = TRUE)
```

# Background

Statistical models can be built to serve different research purposes, such as **explanation** or **prediction**.

-   When one's purpose is **explanation**, the researcher writes a statistical model that tests a specific literature-informed theory (the model may have even been pre-registered).
    -   Typically, the researcher wants the model to have interpretable parameter estimates that provide an understanding of the nature of the relationship between each predictor and the outcome variable
-   When one's purpose is **prediction**, the researcher wants to build a statistical model that does the best at **accurately predicting new data**.
    -   The model does not necessarily need to produce interpretable parameter estimates as long as, when applied to new data, the statistical model does a good job at predicting those values.

You could also be interested in achieving **a mix** of these two goals.

# Machine Learning

Machine learning encompasses a set of techniques that have been developed to **build good predictive models**. That is, the goal of machine learning is to build a statistical model that does the best at predicting new, not-yet-seen data points.

## Common Machine Learning Algorithms

There are many algorithms that can be used to build a machine learning model. You can predict continuous or categorical outcomes - just make sure you choose an algorithm that's appropriate for that type of DV!

-   Linear regression, lm()

-   Logistic regression, LogitBoost()

-   Support vector machines, svm() or svmLinear()

-   Random forests, randomForest()

And there are hundreds more... for a full list: names(getModelInfo())

As the algorithm you use becomes increasingly complex, the parameter estimates may become less and less interpretable, were you to examine them.

# Overfitting

A very important consideration in the building of predictive models is overfitting. **Overfitting is the problem of fitting a model too well to the variability of a specific data set. This could result in the model having poor predictive power when applied to a new data set.**

When building a statistical model, you want to meet a balance between finding a model that seems to fit the data well, **without accounting for so much error that you overfit the model**.

# Training Set vs Testing Set

One of the main innovations of machine learning is that it requires you to build your statistical model using a separate set of data than the one used to test the accuracy of your model.

The original data is split into a **training set** that is used to build and optimize the statistical model, and a **testing set** that is used to test the resulting model's predictive accuracy (i.e., ability to accurately predict new, not-yet-seen data).

You typically want more data in the training set since this data is being used to build your model. Some proportions that are used typically are 70% training, 30% testing, but you can adjust these.

# Cross-Validation

Recall that the statistical model is built using **only the data in the training set**. To avoid overfitting, we do not want the model we build to fit too specifically to the unique error in the training set.

Cross-validation is a technique for splitting the training set multiple times into **mini training-testing sets**. The model is fit using each of the mini training sets and then its predictive accuracy is measured in each of the mini testing sets. **The results are averaged across each of these models**.

Here are a couple of cross-validation approaches:

-   **leave-one-out cross-validation**: Using the training set, fit the statistical model to n-1 observations. Then, test the predictive accuracy of the model on the single observation that's been left out. Repeat this n times.

-   **k-fold cross-validation**: Randomly split the training set into **k** chunks (aka, folds) of roughly equal size. One of the chunks is treated as the testing set. The remaining k-1 chunks are treated as the training set for building the model, and the predictive accuracy of the model is examined examined using the testing set. This process is repeated k times, with each time, a different chunk being treated as the testing set. Typically, people use values of k = 5 or k = 10 (but it can be any value up to n).

# Tuning

In machine learning, there are two types of parameters: model parameters and hyperparameters.

-   Model parameters are estimated from the data (e.g., model coefficients in linear regression)

-   Hyperparameters are values specific to some algorithms that can be "tuned" by the researcher. Different algorithms have different hyperparameters. You won't know the best values of hyperparameters prior to training a model. You have to rely on rules of thumb or try to find the best values through trial and error. This is the tuning process.

# Example: Steps of Machine Learning

We will focus on examples that use a continuous outcome variable (but you can certainly perform machine learning for categorical outcomes, too. Just make sure you use an appropriate algorithm!).

The data set we will use for this example is called `Prestige` from the `car` package. This data set contains various features (aka, variables) across a variety of occupations.

```{r}
head(Prestige) 
str(Prestige)

# see a description of the variables
# ?Prestige

# Not ideal way of handling missing data, but using for sake of today's example
Prestige <- na.omit(Prestige)
```

Let's see if we can build a model that does well at predicting **income** from the rest of the variables in the data set.

## Data Splitting

First, split the original data into a **training set** and a **testing set**. Remember, the statistical model is built using only the training set.

`createDataPartition` is used to split the original data into a training set and a testing set. The inputs include:

-   **y** = the outcome variable

-   **times** = the number of times you want to split the data

-   **p** = the percentage of the data that goes into the training set

-   **list** = FALSE gives the results in a matrix with the row numbers of the partition that you can pass back into the training data to effectively split it

```{r}
set.seed(50) 

# Split the original data set into a training set and a testing set
partition_data <- createDataPartition(Prestige$income, times = 1, p = .7, list = FALSE)

training.set <- Prestige[partition_data, ] # Training set
testing.set <- Prestige[-partition_data, ] # Testing set
```

## Building Model using Training Set

`caret` is an R package that consolidates all of the many various machine learning algorithms into an easy-to-use interface.

The `train` function is used for model training (aka, building a model). It uses the following inputs:

-   **y** = the outcome variable
    -   **y \~.** means predict y from all other variables in the dataset
-   **method** = the machine learning algorithm you want to use
-   **trControl** = the cross-validation method
-   **tuneGrid** = a data frame of the hyperparameters that you want to be evaluated during model training
-   **preProc** = any pre-processing adjustments that you want done on the predictor data

Let's build the model using an algorithm we're familiar with: linear regression model

```{r}
set.seed(21)

ml_linear_model <- train(income ~. , 
                         data = training.set, 
                         method = "lm", 
                         trControl = trainControl(method="cv", number = 5), # k-folds CV with k=5 
                        #tuneGrid = , # leaving empty because there's no hyperparameters to tune in this case
                          preProc = c("center")) 


ml_linear_model
ml_linear_model$results
```

How accurately does this model predict incomes in the testing set?

```{r}
linear_predicted <- predict(ml_linear_model, testing.set) 

# Overall accuracy assessment
postResample(linear_predicted, testing.set$income)
```

The model we built accounts for 68% of the variability in incomes in the testing set.

# Comparing Models

One of the advantages of this approach is you can easily run and compare models fit using different types of algorithms. Another commonly used machine learning model is called a Support Vector Machine (SVM).

## Support Vector Machine

Unlike linear regression, SVMs **do** have a hyperparameter that can be tuned. You can see what it's called by passing the model algorithm, called `svmLinear`, to the `getModelInfo` function.

```{r}
getModelInfo("svmLinear")$svmLinear$parameters 
```

Then, we can choose "tuning" values to try out for this hyperparameter `expand.grid` is the function that allows you to specify tuning values.

```{r}
tuning_values <- expand.grid(C = c(0.001, 0.01, 0.1, 1, 10, 100)) 
```

Now, let's build the model using the training set.

```{r}
set.seed(42)

ml_svm_model <- train(income ~. , 
                      data = training.set, 
                      method = "svmLinear", 
                      trControl = trainControl(method="cv", number = 5), 
                      tuneGrid = tuning_values, 
                      preProc = c("center"))

ml_svm_model
ml_svm_model$results
```

How accurately does this model predict incomes in the testing set?

```{r}
svm_predicted <- predict(ml_svm_model, testing.set) 

# Overall accuracy assessment
postResample(svm_predicted, testing.set$income)
```

The model we built accounts for 71% of the variability in incomes in the testing set.

## Random Forest

Let's try one more algorithm called Random Forest.

Find the hyperparameters that can be tuned.

```{r}
getModelInfo("rf")$rf$parameters
tuning_values <- expand.grid(mtry = c(2, 3, 4, 5))
```

Build the model using the training set.

```{r}
set.seed(47)

ml_rf_model <- train(income ~. , 
                  data = training.set, 
                  method = "rf", 
                  trControl = trainControl(method="cv", number = 5), 
                  tuneGrid = tuning_values,
                  preProc = c("center"))

ml_rf_model
ml_rf_model$results
```

How accurately does this model predict incomes in the testing set?

```{r}
# Testing predictive ability of model in testing set
rf_predicted <- predict(ml_rf_model, testing.set)

postResample(rf_predicted, testing.set$income) 
```

The model we built accounts for 77% of the variability in incomes in the testing set.

## Which model did the best?

```{r}
RMSE_Testing <- c(postResample(linear_predicted, testing.set$income)[1], postResample(svm_predicted, testing.set$income)[1], postResample(rf_predicted, testing.set$income)[1])

RSq_Testing <- c(postResample(linear_predicted, testing.set$income)[2], postResample(svm_predicted, testing.set$income)[2], postResample(rf_predicted, testing.set$income)[2])

Algorithm <- c("Linear Regression", "Support Vector Machine", "Random Forest")

compare_models <- data.frame(cbind(Algorithm, RMSE_Testing, RSq_Testing))

compare_models
```

# Papaja

`Papaja` is a package that allows you to produce APA-formatted reporting, tables, and graphs directly from your .rmd file.

We need to install a version of LaTeX that works with R (I chose `tinytex`), `papaja`, and `corx`.

```{r}
# tinytex::install_tinytex()
library(tinytex) # LaTeX

# install.packages("papaja")
library(papaja) # for producing APA-style output

# install.packages("citr")
library(citr) # for citing works in-text
```

# Inserting Object Values in Text

To insert a value stored within an object into text in Rmarkdown, you can write it as "r name_of_object" (replacing the quotes with back ticks; i.e., \`).

For example, if you stored the value `5` in an object called `x` you would write "r x" (replacing the quotes with back ticks) to have the number 5 inserted in the text when you knit the Rmarkdown document.

```{r}
M <- 2.37
SD <- 1.22
```

The mean of the sample was `r M` and the standard deviation was `r SD`.

# Tables

## Descriptive Statistics

First, let's get descriptive statistics for a set of variables. We can go ahead and keep using the `Prestige` data set.

```{r}
descriptives <- describe(Prestige) %>%
  select(n, mean, sd) 

colnames(descriptives) <- c("N", "Mean", "SD")

descriptives # raw table
```

Next, we format the numerical columns using `printnum()` function from `{papaja}` for this. In essence, this puts everything to 2 decimal points (you can change this), pads with zeroes when necessary, and turns those to strings so they print correctly. "-1" indicates that the first column is non-numerical so does not apply.
```{r}
descriptives <- descriptives %>%
  mutate_if(is.numeric, printnum)
```


To convert this table into APA-style, we can use the `apa_table()` function. **Make sure** in the code chunk options you set `results = "asis"` for it to render correctly. We're also going to give the code chunk a label:

```{r descrips-tab, results = "asis"}
apa_table(descriptives, 
          caption = "Descriptive Statistics", 
          note    = "This table was created with apa_table().")
```


## Regression Results

As an example, let's run a regression predicting income from the other variables in the `Prestige` data set. 
```{r}
head(Prestige)

regr_model <- lm(income ~ ., data = Prestige)
summary(regr_model)
```


The `apa_print()` function can be used to convert results from a statistical analysis into an APA-results object that can then be passed to `apa_table()`.
```{r}
apa_regr_model <- apa_print(regr_model)
```

The table is stored in the `$table` component of this object. Pass it to the `apa_table()` function. **Make sure** in the code chunk options you set `results = "asis"` for it to render correctly. We're also going to give the code chunk a label:
```{r reg-tbl-1, results = "asis"}
apa_regr_model$table %>% 
  apa_table(caption = "Regression Model Predicting Income from Education, Women, Prestige, Census and Type",
            note    = "* p < 0.05; ** p < 0.01; *** p < 0.001")
```



# Figures 

You can convert `ggplot` figures into APA-style by adding `theme_apa()`.

```{r}
ggplot(Prestige, aes(x = prestige, y = income)) +
  geom_point(alpha   = .5) +
  geom_smooth(method = "lm") +
  labs(x = "Prestige",
       y = "Income") +
  theme_apa()
```



# In-text summaries of tests

One of the most useful aspects of writing manuscripts in .Rmd and `{papaja}` is that we can report our results directly from the statistical models we run in R, reducing the likelihood of copy-and-pasting errors.

We can use `apa_print()` from `{papaja}` to report the results of our statistical tests in text. Abbreviated results can be called with `$estimate` from the `apa_print()` object. For example, `apa_regr_model$estimate$prestige` will print a nicely formatted printout of the slope for prestige on income:

> Prestige was a significant, positive predictor of income while controlling for education, percentage of women, census, and position type, `r apa_regr_model$estimate$prestige`.


We can instead get the full result (including *t*, *df*, and *p*) by calling the `$full_result` object, so `apa_regr_model$full_result$prestige` will be rendered like so:

> Prestige was a significant, positive predictor of income while controlling for education, percentage of women, census, and position type, `r apa_regr_model$full_result$prestige`.

You can also extract the R-squared value and its corresponding F-test:

> The model accounted for a significant amount of variability in people's income, (`r apa_regr_model$full_result$modelfit$r2`).
