ANOVA (Analysis of Variance): A Comprehensive Guide with Python Implementation

    
      import numpy as np
      import pandas as pd
      import matplotlib.pyplot as plt
      import seaborn as sns
      from scipy import stats
      import statsmodels.api as sm
      from statsmodels.formula.api import ols

      # Set a random seed for reproducibility
      np.random.seed(42)

      # Generate synthetic data for three treatment groups
      group_a = np.random.normal(loc=5, scale=1.5, size=30)
      group_b = np.random.normal(loc=6, scale=1.5, size=30)
      group_c = np.random.normal(loc=7, scale=1.5, size=30)

      # Create a DataFrame
      data = pd.DataFrame({
          'value': np.concatenate([group_a, group_b, group_c]),
          'group': np.repeat(['A', 'B', 'C'], repeats=30)
      })

      # Display the first few rows
      print(data.head())
    
  

    
      # Create a box plot
      plt.figure(figsize=(10, 6))
      sns.boxplot(x='group', y='value', data=data)
      plt.title('Box Plot of Values by Group')
      plt.xlabel('Group')
      plt.ylabel('Value')
      plt.grid(True, linestyle='--', alpha=0.7)
      plt.show()

      # Create a violin plot for more detailed distribution visualization
      plt.figure(figsize=(10, 6))
      sns.violinplot(x='group', y='value', data=data)
      plt.title('Violin Plot of Values by Group')
      plt.xlabel('Group')
      plt.ylabel('Value')
      plt.grid(True, linestyle='--', alpha=0.7)
      plt.show()
    
  

    
      # 1. Check for normality within each group using Shapiro-Wilk test
      for group_name, group_data in data.groupby('group'):
          stat, p_value = stats.shapiro(group_data['value'])
          print(f'Group {group_name} Shapiro-Wilk test: p-value = {p_value:.4f}')

      # 2. Check for homogeneity of variance using Levene's test
      stat, p_value = stats.levene(group_a, group_b, group_c)
      print(f'Levene\'s test for homogeneity of variance: p-value = {p_value:.4f}')

      # 3. Create Q-Q plots to visually check normality
      fig, axes = plt.subplots(1, 3, figsize=(15, 5))
      for i, (group_name, group_data) in enumerate(data.groupby('group')):
          stats.probplot(group_data['value'], plot=axes[i])
          axes[i].set_title(f'Q-Q Plot for Group {group_name}')
      plt.tight_layout()
      plt.show()
    
  

    
      # Perform one-way ANOVA using SciPy
      f_stat, p_value = stats.f_oneway(group_a, group_b, group_c)
      print(f'F-statistic: {f_stat:.4f}')
      print(f'p-value: {p_value:.4f}')

      # Interpret the results
      alpha = 0.05
      if p_value < alpha:
          print("Reject the null hypothesis: There are significant differences between groups.")
      else:
          print("Fail to reject the null hypothesis: No significant differences between groups.")
    
  

    
      # Perform ANOVA using statsmodels
      model = ols('value ~ C(group)', data=data).fit()
      anova_table = sm.stats.anova_lm(model, typ=2)
      print("\nANOVA Table:")
      print(anova_table)
    
  

    
      # Perform Tukey's HSD test
      from statsmodels.stats.multicomp import pairwise_tukeyhsd

      # Tukey's HSD test
      tukey_results = pairwise_tukeyhsd(data['value'], data['group'], alpha=0.05)
      print("\nTukey's HSD Test:")
      print(tukey_results)

      # Visualize the post-hoc results
      plt.figure(figsize=(10, 6))
      tukey_results.plot_simultaneous()
      plt.title("Tukey's HSD Test for Multiple Comparisons")
      plt.grid(True, linestyle='--', alpha=0.7)
      plt.show()
    
  

    
      # Calculate eta squared (effect size)
      def eta_squared(anova_table):
          sum_sq = anova_table['sum_sq']
          eta_sq = sum_sq.iloc[0] / sum_sq.sum()
          return eta_sq

      # Calculate effect size
      effect_size = eta_squared(anova_table)
      print(f"\nEta squared (effect size): {effect_size:.4f}")

      # Interpret the effect size
      if effect_size < 0.01:
          interpretation = "very small"
      elif effect_size < 0.06:
          interpretation = "small"
      elif effect_size < 0.14:
          interpretation = "medium"
      else:
          interpretation = "large"

      print(f"This indicates a {interpretation} effect size.")
    
  

    
      # Set seed for reproducibility
      np.random.seed(123)

      # Generate data for three investment strategies with realistic returns
      # Value investing: Mean return 8%, higher volatility
      value_returns = np.random.normal(loc=8, scale=15, size=30)

      # Growth investing: Mean return 11%, even higher volatility
      growth_returns = np.random.normal(loc=11, scale=18, size=30)

      # Index investing: Mean return 7%, lower volatility
      index_returns = np.random.normal(loc=7, scale=10, size=30)

      # Create a DataFrame
      investment_data = pd.DataFrame({
          'annual_return': np.concatenate([value_returns, growth_returns, index_returns]),
          'strategy': np.repeat(['Value', 'Growth', 'Index'], repeats=30)
      })

      # Display summary statistics
      print("\nAnnual Returns Summary by Investment Strategy:")
      print(investment_data.groupby('strategy')['annual_return'].describe())

      # Visualize the data
      plt.figure(figsize=(10, 6))
      sns.boxplot(x='strategy', y='annual_return', data=investment_data)
      plt.title('Annual Returns by Investment Strategy')
      plt.xlabel('Investment Strategy')
      plt.ylabel('Annual Return (%)')
      plt.axhline(y=0, color='r', linestyle='-', alpha=0.3)  # Zero return reference line
      plt.grid(True, linestyle='--', alpha=0.7)
      plt.show()

      # Perform one-way ANOVA
      f_stat, p_value = stats.f_oneway(value_returns, growth_returns, index_returns)
      print(f'\nANOVA Results:')
      print(f'F-statistic: {f_stat:.4f}')
      print(f'p-value: {p_value:.8f}')

      # Check if ANOVA is significant
      alpha = 0.05
      if p_value < alpha:
          print("Result: There are significant differences in returns between investment strategies.")

          # Perform post-hoc test to identify which strategies differ
          tukey_results = pairwise_tukeyhsd(investment_data['annual_return'], investment_data['strategy'], alpha=0.05)
          print("\nTukey's HSD Test:")
          print(tukey_results)

          # Calculate risk-adjusted returns (Sharpe ratio simplified)
          # Assuming risk-free rate of 2%
          risk_free = 2
          investment_data['excess_return'] = investment_data['annual_return'] - risk_free

          # Calculate Sharpe ratios for each strategy
          sharpe_ratios = {}
          for strategy, data in investment_data.groupby('strategy'):
              mean_excess_return = data['excess_return'].mean()
              std_dev = data['excess_return'].std()
              sharpe_ratio = mean_excess_return / std_dev if std_dev > 0 else 0
              sharpe_ratios[strategy] = sharpe_ratio

          print("\nRisk-Adjusted Performance (Sharpe Ratio):")
          for strategy, ratio in sharpe_ratios.items():
              print(f"{strategy}: {ratio:.4f}")
      else:
          print("Result: No significant differences in returns between investment strategies.")