# selecting_tests.R
# Required packages (run once if you don't have them)
# install.packages(c("car","rmarkdown","knitr"))

library(car)       # for leveneTest
# library(rmarkdown) # not needed to load to render, but useful if producing docx

# ---------------------------
# Scenario 1: independent groups (revision class vs no)
# ---------------------------
scores <- data.frame(
  group = rep(c("yes", "no"), each = 15),
  score = c(85, 80, 78, 83, 81, 79, 84, 82, 87, 83, 85, 88, 80, 86, 84,
            72, 70, 69, 73, 75, 68, 71, 74, 70, 69, 72, 71, 68, 73, 70)
)

# Look
str(scores); aggregate(score ~ group, data = scores, FUN = function(x) c(mean=mean(x), sd=sd(x)))

# Assumptions
by(scores$score, scores$group, shapiro.test)    # normality by group
leveneTest(score ~ group, data = scores)        # homogeneity of variances

# Test: Welch two-sample t-test (safer than equal-variance t-test)
t1 <- t.test(score ~ group, data = scores, var.equal = FALSE)
t1
# If non-normal or you prefer non-parametric:
wilcox.test(score ~ group, data = scores)

# ---------------------------
# Scenario 2: paired data (blood pressure before vs after)
# ---------------------------
bp <- data.frame(
  before = c(140, 145, 138, 142, 141, 139, 144, 143, 140, 137,
             141, 142, 139, 143, 138, 140, 144, 146, 139, 141),
  after  = c(130, 132, 129, 134, 133, 130, 131, 132, 129, 128,
             131, 132, 130, 133, 129, 130, 134, 135, 130, 132)
)

# Check differences
diffs <- bp$before - bp$after
shapiro.test(diffs)   # normality of differences

# Paired t-test
t2 <- t.test(bp$before, bp$after, paired = TRUE)
t2
# If differences non-normal: wilcox.test(bp$before, bp$after, paired=TRUE)

# ---------------------------
# Scenario 3: two independent groups with similar means (supplement groups A vs B)
# ---------------------------
heart_data <- data.frame(
  group = rep(c("A", "B"), each = 20),
  bpm = c(72, 70, 71, 73, 74, 72, 70, 72, 73, 71,
          70, 72, 74, 72, 71, 73, 72, 72, 73, 71,
          72, 72, 73, 71, 72, 70, 72, 74, 73, 71,
          70, 71, 72, 73, 72, 70, 72, 71, 73, 72)
)

aggregate(bpm ~ group, data = heart_data, FUN = function(x) c(mean=mean(x), sd=sd(x)))
by(heart_data$bpm, heart_data$group, shapiro.test)
leveneTest(bpm ~ group, data = heart_data)

# Independent t-test (Welch)
t3 <- t.test(bpm ~ group, data = heart_data, var.equal = FALSE)
t3
# Non-parametric alternative (if needed)
wilcox.test(bpm ~ group, data = heart_data)

# ---------------------------
# Scenario 4: one-way ANOVA (sleep hours across age groups)
# ---------------------------
sleep_data <- data.frame(
  age_group = rep(c("young", "middle", "senior"), each = 15),
  hours = c(8.5, 8.2, 8.1, 7.8, 8.0, 8.3, 7.9, 8.4, 8.2, 8.1, 8.0, 7.7, 7.9, 8.3, 8.0,
            6.5, 6.3, 6.6, 6.4, 6.2, 6.5, 6.6, 6.4, 6.3, 6.1, 6.4, 6.2, 6.3, 6.5, 6.4,
            5.4, 5.2, 5.5, 5.3, 5.1, 5.3, 5.4, 5.2, 5.3, 5.5, 5.2, 5.3, 5.4, 5.1, 5.3)
)

aggregate(hours ~ age_group, data = sleep_data, FUN = function(x) c(mean=mean(x), sd=sd(x)))
# Assumptions
by(sleep_data$hours, sleep_data$age_group, shapiro.test)
leveneTest(hours ~ age_group, data = sleep_data)

# One-way ANOVA
anova_model <- aov(hours ~ age_group, data = sleep_data)
summary(anova_model)
TukeyHSD(anova_model)   # post-hoc pairwise comparisons

# If ANOVA assumptions badly violated: kruskal.test(hours ~ age_group, data = sleep_data)

# ---------------------------
# Scenario 5: correlation (study hours vs score)
# ---------------------------
study <- data.frame(
  hours = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
  score = c(55, 60, 65, 70, 74, 78, 82, 85, 88, 90)
)

plot(study$hours, study$score, main = "Study hours vs score", xlab = "Hours", ylab = "Score")
cor.test(study$hours, study$score, method = "pearson")

# ---------------------------
# Scenario 6: chi-square (gender vs preference for remote work)
# ---------------------------
gender <- c(rep("Male", 20), rep("Female", 20))
preference <- c("Remote", "Office", "Remote", "Remote", "Remote", "Office", "Remote", "Remote", "Office", "Remote",
                "Office", "Remote", "Remote", "Remote", "Remote", "Office", "Remote", "Office", "Remote", "Remote",
                "Office", "Office", "Office", "Remote", "Office", "Remote", "Office", "Office", "Remote", "Office",
                "Remote", "Remote", "Office", "Remote", "Remote", "Office", "Office", "Remote", "Remote", "Office")

table_data <- table(gender, preference)
table_data
chisq.test(table_data)   # if expected counts <5 in many cells, use fisher.test(table_data)

# ---------------------------
# Scenario 7: one-way ANOVA (work hours across departments)
# ---------------------------
work_hours <- data.frame(
  department = rep(c("HR", "Finance", "IT"), each = 15),
  hours = c(
    38, 39, 40, 41, 39, 40, 38, 40, 39, 41, 40, 39, 38, 40, 39,
    39, 40, 41, 40, 39, 40, 38, 41, 40, 39, 40, 40, 39, 40, 39,
    40, 41, 39, 40, 40, 38, 39, 40, 41, 40, 39, 39, 40, 40, 40)
)

aggregate(hours ~ department, data = work_hours, FUN = function(x) c(mean=mean(x), sd=sd(x)))
by(work_hours$hours, work_hours$department, shapiro.test)
leveneTest(hours ~ department, data = work_hours)

anova_wh <- aov(hours ~ department, data = work_hours)
summary(anova_wh)
# if nonparametric alternative: kruskal.test(hours ~ department, data = work_hours)

# ---------------------------
# Save results suggestion: create a markdown / Word document
# ---------------------------
# Minimal RMarkdown file (save as "results_report.Rmd") and then knit to Word:
rmd_text <- '
---
title: "Statistical test results"
output: word_document
---

```{r, echo=TRUE}
# source the script with all objects (optional)
source("selecting_tests.R")