Day 2: Significant testing & effect size calculating

x <-  c(2,3,5,9,8,11,6)
count <- 0
for (val in x) {
  if(val > 5) count = count +1
}
print(count)

## ---------------------------
## Simple nonparametric bootstrap demo
## ---------------------------

set.seed(5) # Make results reproducible (same random draws every run)

n_iter  <- 1000 # Number of bootstrap resamples
s_size  <- 20 # Size of each resample (should match original sample size)

# Simulate one sample
Treatment_A <- rnorm(s_size, mean = 4, sd = 2)  # 20 draws from normal distribution N(4, 2^2)
# Treatment_A <- c(Treatment_A, 30)             # <- add an outlier and re-run

# Preallocate storage for speed & clarity
mean_values <- vector() # define vector to store mean values of resampled datasets
sd_values   <- vector() # define vector to store standard deviation of resampled datasets

for (k in 1:n_iter) {
  # Draw a bootstrap sample: sample WITH replacement from the original sample
  data_tmp <- sample(Treatment_A, s_size, replace = TRUE) 
  
  # Store the statistic(s) of interest from this resample
  mean_values[k] <- mean(data_tmp) # the resampled MEAN
  sd_values[k]   <- sd(data_tmp) # the resampled SD
}

# --- Bootstrap summaries ---

# The bootstrap estimate of the mean (often close to the original sample mean)
bootstrapped_mean <- mean(mean_values)

# The average of resampled SDs (typical *within-sample* spread), NOT the SE of the mean
bootstrapped_sd   <- mean(sd_values)

# Original (non-bootstrapped) summaries from the observed sample
origi_mean        <- mean(Treatment_A)
origi_sd          <- sd(Treatment_A)

# A simple 95% percentile CI for the mean from the bootstrap distribution:
bootstrap_ci_mean <- quantile(mean_values, c(0.025, 0.975), names = FALSE)

# Collect key outputs
c(original_mean = origi_mean,
  boot_mean     = bootstrapped_mean,
  original_sd   = origi_sd,
  boot_sd   = bootstrapped_sd,
  CI_mean_L     = bootstrap_ci_mean[1],
  CI_mean_U     = bootstrap_ci_mean[2])

# ------------------------------------------------------------
# bootstrap_effects()
# Purpose:
#   Compute the observed mean difference (y - x) and Cohen's d
#   for two independent groups, and get 95% bootstrap CIs for
#   each effect size using a simple nonparametric (percentile)
#   bootstrap (resample within groups with replacement).
#
# Inputs:
#   x, y    : numeric vectors (independent samples; NAs allowed)
#   n_iter  : number of bootstrap resamples (default 1000)
#   seed    : RNG seed for reproducibility
#
# Outputs (list):
#   $mean_difference : c(estimate, CI_L, CI_U)
#   $cohens_d        : c(estimate, CI_L, CI_U)
#
# Notes:
#   - Cohen's d uses the pooled SD (assumes equal variances).
#   - CIs use the percentile method (simple and common for teaching).
#   - For small n, consider Hedges' g (bias-corrected d).
# ------------------------------------------------------------

bootstrap_effects <- function(x, y, n_iter = 1000, seed = 1) {
  set.seed(seed)
  x <- x[!is.na(x)]; y <- y[!is.na(y)]  #Drop missing values so resampling doesn't pull NAs
  n1 <- length(x); n2 <- length(y) # Sample sizes

  # ---- Observed (non-bootstrapped) effect sizes ----
  # Mean difference (y - x): positive => y has higher mean
  md_obs <- mean(y) - mean(x)
  
  # Pooled standard deviation (Cohen's d standardizer)
  # var() in R uses (n-1) denominator, so this matches the pooled-variance formula.
  sp     <- sqrt( ((n1-1)*var(x) + (n2-1)*var(y)) / (n1 + n2 - 2) )
  
  # Cohen's d (standardized mean difference)
  d_obs  <- (mean(y) - mean(x)) / sp
  
  # ---- Bootstrap loop ----
  # Allocate vectors to store bootstrap replicates
  md_boot <- numeric(n_iter); d_boot <- numeric(n_iter)
  
  for (k in 1:n_iter) {
    # Resample with replacement within each group (nonparametric bootstrap)
    x_k <- sample(x, n1, replace = TRUE)
    y_k <- sample(y, n2, replace = TRUE)
    
    # Effect size for this bootstrap sample
    md_boot[k] <- mean(y_k) - mean(x_k)
    
    # Recompute pooled SD from the bootstrap resamples
    sp_k <- sqrt( ((n1-1)*var(x_k) + (n2-1)*var(y_k)) / (n1 + n2 - 2) )
    d_boot[k]  <- (mean(y_k) - mean(x_k)) / sp_k
  }
  
  # ---- 95% percentile confidence intervals ----
  md_ci <- quantile(md_boot, c(.025,.975), names = FALSE)
  d_ci  <- quantile(d_boot,  c(.025,.975), names = FALSE)

  # ---- Return a tidy list ----
  output_list <- list(
    mean_difference = c(estimate = md_obs, CI_L = md_ci[1], CI_U = md_ci[2]),
    cohens_d        = c(estimate = d_obs,  CI_L = d_ci[1],  CI_U = d_ci[2])
  )
  return(output_list)
}

# Example:
# bootstrap_effects(control_WSA, drought_WSA, n_iter = 5000, seed = 123)