Skip to content

Instantly share code, notes, and snippets.

@graebnerc
Created June 7, 2022 15:33
Show Gist options
  • Save graebnerc/4a9e5bb95459d8ada31f43141976efe8 to your computer and use it in GitHub Desktop.
Save graebnerc/4a9e5bb95459d8ada31f43141976efe8 to your computer and use it in GitHub Desktop.
Lecture notes and solutions to the exercises of session 13 on multiple linear regression
Lecture notes and solutions to the exercises of session 13 on multiple linear regression
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)
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)
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)
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"]]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment