---
title: "Introdction to GLM"
output: 
  html_document:
    toc: true
    toc_float: true
    code_folding: show
---

## Learning Objectives

### Basic

1. Define the [components](#glm-components) of the GLM
2. [Simulate data](#sim-glm) using GLM equations
3. Identify the model parameters that correspond to the data-generation parameters
4. Understand and plot [residuals](#residuals)
5. [Predict new values](#predict) using the model
6. Explain the differences among [coding schemes](#coding-schemes) 

### Intermediate

7. Demonstrate the [relationships](#test-rels) among two-sample t-test, one-way ANOVA, and linear regression
8. Given data and a GLM, [generate a decomposition matrix](#decomp) and calculate sums of squares, mean squares, and F ratios for a one-way ANOVA

## Resources

* [Jeff Miller and Patricia Haden, Statistical Analysis with the Linear Model (free online textbook)](http://www.otago.ac.nz/psychology/otago039309.pdf)
* [lecture slides introducing the General Linear Model](slides/08_glm_slides.pdf)
* [GLM shiny app](http://shiny.psy.gla.ac.uk/Dale/GLM/)
* [F distribution](http://shiny.psy.gla.ac.uk/Dale/fdist)

## Setup

You'll need the following packages. 

```{r, results = 'hide', warning = FALSE, message = FALSE}

knitr::opts_chunk$set(echo = TRUE, error = TRUE)

# libraries needed for these examples
library(tidyverse)
library(broom)

set.seed(30250) # makes sure random numbers are reproducible
```


## GLM

### What is the GLM?

The General Linear Model (GLM) a general mathematical framework for expressing relationships among variables that can express or test linear relationships between a numerical **dependent variable** and any combination of categorical or continuous **independent variables**.

### Components {#glm-components}

There are some mathematical conventions that you need to learn to understand the equations representing linear models. Once you understand those, learning about the GLM will get much easier.


| Component of GLM | Notation                      |
|------------------|-------------------------------|
| Dependent Variable (DV) | $Y$                    |
| Grand Average    | $\mu$ (the Greek letter "mu") |
| Main Effects     | $A, B, C, \ldots$             |
| Interactions     | $AB, AC, BC, ABC, \ldots$     |
| Random Error     | $S(Group)$                    |


The linear equation predicts the dependent variable ($Y$) as the sum of the grand average value of $Y$ ($\mu$, also called the intercept), the main effects of all the predictor variables ($A+B+C+ \ldots$), the interactions among all the predictor variables ($AB, AC, BC, ABC, \ldots$), and some random error ($S(Group)$). The equation for a model with two predictor variables ($A$ and $B$) and their interaction ($AB$) is written like this:

$Y$ ~ $\mu+A+B+AB+S(Group)$

Don't worry if this doesn't make sense until we walk through a concrete example.

### Simulating data from GLM {#sim-glm}

A good way to learn about linear models is to simulate data where you know exactly how the variables are related, and then analyse this simulated data to see where the parameters show up in the analysis.

We'll start with a very simple linear model that just has a single categorical factor with two levels. Let's say we're predicting reaction times for congruent and incongruent trials in a Stroop task for a single participant. Average reaction time (`mu`) is 800ms, and is 50ms faster for congruent than incongruent trials (`effect`). 

A **factor** is a categorical variable that is used to divide subjects into groups, usually to draw some comparison. Factors are composed of different **levels**. Do not confuse factors with levels!

You need to represent categorical factors with numbers. The numbers, or **coding** you choose will affect the numbers you get out of the analysis and how you need to interpret them. Here, we will **effect code** the trial types so that congruent trials are coded as +0.5, and incongruent trials are coded as -0.5.

A person won't always respond exactly the same way. They might be a little faster on some trials than others, due to random fluctuations in attention, learning about the task, or fatigue. So we can add an **error term** to each trial. We can't know how much any specific trial will differ, but we can characterise the distribution of how much trials differ from average and then sample from this distribution. 

Here, we'll assume the error term is sampled from a normal distribution with a standard deviation of 100 ms (the mean of the error term distribution is always 0). We'll also sample 100 trials of each type, so we can see a range of variation.

So first create variables for all of the parameters that describe your data.

```{r}
n_per_grp <- NULL
mu <- NULL # average RT
effect <- NULL # average difference between congruent and incongruent trials
error_sd <- NULL # standard deviation of the error term
trial_types <- NULL # effect code
```

Then simulate the data by creating a data table with a row for each trial and columns for the trial type and the error term (random numbers samples from a normal distribution with the SD specified by `error_sd`). For categorical variables, include both a column with the text labels (`trial_type`) and another column with the coded version (`trial_type.e`) to make it easier to check what the codings mean and to use for graphing. Calculate the dependent variable (`RT`) as the sum of the grand mean (`mu`), the coefficient (`effect`) multiplied by the effect-coded predictor variable (`trial_type.e`), and the error term.

```{r}
dat <- NULL
```


The `!!!` (triple bang) in the code `recode(trial_type, !!!trial_types)` is a way to expand the vector `trial_types <- c("congruent" = 0.5, "incongruent" = -0.5)`. It's equivalent to `recode(trial_type, "congruent" = 0.5, "incongruent" = -0.5)`. This pattern avoids making mistakes with recoding because there is only one place where you set up the category to code mapping (in the `trial_types` vector).

Last but not least, always plot simulated data to make sure it looks like you expect.

```{r plot-sim, fig.cap = "Simulated Data"}
ggplot(dat)
```

### Linear Regression

Now we can analyse the data we simulated using the function `lm()`. It takes the formula as the first argument. This is the same as the data-generating equation, but you can omit the error term (this is implied), and takes the data table as the second argument. Use the `summary()` function to see the statistical summary.

```{r}
my_lm <- NULL

summary(my_lm)
```

Notice how the **estimate** for the `(Intercept)` is close to the value we set for `mu` and the estimate for `trial_type.e` is close to the value we set for `effect`.

Change the values of `mu` and `effect`, resimulate the data, and re-run the linear model. What happens to the estimates?

### Residuals {#residuals}

You can use the `residuals()` function to extract the error term for each each data point. This is the DV values, minus the estimates for the intercept and trial type. We'll make a density plot of the residuals below and compare it to the normal distribution we used for the error term.

```{r res-density-plot, fig.cap="Model residuals should be approximately normally distributed for each group"}

res <- NULL

ggplot(dat) + 
  # creates a normal distribution curve
  stat_function(aes(0), color = "grey60",
                fun = dnorm, n = 101,
                args = list(mean = 0, sd = error_sd)) +
  geom_density(aes(res, color = trial_type))
  
```


You can also compare the model residuals to the simulated error values. If the model is accurate, they should be almost identical. If the intercept estimate is slightly off, the points will be slightly above or below the black line. If the estimate for the effect of trial type is slightly off, there will be a small, systematic difference between residuals for congruent and incongruent trials.

```{r res-err-plot, fig.cap="Model residuals should be very similar to the simulated error"}

ggplot(dat) +
  geom_abline(slope = 1) +
  geom_point(aes(error, res, color = trial_type)) +
  ylab("Model Residuals") +
  xlab("Simulated Error")
```


What happens to the residuals if you fit a model that ignores trial type (e.g., `lm(Y ~ 1, data = dat)`)?


### Predict New Values {#predict}

You can use the estimates from your model to predict new data points, given values for the model parameters. For this simple example, we just need to now the trial type to make a prediction.

For congruent trials, you would predict that a new data point would be equal to the intercept estimate plus the trial type estimate multiplied by 0.5 (the effect code for congruent trials).

```{r}
int_est <- NULL
tt_est  <- NULL
tt_code <- NULL

new_congruent_RT <- NULL

new_congruent_RT
```

You can also use the `predict()` function to do this more easily. The second argument is a data table with columns for the factors in the model and rows with the values that you want to use for the prediction.

```{r}
newdata <- NULL
predict(my_lm, newdata)
```


If you look up this function using `?predict`, you will see that "The function invokes particular methods which depend on the class of the first argument."

What this means is that `predict()` works differently depending on whether you're predicting from the output of `lm()` or other analysis functions. You can search for help on the lm version with `?predict.lm`.



### Coding Categorical Variables {#coding-schemes}

In the example above, we used **effect coding** for trial type. You can also use **sum coding**, which assigns +1 and  -1 to the levels instead of +0.5 and -0.5.  More commonly, you might want to use **treatment coding**, which assigns 0 to one level (usually a baseline or control condition) and 1 to the other level (usually a treatment or experimental condition).

Here we will add sum-coded and treatment-coded versions of `trial_type` to the dataset using the `recode()` function.

```{r}

dat <- NULL

```

If you define named vectors with your levels and coding, you can use them with the `recode()` function if you expand them using `!!!`.

```{r}

tt_sum <- NULL
tt_tr  <- NULL

dat <- NULL
)
```

Here are the coefficients for the effect-coded version. They should be the same as those from the last analysis.

```{r}
lm(RT ~ trial_type.e, data = dat)$coefficients
```

Here are the coefficients for the sum-coded version. This give the same results as effect coding, except the estimate for the categorical factor will be exactly half as large, as it represents the difference between each trial type and the hypothetical condition of 0 (the overall mean RT), rather than the difference between the two trial types.

```{r}
lm(RT ~ trial_type.sum, data = dat)$coefficients
```

Here are the coefficients for the treatment-coded version. The estimate for the categorical factor will be the same as in the effect-coded version, but the intercept will decrease. It will be equal to the intercept minus the estimate for trial type from the sum-coded version.

```{r}
lm(RT ~ trial_type.tr, data = dat)$coefficients
```


## Relationships among tests {#test-rels}

### T-test

The t-test is just a special, limited example of a general linear model.

```{r}
t.test(RT ~ trial_type.e, data = dat, var.equal = TRUE)
```

What happens when you use other codings for trial type in the t-test above? Which coding maps onto the results of the t-test best?

```{r}

```


### ANOVA

ANOVA is also a special, limited version of the linear model.

```{r}
my_aov <- NULL

summary(my_aov, intercept = TRUE)
```

The easiest way to get parameters out of an analysis is to use the `broom::tidy()` function. This returns a tidy table that you can extract numbers of interest from. Here, we just want to get the F-value for the effect of trial_type. Compare the square root of this value to the t-value from the t-tests above.

```{r}
f <- broom::tidy(my_aov)$statistic[1]
sqrt(f)
```

## Understanding ANOVA

We'll walk through an example of a one-way ANOVA with the following equation:

$Y_{ij} = \mu + A_i + S(A)_{ij}$

This means that each data point ($Y_{ij}$) is predicted to be the sum of the grand mean ($\mu$), plus the effect of factor A ($A_i$), plus some residual error ($S(A)_{ij}$).

### Means, Variability, and Deviation Scores

Let's create a simple simulation function so you can quickly create a two-sample dataset with specified Ns, means, and SDs.

```{r}
two_sample <- function() {
  
}
```

Now we will use `two_sample()` to create a dataset `dat` with N=5 per group, means of -2 and +2, and SDs of 1 and 1 (yes, this is an effect size of d = 4).

```{r}
dat <- NULL
```

You can calculate how each data point (`Y`) deviates from the overall sample mean ($\hat{\mu}$), which is represented by the horizontal grey line below and the deviations are the vertical grey lines. You can also calculate how different each point is from its group-specific mean ($\hat{A_i}$), which are represented by the horizontal coloured lines below and the deviations are the coloured vertical lines.

```{r deviation, echo = FALSE, fig.cap="Deviations of each data point (Y) from the overall and group means"}

group_by(dat, grp) %>% 
  mutate(Ymean = mean(Y)) %>%
  ungroup() %>%
  mutate(i = row_number()) %>%
  ggplot(aes(i, Y, color = grp)) +
    geom_hline(yintercept = mean(dat$Y), color = "grey40") +
    geom_linerange(aes(x = i-0.025, ymin = Y, ymax = mean(dat$Y)), color = "grey40") +
    geom_linerange(aes(x = i-0.025, ymin = mean(dat$Y), ymax = Y), color = "grey40") +
    geom_hline(aes(yintercept = Ymean, color = grp)) +
    geom_linerange(aes(x = i+0.025, ymin = Y, ymax = Ymean, group = i, color = grp)) +
    geom_linerange(aes(x = i+0.025, ymin = Ymean, ymax = Y, group = i, color = grp)) +
    geom_point(size = 3) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
  xlab("")
  
```

You can use these deviations to calculate variability between groups and within groups. ANOVA tests whether the variability between groups is larger than that within groups, accounting for the number of groups and observations.

### Decomposition matrices {#decomp}

We can use the estimation equations for a one-factor ANOVA to calculate the model components.

* `mu` is the overall mean 
* `a` is how different each group mean is from the overall mean
* `err` is residual error, calculated by subtracting `mu` and `a` from `Y`

This produces a *decomposition matrix*, a table with columns for `Y`, `mu`, `a`, and `err`.

```{r}
decomp <- NULL
```

`r knitr::kable(decomp)`

Calculate sums of squares for `mu`, `a`, and `err`.

```{r}
SS <- NULL

SS
```


If you've done everything right, `SS$mu + SS$a + SS$err` should equal the sum of squares for Y.

```{r}
SS_Y <- sum(decomp$Y^2)
all.equal(SS_Y, SS$mu + SS$a + SS$err)
```

Divide each sum of squares by its corresponding degrees of freedom (df) to calculate mean squares. The df for `mu` is 1, the df for factor `a` is `K-1` (K is the number of groups), and the df for `err` is `N - K` (N is the number of observations).

```{r}
K <- NULL
N <- NULL
df <- NULL
MS <- NULL

MS
```


Then calculate an F-ratio for `mu` and `a` by dividing their mean squares by the error term mean square. Get the p-values that correspond to these F-values using the `pf()` function.

```{r}
F_mu <- NULL
F_a  <- NULL
p_mu <- NULL
p_a  <- NULL
```

Put everything into a data frame to display it in the same way as the ANOVA summary function.

```{r}
my_calcs <- data.frame(
  term = c("Intercept", "grp", "Residuals"),
  Df = df,
  SS = c(SS$mu, SS$a, SS$err),
  MS = c(MS$mu, MS$a, MS$err),
  F = c(F_mu, F_a, NA),
  p = c(p_mu, p_a, NA)
)

my_calcs
```


Now run a one-way ANOVA on your results and compare it to what you obtained in your calculations.

```{r}
aov(Y ~ grp, data = dat) %>% summary(intercept = TRUE)
```

Using the code above, write your own function that takes a table of data and returns the ANOVA results table like above.

