```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, rows.print = 10)

# suppress scientific notation
options(scipen = 999)
```

You can download the Rmd file [here](lab-3.Rmd) to follow along.

# Purpose

Today we will review univariate regression, discuss how to summarize and visualize uncertainty in regression models, and discuss how to estimate regression coefficients using matrix algebra. At the end, we will introduce the General Linear Model and demonstrate how GLM can be used to understand all of the statistical tests we have learned so far (*t*-tests, ANOVA, correlations, regressions) within one unifying framework. 

For today's lab, you will need to load the following libraries:

```{r, message = FALSE, warning = FALSE}
library(tidyverse) # for plotting and data wrangling
library(rio) # for importing data
library(broom) # for cleaning
library(sjPlot) # for plotting
library(ggpubr) # for plotting
library(carData) # for Guyer dataset
```

***

# Visualizing uncertainty

For today's lab, we are going to continue using the same dataset from previous labs that contains variables measuring gender, self-reported conscientiousness (from the BFI), and self-reported physical health.

## Data and review

First, let's load in the data.

```{r}
health <- import("https://raw.githubusercontent.com/uopsych/psy612/master/labs/lab-3/data/consc_health.csv")
```

Our model can be described with the following regression equation:

$$health_i = b_0 + b_1consc_i + e_i$$

We can run the model in R using the `lm()` function.

```{r}
model <- lm(sr_health ~ consc, data = health)
summary(model)
```

> **Question:** What do the intercept and slope mean? What do the *p*-values tell us? 

## Confidence intervals

