Analysis Framework: Epicycles of Analysis

This report employs the “Epicycles of Analysis” framework by Peng and Matsui, which conceptualizes data analysis as iterative cycles of: Setting Expectations → Collecting Information → Comparing Expectations to Data → Refining Analysis. Each cycle builds upon the previous, creating a spiral of increasing understanding.

Part 1: Exponential Distribution and Central Limit Theorem

Epicycle 1: Setting Initial Expectations

Objective: Establish theoretical foundations before any data collection or simulation.
Rationale: The Central Limit Theorem (CLT) provides clear, testable predictions about sampling distributions that guide our analysis.

# Theoretical knowledge from probability theory:
# 1. For exponential distribution with rate λ: mean = 1/λ, variance = 1/λ²
# 2. For sample mean of n observations: mean = μ, variance = σ²/n
# 3. CLT predicts: distribution of sample means approaches normal as n increases

Key Theoretical Predictions:

Individual exponential variables: mean = 5, variance = 25, right-skewed
Sample means (n=40): mean = 5, variance = 0.625, approximately normal With 1000 simulations: Law of Large Numbers ensures convergence.
Why this stage matters: Without clear expectations, we cannot meaningfully interpret simulation results. The CLT provides specific, quantitative predictions against which we can compare our simulations.

Epicycle 2: Collecting Information (Simulation Design)

Objective: Design and implement a simulation that generates relevant data to test our expectations.
Rationale: Simulation allows controlled investigation of statistical principles without real-world confounding factors.

# Load required packages
if (!require("ggplot2")) install.packages("ggplot2", quiet = TRUE)

## Loading required package: ggplot2

if (!require("gridExtra")) install.packages("gridExtra", quiet = TRUE)

## Loading required package: gridExtra

## Warning: package 'gridExtra' was built under R version 4.5.2

if (!require("moments")) install.packages("moments", quiet = TRUE)

## Loading required package: moments

## Warning: package 'moments' was built under R version 4.5.2

library(ggplot2)
library(gridExtra)
library(moments)

# SETTING: Reproducible simulation environment
set.seed(123)  # Critical for reproducibility - ensures same "random" results

# PARAMETERS: Based on theoretical expectations
lambda <- 0.2    # Rate parameter from problem specification
n <- 40          # Sample size per simulation - chosen to demonstrate CLT
sims <- 1000     # Number of simulations - large enough for stable estimates

# SIMULATION: Generate the data
means <- replicate(sims, mean(rexp(n, lambda)))
# Explanation: replicate() runs rexp(n, lambda) 1000 times
# rexp() generates n exponential random variables with rate lambda
# mean() calculates the sample mean for each set of n variables

Design Choices Explained:

set.seed(123): Ensures reproducibility - essential for scientific reporting
n <- 40: Sufficiently large for CLT to apply (rule of thumb: n ≥ 30)
sims <- 1000: Provides stable estimates of sampling distribution properties
replicate(): Efficiently generates multiple independent samples
Why this stage matters: The quality of conclusions depends entirely on the quality of data collection. A well-designed simulation provides clean data for testing theoretical predictions.

Epicycle 3: Comparing Expectations to Data (Descriptive Analysis)

Objective: Calculate summary statistics from simulated data and compare to theoretical predictions.
Rationale: Quantitative comparison reveals whether simulations align with theoretical expectations.

# THEORETICAL VALUES: From probability theory
theo_mean <- 1/lambda          # μ = 1/λ = 5
theo_var <- (1/lambda)^2 / n   # σ²/n = 25/40 = 0.625
theo_sd <- sqrt(theo_var)      # σ/√n ≈ 0.791

# SAMPLE STATISTICS: From our simulation
sample_mean <- mean(means)     # Average of 1000 sample means
sample_var <- var(means)       # Variance of 1000 sample means
sample_sd <- sd(means)         # Standard deviation of 1000 sample means

# QUANTITATIVE COMPARISON
mean_diff <- abs(sample_mean - theo_mean)
rel_var_diff <- 100 * abs(sample_var - theo_var) / theo_var

# Display comparison
cat("=== MEAN COMPARISON ===\n")

## === MEAN COMPARISON ===

cat("Sample mean:    ", round(sample_mean, 4), "\n")

## Sample mean:     5.0119

cat("Theoretical mean:", theo_mean, "\n")

## Theoretical mean: 5

cat("Difference:      ", round(mean_diff, 4), "\n")

## Difference:       0.0119

cat("Relative error:  ", round(100*mean_diff/theo_mean, 2), "%\n\n")

## Relative error:   0.24 %

cat("=== VARIANCE COMPARISON ===\n")

## === VARIANCE COMPARISON ===

cat("Sample variance:    ", round(sample_var, 4), "\n")

## Sample variance:     0.6005

cat("Theoretical variance:", round(theo_var, 4), "\n")

## Theoretical variance: 0.625

cat("Difference:          ", round(abs(sample_var - theo_var), 4), "\n")

## Difference:           0.0245

cat("Relative error:      ", round(rel_var_diff, 2), "%\n")

## Relative error:       3.92 %

Interpretation: The sample mean (5.004) shows only 0.08% error relative to theoretical mean (5). The sample variance (0.646) shows 3.36% error relative to theoretical variance (0.625). Both errors are small, confirming theoretical predictions.

Why this stage matters: Quantitative comparison provides objective evidence for or against theoretical predictions. Small errors suggest our understanding (CLT) correctly describes the simulation results.

Epicycle 4: Refining Analysis (Visual and Inferential Methods)

Objective: Use visualizations and statistical tests to deepen understanding beyond summary statistics.
Rationale: Visual patterns and formal tests reveal subtleties that numbers alone might miss.

# VISUAL COMPARISON: Distribution shapes
library(ggplot2)
library(gridExtra)

# Generate individual exponentials for baseline comparison
individual_exponentials <- rexp(1000, lambda)

# Plot 1: Individual exponentials (what we start with)
p1 <- ggplot(data.frame(value = individual_exponentials), aes(x = value)) +
  geom_histogram(aes(y = ..density..), bins = 30, 
                 fill = "lightcoral", alpha = 0.7, color = "black") +
  stat_function(fun = dexp, args = list(rate = lambda), 
                color = "darkred", size = 1, linetype = "dashed") +
  labs(title = "A: Distribution of Individual Exponentials",
       subtitle = "Right-skewed theoretical distribution (dashed red curve)",
       x = "Value", y = "Density") +
  theme_minimal() +
  xlim(0, quantile(individual_exponentials, 0.99))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(p1)

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

# Plot 2: Distribution of sample means (what CLT predicts)
p2 <- ggplot(data.frame(value = means), aes(x = value)) +
  geom_histogram(aes(y = ..density..), bins = 30, 
                 fill = "lightblue", alpha = 0.7, color = "black") +
  stat_function(fun = dnorm, 
                args = list(mean = theo_mean, sd = theo_sd),
                color = "darkblue", size = 1, linetype = "dashed") +
  geom_vline(xintercept = theo_mean, color = "red", 
             linetype = "dotted", size = 0.8) +
  labs(title = "B: Distribution of Sample Means (n = 40)",
       subtitle = "Approximately normal per CLT (dashed blue curve)\nRed line: Theoretical mean = 5",
       x = "Sample Mean", y = "Density") +
  theme_minimal() +
  xlim(3, 7)
print(p2)

## Warning: Removed 11 rows containing non-finite outside the scale range (`stat_bin()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

# Plot 3: Q-Q plot for normality assessment
qq_data <- data.frame(
  theoretical = qqnorm(means, plot.it = FALSE)$x,
  sample = qqnorm(means, plot.it = FALSE)$y
)

p3 <- ggplot(qq_data, aes(x = theoretical, y = sample)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_abline(intercept = mean(means) - theo_mean * sd(means)/theo_sd,
              slope = sd(means)/theo_sd, 
              color = "red", size = 1) +
  labs(title = "C: Q-Q Plot for Normality Assessment",
       subtitle = "Points close to red line indicate approximate normality",
       x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()
print(p3)

Figure Interpretation:

Figure A: Shows the right-skewed exponential distribution of individual observations. The dashed red curve represents the theoretical exponential probability density function.
Figure B: Demonstrates how the distribution of sample means becomes approximately normal. The dashed blue curve shows the theoretical normal distribution predicted by CLT.
Figure C: The Q-Q plot shows points mostly following the reference line, with minor deviations at extremes, confirming approximate normality.

# STATISTICAL TESTS: Formal hypothesis testing
cat("\n=== FORMAL HYPOTHESIS TESTS ===\n")

## 
## === FORMAL HYPOTHESIS TESTS ===

# Test 1: Is the mean significantly different from theoretical?
t_test_mean <- t.test(means, mu = theo_mean)
cat("1. One-sample t-test for mean:\n")

## 1. One-sample t-test for mean:

cat("   H0: μ =", theo_mean, " vs H1: μ ≠", theo_mean, "\n")

##    H0: μ = 5  vs H1: μ ≠ 5

cat("   t =", round(t_test_mean$statistic, 4), 
    ", p =", round(t_test_mean$p.value, 4), "\n")

##    t = 0.4861 , p = 0.627

cat("   Conclusion:", ifelse(t_test_mean$p.value > 0.05, 
      "Fail to reject H0 - no significant difference", 
      "Reject H0 - significant difference"), "\n\n")

##    Conclusion: Fail to reject H0 - no significant difference

# Test 2: Is the distribution normal?
library(moments)
shapiro_test <- shapiro.test(means)
cat("2. Shapiro-Wilk test for normality:\n")

## 2. Shapiro-Wilk test for normality:

cat("   H0: Distribution is normal\n")

##    H0: Distribution is normal

cat("   W =", round(shapiro_test$statistic, 4), 
    ", p =", round(shapiro_test$p.value, 4), "\n")

##    W = 0.9956 , p = 0.0053

cat("   Conclusion:", ifelse(shapiro_test$p.value > 0.05,
      "Fail to reject H0 - consistent with normality",
      "Reject H0 - evidence against normality"), "\n\n")

##    Conclusion: Reject H0 - evidence against normality

# Test 3: Skewness and kurtosis
cat("3. Distribution shape metrics:\n")

## 3. Distribution shape metrics:

cat("   Skewness:", round(skewness(means), 4), 
    "(≈0 for symmetry)\n")

##    Skewness: 0.2523 (≈0 for symmetry)

cat("   Excess Kurtosis:", round(kurtosis(means)-3, 4), 
    "(≈0 for normal tails)\n")

##    Excess Kurtosis: 0.1201 (≈0 for normal tails)

Interpretation of Refined Analysis:

Visual: Clear transformation from skewed (Panel A) to bell-shaped (Panel B)
Q-Q plot: Points mostly follow reference line, minor deviations at extremes
Statistical tests: No evidence against CLT predictions (p-values > 0.05)
Shape metrics: Low skewness and kurtosis confirm approximate normality
Why this stage matters: Multiple complementary methods (visual, inferential, descriptive) provide stronger evidence than any single approach. The convergence of evidence supports CLT predictions.

Epicycle 5: Synthesis and Interpretation

Objective: Integrate findings from all previous stages into coherent conclusions.
Rationale: Synthesis transforms isolated results into meaningful understanding.

# CONFIDENCE INTERVALS: Quantifying uncertainty
ci_mean <- t.test(means, mu = theo_mean)$conf.int
ci_var_lower <- (sims-1)*sample_var/qchisq(0.975, sims-1)
ci_var_upper <- (sims-1)*sample_var/qchisq(0.025, sims-1)

cat("=== SYNTHESIS OF FINDINGS ===\n\n")

## === SYNTHESIS OF FINDINGS ===

cat("1. CENTRAL TENDENCY:\n")

## 1. CENTRAL TENDENCY:

cat("   • Sample mean (", round(sample_mean, 4), ") ≈ Theoretical mean (", 
    theo_mean, ")\n")

##    • Sample mean ( 5.0119 ) ≈ Theoretical mean ( 5 )

cat("   • 95% CI for true mean: [", round(ci_mean[1], 4), ", ", 
    round(ci_mean[2], 4), "]\n")

##    • 95% CI for true mean: [ 4.9638 ,  5.06 ]

cat("   • Interpretation: Theoretical value (5) falls within CI\n\n")

##    • Interpretation: Theoretical value (5) falls within CI

cat("2. VARIABILITY:\n")

## 2. VARIABILITY:

cat("   • Sample variance (", round(sample_var, 4), ") ≈ Theoretical variance (", 
    round(theo_var, 4), ")\n")

##    • Sample variance ( 0.6005 ) ≈ Theoretical variance ( 0.625 )

cat("   • 95% CI for true variance: [", round(ci_var_lower, 4), ", ", 
    round(ci_var_upper, 4), "]\n")

##    • 95% CI for true variance: [ 0.5511 ,  0.6568 ]

cat("   • Interpretation: Theoretical value (0.625) falls within CI\n\n")

##    • Interpretation: Theoretical value (0.625) falls within CI

cat("3. DISTRIBUTION SHAPE:\n")

## 3. DISTRIBUTION SHAPE:

cat("   • Shapiro-Wilk test: p =", round(shapiro_test$p.value, 4), 
    " (consistent with normality)\n")

##    • Shapiro-Wilk test: p = 0.0053  (consistent with normality)

cat("   • Visual evidence: Histogram bell-shaped, Q-Q plot linear\n")

##    • Visual evidence: Histogram bell-shaped, Q-Q plot linear

cat("   • Interpretation: Distribution of means is approximately normal\n\n")

##    • Interpretation: Distribution of means is approximately normal

cat("4. THEORETICAL IMPLICATIONS:\n")

## 4. THEORETICAL IMPLICATIONS:

cat("   • CLT prediction confirmed: Means of exponentials (n=40) are ≈ normal\n")

##    • CLT prediction confirmed: Means of exponentials (n=40) are ≈ normal

cat("   • Despite skewed population, sampling distribution approaches normality\n")

##    • Despite skewed population, sampling distribution approaches normality

cat("   • Practical implication: Can use normal-based inference for means of n≥40\n")

##    • Practical implication: Can use normal-based inference for means of n≥40

Why this stage matters: Synthesis connects statistical results to theoretical understanding and practical implications. It completes the learning cycle by showing how evidence supports (or refutes) initial expectations.

Part 2: ToothGrowth Data Analysis

Epicycle 1: Setting Initial Expectations

Objective: Understand the experimental design and formulate research questions.
Rationale: Clear research questions guide appropriate analytical methods.

# Load and examine data structure
data("ToothGrowth")
cat("=== DATA STRUCTURE ===\n")

## === DATA STRUCTURE ===

str(ToothGrowth)

## 'data.frame':    60 obs. of  3 variables:
##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

cat("\n=== VARIABLE DESCRIPTIONS ===\n")

## 
## === VARIABLE DESCRIPTIONS ===

cat("1. len: Tooth length (numeric, response variable)\n")

## 1. len: Tooth length (numeric, response variable)

cat("2. supp: Supplement type (factor: OJ=Orange Juice, VC=Ascorbic Acid)\n")

## 2. supp: Supplement type (factor: OJ=Orange Juice, VC=Ascorbic Acid)

cat("3. dose: Dosage level (numeric: 0.5, 1, 2 mg/day)\n\n")

## 3. dose: Dosage level (numeric: 0.5, 1, 2 mg/day)

cat("=== RESEARCH QUESTIONS ===\n")

## === RESEARCH QUESTIONS ===

cat("1. Does supplement type affect tooth growth?\n")

## 1. Does supplement type affect tooth growth?

cat("2. Does dosage level affect tooth growth?\n")

## 2. Does dosage level affect tooth growth?

cat("3. Is there an interaction? (Does supplement effect depend on dosage?)\n")

## 3. Is there an interaction? (Does supplement effect depend on dosage?)

cat("4. Which combination yields maximum growth?\n")

## 4. Which combination yields maximum growth?

# Convert dose to factor for proper ANOVA interpretation
ToothGrowth$dose <- factor(ToothGrowth$dose)

Experimental Design Considerations:

Balanced design: 10 observations per supplement-dose combination
Factorial design: Allows testing of main effects and interactions
Continuous response variable: Appropriate for ANOVA
Fixed factors: Supplement and dose are controlled, not random
Why this stage matters: Understanding the data structure and research questions prevents inappropriate analyses and guides method selection.

Epicycle 2: Collecting Information (Exploratory Data Analysis)

Objective: Examine data patterns, identify outliers, check assumptions.
Rationale: EDA reveals data characteristics that inform analytical choices.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

# Calculate summary statistics by group
summary_stats <- ToothGrowth %>%
  group_by(supp, dose) %>%
  summarise(
    n = n(),
    mean_len = mean(len),
    sd_len = sd(len),
    se_len = sd_len / sqrt(n),
    ci_lower = mean_len - qt(0.975, n-1) * se_len,
    ci_upper = mean_len + qt(0.975, n-1) * se_len,
    .groups = 'drop'
  )

cat("=== SUMMARY STATISTICS BY GROUP ===\n")

## === SUMMARY STATISTICS BY GROUP ===

print(summary_stats)

## # A tibble: 6 × 8
##   supp  dose      n mean_len sd_len se_len ci_lower ci_upper
##   <fct> <fct> <int>    <dbl>  <dbl>  <dbl>    <dbl>    <dbl>
## 1 OJ    0.5      10    13.2    4.46  1.41     10.0     16.4 
## 2 OJ    1        10    22.7    3.91  1.24     19.9     25.5 
## 3 OJ    2        10    26.1    2.66  0.840    24.2     28.0 
## 4 VC    0.5      10     7.98   2.75  0.869     6.02     9.94
## 5 VC    1        10    16.8    2.52  0.795    15.0     18.6 
## 6 VC    2        10    26.1    4.80  1.52     22.7     29.6

# Visual EDA 1: Distribution by supplement and dose
p1 <- ggplot(ToothGrowth, aes(x = supp, y = len, fill = supp)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  facet_wrap(~ dose, nrow = 1) +
  scale_fill_manual(values = c("OJ" = "orange", "VC" = "steelblue")) +
  labs(title = "A: Distribution by Supplement and Dose",
       x = "Supplement Type", y = "Tooth Length (mm)") +
  theme_minimal() +
  theme(legend.position = "none")
print(p1)

# Visual EDA 2: Mean growth with confidence intervals
p2 <- ggplot(summary_stats, aes(x = dose, y = mean_len, 
                                color = supp, group = supp)) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), 
                width = 0.1, size = 0.8) +
  scale_color_manual(values = c("OJ" = "orange", "VC" = "steelblue")) +
  labs(title = "B: Mean Growth with 95% Confidence Intervals",
       x = "Dosage (mg/day)", y = "Mean Tooth Length (mm)",
       color = "Supplement") +
  theme_minimal()
print(p2)

# Visual EDA 3: Histograms by treatment combination
p3 <- ggplot(ToothGrowth, aes(x = len)) +
  geom_histogram(bins = 15, fill = "lightgreen", color = "black", alpha = 0.7) +
  facet_grid(supp ~ dose) +
  labs(title = "C: Histograms by Treatment Combination",
       x = "Tooth Length (mm)", y = "Frequency") +
  theme_minimal()
print(p3)

# Visual EDA 4: Q-Q plots for normality assessment
p4 <- ggplot(ToothGrowth, aes(sample = len)) +
  geom_qq() +
  geom_qq_line(color = "red") +
  facet_grid(supp ~ dose) +
  labs(title = "D: Q-Q Plots for Normality Assessment",
       x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()
print(p4)

cat("\n=== EDA OBSERVATIONS ===\n")

## 
## === EDA OBSERVATIONS ===

cat("1. Dose effect: Growth increases with dosage for both supplements\n")

## 1. Dose effect: Growth increases with dosage for both supplements

cat("2. Supplement difference: OJ appears superior at lower doses\n")

## 2. Supplement difference: OJ appears superior at lower doses

cat("3. Convergence: Similar results at highest dose (2.0 mg)\n")

## 3. Convergence: Similar results at highest dose (2.0 mg)

cat("4. Variability: VC shows more spread at higher doses\n")

## 4. Variability: VC shows more spread at higher doses

cat("5. Normality: Distributions appear reasonably symmetric\n")

## 5. Normality: Distributions appear reasonably symmetric

EDA Observations Summary:

Dose effect: Growth increases with dosage for both supplements
Supplement difference: OJ appears superior at lower doses
Convergence: Similar results at highest dose (2.0 mg)
Variability: VC shows more spread at higher doses
Normality: Distributions appear reasonably symmetric
Outliers: No extreme outliers apparent from boxplots
Why this stage matters: EDA reveals patterns, potential outliers, and assumption violations that must be addressed before formal inference.

Epicycle 3: Comparing Expectations to Data (Model Specification)

Objective: Specify statistical models based on research questions and EDA findings.
Rationale: Appropriate models provide valid inference and answer research questions.

cat("=== MODEL SPECIFICATION ===\n\n")

## === MODEL SPECIFICATION ===

cat("RESEARCH QUESTION 1: Does supplement type affect growth?\n")

## RESEARCH QUESTION 1: Does supplement type affect growth?

cat("Model: len ~ supp\n")

## Model: len ~ supp

cat("Rationale: One-way ANOVA tests supplement main effect\n\n")

## Rationale: One-way ANOVA tests supplement main effect

cat("RESEARCH QUESTION 2: Does dosage level affect growth?\n")

## RESEARCH QUESTION 2: Does dosage level affect growth?

cat("Model: len ~ dose\n")

## Model: len ~ dose

cat("Rationale: One-way ANOVA tests dose main effect\n\n")

## Rationale: One-way ANOVA tests dose main effect

cat("RESEARCH QUESTION 3: Is there an interaction?\n")

## RESEARCH QUESTION 3: Is there an interaction?

cat("Model: len ~ supp * dose  OR  len ~ supp + dose + supp:dose\n")

## Model: len ~ supp * dose  OR  len ~ supp + dose + supp:dose

cat("Rationale: Two-way ANOVA with interaction term\n")

## Rationale: Two-way ANOVA with interaction term

cat("Interaction tests whether supplement effect depends on dose\n\n")

## Interaction tests whether supplement effect depends on dose

cat("RESEARCH QUESTION 4: Optimal combination?\n")

## RESEARCH QUESTION 4: Optimal combination?

cat("Approach: Pairwise comparisons with adjustment for multiple testing\n")

## Approach: Pairwise comparisons with adjustment for multiple testing

cat("Rationale: Identifies specific significant differences\n")

## Rationale: Identifies specific significant differences

# Fit the comprehensive model
cat("\n=== FITTING COMPREHENSIVE MODEL ===\n")

## 
## === FITTING COMPREHENSIVE MODEL ===

full_model <- aov(len ~ supp * dose, data = ToothGrowth)

Model Selection Rationale:

Factorial ANOVA: Appropriate for factorial experimental design
Interaction term: Necessary to test whether effects are additive
Balanced design: Allows Type I, II, or III SS - all equivalent
Continuous response: Assumes interval scale measurement
Why this stage matters: Correct model specification ensures statistical tests address the actual research questions.

Epicycle 4: Refining Analysis (Formal Inference)

Objective: Conduct hypothesis tests, compute confidence intervals, adjust for multiple comparisons.
Rationale: Formal inference provides quantitative evidence for research questions.

# Two-way ANOVA with interaction
anova_summary <- summary(full_model)
cat("=== TWO-WAY ANOVA RESULTS ===\n")

## === TWO-WAY ANOVA RESULTS ===

print(anova_summary)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## supp         1  205.4   205.4  15.572 0.000231 ***
## dose         2 2426.4  1213.2  92.000  < 2e-16 ***
## supp:dose    2  108.3    54.2   4.107 0.021860 *  
## Residuals   54  712.1    13.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Extract key values for interpretation
supp_F <- anova_summary[[1]][1, "F value"]
supp_p <- anova_summary[[1]][1, "Pr(>F)"]
dose_F <- anova_summary[[1]][2, "F value"]
dose_p <- anova_summary[[1]][2, "Pr(>F)"]
inter_F <- anova_summary[[1]][3, "F value"]
inter_p <- anova_summary[[1]][3, "Pr(>F)"]

cat("\n=== INTERPRETATION OF ANOVA ===\n")

## 
## === INTERPRETATION OF ANOVA ===

cat("1. Supplement effect: F =", round(supp_F, 2), 
    ", p =", format(supp_p, scientific = TRUE, digits = 3), "\n")

## 1. Supplement effect: F = 15.57 , p = 2.31e-04

cat("   Conclusion:", ifelse(supp_p < 0.05, 
      "Significant - supplement type affects growth",
      "Not significant - no evidence supplement affects growth"), "\n\n")

##    Conclusion: Significant - supplement type affects growth

cat("2. Dose effect: F =", round(dose_F, 2), 
    ", p =", format(dose_p, scientific = TRUE, digits = 3), "\n")

## 2. Dose effect: F = 92 , p = 4.05e-18

cat("   Conclusion:", ifelse(dose_p < 0.05,
      "Significant - dosage level affects growth",
      "Not significant - no evidence dosage affects growth"), "\n\n")

##    Conclusion: Significant - dosage level affects growth

cat("3. Interaction effect: F =", round(inter_F, 2), 
    ", p =", round(inter_p, 4), "\n")

## 3. Interaction effect: F = 4.11 , p = 0.0219

cat("   Conclusion:", ifelse(inter_p < 0.05,
      "Significant - supplement effect depends on dosage",
      "Not significant - supplement effect consistent across doses"), "\n")

##    Conclusion: Significant - supplement effect depends on dosage

# Post-hoc analysis due to significant interaction
cat("\n=== POST-HOC ANALYSIS (Pairwise Comparisons) ===\n")

## 
## === POST-HOC ANALYSIS (Pairwise Comparisons) ===

cat("Because interaction is significant, analyze each dose separately\n\n")

## Because interaction is significant, analyze each dose separately

pairwise_results <- data.frame()
for(dose_val in levels(ToothGrowth$dose)) {
  subset_data <- ToothGrowth[ToothGrowth$dose == dose_val, ]
  t_test_result <- t.test(len ~ supp, data = subset_data, var.equal = TRUE)
  
  means <- tapply(subset_data$len, subset_data$supp, mean)
  diff_means <- diff(means)  # OJ - VC
  
  pairwise_results <- rbind(pairwise_results,
    data.frame(
      Dose = dose_val,
      Mean_OJ = round(means["OJ"], 2),
      Mean_VC = round(means["VC"], 2),
      Difference = round(diff_means, 2),
      Percent_Increase = round(100 * diff_means / means["VC"], 1),
      CI_95 = paste0("[", round(t_test_result$conf.int[1], 2), 
                    ", ", round(t_test_result$conf.int[2], 2), "]"),
      t_stat = round(t_test_result$statistic, 3),
      p_value = round(t_test_result$p.value, 4),
      Significance = ifelse(t_test_result$p.value < 0.05, "Yes", "No")
    ))
}

print(pairwise_results)

##     Dose Mean_OJ Mean_VC Difference Percent_Increase         CI_95 t_stat
## OJ   0.5   13.23    7.98      -5.25            -65.8  [1.77, 8.73]  3.170
## OJ1    1   22.70   16.77      -5.93            -35.4  [2.84, 9.02]  4.033
## OJ2    2   26.06   26.14       0.08              0.3 [-3.72, 3.56] -0.046
##     p_value Significance
## OJ   0.0053          Yes
## OJ1  0.0008          Yes
## OJ2  0.9637           No

# Visualization of pairwise results
p_posthoc <- ggplot(pairwise_results, aes(x = Dose, y = Difference, 
                                          color = Significance)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = as.numeric(gsub("\\[([^,]+),.*", "\\1", CI_95)),
                    ymax = as.numeric(gsub(".*, ([^]]+)\\]", "\\1", CI_95))),
                width = 0.1, size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_text(aes(label = paste0(round(Difference, 1), " mm\n", 
                               Percent_Increase, "%")), 
            vjust = -0.5, size = 3) +
  scale_color_manual(values = c("Yes" = "steelblue", "No" = "gray")) +
  labs(title = "Mean Differences: OJ minus VC with 95% Confidence Intervals",
       subtitle = "Labels show absolute difference and percentage increase",
       x = "Dosage (mg/day)", y = "Difference in Mean Length (mm)",
       color = "Statistically\nSignificant?") +
  theme_minimal() +
  ylim(min(pairwise_results$Difference) - 2, 
       max(pairwise_results$Difference) + 2)

print(p_posthoc)

# Dose effects within each supplement
cat("\n=== DOSE EFFECTS WITHIN EACH SUPPLEMENT ===\n")

## 
## === DOSE EFFECTS WITHIN EACH SUPPLEMENT ===

library(multcomp)

## Warning: package 'multcomp' was built under R version 4.5.2

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Warning: package 'TH.data' was built under R version 4.5.2

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

# For OJ
oj_data <- ToothGrowth[ToothGrowth$supp == "OJ", ]
oj_anova <- aov(len ~ dose, data = oj_data)
oj_tukey <- TukeyHSD(oj_anova)

# For VC
vc_data <- ToothGrowth[ToothGrowth$supp == "VC", ]
vc_anova <- aov(len ~ dose, data = vc_data)
vc_tukey <- TukeyHSD(vc_anova)

cat("Orange Juice (OJ) dose comparisons:\n")

## Orange Juice (OJ) dose comparisons:

print(oj_tukey)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = len ~ dose, data = oj_data)
## 
## $dose
##        diff        lwr       upr     p adj
## 1-0.5  9.47  5.3096046 13.630395 0.0000158
## 2-0.5 12.83  8.6696046 16.990395 0.0000001
## 2-1    3.36 -0.8003954  7.520395 0.1309258

cat("\nAscorbic Acid (VC) dose comparisons:\n")

## 
## Ascorbic Acid (VC) dose comparisons:

print(vc_tukey)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = len ~ dose, data = vc_data)
## 
## $dose
##        diff       lwr      upr    p adj
## 1-0.5  8.79  4.901765 12.67824 1.75e-05
## 2-0.5 18.16 14.271765 22.04824 0.00e+00
## 2-1    9.37  5.481765 13.25824 6.60e-06

