Week 10 Data Dive

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(patchwork)
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)
theme_set(theme_minimal())

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

Select an interesting binary column of data, or one which can be reasonably converted into a binary variable 

For this chocolate CSV I have the columns of: ref, company_manufacturer, company_location, review_date, country_of_bean_origin, specific_bean_origin_or_bar_name, cocoa_percent, ingredients, most_memorable_characteristics, and rating. It seems that the only column of data that can reasonably be converted into a binary variable is rating, I could make ratings 0-2 equate to 0 and 3-5 equate to 1.

chocolate$rating<-ifelse(chocolate$rating >= 2, 1, 0)

Build a logistic regression model for this variable, using between 1-4 explanatory variables

For my data set there really is only one numerical explanatory variable I could use which is cocoa_percent which I will be using.

library(tidyverse)
library(broom)


data <- read.csv("~/Desktop/chocolate.csv")

threshold <- 2.0
data$binary_rating <- ifelse(data$rating > threshold, 1, 0)

data <- na.omit(data)

logistics_model <- glm(binary_rating ~ cocoa_percent,
                   data = data, family = binomial)

summary(logistics_model)
## 
## Call:
## glm(formula = binary_rating ~ cocoa_percent, family = binomial, 
##     data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     13.532      1.178  11.484   <2e-16 ***
## cocoa_percent  -12.941      1.502  -8.617   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.40  on 2529  degrees of freedom
## Residual deviance: 433.42  on 2528  degrees of freedom
## AIC: 437.42
## 
## Number of Fisher Scoring iterations: 7
conf_intervals <- confint(logistics_model)
## Waiting for profiling to be done...
print(conf_intervals)
##                   2.5 %   97.5 %
## (Intercept)    11.21207 15.85562
## cocoa_percent -15.87331 -9.95024
library(ggplot2)

data$predicted_prob <- predict(logistics_model, type = "response")

ggplot(data, aes(x = cocoa_percent, y = predicted_prob)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  labs(x = "Cocoa Percent", y = "Predicted Probability of Binary Outcome") +
  ggtitle("Logistic Regression: Predicted Probability vs. Cocoa Percent")
## `geom_smooth()` using formula = 'y ~ x'
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

The coefficients intercept (13.532) represents the log odds of the binary variable (rating) when the value of the predictor variable (cocoa_percent) is zero. The cocoa_percent coefficient (-12.941) represents the change in the log odds of the binary outcome associated with a one-unit increase in the cocoa_percent. It indicates that as cocoa_percent increases by one unit, the log odds of the binary outcome decreases by approximately 12.941 units. This makes sense as larger percents of cocoa in chocolate bars could have a more bitter taste which could result in lower ratings. Both of the coefficients have small p-values that means both coefficients are significant predictors of the binary outcome.

The confidence intervals for both coefficients, the 95% confidence intervals do not include zero, indicating that they are significantly different from zero. Overall, the logistic regression model suggests a strong negative relationship between cocoa_percent and the probability of the binary outcome rating. As cocoa_percent increases, the log odds of the binary outcome decrease significantly.

Consider a transformation for any explanatory variable, and illustrate why you need the transformation (or why you do not)

I decided to use the explanatory variable of cocoa_percent because there are no other numeric explanatory variables other than review_date which has no real influence on rating. Here I plotted the binary_rating with cocoa_percent.

# Scatterplot of cocoa_percent vs. binary_rating
ggplot(data, aes(x = cocoa_percent, y = binary_rating)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue")+
  labs(x = "Cocoa Percent", y = "Binary Rating") +
  ggtitle("Scatterplot of Cocoa Percent vs. Binary Rating")
## `geom_smooth()` using formula = 'y ~ x'

From this plot we see a negative linear relationship that hits some data points so I take it as the line will not be perfect as we are plotting cocoa_percent against a binary variable that will only produce us two y-values (0 or 1) so we will always get two parallel lines for this. I assume from this I do not need a transformation, but just to give it a shot I did the following: Log Transformation, Square Root Transformation, and Inverse Transformation.

Log Transformation:

data$log_cocoa_percent <- log(data$cocoa_percent + 1)  # Added 1 to avoid taking log of zero


ggplot(data, aes(x = log_cocoa_percent, y = binary_rating)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue")+
  labs(x = "Log Cocoa Percent", y = "Binary Rating") +
  ggtitle("Scatterplot of Log Cocoa Percent vs. Binary Rating")
## `geom_smooth()` using formula = 'y ~ x'

Again I see a negative linear relationship which makes sense as cocoa_percent increases rating decreases as the chocolate becomes bitter at higher cocoa percentages.

Square Root Transformation:

data$sqrt_cocoa_percent <- sqrt(data$cocoa_percent)


ggplot(data, aes(x = sqrt_cocoa_percent, y = binary_rating)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue")+
  labs(x = "Square Root Cocoa Percent", y = "Binary Rating") +
  ggtitle("Scatterplot of Square Root Cocoa Percent vs. Binary Rating")
## `geom_smooth()` using formula = 'y ~ x'

Again I see a negative linear relationship which makes sense as cocoa_percent increases rating decreases as the chocolate becomes bitter at higher cocoa percentages.

Inverse Transformation:

data$inv_cocoa_percent <- 1 / (data$cocoa_percent + 1)  # Adding 1 to avoid division by zero


ggplot(data, aes(x = inv_cocoa_percent, y = binary_rating)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue")+
  labs(x = "Inverse Cocoa Percent", y = "Binary Rating") +
  ggtitle("Scatterplot of Inverse Cocoa Percent vs. Binary Rating")
## `geom_smooth()` using formula = 'y ~ x'

I see a positive linear relationship which makes sense as this is the INVERSE plot but it still suggests as cocoa_percent increases rating decreases as the chocolate becomes bitter at higher cocoa percentages.

I decided to plot cocoa_percent against rating to see if it would give us a decreasing linear relationship which it did.

ggplot(data, aes(x = cocoa_percent, y = rating)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue")+
  labs(x = "Cocoa Percent", y = "Rating") +
  ggtitle("Scatterplot of Cocoa Percent vs. Rating")
## `geom_smooth()` using formula = 'y ~ x'