Week 12 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(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(tsibble)
## 
## Attaching package: 'tsibble'
## 
## The following object is masked from 'package:zoo':
## 
##     index
## 
## The following object is masked from 'package:lubridate':
## 
##     interval
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
theme_set(theme_minimal())
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

Select a column of your data that encodes time (e.g., “date”, “time stamp”, “year”, etc.). Convert this into a Date in R.

My chocolate csv has the review_date in the format of year, so I decided to make up a month and day and just keep the year.

chocolate$review_date <- as.Date(paste(chocolate$review_date, "-01-01", sep = ""), format = "%Y-%m-%d")

Choose a column of data to analyze over time. This should be a “response-like” variable that is of particular interest.

I decided to choose chocolate rating as the “response-like” variable to see how chocolate ratings have changed over the years!

average_rating <- chocolate %>%
  group_by(review_date) %>%
  summarize(avg_rating = mean(rating, na.rm = TRUE))

ggplot(average_rating, aes(x = review_date, y = avg_rating)) +
  geom_line() +
  labs(x = "Review Date", y = "Average Rating", title = "Average Chocolate Rating Over Time")

Create a tsibble object of just the date and response variable. Then, plot your data over time. Consider different windows of time.

library(tsibble)
library(dplyr)
library(ggplot2)


average_ratings <- chocolate %>%
  group_by(year = lubridate::year(review_date)) %>%
  summarize(avg_rating = mean(rating, na.rm = TRUE))


chocolate_tsibble <- tsibble(
  review_date = as.Date(paste(average_ratings$year, "-01-01", sep = ""), format = "%Y-%m-%d"),
  avg_rating = average_ratings$avg_rating,
  index = review_date)


ggplot(chocolate_tsibble, aes(x = review_date, y = avg_rating)) +
  geom_line() +
  labs(x = "Review Date", y = "Average Rating", title = "Average Chocolate Rating Over Time")

  • What stands out immediately?

From creating this tsibble object of just the date and response variable, I ran into the issue that due to my data set only giving me the review_date in the format of year, I converted the dates to all have 01-01-year format. This then caused the issue that there are several reviews with the exact same date. Therefore, in order to resolve this I decided to take the average rating of the duplicate dates and plot that against the review_date year.

Use linear regression to detect any upwards or downwards trends.

lm_model <- lm(avg_rating ~ review_date, data = chocolate_tsibble)
summary(lm_model)
## 
## Call:
## lm(formula = avg_rating ~ review_date, data = chocolate_tsibble)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13341 -0.04169  0.01386  0.02550  0.12549 
## 
## Coefficients:
##               Estimate Std. Error t value       Pr(>|t|)    
## (Intercept) 2.49063611 0.16712627  14.903 0.000000000555 ***
## review_date 0.00004344 0.00001046   4.153       0.000976 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07045 on 14 degrees of freedom
## Multiple R-squared:  0.552,  Adjusted R-squared:   0.52 
## F-statistic: 17.25 on 1 and 14 DF,  p-value: 0.0009757
  • Do you need to subset the data for multiple trends?

  • How strong are these trends?

From this I was able to see that the coefficient for the review date is meaningful as the p-value is 0.0009757 which is really small indicating that the linear regression is something here. This suggests that there is evidence to conclude that there is a linear trend in the average rating over time. This makes much more sense as one would expect as technology becomes more advanced the chocolate becomes better!

Use smoothing to detect at least one season in your data, and interpret your results.

chocolate_tsibble$avg_rating_ma <- zoo::rollmean(chocolate_tsibble$avg_rating, k = 3, fill = NA)

ggplot(chocolate_tsibble, aes(x = review_date)) +
  geom_line(aes(y = avg_rating), color = "blue", alpha = 0.7) +
  geom_line(aes(y = avg_rating_ma), color = "red") +
  labs(x = "Review Date", y = "Average Rating", title = "Average Chocolate Rating Over Time (with Moving Average)")
## Warning: Removed 2 rows containing missing values (`geom_line()`).

I did this to try a different approach to detect seasonality in my chocolate rating over time data. I decided to calculate a rolling mean of the average rating variable in the chocolate_tsibble. The k=3 parameter specified the width or the number of observations to include in each calculation of the rolling mean, so in my case the k=3 indicates the rolling mean was calculated using the current observation and the two preceding observations. There is a difference in trend in the average rating about once every ~3 years, meaning the first 3 years there was a decrease in the average rating, then the next 3 years there was an increase, and the next 3 years there was a decrease.

  • Can you illustrate the seasonality using ACF or PACF?
library(tsibble)
library(ggplot2)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
acf_result <- Acf(chocolate_tsibble$avg_rating, lag.max = 24, plot = FALSE)
pacf_result <- Pacf(chocolate_tsibble$avg_rating, lag.max = 24, plot = FALSE)

autoplot(acf_result, main = "ACF Plot")
## Warning in ggplot2::geom_segment(lineend = "butt", ...): Ignoring unknown
## parameters: `main`

autoplot(pacf_result, main = "PACF Plot")
## Warning in ggplot2::geom_segment(lineend = "butt", ...): Ignoring unknown
## parameters: `main`

In this section I attempted to illustrate the seasonality using autocorrelation function (ACF) and partial autocorrelation function (PACF). These plots do not cross the dotted-blue lines, which indicate the bounds of the 95% confidence interval, therefore there is no evidence of any correlation with previous lags. There are no obvious, repeating patterns in either graph that indicate seasonality. This lack of seasonality is associated with the available data from the csv, the review date is listed in years which results in a smooth view of the avg_rating over time from the start. I think the lack of seasonality makes sense in this data set because we are given a general average rating of the chocolate over the years. If we had average rating on specific batches of chocolates and the dates of the cocoa beans used to produce that chocolate, then maybe we could find whether there’s a seasonality between chocolate batches (their average ratings) over time.