# Load tidyverse
library(tidyverse)
library(janitor)
library(kableExtra)
# Set seed for reproducibility
set.seed(42)

# Simulate data
ufo_data <- tibble(
  alien_species = sample(c("Martian", "Venusian", "Zeta Reticulan", "Reptilian"),
                         size = 100, replace = TRUE, prob = c(0.4, 0.2, 0.3, 0.1)),
  abduction_method = sample(c("Tractor Beam", "Teleportation", "Sleep Gas", "Invisibility Cloak"),
                            size = 100, replace = TRUE)
)

# View a few rows
# ufo_data %>% head()

# Create joint frequency table
ufo_table <- ufo_data |> 
  tabyl(alien_species, abduction_method)


ufo_table %>%
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> column_spec(1, bold = TRUE)

ufo_table |> adorn_totals(where = c("row","col"))|> 
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

ufo_table |> 
  adorn_totals(where = c("row","col"))|> 
  adorn_percentages(denominator = "all") |>
  adorn_rounding(digits = 2) |>
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

  ufo_table |> 
  adorn_totals(where = c("row","col"))|> 
  adorn_percentages(denominator = "row") |>
  adorn_rounding(digits = 2) |>
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

  ufo_table |> 
  adorn_totals(where = c("row","col"))|> 
  adorn_percentages(denominator = "col") |>
  adorn_rounding(digits = 2) |>
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

library(tidyverse)
library(janitor)

# Step 1: Simulate data
set.seed(42)
ufo_data <- tibble(
  alien_species = sample(c("Martian", "Venusian", "Zeta Reticulan", "Reptilian"),
                         size = 100, replace = TRUE, prob = c(0.4, 0.2, 0.3, 0.1)),
  abduction_method = sample(c("Tractor Beam", "Teleportation", "Sleep Gas", "Invisibility Cloak"),
                            size = 100, replace = TRUE)
)

# Step 2: Create observed tabyl with totals
observed_tab <- ufo_data %>%
  tabyl(alien_species, abduction_method) %>%
  adorn_totals("both")

# Step 3: Extract row and column totals
row_totals <- observed_tab %>%
  filter(alien_species != "Total") %>%
  pull(Total)

col_totals <- observed_tab %>%
  filter(alien_species == "Total") %>%
  select(-alien_species, -Total) %>%
  as.numeric()

grand_total <- sum(row_totals)

# Step 4: Compute expected frequencies under independence
expected_matrix <- outer(row_totals, col_totals) / grand_total
dimnames(expected_matrix) <- list(
  observed_tab$alien_species[observed_tab$alien_species != "Total"],
  names(observed_tab)[2:(ncol(observed_tab) - 1)]
)

# Step 5: Convert expected matrix to tabyl-style tibble
expected_tab <- as_tibble(expected_matrix, rownames = "alien_species") %>%
  adorn_totals("both") 


# Print expected frequencies under independence
expected_tab |> adorn_percentages("row") |>  
  kbl(caption = "Abduction Method by Alien Species", align = "c") |> 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

expected_tab |> adorn_percentages("col") |>  
  kbl(caption = "Abduction Method by Alien Species", align = "c") |> 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

expected_tab |> adorn_percentages("all") |>  
  kbl(caption = "Abduction Method by Alien Species", align = "c") |> 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

ufo_table |> 
  adorn_totals(where = c("row","col"))|> 
  adorn_percentages(denominator = "all") |>
  adorn_rounding(digits = 2) |>
  kbl(caption = "Abduction Method by Alien Species", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center") |> 
  column_spec(1, bold = TRUE) |> 
  column_spec(6, color = "dodgerblue") |> 
  row_spec(5, color = "forestgreen") 

chisq.test(table(ufo_data), correct = FALSE) |> 
  broom::tidy() |> 
  select(statistic, parameter, p.value) |> 
  mutate(across(everything(), ~ round(.x, 2))) |> kbl()

library(tidyverse)

set.seed(123)

# 1. Strong Positive: Ice cream vs shark attacks
df1 <- tibble(
  ice_cream_sales = seq(100, 1000, length.out = 50),
  shark_attacks = ice_cream_sales * 0.05 + rnorm(50, 0, 10),
  pair = "Ice Cream vs Shark Attacks"
)

df1 |> 
  ggplot(aes(x = ice_cream_sales, y = shark_attacks)) +
  geom_point(color="indianred") +
  # geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(x = "Ice Cream Sales", y = "Shark Attacks") + #title = "Ice Cream Sales vs Shark Attacks",
  theme_minimal()

df1 |> 
  mutate(
    # Sign of contribution to correlation (+ if both above or both below mean, - otherwise)
    sign_index = ifelse(
      (ice_cream_sales > mean(ice_cream_sales) & shark_attacks > mean(shark_attacks)) |
      (ice_cream_sales < mean(ice_cream_sales) & shark_attacks < mean(shark_attacks)),
      "+", "-"
    ),
    # Size index based on distance from means (scaled absolute contribution)
    size_index = abs((ice_cream_sales - mean(ice_cream_sales)) * 
                     (shark_attacks - mean(shark_attacks)))
  ) |>
  ggplot(aes(x = ice_cream_sales, y = shark_attacks, color = sign_index)) +
  geom_text(aes(label = sign_index, size = size_index), alpha = 0.8) +
  geom_vline(xintercept = mean(df1$ice_cream_sales), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df1$shark_attacks), color = "forestgreen",size=1.5,alpha=.5) +
  labs(
    # title = "Contribution Signs to Correlation",
    # subtitle = "Positive ( + ) when above/below both means, Negative ( - ) otherwise",
    x = "Ice Cream Sales",
    y = "Shark Attacks"
  ) +
  scale_size(range = c(3, 8)) +
  theme_minimal() +
  theme(legend.position = "none")

