```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      eval = TRUE,
                      size = "tiny",
                      message = FALSE,
                      warning = FALSE,
                      fig.align = "center",
                      out.width = "80%")

def.chunk.hook  <- knitr::knit_hooks$get("chunk")
knitr::knit_hooks$set(chunk = function(x, options) {
  x <- def.chunk.hook(x, options)
  ifelse(options$size != "normalsize", paste0("\n \\", options$size,"\n\n", x, "\n\n \\normalsize"), x)
})
```

```{r packages, include=FALSE}
devtools::load_all()
library(tidyverse)
library(kableExtra)
library(patchwork)
library(here)
```

```{r functions, include = FALSE}
funs <- get_funs(here("R", "utils.R"))
```

## Outline

<!-- \tableofcontents[hideallsubsections] -->
\scriptsize
\tableofcontents

# Beyond the Gaussian distribution

## Quick recap about Gaussian distribution

- The Gaussian distribution is part of the Exponential family
- It is defined with mean ($\mu$) and the standard deviation ($\sigma$) that are independent
- It is symmetric with the same value for mean, mode and median
- The support is $[- \infty, + \infty]$

The Probability Density Function (PDF) is:

$$
f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x - \mu}{\sigma})^2}
$$

## Quick recap about Gaussian distribution

```{r, echo = FALSE}
ggnorm(0, 1) + 
    mytheme() +
    ggtitle(latex2exp::TeX("$\\mu = 0$, $\\sigma = 1$"))
```

\begin{center}
But not always gaussian-like variables!
\end{center}

## Quick recap about Gaussian distribution

In fact, in Psychology, variables do not always satisfy the properties of the Gaussian distribution. For example:

- Reaction times
- Accuracy
- Percentages or proportions
- Discrete counts
- Likert scales
- ...

## Reaction times

Measuring **reaction times during a cognitive task**. Non-negative and probably skewed data.

```{r, echo = FALSE}
dat <- data.frame(
    x = rgamma(1e5, 9, scale = 0.5)*100
) 

dat |> 
    ggplot(aes(x = x)) +
    geom_histogram(fill = "lightblue",
                   color = "black") +
    xlab("Reaction Times (ms)") +
    ylab("Count") +
    mytheme()
```

## Binary outcomes

Counting the number of people passing the exam out of the total. Discrete and non-negative. A series of binary (i.e., *bernoulli*) experiments.

```{r, echo = FALSE}
dat <- data.frame(y = c(70, 30), x = c("Passed", "Failed"))

dat |> 
    ggplot(aes(x = x, y = y)) +
    geom_col(color = "black",
             fill = "lightblue") +
    ylab("%") +
    ylim(c(0, 100)) +
    mytheme() +
    theme(axis.title.x = element_blank()) +
    ggtitle("Statistics Final Exam (n = 100)")
```

## Counts

Counting the number of new hospitalized patients during one month in different cities. Discrete and non-negative values.

```{r, echo = FALSE}
dat <- data.frame(x = rpois(1e5, 15))

dat |> 
    ggplot(aes(x = x)) +
    geom_bar(fill = "lightblue",
             color = "black") +
    mytheme() +
    scale_x_continuous(breaks = seq(0, 50, 5)) +
    xlab("Number of new patients during one month") +
    ylab("Count") +
    mytheme()
```

## Question...

`r question("Should we use a linear model for these variables?")`

## Should we use a linear model for these variables?

Let's try to fit a linear model on the probability of passing the exam ($N = 50$) as a function of the hours of study:

::: columns
:::: column
```{r, echo = FALSE}
# y = number of exercises solved in 1 semester
# x = percentage of attended lectures

n <- 50
x <- round(runif(n, 0, 100))
b0 <- 0.01
b1 <- 0.08
y <- rbinom(length(x), 1, plogis(qlogis(b0) + b1*x))

dat <- data.frame(student = 1:n, study_hours = x, passing = y)

dat |> 
    trim_df() |> 
    mytab()
```
::::
:::: column
```{r, echo = FALSE}
n <- length(dat$passing)
data.frame(n = n, npassing = sum(dat$passing), nfailing = n - sum(dat$passing), ppassing = mean(dat$passing)) |> mytab()
```

::::
:::

## Should we use a linear model for these variables?

Let's plot the data:

```{r, echo = FALSE}
exam_plot <- dat |>
    ggplot(aes(x = study_hours, y = passing)) +
    geom_hline(yintercept = 0, linetype = "dashed", col = "red") +
    geom_point(size = 3,
               alpha = 0.5,
               position = position_jitter(height = 0.03)) +
    mytheme() +
    xlab("Hours of study") +
    ylab("Passing the exam")
exam_plot
```

## Should we use a linear model for these variables?

Let's fit a linear model `passing ~ study_hours` using `lm`:

```{r, echo = FALSE}
exam_plot +
    geom_smooth(method = "lm", se = F)
```

`r center(color(bold("Do you see something strange?"), "myblue"))`

## Should we use a linear model for these variables?

A little **spoiler**, the relationship should be probably like this:

```{r, echo = FALSE}
exam_plot +
    stat_smooth(method = "glm", 
                se = FALSE,
                method.args = list(family = binomial))
```

## Should we use a linear model for these variables?

Another example, the number of solved exercises in a semester as a function of the number of attended lectures ($N = 100$):

::: columns
:::: column
```{r, echo = FALSE}
# y = number of exercises solved in 1 semester
# x = percentage of attended lectures

n <- 100
x <- round(runif(n, 0, 63))
y <- rpois(n, exp(0 + 0.05*x))

dat <- data.frame(student = 1:n, attended_lectures = x, nsolved = y)

dat |> 
    trim_df() |> 
    mytab()
```
::::
:::: column
```{r echo = FALSE}
dat |> 
    ggplot(aes(x = nsolved)) +
    geom_bar() +
    mytheme()
```
::::
:::

## Should we use a linear model for these variables?

```{r, echo = FALSE}
exam_plot <- dat |> 
    ggplot(aes(x = attended_lectures, y = nsolved)) +
    geom_hline(yintercept = 0, linetype = "dashed", col = "red") +
    geom_point(size = 3) +
    mytheme() +
    xlab("Number of attended lectures") +
    ylab("Solved exercises")
exam_plot
```

## Should we use a linear model for these variables?

Again, fitting the linear model seems partially appropriate but there are some problems:

```{r, echo = FALSE}
exam_plot +
    geom_smooth(method = "lm", se = FALSE)
```

## Should we use a linear model for these variables?

Also the residuals are quite problematic:

```{r, echo = FALSE}
fit <- lm(nsolved ~ attended_lectures, data = dat)

dfit <- data.frame(
    fitted = fitted(fit),
    residuals = residuals(fit)
)

qqn <- dfit |> 
    ggplot(aes(sample = residuals)) + 
    stat_qq() + 
    stat_qq_line() +
    xlab("Theoretical Quantiles") +
    ylab("Residuals") +
    mytheme()

res_fit <- dfit |> 
    ggplot(aes(x = fitted, y = residuals)) +
    geom_point() +
    mytheme() +
    ylab("Residuals") +
    xlab("Fitted") +
    geom_smooth(se = FALSE, color = "red")
qqn + res_fit
```

## Should we use a linear model for these variables?

Another little spoiler, the model should consider both the support of the `y` variable and the non-linear pattern. Probably something like this:

```{r, echo = FALSE}
exam_plot +
    stat_smooth(method = "glm", 
                se = FALSE,
                method.args = list(family = poisson))
```

## So what?

Both linear models somehow capture the expected relationship but there are serious fitting problems:

- impossible predictions
- poor fitting for non-linear patterns

As a general rule in ~~life~~ statistics:

`r quote("\\textbf{All models are wrong, some are useful.}", "George Box")`

## We need a new class of models...

- We need that our model take into account the **features of our response variable**
- We need a model that, **with appropriate transformation**, keep **properties of standard linear models**
- We need a model that is **closer to the true data generation process**

\begin{center}
\textcolor{myblue}{Let's switch to Generalized Linear Models!}
\end{center}

# Generalized Linear Models

## Main references

For a detailed introduction about GLMs

- Chapters: 1 (intro), 4 (GLM fitting), 5 (GLM for binary data)

```{r, echo = FALSE, out.width="30%"}
knitr::include_graphics("img/agresti2015-foundations-lm-glm.jpg")
```

## Main references

For a basic and well written introduction about GLM, especially the Binomial GLM

- Chapters: 3 (intro GLMs), 4-5 (Binomial Logistic Regression)

```{r, echo = FALSE, out.width="50%"}
knitr::include_graphics("img/agresti2019-intro-to-categorical.jpg")
```

## Main references

Great resource for interpreting Binomial GLM parameters:

- Chapters: 13-14 (Binomial Logistic GLM), 15 (Poisson and others GLMs)
```{r, echo = FALSE, out.width="30%"}
knitr::include_graphics("img/gelman2020-reg-and-other-stories.jpg")
```

## Main references

Detailed GLMs book. Very useful especially for the diagnostic part:

- Chapters: 8 (intro), 9 (Binomial GLM), 10 (Poisson GLM and overdispersion)

```{r, echo = FALSE, out.width="30%"}
knitr::include_graphics("img/dunn2018-glm.jpg")
```

## Main references

The holy book :)

