Lecture notes and solutions to the exercises of session 13 on multiple linear regression

#T12: lecture notes

Lecture notes and solutions to the exercises of session 13 on multiple linear regression

T13-CatVar-MultPlots.R

here::i_am("R/T13-CatVar-MultPlots.R")
library(here)
library(DataScienceExercises)
library(ggplot2)
library(dplyr)
library(icaeDesign)

journal_data <- DataScienceExercises::econjournals %>%
  dplyr::filter(papers>10, sub_price<5000)

# Parallel slopes model----------------
journal_linmod_pslopes <- journal_data %>%
  ggplot(data = ., aes(x=pages_py, y=sub_price, color=publisher_type)) +
  geom_point() +
  geom_parallel_slopes(se = FALSE) +
  labs(y="Subscription price", x = "Page length") +
  scale_color_brewer(palette = "Set1", direction = -1) +
  coord_cartesian(xlim = c(0, 2800)) +
  theme_icae() + guides(col=guide_legend(title = "Publisher type: ")) +
  theme(legend.position = "bottom", legend.title = element_text())
journal_linmod_pslopes

ggsave(plot = journal_linmod_pslopes,
       filename = here("output/T13-ParSlopes-Model.pdf"),
       width = 5.5, height = 3)

# Interaction model--------------------
journal_linmod_interact <- journal_data %>%
  ggplot(data = ., aes(x=pages_py, y=sub_price, color=publisher_type)) +
  geom_point() +
  stat_smooth(method="lm", fullrange=TRUE, se=FALSE) +
  labs(y="Subscription price", x = "Page length") +
  scale_color_brewer(palette = "Set1", direction = -1) +
  coord_cartesian(xlim = c(0, 2800)) +
  theme_icae() + guides(col=guide_legend(title = "Publisher type: ")) +
  theme(legend.position = "bottom", legend.title = element_text())
journal_linmod_interact

ggsave(plot = journal_linmod_interact,
       filename = here("output/T13-Interaction-Model.pdf"),
       width = 5.5, height = 3)

journal_linmod_interact_zoom <- journal_linmod_interact +
  coord_cartesian(xlim = c(-10, 500), ylim = c(-10, 500)) +
  stat_smooth(method="lm", fullrange=TRUE, se=FALSE) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  theme(legend.position = "none")
journal_linmod_interact_zoom

ggsave(plot = journal_linmod_interact_zoom,
       filename = here("output/T13-Interaction-Model-Zoomed.pdf"),
       width = 2.5, height = 4)

T13-CatVar-ViolinPlot.R

here::i_am("R/T13-CatVar-ViolinPlot.R")
library(here)
library(ggplot2)
library(DataScienceExercises)
library(icaeDesign)

life_exp <- DataScienceExercises::gdplifexp2007 %>%
  dplyr::select(continent, lifeExp, gdpPercap) %>%
  dplyr::mutate(continent=factor(continent))

life_exp_categories_plot <- ggplot(
  data = life_exp,
  aes(x=continent, color=continent, fill=continent, y=lifeExp)
) +
  geom_point(alpha=0.25, shape=21) +
  geom_point(
    data = life_exp_categories,
    mapping = aes(x=continent, fill=continent, y=lifeExp_mean),
    shape = 13, alpha=1.0, size=4) +
  labs(y="Life Expct.") +
  theme_icae() +
  theme(
    axis.title.x = element_blank(),
    axis.text = element_text(size=8),
    axis.title.y = element_text(size=10),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank()
  )

ggsave(plot = life_exp_categories_plot,
       filename = here("output/lifeExp_catReg-plot.pdf"),
       width = 5.5, height = 2.5)

T13-Exercise-1.R

library(DataScienceExercises)
library(texreg)
# This solution explores different specifications, not only the one required
#  by the exercise on slide 15
beer_data <- DataScienceExercises::beer

mult_model_inc <- lm(
  formula = consumption ~ income,
  data = beer_data)

mult_model_inc_price <- lm(
  formula = consumption ~ income + price,
  data = beer_data)

mult_model_inc_price_liq <- lm(
  formula = consumption ~ income + price + price_liquor,
  data = beer_data)

mult_model_full <- lm(
  formula = consumption ~ income + price + price_liquor + price_other,
  data = beer_data)

texreg::screenreg(
  l = list(mult_model_inc,
           mult_model_inc_price,
           mult_model_inc_price_liq,
           mult_model_full),
  digits = 5)

T13-Script.R

library(DataScienceExercises)
library(dplyr)
library(skimr)
library(moderndive)
library(texreg)
beer_data <- DataScienceExercises::beer

# Slide 14: different regression models----------

simple_model_inc <- lm(
  formula = consumption ~ income,
  data = beer_data)

simple_model_price <- lm(
  formula = consumption ~ income,
  data = beer_data)

mult_model_inc_price <- lm(
  formula = consumption ~ income + price,
  data = beer_data)

texreg::screenreg(
  l = list(
    simple_model_inc,
    simple_model_price,
    mult_model_inc_price),
  digits = 5)

# Exploratory analysis on slide 18
life_exp <- DataScienceExercises::gdplifexp2007 %>%
  dplyr::select(continent, lifeExp, gdpPercap) %>%
  dplyr::mutate(continent=factor(continent))
skimr::skim(life_exp)

# Illustrating regressions with categorical LHS-------------

cont_linmod <- lm(lifeExp~continent, data = life_exp)
get_regression_table(cont_linmod)

africa_mean_exp <- life_exp %>%
  dplyr::filter(continent=="Africa") %>%
  dplyr::summarise(lifeExp_mean=mean(lifeExp)) %>%
  dplyr::pull(lifeExp_mean)

life_exp_categories <- life_exp %>%
  dplyr::group_by(continent) %>%
  dplyr::summarise(
    lifeExp_mean=mean(lifeExp),
    diff_africa=lifeExp_mean-africa_mean_exp
  )
cont_linmod

# For the violin plot see T13-CatVar-ViolinPlot.R

# Estimations for multiple regression part-----------------
# For the plots see T13-CatVar-MultPlots.R

journal_data <- DataScienceExercises::econjournals %>%
  dplyr::filter(papers>10, sub_price<5000)

# The parallel slopes model:
journal_linmod_psm <- lm(
  formula = sub_price~pages_py+publisher_type,
  data = journal_data)
get_regression_table(journal_linmod_psm)

# The interaction model:

journal_linmod_intct <- lm(sub_price~pages_py*publisher_type, data = journal_data)
get_regression_table(journal_linmod_intct)

# Model selection using R-squared----------------

summary(journal_linmod_intct)[["r.squared"]]
summary(journal_linmod_intct)[["adj.r.squared"]]

summary(journal_linmod_psm)[["r.squared"]]
summary(journal_linmod_psm)[["adj.r.squared"]]