Interpretation of Formal Inference Results:

ANOVA Results: All three effects (supplement, dose, interaction) are statistically significant (p < 0.05).

Pairwise Comparisons:

At 0.5 mg: OJ produces 5.25 mm greater growth than VC (p = 0.006)

At 1.0 mg: OJ produces 5.93 mm greater growth than VC (p = 0.001)
At 2.0 mg: No significant difference (p = 0.963)

Dose Effects Within Supplements:

For OJ: Significant increases from 0.5 to 1.0 mg and from 0.5 to 2.0 mg
For VC: All dose comparisons show significant differences
Why this stage matters: Formal inference provides objective, quantitative answers to research questions with controlled error rates.

Epicycle 5: Assumption Checking and Model Diagnostics

Objective: Verify that model assumptions are satisfied for valid inference.
Rationale: Violated assumptions can lead to incorrect conclusions.

library(car)

## Warning: package 'car' was built under R version 4.5.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.5.2

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

cat("=== MODEL ASSUMPTION CHECKS ===\n\n")

## === MODEL ASSUMPTION CHECKS ===

# 1. Normality of residuals
shapiro_test <- shapiro.test(residuals(full_model))
cat("1. NORMALITY OF RESIDUALS (Shapiro-Wilk test):\n")

## 1. NORMALITY OF RESIDUALS (Shapiro-Wilk test):

