xbe.at

Statistics and maths / Probability and distributions

Correlation Pitfalls Avoiding Misleading Statistics

11 min read, 14 sections

Contents
  1. Understanding Correlation Pitfalls
  2. Anscombe's Quartet Demonstration
  3. Outlier Impact on Correlation
  4. Robust Correlation Measures
  5. Non-linear Relationships Detection
  6. Correlation Matrix Visualization
  7. Bootstrap Correlation Analysis
  8. Correlation vs Causation Demonstration
  9. Partial Correlation Analysis
  10. Identifying Spurious Correlations
  11. Robust Correlation Through Data Transformation
  12. Time Series Correlation Analysis
  13. Machine Learning Impact of Correlation
  14. Additional Resources

Understanding Correlation Pitfalls

Correlation coefficients can be deceptive when analyzing relationships between variables. While a correlation value provides a single numeric measure of association, it can mask underlying patterns, non-linear relationships, and be significantly influenced by outliers in the dataset.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

# Generate base dataset
np.random.seed(42)
x = np.linspace(0, 10, 100)
y = x + np.random.normal(0, 1, 100)

# Calculate correlation
correlation = stats.pearsonr(x, y)[0]

# Plotting
plt.figure(figsize=(10, 6))
plt.scatter(x, y, alpha=0.5)
plt.title(f'Linear Relationship\nPearson Correlation: {correlation:.3f}')
plt.xlabel('X variable')
plt.ylabel('Y variable')
plt.show()

Anscombe's Quartet Demonstration

Anscombe's Quartet is a powerful illustration of why visualization is crucial in data analysis. This famous dataset comprises four distinct x-y pairs that have nearly identical statistical properties but vastly different distributions when plotted.

import pandas as pd

