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 errors
library(ggthemes)
library(ggrepel)
library(broom)
library(lindia)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
options(scipen = 6)
library(tidyverse)
library(dplyr)
chocolate <- 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.000
Build a linear(or generalized linear) model:
model <- lm(rating ~ cocoa_percent, data = chocolate)
summary(model)
##
## 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-13
Use tools from previous weeks to diagnose the model
Residual Analysis: I decided to plot the residuals versus fitted values. The x-axis represents the fitted (predicted) values from the linear regression model, the y-axis is the residuals (the difference between the observed and predicted values). This plot allows us to assess the linearity (through the distribution of residuals along the x-axis), which we observe the values to hover over our horizontal 0 line, creating a distributed roughly symmetrically about the origin, which signifies that the linear regression model is a good fit for the data.
residuals <- residuals(model)
fitted_values <- fitted(model)
plot(fitted_values, residuals,
xlab = "Fitted Values", ylab = "Residuals",
main = "Residuals vs. Fitted Values",
pch = 16, col = "blue")
Cook’s Distance: I used this to identify influential data points that might have a large effect on the model. Cook’s distance measures how much the model’s coefficients change when an observation is removed from the data set. From this plot we can see that there is just one point that is largely over the value of 0.08 meaning if this point was removed from the data set it would result in a substantial change in the regression coefficients.
cooksd <- cooks.distance(model)
plot(cooksd, pch = "o", main = "Cook's distance plot")
abline(h = 4/length(cooksd), col = "red")
Heteroscedasticity: I decided to use a test for
homoscedasticity of residuals.
The studentized Breusch-Pagan test is a statistical test used to assess
the presence of heteroscedasticity in a linear regression model.
Heteroscedasticity is the situation where the variance of the
errors(residuals) in a regression model is not constant across all
levels of the predictor variables. The null hypothesis of this test is
that the variance of the residuals is constant across all levels of the
predictor variable(s), indicating homoscedasticity. The alternative
hypothesis is that the variance of the residuals is not constant,
indicating heteroscedasticity.
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
bptest(model)
##
## studentized Breusch-Pagan test
##
## data: model
## BP = 8.52, df = 1, p-value = 0.003513
Based on the p-value, there is evidence to suggest that the variance of the residuals is not constant across all levels of the predictor variable(s), indicating the presence of heteroscedasticity in the linear regression model. When heteroscedasticity is present, it can lead to biased and inefficient estimates of the regression coefficients.
Interpret at least one of the coefficients
summary(model)
##
## 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-13
I will be looking at the cocoa_percent coefficient which is -1.1630. This means that on average for every one unit increase in the cocoa percentage of a chocolate product, the predicted rating decreases by 1.1630 units which makes sense as a more cocoa concentrated chocolate results in a more bitter taste, decreasing its rating. The coefficient is statistically significant by its very small p-value (1.22e-13) which suggests the cocoa percentage has a significant impact on the rating of chocolate products.
The intercept term (4.0295) is the estimated rating of a chocolate product when the cocoa percentage is zero. This doesn’t provide much insight as it doesn’t seem practical for a chocolate product to have zero cocoa percentage.
The residual standard error (0.4406) represents the deviation of the observed ratings from the predicted ratings by the model, it indicates how well the model fits the data.
The multiple R-squared value (0.02152) indicates the proportion of variance in the ratings, so only 2.15% of the variance in ratings can be explained by the cocoa percentage alone.
The adjusted R-squared value (0.02113) adjusts for the number of predictors in the model.
Overall, the coefficient for cocoa_percent suggests that cocoa percentage has a statistically significant negative effect on the rating of chocolate products. However, the low R-squared value indicates that the model explains only a small proportion of the variance in ratings, suggesting that other factors not included in the model may also influence chocolate ratings.