Week 9
2024-03-19
Week 9 Data Dive
Simple linear regression model including 1-3 more variables
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorschocolate <- read.csv("~/Desktop/chocolate.csv")
library(readr)
summary(chocolate)## ref company_manufacturer company_location review_date
## Min. : 5 Length:2530 Length:2530 Min. :2006
## 1st Qu.: 802 Class :character Class :character 1st Qu.:2012
## Median :1454 Mode :character Mode :character Median :2015
## Mean :1430 Mean :2014
## 3rd Qu.:2079 3rd Qu.:2018
## Max. :2712 Max. :2021
## country_of_bean_origin specific_bean_origin_or_bar_name cocoa_percent
## Length:2530 Length:2530 Min. :0.4200
## Class :character Class :character 1st Qu.:0.7000
## Mode :character Mode :character Median :0.7000
## Mean :0.7164
## 3rd Qu.:0.7400
## Max. :1.0000
## ingredients most_memorable_characteristics rating
## Length:2530 Length:2530 Min. :1.000
## Class :character Class :character 1st Qu.:3.000
## Mode :character Mode :character Median :3.250
## Mean :3.196
## 3rd Qu.:3.500
## Max. :4.000library(dplyr)
library(tidyverse)
library(ggthemes)
library(ggrepel)
library(AmesHousing)
library(boot)
library(broom)
library(lindia)
# remove scientific notation
options(scipen = 6)
# default theme, unless otherwise noted
theme_set(theme_minimal())For last weeks data dive, I included the variables of cocoa_percent, rating, and country_of_bean_origin. This week I decided to add company_manufacturer as the rating of chocolates could be affected by the quality of the company manufacturer.
model1 <- lm(rating ~ cocoa_percent, chocolate) |>
summary()
print(model1)##
## Call:
## lm(formula = rating ~ cocoa_percent, data = chocolate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.21541 -0.23867 0.03459 0.28459 0.99393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0295 0.1121 35.949 < 2e-16 ***
## cocoa_percent -1.1630 0.1560 -7.456 1.22e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4406 on 2528 degrees of freedom
## Multiple R-squared: 0.02152, Adjusted R-squared: 0.02113
## F-statistic: 55.59 on 1 and 2528 DF, p-value: 1.218e-13model2 <- lm(rating ~ cocoa_percent+ country_of_bean_origin, chocolate) |>
summary()model3 <- lm(rating ~ cocoa_percent+ country_of_bean_origin+company_manufacturer, chocolate) |>
summary()model1 <- lm(rating ~ cocoa_percent, chocolate) |>
summary()
model2 <- lm(rating ~ cocoa_percent+ country_of_bean_origin, chocolate) |>
summary()
model3 <- lm(rating ~ cocoa_percent+ country_of_bean_origin+company_manufacturer, chocolate) |>
summary()
df_plot <- data.frame(p = c(2, 3, 4),
r_squared = c(model1$r.squared,
model2$r.squared,
model3$r.squared),
adj_r_squared = c(model1$adj.r.squared,
model2$adj.r.squared,
model3$adj.r.squared)) |>
mutate(r_perc_increase = 100 * (r_squared - model1$r.squared) / model1$r.squared,
r_adj_perc_increase = 100 * (adj_r_squared - model1$adj.r.squared) / model1$adj.r.squared,
diff = r_perc_increase - r_adj_perc_increase)
ggplot(data=df_plot) +
geom_line(mapping = aes(x=p, y = r_perc_increase,
color = 'R Squared Improvement')) +
geom_line(mapping = aes(x=p, y = r_adj_perc_increase,
color = 'Adjusted R Squared Improvement')) +
geom_text_repel(mapping = aes(x=p,
y = r_adj_perc_increase,
label = paste("diff = ", round(diff, 1))),
color = "darkred",
nudge_x = 0.3) +
labs(color='',
x='Number of Coefficients (p)',
y='% Increase in R-Squared',
title = "Improving on a Model with 2 Coefficients")+
scale_color_brewer(palette = 'Dark2') +
theme(legend.position="bottom")
From the graph above there is a huge concern that the % increase in R-squared is from 0-2000%.
model <- lm(rating ~ cocoa_percent + country_of_bean_origin + company_manufacturer, data = chocolate)
tidy(model, conf.int = TRUE)## # A tibble: 642 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 4.63 0.431 10.7 3.63e-26 3.78 5.47
## 2 cocoa_percent -1.40 0.167 -8.39 9.50e-17 -1.73 -1.08
## 3 country_of_bean_ori… -0.128 0.391 -0.327 7.43e- 1 -0.895 0.639
## 4 country_of_bean_ori… -0.357 0.390 -0.916 3.60e- 1 -1.12 0.407
## 5 country_of_bean_ori… -0.317 0.391 -0.812 4.17e- 1 -1.08 0.449
## 6 country_of_bean_ori… -0.244 0.392 -0.623 5.33e- 1 -1.01 0.525
## 7 country_of_bean_ori… -0.721 0.546 -1.32 1.87e- 1 -1.79 0.350
## 8 country_of_bean_ori… -0.610 0.454 -1.34 1.79e- 1 -1.50 0.280
## 9 country_of_bean_ori… -0.230 0.549 -0.419 6.75e- 1 -1.31 0.847
## 10 country_of_bean_ori… -0.224 0.393 -0.572 5.68e- 1 -0.994 0.545
## # ℹ 632 more rowsDiagnostic Plots
For each plot, I will point out any indications of issues with the model.
- Residuals vs. \(\hat{y}\) values
- Residuals vs. \(x\) values
- Residual Histogram
- QQ-Plot
- Cook’s D by Observation
Residuals vs. Fitted Values
The first plot we will look at is the residual values vs. the fitted values:
gg_resfitted(model) +
geom_smooth(se=FALSE)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
The fitted values display the chocolate rating, and the plot displays that the variance in rating is huge for ratings 2-4.
Residuals vs. X Values
plots <- gg_resX(model, plot.all = FALSE)
# for each variable of interest ...
plots$cocoa_percent+ theme(axis.text.x = element_text(angle = 90, hjust = 1))
plots$country_of_bean_origin+theme(axis.text.x = element_text(angle = 90, hjust = 1))
plots$company_manufacturer+theme(axis.text.x = element_text(angle = 90, hjust = 1))
What we want is for a constant variance of residuals across values of \(X\) variables. The plot for cocoa_percent follows the a pattern in the cocoa_percent from 0.6-0.8.
Residual Histogram
gg_reshist(model)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There is a concern as the tails are not equally distributed, there is a “shorter” tail on the left, and a longer one on the right. This is not severe, but it implies that some aspect of the phenomenon is missing from the model such as we may need more/different explanatory variables.
QQ-Plots
# # the normal QQ plot
gg_qqplot(model)
Cook’s D by Observation
Here, we can see a few rows of data have a high influence on the model. Each number represents the row in the data. Rows 9, 134, 177, and 254 may be troublesome, so we can investigate them by plotting various variables against the response variable.
For example:
# Fit your linear regression model
model <- lm(rating ~ cocoa_percent + country_of_bean_origin + company_manufacturer, data = chocolate)
# Compute Cook's distances
cooksd <- cooks.distance(model)
# Create a data frame with observation index and Cook's distances
cooksd_data <- data.frame(index = 1:length(cooksd), cooksd = cooksd)
# Plot Cook's distances by observation
library(ggplot2)
ggplot(cooksd_data, aes(x = index, y = cooksd)) +
geom_point() +
labs(title = "Cook's Distance by Observation",
x = "Observation Index",
y = "Cook's Distance") +
theme_minimal()## Warning: Removed 117 rows containing missing values (`geom_point()`).
Clearly, some of these rows have a stronger influence when we
consider just the lot_area, but not all of them. We may
want to find where “highly influential rows” like these intersect across
the different explanatory variables of the model. I.e., which rows are
very highly influential among most (or all) explanatory variables.