- Chapters: 14 and 15

```{r, echo = FALSE, out.width="30%"}
knitr::include_graphics("img/fox2015-applied-glm.jpg")
```

## Main references

Another good reference...

- Chapters: 8

```{r, echo = FALSE, out.width="30%"}
knitr::include_graphics("img/faraway2016-extending-glm.jpg")
```

## General idea

- models that assume **distributions other than the normal distributions**
- models that considers **non-linear relationships**
- models that allow **heteroscedasticity**

## Recipe for a GLM

- **Random Component**
- **Systematic Component**
- **Link Function**

## Random Component

The **random component** of a GLM identify the response variable $Y$ and the appropriate probability distribution. For example for a numerical and continuous variable we could use a Normal distribution (i.e., a standard linear model). For a discrete variable representing counts of events we could use a Poisson distribution, etc.

```{r, echo = FALSE, out.width="50%"}
par(mfrow = c(1,3))

curve(dnorm(x), -4, 4, main = "Normal", ylab = "density(x)", cex.lab = 1.5, cex.main = 1.5)
plot(0:10, dbinom(0:10, 10, 0.5), type = "h", main = "Binomial", ylab = "density(x)",
     cex.lab = 1.5, cex.main = 1.5, xlab = "x")
points(0:10, dbinom(0:10, 10, 0.5), pch = 19)
plot(0:20, dpois(0:20, 8), type = "h", main = "Poisson", ylab = "density(x)",
     cex.lab = 1.5, cex.main = 1.5,
     xlab = "x")
```

## Systematic Component

The **systematic component** or *linear predictor* ($\eta$) of a GLM is the combination of explanatory variables i.e. $\beta_0 + \beta_1x_1 + ... + \beta_px_p$.

\begin{align*}
\eta = \beta_0 + \beta_1x_1 + ... + \beta_px_p
\end{align*}

When the **link function** (see next slide) is used, the relationship between $\eta$ and the expected value $\mu$ of the **random component** is linear (as in standard linear models)

## Link Function

The **link function** $g(\mu)$ is an **invertible** function that connects the expected value (i.e., the mean $\mu$) of the probability distribution (i.e., the random component) with the *linear combination* of predictors $g(\mu) = \beta_0 + \beta_1x_1 + ... + \beta_px_p$. The inverse of the link function $g^{-1}$ map the linear predictor ($\eta$) into the original scale.

\begin{align*}
g(\mu) = \beta_0 + \beta_1x_1 + ... + \beta_px_p \\
\mu = g^{-1}(\beta_0 + \beta_1x_1 + ... + \beta_px_p) 
\end{align*}

Thus, the relationship between $\mu$ and $\eta$ is linear only when the **link function ** is applied i.e. $g(\mu) = \eta$. 

## Link function

The simplest **link function** is the **identity link** where $g(\mu) = \mu$ and correspond to the standard linear model. In fact, the linear regression is just a GLM with a **Gaussian random component** and the **identity** link function.

There are multiple **random components** and **link functions** for example with a 0/1 binary variable the usual choice is using a **Binomial** random component and the **logit** link function.

```{r, echo = FALSE}
fam <- c("gaussian", "binomial", "binomial", "poisson")
link <- c("identity", "logit", "probit", "log")
range <- c("$(-\\infty,+\\infty)$", "$\\frac{0, 1, ..., n_{i}}{n_{i}}$",
           "$\\frac{0, 1, ..., n_{i}}{n_{i}}$",
           "$0, 1, 2, ...$")

linktab <- data.frame(Family = fam, Link = link, Range = range)
linktab |> collapse_rows_df(Family) |> mytab()
```

# Relevant distributions

## Binomial distribution

The probability of having $k$ success (e.g., 0, 1, 2, etc.) out of $n$ trials with a probability of success $p$ is:

$$
f(n, k, p) = Pr(X = k) = \binom{n}{k} p^k(1 - p)^{n - k}
$$

The $np$ is the mean of the binomial distribution and $np(1 - p)$ is the variance.

## Bernoulli distribution

The **binomial** distribution is just a repetition of $k$ **Bernoulli** trials. A single Bernoulli trial is:

\begin{align*}
f(x, p) = p^x (1 - p)^{1 - x} \\
x \in \{0, 1\}
\end{align*}

The mean is $p$ and the variance is $p(1 - p)$

## Bernoulli and Binomial

The simplest situation for a Bernoulli trial is a coin flip. In R:

::: columns
:::: column

```{r}
n <- 1
p <- 0.7
rbinom(1, n, p) # a single bernoulli trial

n <- 10
rbinom(10, 1, p) # n bernoulli trials

rbinom(1, n, p) # binomial version
```

::::
:::: column
```{r, echo = FALSE, out.width="90%"}
n <- 30
p <- 0.7
dat <- data.frame(k = 0:n)
dat$y <- dbinom(dat$k, n, p)

dat |> 
    ggplot(aes(x = k, y = y)) +
    geom_point() +
    geom_segment(aes(x = k, xend = k, y = 0, yend = y)) +
    mytheme() +
    ylab("dbinom(x, n, p)") +
    ggtitle(latex2exp::TeX("$n = 30$, $p = 0.7$"))

```
::::
:::

## Bernoulli and Binomial

The Bernoulli and the Binomial distributions are used as **random components** when we have the dependent variable assuming 2 values (e.g., *correct* and *incorrect*) and we have the total number of trials:

- Accuracy on a cognitive task
- Patients recovered or not after a treatment
- People passing or not an exam

## Poisson distribution

The number of events $k$ during a fixed time interval (e.g., number of new user on a website in 1 week) is:

\begin{align*}
f(k,\lambda) = Pr(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}
\end{align*}

Where $k$ is the number of occurrences ($k = 0, 1, 2, ...$), $e$ is Euler's number ($e = 2.71828...$) and $!$ is the factorial function. The mean and the variance of the Poisson distribution is $\lambda$

## Poisson distribution

```{r, echo = FALSE}
lambda <- c(5, 10, 15)
dat <- expand.grid(x = 0:100, lambda = lambda)
dat$y <- dpois(dat$x, dat$lambda)

dat |> 
    filter(x < 30) |> 
    ggplot(aes(x = x, y = y,
               color = factor(lambda),
               group = factor(lambda))) +
    geom_point(size = 3) +
    geom_line() +
    mytheme() +
    ylab(latex2exp::TeX("dpois(x, $\\lambda$)")) +
    scale_color_discrete(name = latex2exp::TeX("$\\lambda$")) +
    theme(legend.position = "bottom")
```

As $\lambda$ increases, the distribution is well approximated by a Gaussian distribution, but the Poisson is discrete.

# Data simulation `#extra`

## Data simulation `#extra`

- During the course we will try to simulate some data. Simulating data is an amazing education tool to understand a statistical model.
- By simulating from a **generative model** we are doing a so-called **Monte Carlo Simulations** [@Gentle2009-cj]

## Data simulation `#extra`

In R there are multiple functions to generate data from probability distributions:

```{r, echo = FALSE}
distr <- c(dnorm, dbinom)
args <- lapply(distr, formals)

distr <- c("norm", "pois", "binom")
what <- c(r = "Generate random numbers", 
          d = "Compute the density", 
          q = "Return the quantile given a cumulative proability",
          p = "Return the cumulative probability given a quantile")

distr_tab <- tidyr::expand_grid(distr, what)

distr_tab$fun <- sprintf("%s", code(names(distr_tab$what)))
names(distr_tab) <- c("Distribution", "Action", "Function")

distr_tab |> 
    select(Function, Distribution, Action) |> 
    arrange(Function) |> 
    mytab() |> 
    collapse_rows(columns = c(1,3))
```

## Data simulation `#extra`

::: columns

:::: column
```{r, size = "tiny"}
n <- 1e5 # number of experiments
nt <- 100 # number of subjects
p <- 0.7 # probability of success
nc <- rbinom(n, nt, p)
```

```{r, echo = FALSE}
hist(nc/nt)
```
::::

:::: column
```{r, size = "tiny"}
n <- 1e5 # number of subjects
lambda <- 30 # mean/variance
y <- rpois(n, lambda)
```

```{r, echo = FALSE}
hist(y)
```
::::
:::

## References

\footnotesize
