# Setup
if (!require("dplyr")) install.packages("dplyr")
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
if (!require("gt")) install.packages("gt")
## Loading required package: gt
if (!require("gtExtras")) install.packages("gtExtras")
## Loading required package: gtExtras
library(dplyr)
library(ggplot2)
library(gt)
library(gtExtras)

options(scipen = 999)
# Load Data
mydata <- read.csv("SandwichAd.csv") #Edit: TYPE YOUR DATA FILE'S NAME

mydata <- mydata %>%
  mutate(
    DV = Calories, #Edit: TYPE YOUR DEPENDENT VARIABLE'S NAME
    IV = Actor  #Edit: TYPE YOUR INDEPENDENT VARIABLE'S NAME
  )
# Histograms of DV per IV group with group means
# Calculate group means
group_means <- mydata %>%
  group_by(IV) %>%
  summarise(
    mean_DV = mean(DV, na.rm = TRUE),
    .groups = "drop"
  )

Graphic <- ggplot(mydata, aes(x = DV)) +
  geom_histogram(
    binwidth = diff(range(mydata$DV, na.rm = TRUE)) / 30,
    color = "black",
    fill = "#1f78b4",
    alpha = 0.7
  ) +
  geom_vline(
    data = group_means,
    aes(xintercept = mean_DV),
    color = "red",
    linetype = "dashed",
    linewidth = 1
  ) +
  facet_grid(IV ~ .) +
  labs(
    title = "DV Distributions by IV Group",
    x = "Dependent Variable",
    y = "Count"
  ) +
  theme_minimal()

# Show graphic
Graphic

# Descriptive Statistics
Descriptives <- mydata %>%
  group_by(IV) %>%
  summarise(
    count = n(),
    mean = mean(DV, na.rm = TRUE),
    sd = sd(DV, na.rm = TRUE),
    min = min(DV, na.rm = TRUE),
    max = max(DV, na.rm = TRUE),
    .groups = "drop"
  )

DescriptivesTable <- Descriptives %>%
  gt() %>%
  fmt_number(
    columns = c(mean, sd, min, max),
    decimals = 2
  ) %>%
  tab_header(
    title = "Descriptive Statistics",
    subtitle = "Summary statistics by group"
  ) %>%
  cols_label(
    IV = "Group",
    count = "N",
    mean = "Mean",
    sd = "Standard Deviation",
    min = "Minimum",
    max = "Maximum"
  )
# Normality Check (Shapiro-Wilk)
shapiro_summary <- mydata %>%
  group_by(IV) %>%
  summarise(
    W_statistic = shapiro.test(DV)$statistic,
    p_value = shapiro.test(DV)$p.value,
    .groups = "drop"
  )

ShapiroTable <- shapiro_summary %>%
  gt() %>%
  fmt_number(
    columns = c(W_statistic, p_value),
    decimals = 3
  ) %>%
  tab_header(
    title = "Shapiro-Wilk Normality Test Results",
    subtitle = "Tests of normality within each group"
  ) %>%
  cols_label(
    IV = "Group",
    W_statistic = "W Statistic",
    p_value = "p-value"
  ) %>%
  tab_source_note(
    source_note = md(
      "NOTE: If one or both of the p-values are less than .05 and if the groups are smaller than 40, the t-test's normality assumption may be violated, and you should use the Wilcoxon rank-sum test results instead of the t-test results."
    )
  )
# Inferential Tests
# Run Welch's t-test (default for two groups)
t_res <- t.test(
  DV ~ IV,
  data = mydata,
  var.equal = FALSE
)

# Run Wilcoxon rank-sum test
wilcox_res <- wilcox.test(
  DV ~ IV,
  data = mydata
)

# Create a tidy summary of t-test results
t_summary <- tibble(
  Group1 = levels(as.factor(mydata$IV))[1],
  Group2 = levels(as.factor(mydata$IV))[2],
  Mean1 = t_res$estimate[1],
  Mean2 = t_res$estimate[2],
  t = t_res$statistic,
  df = t_res$parameter,
  p = t_res$p.value,
  CI_low = t_res$conf.int[1],
  CI_high = t_res$conf.int[2]
)

# Create a tidy summary of Wilcoxon results
wilcox_summary <- tibble(
  Group1 = levels(as.factor(mydata$IV))[1],
  Group2 = levels(as.factor(mydata$IV))[2],
  W = wilcox_res$statistic,
  p = wilcox_res$p.value
)
# Present t-Test Results as gt Table
Table <- t_summary %>%
  gt() %>%
  fmt_number(
    columns = c(Mean1, Mean2, CI_low, CI_high),
    decimals = 2
  ) %>%
  fmt_number(
    columns = c(t, df, p),
    decimals = 3
  ) %>%
  tab_header(
    title = "Independent Samples t-Test Results",
    subtitle = "Welch's t-test (unequal variances assumed)"
  ) %>%
  cols_label(
    Group1 = "Group 1",
    Group2 = "Group 2",
    Mean1 = "Mean (Group 1)",
    Mean2 = "Mean (Group 2)",
    t = "t Statistic",
    df = "Degrees of Freedom",
    p = "p-value",
    CI_low = "95% CI (Lower)",
    CI_high = "95% CI (Upper)"
  )
# Present Wilcoxon Results as gt Table
WilcoxonTable <- wilcox_summary %>%
  gt() %>%
  fmt_number(
    columns = c(W, p),
    decimals = 3
  ) %>%
  tab_header(
    title = "Wilcoxon Rank-Sum Test Results",
    subtitle = "Nonparametric comparison of group distributions"
  ) %>%
  cols_label(
    Group1 = "Group 1",
    Group2 = "Group 2",
    W = "W Statistic",
    p = "p-value"
  )
# Visual Output
# Histogram
Graphic

# Descriptive statistics table
DescriptivesTable
Descriptive Statistics
Summary statistics by group
Group N Mean Standard Deviation Minimum Maximum
Buff 50 584.26 128.03 316.00 827.00
Dad bod 50 656.28 141.25 358.00 946.00
 
# Shapiro-Wilk table
ShapiroTable
Shapiro-Wilk Normality Test Results
Tests of normality within each group
Group W Statistic p-value
Buff 0.976 0.385
Dad bod 0.985 0.790
NOTE: If one or both of the p-values are less than .05 and if the groups are smaller than 40, the t-test’s normality assumption may be violated, and you should use the Wilcoxon rank-sum test results instead of the t-test results. 
# Wilcoxon rank-sum table
WilcoxonTable
Wilcoxon Rank-Sum Test Results
Nonparametric comparison of group distributions
Group 1 Group 2 W Statistic p-value
Buff Dad bod 882.000 0.011
 
# Welch's t-test table
Table
Independent Samples t-Test Results
Welch's t-test (unequal variances assumed)
Group 1 Group 2 Mean (Group 1) Mean (Group 2) t Statistic Degrees of Freedom p-value 95% CI (Lower) 95% CI (Upper)
Buff Dad bod 584.26 656.28 −2.671 97.069 0.009 −125.53 −18.51