Instantly share code, notes, and snippets.

graebnerc/#T12: lecture notes

Created June 7, 2022 15:33
Show Gist options
• 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 Lecture notes and solutions to the exercises of session 13 on multiple linear regression
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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"]]