# Create Anscombe's Quartet
anscombe = {
    'x1': [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    'y1': [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68],
    'x2': [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    'y2': [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74],
    'x3': [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
    'y3': [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73],
    'x4': [8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8],
    'y4': [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]
}

df = pd.DataFrame(anscombe)

# Plot all four datasets
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
datasets = [(df['x1'], df['y1']), (df['x2'], df['y2']), 
           (df['x3'], df['y3']), (df['x4'], df['y4'])]

for i, (ax, (x, y)) in enumerate(zip(axes.flat, datasets)):
    ax.scatter(x, y)
    ax.set_title(f'Dataset {i+1}\nr = {stats.pearsonr(x, y)[0]:.3f}')
    
plt.tight_layout()
plt.show()

Outlier Impact on Correlation

Outliers can drastically alter correlation coefficients, leading to misleading interpretations. A single extreme point can artificially inflate or deflate the correlation, making it crucial to identify and investigate outliers before drawing conclusions.

# Generate clean dataset
np.random.seed(42)
x_clean = np.linspace(0, 10, 50)
y_clean = 2 * x_clean + np.random.normal(0, 1, 50)

# Add outliers
x_outliers = np.append(x_clean, [0, 10])
y_outliers = np.append(y_clean, [20, 0])

# Calculate correlations
corr_clean = stats.pearsonr(x_clean, y_clean)[0]
corr_outliers = stats.pearsonr(x_outliers, y_outliers)[0]

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

ax1.scatter(x_clean, y_clean)
ax1.set_title(f'Without Outliers\nr = {corr_clean:.3f}')

ax2.scatter(x_outliers, y_outliers)
ax2.scatter([0, 10], [20, 0], color='red', s=100, label='Outliers')
ax2.set_title(f'With Outliers\nr = {corr_outliers:.3f}')
ax2.legend()

plt.show()

Robust Correlation Measures

Standard Pearson correlation is sensitive to outliers, but alternative robust correlation measures like Spearman's rank correlation and Kendall's tau can provide more reliable assessments of relationships in the presence of outliers or non-linear patterns.

import numpy as np
from scipy.stats import pearsonr, spearmanr, kendalltau

# Generate data with outliers
np.random.seed(42)
x = np.concatenate([np.random.normal(0, 1, 98), [10, -10]])  # 98 normal points + 2 outliers
y = np.concatenate([x[:98] + np.random.normal(0, 0.5, 98), [-8, 8]])  # Related to x with noise + outliers

# Calculate different correlation measures
pearson_corr, _ = pearsonr(x, y)
spearman_corr, _ = spearmanr(x, y)
kendall_corr, _ = kendalltau(x, y)

# Plotting
plt.figure(figsize=(12, 6))
plt.scatter(x, y, alpha=0.5)
plt.title(f'Correlation Measures Comparison\n' 
         f'Pearson: {pearson_corr:.3f}, Spearman: {spearman_corr:.3f}, '
         f'Kendall: {kendall_corr:.3f}')
plt.xlabel('X variable')
plt.ylabel('Y variable')
plt.show()

Non-linear Relationships Detection

Traditional correlation measures may fail to capture non-linear relationships between variables. This example demonstrates how correlation coefficients can be misleading when analyzing data with clear non-linear patterns.

# Generate non-linear relationship data
x = np.linspace(-5, 5, 100)
y = x**2 + np.random.normal(0, 1, 100)

# Calculate correlations
pearson = pearsonr(x, y)[0]
spearman = spearmanr(x, y)[0]

# Create visualization
plt.figure(figsize=(12, 6))
plt.subplot(121)
plt.scatter(x, y, alpha=0.5)
plt.title(f'Quadratic Relationship\nPearson r: {pearson:.3f}')

# Add polynomial fit
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
plt.plot(x, p(x), 'r-', alpha=0.8)

# Add residual plot
plt.subplot(122)
residuals = y - p(x)
plt.scatter(x, residuals, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.title('Residual Plot')

plt.tight_layout()
plt.show()

Correlation Matrix Visualization

Understanding correlations between multiple variables requires effective visualization techniques. A correlation matrix heatmap can reveal patterns and relationships, but it's essential to validate findings through scatter plots and statistical tests.

# Generate correlated data
np.random.seed(42)
n_samples = 1000
n_features = 5

# Create correlation structure
true_corr = np.array([
    [1.0, 0.8, -0.6, 0.1, 0.3],
    [0.8, 1.0, -0.5, 0.2, 0.4],
    [-0.6, -0.5, 1.0, -0.3, -0.2],
    [0.1, 0.2, -0.3, 1.0, 0.6],
    [0.3, 0.4, -0.2, 0.6, 1.0]
])

# Generate multivariate normal data
data = np.random.multivariate_normal(
    mean=np.zeros(n_features),
    cov=true_corr,
    size=n_samples
)

# Create DataFrame
df = pd.DataFrame(data, columns=[f'Var{i+1}' for i in range(n_features)])

# Plot correlation matrix
plt.figure(figsize=(10, 8))
sns.heatmap(df.corr(), annot=True, cmap='RdBu', vmin=-1, vmax=1, center=0)
plt.title('Correlation Matrix Heatmap')
plt.show()

Bootstrap Correlation Analysis

Bootstrapping provides a robust way to assess the stability of correlation coefficients by resampling the data multiple times. This technique helps quantify uncertainty in correlation estimates.

def bootstrap_correlation(x, y, n_bootstrap=1000):
    correlations = []
    n_samples = len(x)
    
    for _ in range(n_bootstrap):
        # Random sampling with replacement
        indices = np.random.randint(0, n_samples, n_samples)
        corr = pearsonr(x[indices], y[indices])[0]
        correlations.append(corr)
    
    return np.array(correlations)

# Generate example data
np.random.seed(42)
x = np.random.normal(0, 1, 100)
y = 0.5 * x + np.random.normal(0, 0.5, 100)

# Perform bootstrap analysis
bootstrap_corrs = bootstrap_correlation(x, y)

# Plot results
plt.figure(figsize=(12, 5))
plt.subplot(121)
plt.scatter(x, y, alpha=0.5)
plt.title(f'Original Data\nCorrelation: {pearsonr(x, y)[0]:.3f}')

plt.subplot(122)
plt.hist(bootstrap_corrs, bins=50, density=True)
plt.axvline(np.mean(bootstrap_corrs), color='r', linestyle='--',
            label=f'Mean: {np.mean(bootstrap_corrs):.3f}')
plt.title('Bootstrap Correlation Distribution')
plt.legend()

plt.tight_layout()
plt.show()

Correlation vs Causation Demonstration

Statistical correlation between variables does not imply causation. This example demonstrates how two completely unrelated variables can show strong correlation due to a common underlying factor or pure coincidence.

# Generate time series data
np.random.seed(42)
time_points = 100
t = np.linspace(0, 10, time_points)

# Create two variables that follow similar seasonal patterns
seasonal_component = np.sin(t) + np.cos(t * 0.5)
var1 = seasonal_component + np.random.normal(0, 0.2, time_points)
var2 = seasonal_component + np.random.normal(0, 0.2, time_points)

# Calculate correlation
correlation = pearsonr(var1, var2)[0]

# Plotting
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# Time series plot
ax1.plot(t, var1, label='Variable 1')
ax1.plot(t, var2, label='Variable 2')
ax1.set_title('Time Series of Two Unrelated Variables')
ax1.legend()

# Scatter plot
ax2.scatter(var1, var2, alpha=0.5)
ax2.set_title(f'Correlation Plot (r = {correlation:.3f})')
ax2.set_xlabel('Variable 1')
ax2.set_ylabel('Variable 2')

plt.tight_layout()
plt.show()

Partial Correlation Analysis

Partial correlation helps identify the true relationship between two variables while controlling for the effects of other variables. This technique is crucial for understanding complex multivariate relationships.

def partial_correlation(data, x, y, controlling):
    # Calculate residuals for x
    model_x = np.polynomial.polynomial.polyfit(controlling, x, 1)
    residuals_x = x - np.polynomial.polynomial.polyval(controlling, model_x)
    
    # Calculate residuals for y
    model_y = np.polynomial.polynomial.polyfit(controlling, y, 1)
    residuals_y = y - np.polynomial.polynomial.polyval(controlling, model_y)
    
    # Calculate partial correlation
    return pearsonr(residuals_x, residuals_y)[0]

# Generate data
np.random.seed(42)
n = 100
z = np.random.normal(0, 1, n)  # Confounding variable
x = 0.7 * z + np.random.normal(0, 0.3, n)
y = 0.6 * z + np.random.normal(0, 0.3, n)

# Calculate correlations
direct_corr = pearsonr(x, y)[0]
partial_corr = partial_correlation(None, x, y, z)

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

ax1.scatter(x, y, alpha=0.5)
ax1.set_title(f'Direct Correlation\nr = {direct_corr:.3f}')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')

# Residual plot after controlling for Z
ax2.scatter(residuals_x, residuals_y, alpha=0.5)
ax2.set_title(f'Partial Correlation (controlling for Z)\nr = {partial_corr:.3f}')
ax2.set_xlabel('X residuals')
ax2.set_ylabel('Y residuals')

plt.tight_layout()
plt.show()

Identifying Spurious Correlations

Spurious correlations can arise from various sources including sampling bias, hidden variables, or coincidental patterns. This example demonstrates how to detect and analyze potentially misleading correlations.

# Generate data with apparent correlation due to time trend
np.random.seed(42)
n_points = 100
time = np.linspace(0, 10, n_points)

# Two increasing trends with different rates
trend1 = 0.5 * time + np.random.normal(0, 0.3, n_points)
trend2 = 0.3 * time + np.random.normal(0, 0.2, n_points)

# Calculate correlations
raw_corr = pearsonr(trend1, trend2)[0]
detrended1 = trend1 - np.polyfit(time, trend1, 1)[0] * time
detrended2 = trend2 - np.polyfit(time, trend2, 1)[0] * time
detrended_corr = pearsonr(detrended1, detrended2)[0]

# Plotting
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Original time series
ax1.plot(time, trend1, label='Trend 1')
ax1.plot(time, trend2, label='Trend 2')
ax1.set_title('Original Time Series')
ax1.legend()

# Original correlation
ax2.scatter(trend1, trend2, alpha=0.5)
ax2.set_title(f'Raw Correlation\nr = {raw_corr:.3f}')

# Detrended time series
ax3.plot(time, detrended1, label='Detrended 1')
ax3.plot(time, detrended2, label='Detrended 2')
ax3.set_title('Detrended Time Series')
ax3.legend()

# Detrended correlation
ax4.scatter(detrended1, detrended2, alpha=0.5)
ax4.set_title(f'Detrended Correlation\nr = {detrended_corr:.3f}')

plt.tight_layout()
plt.show()

Robust Correlation Through Data Transformation

Data transformations can help reveal true relationships by normalizing distributions and reducing the impact of outliers. This example demonstrates how different transformations affect correlation analysis.

def plot_correlations(x, y, transforms):
    fig, axes = plt.subplots(2, len(transforms), figsize=(15, 8))
    
    for i, (name, func) in enumerate(transforms.items()):
        # Apply transformation
        x_trans = func(x)
        y_trans = func(y)
        
        # Calculate correlation
        corr = pearsonr(x_trans, y_trans)[0]
        
        # Scatter plot
        axes[0, i].scatter(x_trans, y_trans, alpha=0.5)
        axes[0, i].set_title(f'{name}\nr = {corr:.3f}')
        
        # QQ plot
        stats.probplot(x_trans, dist="norm", plot=axes[1, i])
        axes[1, i].set_title(f'Q-Q Plot ({name})')
    
    plt.tight_layout()
    return fig

# Generate skewed data with outliers
np.random.seed(42)
x = np.exp(np.random.normal(0, 1, 1000))
y = x + np.random.normal(0, x/2)

# Define transformations
transforms = {
    'Raw': lambda x: x,
    'Log': lambda x: np.log1p(x),
    'Square Root': lambda x: np.sqrt(x),
    'Box-Cox': lambda x: stats.boxcox(x + 1)[0]
}

# Plot results
plot_correlations(x, y, transforms)
plt.show()

Time Series Correlation Analysis

Time series data requires special consideration due to temporal dependencies. This example shows how to analyze correlations in time series using various lag structures and rolling windows.

def rolling_correlation(x, y, window):
    return pd.Series(x).rolling(window).corr(pd.Series(y))

# Generate time series data
np.random.seed(42)
n_points = 500
t = np.linspace(0, 10, n_points)

# Create two related series with varying correlation
series1 = np.sin(t) + np.random.normal(0, 0.2, n_points)
series2 = np.sin(t + np.pi/4) + np.random.normal(0, 0.2, n_points)

# Calculate rolling correlation
window_size = 50
rolling_corr = rolling_correlation(series1, series2, window_size)

# Plot results
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 12))

# Original series
ax1.plot(t, series1, label='Series 1')
ax1.plot(t, series2, label='Series 2')
ax1.set_title('Original Time Series')
ax1.legend()

# Rolling correlation
ax2.plot(t[window_size-1:], rolling_corr[window_size-1:])
ax2.set_title(f'Rolling Correlation (window = {window_size})')
ax2.axhline(y=0, color='r', linestyle='--')

# Lag plot
max_lag = 20
correlations = [pearsonr(series1[lag:], series2[:-lag])[0] 
                if lag > 0 else pearsonr(series1, series2)[0] 
                for lag in range(max_lag)]

ax3.stem(range(max_lag), correlations)
ax3.set_title('Cross-Correlation by Lag')
ax3.set_xlabel('Lag')
ax3.set_ylabel('Correlation')

plt.tight_layout()
plt.show()

Machine Learning Impact of Correlation

Correlation analysis is crucial in feature selection and dimensionality reduction for machine learning. This example demonstrates how correlation affects model performance and feature importance.

from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# Generate correlated features
X, y = make_classification(n_samples=1000, n_features=5, n_informative=3,
                          n_redundant=2, random_state=42)

# Calculate feature correlations
feature_corr = pd.DataFrame(X).corr()

# Train model and get feature importance
rf = RandomForestClassifier(random_state=42)
rf.fit(X, y)
importance = rf.feature_importances_

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Correlation heatmap
sns.heatmap(feature_corr, annot=True, cmap='RdBu', ax=ax1)
ax1.set_title('Feature Correlation Matrix')

# Feature importance
ax2.bar(range(len(importance)), importance)
ax2.set_title('Random Forest Feature Importance')
ax2.set_xlabel('Feature Index')
ax2.set_ylabel('Importance')

plt.tight_layout()
plt.show()

Additional Resources

More in Probability and distributions

  1. Discrete Probability Distributions in Python7 min

    Discrete distributions are probability distributions that describe random variables with a finite or countably infinite set of possible values. They are...

  2. Bernoulli Distribution Probability Modeling Explained10 min

    The Bernoulli distribution models binary outcomes where a random variable X can take only two possible values: success (1) with probability p or failure (0)...

  3. Empirical Rule and Normal Distribution in Python9 min

    The Empirical Rule, also known as the 68-95-99.7 rule, is a fundamental concept in statistics that describes the distribution of data in a normal...

  4. Bernoulli Distribution Explained with Python8 min

    The Bernoulli distribution is a discrete probability distribution for a binary random variable. It models scenarios with two possible outcomes, often...

  5. Essentials of Probability and Statistics5 min

    Probability is the branch of mathematics that deals with the likelihood of events occurring. It forms the foundation for statistical analysis and...

  6. Beginners Guide to Analysis I in Python:2 min