cat("   H0: Residuals are normally distributed\n")

##    H0: Residuals are normally distributed

cat("   W =", round(shapiro_test$statistic, 4), 
    ", p =", round(shapiro_test$p.value, 4), "\n")

##    W = 0.985 , p = 0.6694

cat("   Conclusion:", ifelse(shapiro_test$p.value > 0.05,
      "Fail to reject H0 - residuals consistent with normality",
      "Reject H0 - evidence against normality"), "\n\n")

##    Conclusion: Fail to reject H0 - residuals consistent with normality

# 2. Homogeneity of variances
levene_test <- leveneTest(len ~ interaction(supp, dose), data = ToothGrowth)
cat("2. HOMOGENEITY OF VARIANCES (Levene's test):\n")

## 2. HOMOGENEITY OF VARIANCES (Levene's test):

cat("   H0: Variances are equal across groups\n")

##    H0: Variances are equal across groups

cat("   F =", round(levene_test$`F value`[1], 3), 
    ", p =", round(levene_test$`Pr(>F)`[1], 4), "\n")

##    F = 1.709 , p = 0.1484

cat("   Conclusion:", ifelse(levene_test$`Pr(>F)`[1] > 0.05,
      "Fail to reject H0 - variances appear homogeneous",
      "Reject H0 - evidence of heterogeneous variances"), "\n\n")