Our `b's` (intercept and slope) are *estimates* from our sample of true population parameters ($\beta$'s). Whenever we calculate an estimate, we should also determine how precise our estimate is. 

Recall the formula for calculating confidence intervals:

$$CI_b = b \pm CV(SE_b)$$

We can calculate the confidence interval around the estimates in R with the function `stats::confint()`. This function takes the model object as the first argument. By default it will give you 95% CI's. 

```{r}
confint(model)
```


>**Question:** What does these 95% CI for the slope of conscientiousness mean in plain English? 


## Confidence bands

In addition to estimating precision around the our coefficients, we can also estimate our precision around each fitted value, $\hat{Y_i}$. These standard errors are generated by `broom::augment()` (and are labeled `.se.fit`).

```{r}
model %>% 
  augment(se_fit = TRUE) %>% 
  select(sr_health, .fitted, .se.fit) %>% 
  head() 
```

If we were to string all of this information together, it would generate a confidence **band** around our regression line. As we've already seen, you can get this confidence band when creating a scatter plot by adding `geom_smooth(method = "lm")`. 

```{r, message = FALSE}
health %>%
  ggplot(aes(x = consc, y = sr_health)) + 
  geom_point() + 
  geom_smooth(method = "lm") + # adds a layer that includes the regression line & 95% confidence band
  labs(x = "Conscientiousness", y = "Self-rated health") +
  theme_minimal()
```


The animation below is an example of a ["Hypothetical Outcomes Plot"](https://github.com/wilkelab/ungeviz){target="_blank"} (HOP) that visually demonstrates what this 95% CI band represents. In essence, this plot shows what the regression line could look like if we were to repeat our experiment over and over (sampling from the same population each time). 

```{r echo=FALSE, warning = FALSE, message = FALSE}
library(ungeviz)
library(gganimate)
library(transformr)
library(gifski)

set.seed(012220)

boots <- bootstrapper(100)

p <- health %>%
  ggplot(aes(x = consc, y = sr_health)) +
  geom_smooth(method = "lm", color = NA) +
  geom_point(alpha = 0.3) +
  geom_smooth(data = boots, method = "lm", fullrange = TRUE, se = FALSE) +
  theme_minimal() +
  labs(x = "Conscientiousness", y = "Self-rated health") +
  transition_states(.draw, 1, 1) +
  enter_fade() +
  exit_fade() +
  ease_aes()

animate(p, fps = 3)

```


## Prediction bands

The `predict()` function allows us to get the fitted `Y` values from all 60 conscientiousness values in our dataset.  Note: you can also get the fitted `Y` values with the `broom::augment()` function like we did earlier in the lab.


```{r}
predict(model)
```

We can use this information to create "prediction bands". First we will generate our fitted values along with a "prediction interval" (lower and upper bound) for each of these values. 

```{r, warning = FALSE}
fitted_interval <- predict(model, interval = "prediction")
head(fitted_interval)
```

Next we'll bind these fitted values (and their intervals) to our original dataset. 

```{r}
new_df <- cbind(health, fitted_interval)
head(new_df)
```

And finally, we'll plot a prediction band on top of the data by adding a `geom_ribbon()` layer.

```{r, message = FALSE}
new_df %>% 
  ggplot(aes(x = consc, y = sr_health)) +
  geom_point() +
  geom_smooth(method = lm, se = TRUE) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "blue", alpha = 0.1) + 
  labs(x = "Conscientiousness", y = "Self-rated health") +
  theme_minimal()
```

>**Question:** What does the prediction band represent? How does it differ from the confidence band? 


## Other visualization tools

sjPlot::plot_model()`

```{r, eval=FALSE}
plot_model(model = model,    # name of model object
           type = "pred",    # show predicted values (i.e. regression line)
           show.data = TRUE, # include data points on plot
           jitter = TRUE)    # add small amount of random variation to  to prevent overlap
```


`ggpubr::ggscatter()`

```{r, message = FALSE}
ggscatter(data = health,              # name of data.frame
          x = "consc",                # IV (must be quoted)
          y = "sr_health",            # DV (must be quoted)
          add = "reg.line",           # add regression line
          xlab = "Conscientiousness", # x-axis label
          ylab = "Self-rated health", # y-axis label
          conf.int = TRUE,            # show 95% confidence band around regression line
          cor.coef = TRUE)            # display correlation coefficient and p-value
```

***

# Regression with matrix algebra

Consider the following matrix algebra equation:

$$\mathbf{Y} = \mathbf{XB}$$
In this equation:
1. `Y` is a nx1 matrix of outcome variable values from the dataset (e.g., health scores) <br>
2. `X` is a nx2 matrix. The first column is filled with `1`s and the second column is the predictor variable values from the dataset (e.g., consciousness scores) <br>
3. `B` is a 2x1 matrix containing our estimates of the intercept and slope. <br> 
Using this equation, we can plug in our data and solve for the intercept and slope estimates.

## Example 

1. Create the `X` matrix (first column 1s, second column consciousness scores)

```{r}
x_mat <- health %>% # start with the original data frame
  mutate(ones = rep(1, nrow(.))) %>% # create a column of 1's to represent the intercept
  select(ones, consc) %>% # select only the column of 1's and X variable
  as.matrix() %>%  # coerce to a matrix
  unname() # get rid of dimnames 

head(x_mat)
```

2. Create the **`Y`** matrix (self-reported health scores)


```{r}
y_mat <- health %>% # start with the original data frame
  select(sr_health) %>% # select just the Y variable
  as.matrix() %>% # coerce to a matrix
  unname() # get rid of dimnames 

head(y_mat)
```


In matrix algebra, you can't divide by a matrix, so you will need to do some more complicated arithmetic to solve for B (the vector of slope and intercept estimates). You can calculate B with the following formula: 

$$\mathbf{B} = (\mathbf{X'X})^{-1} \mathbf{X'Y}$$

```{r}
solve(t(x_mat) %*% x_mat) %*% (t(x_mat) %*% y_mat)
```

Let's check our work:

```{r}
model$coefficients
```

***

# The General Linear Model

We just saw that regression works "under the hood" by solving a matrix algebra equation to get the intercept and slope of the regression line. As Sara mentioned in class, this matrix algebra equation works for *any* type of data we have. This should clue us into the idea that there is some fundamental process going on behind the scenes of all of our models...

The **general linear model (GLM)** is a family of models that assume the relationship between your DV and IV(s) is additive, and that your outcome is normally distributed. In its simplest form, we can think of the general linear model as: 

$$Data = Model + Error $$ 

This provides a unifying framework for all of the statistical tests we have learned (or at least touched on) so far: *t*-tests, correlations, ANOVA, and linear regression. All of these tests are really doing the same math at the end of the day.

To illustrate this, let's look at the relationship between gender and self-reported health from the `health` data set.


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

Let's covert the `gender` variable to `0's` and `1's`. We will talk more about dummy coding in a future lab. 


```{r}
health <-  health %>% 
  mutate(gender = case_when(gender == "male" ~ 0,
                               gender == "female" ~ 1))
head(health)
```

## *t*-test

```{r}
t_test <- t.test(formula = sr_health ~ gender, data = health, var.equal = TRUE)
t_test
```

## Correlation 

```{r}
cor_test <- cor.test(formula = ~ sr_health + gender, data = health)
cor_test
```

## ANOVA

```{r}
anova_test <- aov(formula = sr_health ~ gender, data = health)
summary(anova_test)

```

## Regression

```{r}
regression <- lm(formula = sr_health ~ gender, data = health)
summary(regression)
```

## Comparison

```{r}
t_test$p.value
cor_test$p.value
anova(anova_test)
summary(regression)
```

***

# Minihacks

## Minihack 1: Matrix algebra

Load the following data from Fox and Guyer’s (1978) anonymity and cooperation study. In the study, twenty groups of four participants each played 30 trials of the the prisoner’s dilemma game. The number of cooperative choices (`cooperation`) made by each group were scored out of 120 (i.e., cooperative choices made by 4 participants over 30 trials). The groups either made decisions publicly or privately (`condition`).

```{r}
guyer <- carData::Guyer
```

Use matrix algebra to manually calculate the intercept and slope for the regression equation:

$$cooperation_i = b_0 + b_1condition_i + e_i$$
1. Dummy code the condition (public = 0 and anonymous = 1).

```{r}
#your code here
```

2. Create the `X` matrix.

```{r}
#your code here
```

3. Create the `Y` matrix.

```{r}
#your code here
```

4. Solve for the `B` matrix.

```{r}
#your code here
```


5. Confirm your results using the `lm()` function.

```{r}
#your code here
```

***

## Minihack 2: Prediction and confidence bands

Using the `mtcars` dataset, create a scatterplot of cyl (x-axis) and mpg (y-axis) that includes both a prediction band and a confidence band.

1. Run a regression model.

```{r}
#your code here
```

2. Generate fitted values and their prediction interval using the `predict()` function.
```{r}
#your code here
```

3. Bind the fitted values and the prediction intervals to the original dataset. 

```{r}
#your code here
```

4. Using ggplot(), create a scatterplot that includes both a prediction band and a confidence band.

```{r}
#your code here
```

***

## Minihack 3: Confidence intervals

For this minihack, we will refer back to the example about conscientiousness and health. We used `confint()` to calculate the 95% CI for our regression coefficients. Your job is start with the model output, stored as a list object (see below), and extract the relevant pieces of information to calculate the 95% CI around the intercept and slope.

Note: the formula for a confidence interval is estimate +- cv*se

```{r}
model <- lm(sr_health ~ consc, data = health)
summary_model <- summary(model)
```

1. Extract the estimates.

```{r}
#your code here
```

2. Extract standard errors. Hint: use summary_model$coefficients.

```{r}
#your code here
```

3. Extract the df (It will be the denominator df from F statistic)).

```{r}
#your code here
```

4. Calculate the critical value.
```{r}
#your code here
```

5. Calculate a 95% CI for both estimates.

```{r}
#your code here
```


6. Verify that your answer corresponds to the result from `confint()`. 

```{r}
#your code here
```
