---
title: "Chapter 8: The Linear Model"
author: "Andy Field"
date: "2025-04-03"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(rstanarm)
library(tidyverse)

# load album sales data
album_tib  <- here::here("data/album_sales.csv") |> 
  read_csv()

# or load from the discovr package
# album_tib <- discovr::album_sales

```

## Scatterplot

```{r}
#| label: scattergg

ggplot(album_tib, aes(adverts, sales)) +
  geom_point(colour = "#88CCEE") +
  coord_cartesian(ylim = c(0, 400), xlim = c(0, 2500)) +
  scale_x_continuous(breaks = seq(0, 2500, 500)) +
  scale_y_continuous(breaks = seq(0, 400, 100)) +
  labs(x = "Advertising budget (thousands)", y = "Album sales (thousands)") +
  geom_smooth(method = "lm", colour = "#CC6677", fill = "#CC6677", alpha = 0.2) +
  theme_minimal()

```

## Summary statistics

```{r}
#| label: album_summary

describe_distribution(album_tib, range = FALSE, ci = 0.95, iterations = 500) |> 
  display()

```

## Linear model with one predictor

```{r}
#| label: lm_adverts

album_lm <- lm(sales ~ adverts, data = album_tib)
# text output
summary(album_lm)

# get model fit
model_performance(album_lm) |>
  display()

# test model fit
test_wald(album_lm) |> 
  display()

# get parameter estimates
model_parameters(album_lm) |>
  display(digits = 3)
```

## Linear model with several predictors


```{r}
#| label: lm_full

# plot the variables

GGally::ggscatmat(album_tib, columns = c("adverts", "airplay", "image", "sales")) +
  theme_minimal()


#fit the model
album_full_lm <- lm(sales ~ adverts + airplay + image, data = album_tib, na.action = na.exclude)

# compare models
compare_performance(album_lm, album_full_lm) |> 
      display(digits = 3)

# test models
test_wald(album_lm, album_full_lm) |> 
      display(digits = 3)

# get parameter estimates
model_parameters(album_full_lm) |>
  display(digits = 3)

# get standardized parameter estimates
model_parameters(album_full_lm, standardize = "refit") |>
  display(digits = 3)

```


## Model assumptions

### Multicollinearity

```{r}
#| label: vif

check_collinearity(album_full_lm)
# mean VIF
check_collinearity(album_full_lm)$VIF |> 
  mean()
```

### Assumption checks

```{r}
#| label: check_model

# plot residuals
check_model(album_full_lm)

```



### Casewise diagnostics

```{r}
#| label: resids

# add diagnostic stats to the original data and store
album_full_rsd <-  broom::augment(album_full_lm) |> 
  rowid_to_column(var = "case_no")

# influence statistics
album_full_inf <- influence.measures(album_full_lm)$infmat |> 
  as_tibble() |> 
  rowid_to_column(var = "case_no")

# write our own function to count the proportion of z-scores
 get_cum_percent <- function(var,  cut_off = 1.96){
  ecdf_var <- abs(var) |> ecdf()
  100*(1 - ecdf_var(cut_off))
 }
 
album_full_rsd |> 
  summarize(
    z_1.96 = get_cum_percent(.std.resid),
    z_2.58 = get_cum_percent(.std.resid,cut_off = 2.58),
    z_3.29 = get_cum_percent(.std.resid,cut_off = 3.29)
  )


# filter residuals
check_outliers(album_full_rsd$.resid,
               method = "zscore",
               threshold = list("zscore" = 1.96))

# check Cook's distance
check_outliers(album_full_lm,
               method = "cook",
               threshold = list('cook' = 0.5))

# filter to see any cases with any value of dfbeta > 1
album_full_inf  |> 
  filter(
    if_any(starts_with("dfb"), ~ abs(.) > 1)
  )

# filter leverage values

leverage_thresh <- 3*mean(album_full_inf$hat, na.rm = TRUE) # store threshold   
album_full_inf |>
  filter(
    hat > leverage_thresh | !between(cov.r, 1 - leverage_thresh, 1 + leverage_thresh)
    ) |> 
  select(!starts_with("dfb")) |> 
  display(digits = 3)


```

## Robust models

### Robust estimation of parameters

```{r}
#| label: lmrob

album_full_rob <- robust::lmRob(sales ~ adverts + airplay + image, data = album_tib, na.action = na.exclude)

#text output
summary(album_full_rob)

# table of parameters that renders nicely
model_parameters(album_full_rob) |> 
  display(digits = 3)


# table of bias tests that renders nicely
summary(album_full_rob)$biasTest |> 
  knitr::kable(digits = 3)

```

### Robust standard errors

```{r}
#| label: hcse

# table of parameters that renders nicely
model_parameters(album_full_lm, vcov = "HC4") |> 
  display(digits = 3)
```


### Bootstrap standard errors

```{r}
#| label: boot

model_parameters(album_full_lm, bootstrap = TRUE) |> 
  display(digits = 3)

```


## Bayesian models

```{r}
#| label: bayes

# Bayesian model with default priors and bayes factors
album_dp <- stan_glm(sales ~ adverts + airplay + image, data = album_tib)

model_parameters(model = album_dp,
                 test = c("pd", "bf"),
                 null = 0) |> 
  display()


# plot potential priors

# prior for intercept
ggplot(tibble(x = c(-10, 50)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 20, sd = 7), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-10, 50, 5)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 20, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for intercept') +
  theme_minimal()

# prior for advertising
ggplot(tibble(x = c(-1, 5)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 1, sd = 0.5), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-1, 5, 1)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 1, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for advertising') +
  theme_minimal()

# prior for airplay
ggplot(tibble(x = c(-2, 12)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 5, sd = 2), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-2, 12, 1)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 5, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for airplay') +
  theme_minimal()

# prior for image
ggplot(tibble(x = c(-10, 25)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 8, sd = 5), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-10, 25, 5)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 8, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for image') +
  theme_minimal()

# Bayesian model with informative priors and bayes factors

album_rstn <- stan_glm(sales ~ adverts + airplay + image,
                      data = album_tib,
                      prior = normal(location = c(1, 5, 8), scale = c(0.5, 2, 5)),
                      prior_intercept = normal(location = 20, scale = 7))

model_parameters(album_rstn, test = c("pd", "bf"), null = 0) |> 
  display()
```