##    Conclusion: Fail to reject H0 - variances appear homogeneous

# 3. Independence (design-based check)
cat("3. INDEPENDENCE OF OBSERVATIONS:\n")

## 3. INDEPENDENCE OF OBSERVATIONS:

cat("   Design check: Different guinea pigs in each group\n")

##    Design check: Different guinea pigs in each group

cat("   Conclusion: Independence assumption reasonable based on design\n\n")

##    Conclusion: Independence assumption reasonable based on design

# Residuals vs Fitted plot
plot(full_model, which = 1, main = "Residuals vs Fitted Values",
     pch = 19, col = "steelblue", cex = 1.2)
abline(h = 0, col = "red", lty = 2, lwd = 2)

cat("Plot 1: Residuals vs Fitted Values\n")

## Plot 1: Residuals vs Fitted Values

cat("   Interpretation: No clear pattern in residuals suggests linearity assumption is reasonable.\n")

##    Interpretation: No clear pattern in residuals suggests linearity assumption is reasonable.

cat("   The red line at y=0 shows the reference; residuals should be randomly scattered around it.\n\n")

##    The red line at y=0 shows the reference; residuals should be randomly scattered around it.

# Normal Q-Q plot
plot(full_model, which = 2, main = "Normal Q-Q Plot",
     pch = 19, col = "steelblue", cex = 1.2)

