Statistical Framework

set.seed(42)

# Generate example data
n <- 50
x <- rnorm(n)
y_linear <- x + rnorm(n, sd = 0.5)
y_monotonic <- sign(x) * abs(x)^0.5 + rnorm(n, sd = 0.3)
y_outlier <- y_linear
y_outlier[c(1, 2)] <- c(10, -10)  # Add outliers

# Calculate correlations
results <- data.frame(
  Scenario = c("Linear + Noise", "Monotonic (Non-linear)", "With Outliers"),
  Pearson = c(cor(x, y_linear), cor(x, y_monotonic), cor(x, y_outlier)),
  Spearman = c(cor(x, y_linear, method = "spearman"),
               cor(x, y_monotonic, method = "spearman"),
               cor(x, y_outlier, method = "spearman"))
)

# Visualization
df_plot <- data.frame(
  x = rep(x, 3),
  y = c(y_linear, y_monotonic, y_outlier),
  Scenario = factor(rep(c("Linear + Noise", "Monotonic (Non-linear)", "With Outliers"), 
                        each = n))
)

ggplot(df_plot, aes(x = x, y = y)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red", linetype = "dashed") +
  facet_wrap(~Scenario, scales = "free_y") +
  labs(
    title = "Pearson vs Spearman Correlation",
    subtitle = "Spearman is more robust to outliers and non-linear relationships",
    x = "Variable X",
    y = "Variable Y"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold"))

knitr::kable(results, digits = 3,
             caption = "Comparison of Pearson and Spearman correlations")

# Simulate sampling distributions
set.seed(123)
n_samples <- 1000
sample_size <- 30

simulate_cors <- function(true_rho, n, n_samples) {
  cors <- sapply(1:n_samples, function(i) {
    sigma <- matrix(c(1, true_rho, true_rho, 1), 2, 2)
    data <- MASS::mvrnorm(n, mu = c(0, 0), Sigma = sigma)
    cor(data[,1], data[,2])
  })
  return(cors)
}

# Different true correlations
rhos <- c(0, 0.5, 0.8)
results_list <- lapply(rhos, function(r) {
  cors <- simulate_cors(r, sample_size, n_samples)
  z <- 0.5 * log((1 + cors) / (1 - cors))
  data.frame(
    Correlation = cors,
    FisherZ = z,
    TrueRho = paste("ρ =", r)
  )
})

df_all <- do.call(rbind, results_list)

# Plot correlation distributions
p1 <- ggplot(df_all, aes(x = Correlation, fill = TrueRho)) +
  geom_density(alpha = 0.5) +
  labs(title = "A) Raw Correlation Distributions",
       x = "Sample Correlation", y = "Density") +
  theme_minimal() +
  theme(legend.position = "top")

# Plot Fisher Z distributions
p2 <- ggplot(df_all, aes(x = FisherZ, fill = TrueRho)) +
  geom_density(alpha = 0.5) +
  labs(title = "B) Fisher Z Distributions",
       x = "Fisher's Z", y = "Density") +
  theme_minimal() +
  theme(legend.position = "top")

gridExtra::grid.arrange(p1, p2, ncol = 2)

# Example design matrix
n_cells <- 20
cell_data <- data.frame(
  Cell = paste0("Cell_", 1:n_cells),
  Subject = rep(c("S1", "S2", "S3", "S4"), each = 5),
  Sensitivity = c(rep("Sensitive", 10), rep("Resistant", 10))
)

# Create design matrix
design <- model.matrix(~ Subject + Sensitivity, data = cell_data)
colnames(design) <- c("Intercept", "Subject_S2", "Subject_S3", "Subject_S4", "Sensitive")

knitr::kable(head(design, 10), caption = "Example Design Matrix (first 10 cells)")

# Simulate the benefit of empirical Bayes
set.seed(42)
n_compounds <- 500
true_effects <- rnorm(n_compounds, mean = 0, sd = 0.5)
n_cells_per_group <- 10

# Simulate with different sample sizes
simulate_effect <- function(true_effect, n) {
  sensitive <- rnorm(n, mean = true_effect, sd = 1)
  resistant <- rnorm(n, mean = 0, sd = 1)
  t.test(sensitive, resistant)$statistic
}

t_stats <- sapply(true_effects, simulate_effect, n = n_cells_per_group)

# Simple posterior (shrinkage toward zero)
shrinkage <- 0.7  # Simplified representation
moderated_t <- shrinkage * t_stats

df_ebayes <- data.frame(
  True = true_effects,
  Ordinary = t_stats,
  Moderated = moderated_t
)

ggplot(df_ebayes) +
  geom_point(aes(x = True, y = Ordinary), alpha = 0.3, color = "red") +
  geom_point(aes(x = True, y = Moderated), alpha = 0.3, color = "blue") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  labs(
    title = "Empirical Bayes Moderation Effect",
    subtitle = "Blue (moderated) estimates are shrunk toward zero, reducing false positives",
    x = "True Effect Size",
    y = "Estimated t-statistic"
  ) +
  theme_minimal() +
  annotate("text", x = -1.5, y = 3, label = "Ordinary t-test", color = "red") +
  annotate("text", x = -1.5, y = 2.5, label = "Moderated t-test", color = "blue")

# Demonstrate multiple testing
set.seed(42)
n_tests <- 1679
n_true_effects <- 50

# Generate p-values
true_p <- c(runif(n_true_effects, 0, 0.01),  # True effects
            runif(n_tests - n_true_effects, 0, 1))  # Null

# Apply corrections
p_bonferroni <- p.adjust(true_p, method = "bonferroni")
p_fdr <- p.adjust(true_p, method = "fdr")

results_mult <- data.frame(
  Raw = true_p,
  Bonferroni = p_bonferroni,
  FDR = p_fdr,
  TrueEffect = c(rep("Yes", n_true_effects), 
                  rep("No", n_tests - n_true_effects))
)

# Count discoveries at alpha = 0.05
discoveries <- data.frame(
  Method = c("Raw p-value", "Bonferroni", "FDR (BH)"),
  `Total Discoveries` = c(sum(true_p < 0.05),
                           sum(p_bonferroni < 0.05),
                           sum(p_fdr < 0.05)),
  `True Positives` = c(sum(true_p[1:n_true_effects] < 0.05),
                        sum(p_bonferroni[1:n_true_effects] < 0.05),
                        sum(p_fdr[1:n_true_effects] < 0.05)),
  check.names = FALSE
)
discoveries$`False Positives` <- discoveries$`Total Discoveries` - discoveries$`True Positives`

knitr::kable(discoveries, 
             caption = "Multiple Testing Correction Comparison (α = 0.05)")

# Effect size interpretation guide
effect_sizes <- data.frame(
  logFC = c(-0.5, -0.3, -0.1, 0, 0.1, 0.3, 0.5),
  Interpretation = c("Strong negative", "Moderate negative", "Weak negative",
                     "No effect", "Weak positive", "Moderate positive", "Strong positive"),
  Meaning = c("Much lower connectivity in resistant",
              "Lower connectivity in resistant",
              "Slightly lower in resistant",
              "No difference",
              "Slightly higher in resistant",
              "Higher connectivity in resistant",
              "Much higher connectivity in resistant")
)

knitr::kable(effect_sizes, 
             caption = "Interpretation of logFC values in differential connectivity")

# Power calculation
power_calc <- function(n, effect, alpha = 0.05) {
  se <- sqrt(2/n)
  z_crit <- qnorm(1 - alpha/2)
  power <- pnorm(effect/se - z_crit) + pnorm(-effect/se - z_crit)
  return(power)
}

n_range <- seq(10, 500, by = 10)
effects <- c(0.1, 0.2, 0.3, 0.5)

power_df <- expand.grid(n = n_range, effect = effects)
power_df$power <- mapply(power_calc, power_df$n, power_df$effect)
power_df$effect <- factor(power_df$effect)

ggplot(power_df, aes(x = n, y = power, color = effect)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0.8, linetype = "dashed", color = "gray50") +
  scale_color_brewer(palette = "Set1", name = "Effect Size\n(logFC)") +
  labs(
    title = "Statistical Power vs Sample Size",
    subtitle = "Dashed line indicates 80% power threshold",
    x = "Number of Cells per Group",
    y = "Statistical Power"
  ) +
  theme_minimal() +
  annotate("text", x = 400, y = 0.85, label = "80% power", color = "gray40")

# Example gene coverage visualization
coverage_data <- data.frame(
  Category = c("L1000 Genes (978)", "In Your Data", "Used for Analysis"),
  Count = c(978, 850, 850),
  Percentage = c(100, 86.9, 86.9)
)

ggplot(coverage_data, aes(x = Category, y = Count, fill = Category)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = paste0(Count, " (", round(Percentage, 1), "%)")),
            vjust = -0.5) +
  scale_fill_brewer(palette = "Blues") +
  labs(
    title = "Gene Coverage Quality Control",
    x = "",
    y = "Number of Genes"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  ylim(0, 1100)