df1 |> 
  mutate(
    # Sign of contribution to correlation (+ if both above or both below mean, - otherwise)
    sign_index = ifelse(
      (ice_cream_sales > mean(ice_cream_sales) & shark_attacks > mean(shark_attacks)) |
      (ice_cream_sales < mean(ice_cream_sales) & shark_attacks < mean(shark_attacks)),
      "+", "-"
    ),
    # Size index based on distance from means (scaled absolute contribution)
    size_index = abs((ice_cream_sales - mean(ice_cream_sales)) * 
                     (shark_attacks - mean(shark_attacks)))
  ) |>
  ggplot(aes(x = ice_cream_sales, y = shark_attacks, color = sign_index)) +
  geom_text(aes(label = sign_index, size = size_index), alpha = 0.8) +
  geom_vline(xintercept = mean(df1$ice_cream_sales), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df1$shark_attacks), color = "forestgreen",size=1.5,alpha=.5) +
  labs(
    # title = "Contribution Signs to Correlation",
    # subtitle = "Positive ( + ) when above/below both means, Negative ( - ) otherwise",
    x = "Ice Cream Sales",
    y = "Shark Attacks"
  ) +
  scale_size(range = c(3, 8)) +
  theme_minimal() +
  theme(legend.position = "none")

# 2. Moderate Positive: Beachgoers vs sunburns
df2 <- tibble(
  beachgoers = seq(50, 500, length.out = 50),
  sunburns = beachgoers * 0.045 + rnorm(50, 0, 20),
  pair = "Beachgoers vs Sunburns"
)

df2 |> 
  mutate(
    # Sign of contribution to correlation (+ if both above or both below mean, - otherwise)
    sign_index = ifelse(
      (beachgoers > mean(beachgoers) & sunburns > mean(sunburns)) |
      (beachgoers < mean(beachgoers) & sunburns < mean(sunburns)),
      "+", "-"
    ),
    # Size index based on distance from means (scaled absolute contribution)
    size_index = abs((beachgoers - mean(beachgoers)) * 
                     (sunburns - mean(sunburns)))
  ) |>
  ggplot(aes(x = beachgoers, y = sunburns, color = sign_index)) +
  geom_text(aes(label = sign_index, size = size_index), alpha = 0.8) +
  geom_vline(xintercept = mean(df2$beachgoers), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df2$sunburns), color = "forestgreen",size=1.5,alpha=.5) +
  labs(
    # title = "Contribution Signs to Correlation",
    # subtitle = "Positive ( + ) when above/below both means, Negative ( - ) otherwise",
    x = "beachgoers",
    y = "sunburns"
  ) +
  scale_size(range = c(3, 8)) +
  theme_minimal() +
  theme(legend.position = "none")

# 2. Moderate Positive: Beachgoers vs sunburns
df3 <- tibble(
  ufo_sightings = rnorm(50, 20, 5),
  potato_sales = rnorm(50, 500, 30),
  pair = "UFO Sightings vs Potato Sales"
)


df3 |> 
  mutate(
    # Sign of contribution to correlation (+ if both above or both below mean, - otherwise)
    sign_index = ifelse(
      (ufo_sightings > mean(ufo_sightings) & potato_sales > mean(potato_sales)) |
      (ufo_sightings < mean(ufo_sightings) & potato_sales < mean(potato_sales)),
      "+", "-"
    ),
    # Size index based on distance from means (scaled absolute contribution)
    size_index = abs((ufo_sightings - mean(ufo_sightings)) * 
                     (potato_sales - mean(potato_sales)))
  ) |>
  ggplot(aes(x = ufo_sightings, y = potato_sales, color = sign_index)) +
  geom_text(aes(label = sign_index, size = size_index), alpha = 0.8) +
  geom_vline(xintercept = mean(df3$ufo_sightings), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df3$potato_sales), color = "forestgreen",size=1.5,alpha=.5) +
  labs(
    # title = "Contribution Signs to Correlation",
    # subtitle = "Positive ( + ) when above/below both means, Negative ( - ) otherwise",
    x = "ufo_sightings",
    y = "potato_sales"
  ) +
  scale_size(range = c(3, 8)) +
  theme_minimal() +
  theme(legend.position = "none")

df4 <- tibble(
  doc_hours = seq(0, 10, length.out = 50),
  social_life = 100 - doc_hours * 8 + rnorm(50, 0, 5),
  pair = "Alien Docs vs Social Life"
)