cat("\nPlot 2: Normal Q-Q Plot\n")

## 
## Plot 2: Normal Q-Q Plot

cat("   Interpretation: Points generally follow the reference line, suggesting residuals are approximately normally distributed.\n")

##    Interpretation: Points generally follow the reference line, suggesting residuals are approximately normally distributed.

cat("   Minor deviations at extremes are expected and not concerning given the sample size.\n\n")

##    Minor deviations at extremes are expected and not concerning given the sample size.

# Scale-Location plot
plot(full_model, which = 3, main = "Scale-Location Plot (Homoscedasticity Check)",
     pch = 19, col = "steelblue", cex = 1.2)
abline(h = 0, col = "red", lty = 2, lwd = 2)

cat("Plot 3: Scale-Location Plot\n")

## Plot 3: Scale-Location Plot

cat("   Interpretation: Relatively constant spread of points across fitted values suggests homoscedasticity (equal variances).\n")

##    Interpretation: Relatively constant spread of points across fitted values suggests homoscedasticity (equal variances).

cat("   No clear funnel pattern indicates variance is approximately constant.\n\n")

##    No clear funnel pattern indicates variance is approximately constant.

# Residuals vs Leverage plot
plot(full_model, which = 5, main = "Residuals vs Leverage",
     pch = 19, col = "steelblue", cex = 1.2)
# Add Cook's distance contours
invisible(sapply(c(0.5, 1), function(x) {
  curve(sqrt(x * (1-x) / length(full_model$residuals)), 
        from = 0, to = 1, add = TRUE, col = "red", lty = 2)
}))

cat("Plot 4: Residuals vs Leverage\n")

## Plot 4: Residuals vs Leverage

cat("   Interpretation: No points outside the red Cook's distance contours (0.5 and 1.0) suggests no influential outliers.\n")

##    Interpretation: No points outside the red Cook's distance contours (0.5 and 1.0) suggests no influential outliers.

cat("   Most points have low leverage, indicating balanced experimental design.\n\n")

##    Most points have low leverage, indicating balanced experimental design.

cat("=== DIAGNOSTIC PLOTS SUMMARY ===\n\n")

## === DIAGNOSTIC PLOTS SUMMARY ===

cat("1. RESIDUALS VS FITTED:\n")

## 1. RESIDUALS VS FITTED:

cat("   • Assessment: Satisfactory - no clear pattern\n")

##    • Assessment: Satisfactory - no clear pattern

cat("   • Implication: Linearity assumption appears reasonable\n\n")

##    • Implication: Linearity assumption appears reasonable

cat("2. NORMAL Q-Q PLOT:\n")

## 2. NORMAL Q-Q PLOT:

cat("   • Assessment: Satisfactory - points follow reference line\n")

##    • Assessment: Satisfactory - points follow reference line

cat("   • Implication: Normality assumption appears reasonable\n\n")

##    • Implication: Normality assumption appears reasonable

cat("3. SCALE-LOCATION PLOT:\n")

## 3. SCALE-LOCATION PLOT:

cat("   • Assessment: Satisfactory - relatively constant spread\n")

##    • Assessment: Satisfactory - relatively constant spread

cat("   • Implication: Homoscedasticity assumption appears reasonable\n\n")

##    • Implication: Homoscedasticity assumption appears reasonable

cat("4. RESIDUALS VS LEVERAGE:\n")

## 4. RESIDUALS VS LEVERAGE:

cat("   • Assessment: Satisfactory - no influential outliers\n")

##    • Assessment: Satisfactory - no influential outliers

cat("   • Implication: No single observation unduly influences results\n\n")

##    • Implication: No single observation unduly influences results

cat("OVERALL CONCLUSION: All diagnostic plots suggest ANOVA assumptions are reasonably met.\n")

## OVERALL CONCLUSION: All diagnostic plots suggest ANOVA assumptions are reasonably met.

cat("The statistical inference from the two-way ANOVA is therefore valid.\n")

## The statistical inference from the two-way ANOVA is therefore valid.

Assumption Check Results:

Normality: Shapiro-Wilk p = r round(shapiro_test$p.value, 4) > 0.05 ✓
Equal Variances: Levene’s test p = r round(levene_test$Pr(>F)[1], 4) > 0.05 ✓
Independence: Reasonable based on experimental design ✓
Diagnostic Plots: All show satisfactory patterns ✓
Why this stage matters: Statistical conclusions are only valid if model assumptions hold. Systematic checking prevents drawing false conclusions.

Epicycle 6: Synthesis and Interpretation

Objective: Integrate all analyses into coherent answers to research questions.
Rationale: Synthesis transforms statistical results into meaningful scientific conclusions.

cat("=== COMPREHENSIVE SYNTHESIS ===\n\n")

## === COMPREHENSIVE SYNTHESIS ===

cat("RESEARCH QUESTION 1: Does supplement type affect tooth growth?\n")

## RESEARCH QUESTION 1: Does supplement type affect tooth growth?

cat("Answer: YES, but it depends on dosage.\n")

## Answer: YES, but it depends on dosage.

cat("Evidence: Significant supplement main effect (p =", 
    format(supp_p, scientific = TRUE, digits = 3), ")\n")

## Evidence: Significant supplement main effect (p = 2.31e-04 )

cat("          OJ generally produces greater growth than VC\n\n")

##           OJ generally produces greater growth than VC

cat("RESEARCH QUESTION 2: Does dosage level affect tooth growth?\n")

## RESEARCH QUESTION 2: Does dosage level affect tooth growth?

cat("Answer: YES, strongly.\n")

## Answer: YES, strongly.

cat("Evidence: Highly significant dose effect (p < 2.2e-16)\n")

## Evidence: Highly significant dose effect (p < 2.2e-16)

cat("          Growth increases with dosage for both supplements\n\n")

##           Growth increases with dosage for both supplements

cat("RESEARCH QUESTION 3: Is there an interaction?\n")

## RESEARCH QUESTION 3: Is there an interaction?

cat("Answer: YES, supplement effectiveness depends on dosage.\n")

## Answer: YES, supplement effectiveness depends on dosage.

cat("Evidence: Significant interaction (p =", round(inter_p, 4), ")\n")

## Evidence: Significant interaction (p = 0.0219 )

cat("          OJ advantage diminishes at higher doses\n\n")

##           OJ advantage diminishes at higher doses

cat("RESEARCH QUESTION 4: Which combination yields maximum growth?\n")

## RESEARCH QUESTION 4: Which combination yields maximum growth?

cat("Answer: 2.0 mg/day of either supplement.\n")

## Answer: 2.0 mg/day of either supplement.

cat("Evidence: At 2.0 mg, no significant difference between OJ and VC\n")

## Evidence: At 2.0 mg, no significant difference between OJ and VC

cat("          (p =", pairwise_results$p_value[pairwise_results$Dose == 2], ")\n\n")

##           (p = 0.9637 )

cat("=== PRACTICAL RECOMMENDATIONS ===\n")

## === PRACTICAL RECOMMENDATIONS ===

cat("1. FOR MAXIMUM GROWTH:\n")

## 1. FOR MAXIMUM GROWTH:

cat("   • Use 2.0 mg/day of either supplement\n")

##    • Use 2.0 mg/day of either supplement

cat("   • Expected length: ~26.1 mm\n\n")

##    • Expected length: ~26.1 mm

cat("2. FOR COST-EFFECTIVENESS AT LOWER DOSES:\n")

## 2. FOR COST-EFFECTIVENESS AT LOWER DOSES:

cat("   • Use Orange Juice (OJ)\n")

##    • Use Orange Juice (OJ)

cat("   • At 0.5 mg: OJ yields", 
    pairwise_results$Percent_Increase[pairwise_results$Dose == 0.5], 
    "% greater growth than VC\n")

##    • At 0.5 mg: OJ yields -65.8 % greater growth than VC

cat("   • At 1.0 mg: OJ yields", 
    pairwise_results$Percent_Increase[pairwise_results$Dose == 1.0], 
    "% greater growth than VC\n\n")

##    • At 1.0 mg: OJ yields -35.4 % greater growth than VC

cat("3. BIOLOGICAL INTERPRETATION:\n")

## 3. BIOLOGICAL INTERPRETATION:

cat("   • Dose-response: Vitamin C promotes tooth growth in dose-dependent manner\n")

##    • Dose-response: Vitamin C promotes tooth growth in dose-dependent manner

cat("   • Delivery method: Orange Juice may enhance absorption at lower doses\n")

##    • Delivery method: Orange Juice may enhance absorption at lower doses

cat("   • Saturation: At high doses, both delivery methods equally effective\n\n")

##    • Saturation: At high doses, both delivery methods equally effective

cat("=== STATISTICAL VALIDITY ===\n")

## === STATISTICAL VALIDITY ===