df4 |> 
  mutate(
    # Sign of contribution to correlation (+ if both above or both below mean, - otherwise)
    sign_index = ifelse(
      (doc_hours > mean(doc_hours) & social_life > mean(social_life)) |
      (doc_hours < mean(doc_hours) & social_life < mean(social_life)),
      "+", "-"
    ),
    # Size index based on distance from means (scaled absolute contribution)
    size_index = abs((doc_hours - mean(doc_hours)) * 
                     (social_life - mean(social_life)))
  ) |>
  ggplot(aes(x = doc_hours, y = social_life, color = sign_index)) +
  geom_text(aes(label = sign_index, size = size_index), alpha = 0.8) +
  geom_vline(xintercept = mean(df4$doc_hours), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df4$social_life), color = "forestgreen",size=1.5,alpha=.5) +
  labs(
    # title = "Contribution Signs to Correlation",
    # subtitle = "Positive ( + ) when above/below both means, Negative ( - ) otherwise",
    x = "doc_hours",
    y = "social_life"
  ) +
  scale_size(range = c(3, 8)) +
  theme_minimal() +
  theme(legend.position = "none")

# Simulate non-linear relationship (parabolic)
set.seed(123)
df_nonlinear <- tibble(
  x = seq(-10, 10, length.out = 100),
  y = x^2 + rnorm(100, 0, 5)  # strong but non-linear
)


# Plot
ggplot(df_nonlinear, aes(x = x, y = y)) +
  geom_point(color = "indianred", alpha = 0.7) +
  geom_vline(xintercept = mean(df_nonlinear$x), color = "dodgerblue",size=1.5,alpha=.5) +
  geom_hline(yintercept = mean(df_nonlinear$y), color = "forestgreen",size=1.5,alpha=.5)+
  # geom_smooth(method = "lm", se = FALSE, linetype = "dashed", color = "blue") +  # linear fit
  # geom_smooth(method = "loess", se = FALSE, color = "green") +  # non-linear fit
  labs(
    title = "Example of Non-Linear Correlation",
    # subtitle = "Quadratic relationship: linear fit fails, LOESS succeeds",
    x = "Variable X",
    y = "Variable Y"
  ) +
  theme_minimal()

library(tidyverse)
library(broom)

# Simulate data
set.seed(42)
alien_data <- tibble(
  alien_species = rep(c("Martian", "Venusian", "Zeta Reticulan", "Reptilian"), each = 30),
  laser_power = c(
    rnorm(30, mean = 100, sd = 10),  # Martians like moderate
    rnorm(30, mean = 130, sd = 15),  # Venusians love powerful beams
    rnorm(30, mean = 90, sd = 12),   # Zetas are conservative
    rnorm(30, mean = 110, sd = 8)    # Reptilians are balanced
  )
)

alien_groups=alien_data |> group_by(alien_species) |> 
  summarise(
    mean_laser_power = mean(laser_power),
    size_alien = n()
  ) 



# Visualize
alien_data |> ggplot(aes(x = alien_species, y = laser_power, fill = alien_species)) +
  # geom_violin(trim = FALSE, alpha = 0.7) +
  geom_jitter(aes(color=alien_species),width = 0.05, alpha = 0.5, size = 1.5) +
  geom_boxplot(width = 0.1, outlier.shape = NA, color = "black") +
  labs(
    title = "Laser Power Preferences by Alien Species",
    x = "Alien Species",
    y = "Laser Gun Power Setting"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Run ANOVA
anova_model <- aov(laser_power ~ alien_species, data = alien_data)
anova_summary <- summary(anova_model)

# Calculate eta squared manually
ss_total <- sum((alien_data$laser_power - mean(alien_data$laser_power))^2)
ss_between <- sum(alien_groups$size_alien * (alien_groups$mean_laser_power - mean(alien_data$laser_power))^2)
eta_squared <- ss_between / ss_total

library(tidyverse)
library(broom)
library(glue)

# Simulate data — no meaningful group differences
set.seed(1337)
alien_sleep <- tibble(
  alien_species = rep(c("Martian", "Venusian", "Zeta Reticulan", "Reptilian"), each = 30),
  sleep_hours = rnorm(120, mean = 8, sd = 1.5)  # same distribution for all
)

# Plot distributions
ggplot(alien_sleep, aes(x = alien_species, y = sleep_hours, fill = alien_species)) +
  geom_jitter(aes(color=alien_species),width = 0.05, alpha = 0.5, size = 1.5) +
  geom_boxplot(width = 0.1, outlier.shape = NA, color = "black") +
  labs(
    title = "Sleep Hours by Alien Species (No Real Difference)",
    x = "Alien Species",
    y = "Sleep Hours"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# ANOVA
anova_model <- aov(sleep_hours ~ alien_species, data = alien_sleep)
anova_summary <- summary(anova_model)

# Compute eta squared
ss_total <- sum((alien_sleep$sleep_hours - mean(alien_sleep$sleep_hours))^2)
ss_between <- anova_summary[[1]][["Sum Sq"]][1]
eta_squared <- ss_between / ss_total