cat("1. Assumptions satisfied: Normality (p =", round(shapiro_test$p.value, 4),
    "), Equal variances (p =", round(levene_test$`Pr(>F)`[1], 4), ")\n")

## 1. Assumptions satisfied: Normality (p = 0.6694 ), Equal variances (p = 0.1484 )

cat("2. Appropriate methods: Factorial ANOVA for factorial design\n")

## 2. Appropriate methods: Factorial ANOVA for factorial design

cat("3. Sample size: Adequate (10 per group) for detecting effects\n")

## 3. Sample size: Adequate (10 per group) for detecting effects

cat("4. Multiple comparisons: Addressed through separate analyses by dose\n")

## 4. Multiple comparisons: Addressed through separate analyses by dose

Overall Project Synthesis

Connecting Part 1 and Part 2

cat("=== CONNECTING THEORETICAL AND APPLIED ANALYSES ===\n\n")

## === CONNECTING THEORETICAL AND APPLIED ANALYSES ===

cat("THEORETICAL FOUNDATION (Part 1):\n")

## THEORETICAL FOUNDATION (Part 1):

cat("• Central Limit Theorem: Means become normally distributed as n increases\n")

## • Central Limit Theorem: Means become normally distributed as n increases

cat("• Demonstrated with: Exponential distribution → Normal distribution of means\n")

## • Demonstrated with: Exponential distribution → Normal distribution of means

cat("• Key insight: Can use normal-based inference for means even with non-normal data\n\n")

## • Key insight: Can use normal-based inference for means even with non-normal data

cat("APPLIED ANALYSIS (Part 2):\n")

## APPLIED ANALYSIS (Part 2):

cat("• Relies on CLT: ANOVA assumes normally distributed means\n")

## • Relies on CLT: ANOVA assumes normally distributed means

cat("• Even if individual tooth lengths aren't perfectly normal, means are ≈ normal\n")

## • Even if individual tooth lengths aren't perfectly normal, means are ≈ normal

cat("• Sample size per group (n=10) smaller than Part 1 (n=40), but ANOVA robust\n\n")

## • Sample size per group (n=10) smaller than Part 1 (n=40), but ANOVA robust

cat("METHODOLOGICAL PROGRESSION:\n")

## METHODOLOGICAL PROGRESSION:

cat("1. Simulation demonstrates theoretical principle (CLT)\n")

## 1. Simulation demonstrates theoretical principle (CLT)

cat("2. Real data analysis applies principle (ANOVA relies on CLT)\n")

## 2. Real data analysis applies principle (ANOVA relies on CLT)

cat("3. Both use: Means, variances, confidence intervals, hypothesis tests\n")

## 3. Both use: Means, variances, confidence intervals, hypothesis tests

cat("4. Both require: Assumption checking, appropriate inference methods\n\n")

## 4. Both require: Assumption checking, appropriate inference methods

cat("BROADER IMPLICATIONS:\n")

## BROADER IMPLICATIONS:

cat("• Statistical theory (CLT) enables practical data analysis (ANOVA)\n")

## • Statistical theory (CLT) enables practical data analysis (ANOVA)

cat("• Simulation validates theoretical understanding before applying to real data\n")

## • Simulation validates theoretical understanding before applying to real data

cat("• Iterative analysis (Epicycles) ensures rigorous, reproducible science\n")

## • Iterative analysis (Epicycles) ensures rigorous, reproducible science

Epicycles Summary

This analysis exemplifies the iterative nature of statistical investigation through the Epicycles of Analysis framework:

Epicycle 1: Setting Expectations - Part 1: Theoretical predictions from CLT about mean and variance - Part 2: Research questions about supplement and dose effects - Purpose: Establishes clear targets for analysis

Epicycle 2: Collecting Information - Part 1: Simulation design with controlled parameters - Part 2: Exploratory data analysis of ToothGrowth dataset - Purpose: Gathers relevant, high-quality data

Epicycle 3: Comparing Expectations to Data - Part 1: Quantitative comparison of sample vs theoretical statistics - Part 2: Model specification based on research questions - Purpose: Initial assessment of whether data match expectations

Epicycle 4: Refining Analysis - Part 1: Visualizations, normality tests, confidence intervals - Part 2: Formal ANOVA, pairwise comparisons, interaction analysis - Purpose: Deepens understanding beyond basic comparisons

Epicycle 5: Assumption Checking - Part 1: Independence and distribution checks - Part 2: Normality, homogeneity of variances, diagnostic plots - Purpose: Ensures validity of statistical conclusions

Epicycle 6: Synthesis and Interpretation - Part 1: Connecting simulation results to CLT theory - Part 2: Answering research questions with practical recommendations - Purpose: Transforms statistical results into meaningful knowledge

Key Features of This Epicycles Approach: - Iterative Nature: Insights from later stages can refine earlier stages - Multiple Lines of Evidence: Combines descriptive, visual, and inferential methods - Assumption Awareness: Explicit checking of statistical requirements - Contextual Interpretation: Results interpreted in theoretical and practical contexts - Reproducibility: Clear documentation of all analytical steps

This structured approach ensures rigorous, transparent, and reproducible statistical analysis that moves beyond mere computation to develop genuine understanding of both statistical principles and their practical applications.

Comprehensive Analysis: Investigating the Exponential Distribution, Central Limit Theorem, and ToothGrowth Data

Edwin R. Mutandi

2026-01-30

Analysis Framework: Epicycles of Analysis

Part 1: Exponential Distribution and Central Limit Theorem

Epicycle 1: Setting Initial Expectations

Epicycle 2: Collecting Information (Simulation Design)

Epicycle 3: Comparing Expectations to Data (Descriptive Analysis)

Epicycle 4: Refining Analysis (Visual and Inferential Methods)

Epicycle 5: Synthesis and Interpretation

Part 2: ToothGrowth Data Analysis

Epicycle 1: Setting Initial Expectations

Epicycle 2: Collecting Information (Exploratory Data Analysis)

Epicycle 3: Comparing Expectations to Data (Model Specification)

Epicycle 4: Refining Analysis (Formal Inference)

Epicycle 5: Assumption Checking and Model Diagnostics

Epicycle 6: Synthesis and Interpretation

Overall Project Synthesis

Connecting Part 1 and Part 2

Epicycles Summary