[OK] CHEATSHEET COMPLET - STATSMODELS (ANALYSE STATISTIQUE) PYTHON

Parfait ! J'ai complété le cheatsheet sur statsmodels ! [BRAVO]

**[LISTE] Contenu complet du fichier :**

1. **Introduction** - Vue d'ensemble de statsmodels
2. **Régression linéaire simple** - OLS basique
3. **Régression multiple** - Plusieurs variables
4. **Diagnostics** - Tests et vérifications
5. **Régression logistique** - Variables binaires
6. **GLM** - Modèles linéaires généralisés
7. **ANOVA** - Analyse de variance
8. **ARIMA** - Séries temporelles autorégressives
9. **SARIMA** - Avec saisonnalité
10. **Décomposition** - Tendance, saisonnalité, résidus
11. **Tests de stationnarité** - ADF, KPSS
12. **Tests statistiques** - t-test, chi2, etc.
13. **Tests de proportions** - Comparaisons
14. **Puissance statistique** - Calculs d'échantillon
15. **Régression robuste** - Résistante aux outliers
16. **Régression quantile** - Médiane et quantiles
17. **Modèles mixtes** - Effets aléatoires
18. **WLS** - Régression pondérée
19. **Kaplan-Meier** - Analyse de survie
20. **Modèle de Cox** - Hazard proportionnel
21. **Statistiques descriptives** - Résumés
22. **Corrélations** - Pearson, Spearman, Kendall
23. **Exemple complet régression** - Workflow complet
24. **Exemple complet séries temporelles** - SARIMA end-to-end
25. **Sélection de modèles** - AIC, BIC, ANOVA
26. **Validation croisée** - K-fold
27. **Bonnes pratiques** - Guidelines
28. **Résumé** - Concepts clés
29. **Erreurs courantes** - À éviter
30. **Checklist** - Liste de contrôle
31. **Ressources** - Documentation

Le cheatsheet couvre tous les aspects essentiels de statsmodels pour l'analyse statistique en Python ! [GRAPHIQUE]*

[OK] 1. INTRODUCTION À STATSMODELS
[OK] Bibliothèque pour modélisation statistique et tests
import statsmodels.api as sm
import statsmodels.formula.api as smf
import numpy as np
import pandas as pd

# Statsmodels fournit:
# - Modèles de régression (linéaire, logistique, etc.)
# - Séries temporelles (ARIMA, SARIMA, VAR)
# - Tests statistiques
# - Analyse de variance (ANOVA)
# - Modèles linéaires généralisés (GLM)
# - Modèles mixtes
# - Statistiques descriptives

# Installation
# pip install statsmodels


[OK] 2. RÉGRESSION LINÉAIRE SIMPLE
[OK] Modèle: y = β₀ + β₁x + ε

import numpy as np
import pandas as pd
import statsmodels.api as sm

# Données exemple
np.random.seed(42)
X = np.random.randn(100)
y = 2 + 3 * X + np.random.randn(100) * 0.5

# Méthode 1: Avec sm.OLS (Ordinary Least Squares)
# Ajouter une constante (intercept)
X_with_const = sm.add_constant(X)

# Créer et ajuster le modèle
model = sm.OLS(y, X_with_const)
results = model.fit()

# Résumé complet
print(results.summary())

# Accéder aux résultats
print(f"Coefficients: {results.params}")
print(f"R²: {results.rsquared}")
print(f"R² ajusté: {results.rsquared_adj}")
print(f"P-values: {results.pvalues}")
print(f"Intervalles de confiance:\n{results.conf_int()}")

# Prédictions
X_new = sm.add_constant([0, 1, 2])
predictions = results.predict(X_new)


# Méthode 2: Avec formules (style R)
data = pd.DataFrame({'x': X, 'y': y})
model = smf.ols('y ~ x', data=data)
results = model.fit()
print(results.summary())


[OK] 3. RÉGRESSION LINÉAIRE MULTIPLE
[OK] Modèle: y = β₀ + β₁x₁ + β₂x₂ + ... + βₙxₙ + ε

# Générer des données
np.random.seed(42)
n = 100
X1 = np.random.randn(n)
X2 = np.random.randn(n)
X3 = np.random.randn(n)
y = 2 + 3*X1 - 1.5*X2 + 2*X3 + np.random.randn(n) * 0.5

# Méthode 1: Avec matrices
X = np.column_stack((X1, X2, X3))
X = sm.add_constant(X)

model = sm.OLS(y, X)
results = model.fit()
print(results.summary())

# Méthode 2: Avec formules (plus lisible)
data = pd.DataFrame({
    'y': y,
    'x1': X1,
    'x2': X2,
    'x3': X3
})

model = smf.ols('y ~ x1 + x2 + x3', data=data)
results = model.fit()
print(results.summary())

# Interaction entre variables
model = smf.ols('y ~ x1 + x2 + x1:x2', data=data)  # x1 * x2
results = model.fit()

# Variables catégorielles
data['category'] = np.random.choice(['A', 'B', 'C'], n)
model = smf.ols('y ~ x1 + C(category)', data=data)
results = model.fit()


[OK] 4. DIAGNOSTICS DE RÉGRESSION

# Après avoir ajusté un modèle
results = model.fit()

# 1. Résidus
residuals = results.resid
fitted_values = results.fittedvalues

# 2. Résidus standardisés
standardized_resid = results.resid_pearson

# 3. Influence des observations (Cook's distance)
influence = results.get_influence()
cooks_d = influence.cooks_distance[0]

# 4. Leverage (effet de levier)
leverage = influence.hat_matrix_diag

# 5. Tests sur les résidus
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox

# Test d'hétéroscédasticité (Breusch-Pagan)
lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(results.resid, results.model.exog)
print(f"Test Breusch-Pagan - p-value: {lm_pvalue}")

# Test d'autocorrélation des résidus
lb_test = acorr_ljungbox(results.resid, lags=10)
print(lb_test)

# 6. Test de normalité des résidus
from scipy import stats
statistic, p_value = stats.shapiro(results.resid)
print(f"Test Shapiro-Wilk - p-value: {p_value}")

# 7. VIF (Variance Inflation Factor) - multicolinéarité
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame()
vif_data["Variable"] = data[['x1', 'x2', 'x3']].columns
vif_data["VIF"] = [variance_inflation_factor(X[:, 1:], i) for i in range(3)]
print(vif_data)


[OK] 5. RÉGRESSION LOGISTIQUE
[OK] Pour variables dépendantes binaires (0/1)

# Générer des données binaires
np.random.seed(42)
n = 200
X1 = np.random.randn(n)
X2 = np.random.randn(n)

# Probabilité et variable binaire
z = -1 + 2*X1 + 1.5*X2
prob = 1 / (1 + np.exp(-z))
y = np.random.binomial(1, prob)

# Créer DataFrame
data = pd.DataFrame({'y': y, 'x1': X1, 'x2': X2})

# Modèle logistique
model = smf.logit('y ~ x1 + x2', data=data)
results = model.fit()
print(results.summary())

# Prédictions (probabilités)
predictions_prob = results.predict(data[['x1', 'x2']])

# Prédictions binaires (avec seuil 0.5)
predictions_binary = (predictions_prob > 0.5).astype(int)

# Odds ratios
odds_ratios = np.exp(results.params)
print(f"Odds ratios:\n{odds_ratios}")

# Matrice de confusion
from sklearn.metrics import confusion_matrix, classification_report

cm = confusion_matrix(y, predictions_binary)
print(f"Matrice de confusion:\n{cm}")
print(classification_report(y, predictions_binary))

# Courbe ROC et AUC
from sklearn.metrics import roc_curve, auc
import matplotlib.pyplot as plt

fpr, tpr, thresholds = roc_curve(y, predictions_prob)
roc_auc = auc(fpr, tpr)

plt.figure()
plt.plot(fpr, tpr, label=f'ROC curve (AUC = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], 'k--')
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
plt.legend()
plt.show()


[OK] 6. MODÈLES LINÉAIRES GÉNÉRALISÉS (GLM)

# GLM avec différentes familles de distribution

# 1. Régression de Poisson (pour comptages)
np.random.seed(42)
X = np.random.randn(100, 2)
X = sm.add_constant(X)
y = np.random.poisson(np.exp(0.5 + 0.3*X[:, 1] + 0.2*X[:, 2]))

model = sm.GLM(y, X, family=sm.families.Poisson())
results = model.fit()
print(results.summary())

# 2. Régression Gamma (pour données continues positives)
y_gamma = np.random.gamma(2, 2, 100)
model = sm.GLM(y_gamma, X, family=sm.families.Gamma())
results = model.fit()

# 3. Régression binomiale négative
from statsmodels.discrete.discrete_model import NegativeBinomial

model = NegativeBinomial(y, X)
results = model.fit()

# 4. Familles disponibles
# - Gaussian (normale)
# - Binomial (logistique)
# - Poisson
# - Gamma
# - InverseGaussian
# - NegativeBinomial
# - Tweedie


[OK] 7. ANOVA (ANALYSE DE VARIANCE)

# Générer des données pour ANOVA
np.random.seed(42)
groups = ['A', 'B', 'C']
data = pd.DataFrame({
    'value': np.concatenate([
        np.random.normal(10, 2, 30),  # Groupe A
        np.random.normal(12, 2, 30),  # Groupe B
        np.random.normal(11, 2, 30)   # Groupe C
    ]),
    'group': np.repeat(groups, 30)
})

# ANOVA à un facteur
model = smf.ols('value ~ C(group)', data=data)
results = model.fit()

# Table ANOVA
from statsmodels.stats.anova import anova_lm
anova_table = anova_lm(results, typ=2)
print(anova_table)

# ANOVA à deux facteurs
data['factor2'] = np.tile(['X', 'Y'], 45)
model = smf.ols('value ~ C(group) + C(factor2) + C(group):C(factor2)', data=data)
results = model.fit()
anova_table = anova_lm(results, typ=2)
print(anova_table)

# Tests post-hoc (comparaisons multiples)
from statsmodels.stats.multicomp import pairwise_tukeyhsd

tukey = pairwise_tukeyhsd(data['value'], data['group'], alpha=0.05)
print(tukey)


[OK] 8. SÉRIES TEMPORELLES - ARIMA

import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Générer une série temporelle
np.random.seed(42)
n = 200
dates = pd.date_range('2020-01-01', periods=n, freq='D')
ts = pd.Series(
    np.cumsum(np.random.randn(n)) + 100,
    index=dates
)

# Visualiser ACF et PACF (pour choisir p, q)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ts, lags=20, ax=axes[0])
plot_pacf(ts, lags=20, ax=axes[1])
plt.show()

# Modèle ARIMA(p, d, q)
# p: ordre autorégressif (AR)
# d: degré de différenciation
# q: ordre moyenne mobile (MA)

model = ARIMA(ts, order=(1, 1, 1))
results = model.fit()
print(results.summary())

# Prédictions
forecast = results.forecast(steps=30)
print(forecast)

# Prédictions avec intervalle de confiance
forecast_obj = results.get_forecast(steps=30)
forecast_mean = forecast_obj.predicted_mean
forecast_ci = forecast_obj.conf_int()

# Visualisation
plt.figure(figsize=(12, 6))
plt.plot(ts, label='Données originales')
plt.plot(forecast_mean, label='Prévisions', color='red')
plt.fill_between(forecast_ci.index,
                 forecast_ci.iloc[:, 0],
                 forecast_ci.iloc[:, 1],
                 alpha=0.3)
plt.legend()
plt.show()


[OK] 9. SÉRIES TEMPORELLES - SARIMA

from statsmodels.tsa.statespace.sarimax import SARIMAX

# SARIMA pour séries avec saisonnalité
# SARIMA(p,d,q)(P,D,Q)s
# s: période de saisonnalité (12 pour mensuel, 4 pour trimestriel)

# Générer série avec saisonnalité
n = 120
dates = pd.date_range('2015-01-01', periods=n, freq='M')
seasonal = 10 * np.sin(2 * np.pi * np.arange(n) / 12)
trend = 0.5 * np.arange(n)
noise = np.random.randn(n) * 2
ts = pd.Series(trend + seasonal + noise + 100, index=dates)

# Modèle SARIMA
model = SARIMAX(ts,
                order=(1, 1, 1),           # (p,d,q)
                seasonal_order=(1, 1, 1, 12))  # (P,D,Q,s)
results = model.fit()
print(results.summary())

# Prédictions
forecast = results.forecast(steps=24)

# Diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.show()


[OK] 10. DÉCOMPOSITION DE SÉRIES TEMPORELLES

from statsmodels.tsa.seasonal import seasonal_decompose

# Décomposition additive: Y = Tendance + Saisonnalité + Résidu
decomposition = seasonal_decompose(ts, model='additive', period=12)

# Décomposition multiplicative: Y = Tendance × Saisonnalité × Résidu
# decomposition = seasonal_decompose(ts, model='multiplicative', period=12)

# Visualisation
fig = decomposition.plot()
fig.set_size_inches(12, 8)
plt.show()

# Accéder aux composantes
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid


[OK] 11. TESTS DE STATIONNARITÉ

from statsmodels.tsa.stattools import adfuller, kpss

# Test ADF (Augmented Dickey-Fuller)
# H0: La série a une racine unitaire (non stationnaire)
result = adfuller(ts)
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:', result[4])

if result[1] < 0.05:
    print("Série stationnaire")
else:
    print("Série non stationnaire")

# Test KPSS
# H0: La série est stationnaire
result = kpss(ts)
print('KPSS Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:', result[3])


[OK] 12. TESTS STATISTIQUES

from scipy import stats
import statsmodels.stats.api as sms

# 1. Test t (comparaison de moyennes)
group1 = np.random.normal(10, 2, 30)
group2 = np.random.normal(12, 2, 30)

t_stat, p_value = stats.ttest_ind(group1, group2)
print(f"Test t: t={t_stat:.4f}, p={p_value:.4f}")

# 2. Test t apparié
before = np.random.normal(10, 2, 30)
after = before + np.random.normal(1, 0.5, 30)

t_stat, p_value = stats.ttest_rel(before, after)
print(f"Test t apparié: t={t_stat:.4f}, p={p_value:.4f}")

# 3. Test chi-carré d'indépendance
contingency_table = pd.crosstab(
    pd.Series(['A']*50 + ['B']*50),
    pd.Series(['X']*30 + ['Y']*70)
)
chi2, p_value, dof, expected = stats.chi2_contingency(contingency_table)
print(f"Chi2: {chi2:.4f}, p={p_value:.4f}")

# 4. Test de normalité (Shapiro-Wilk)
data = np.random.normal(0, 1, 100)
stat, p_value = stats.shapiro(data)
print(f"Shapiro-Wilk: W={stat:.4f}, p={p_value:.4f}")

# 5. Test de Jarque-Bera (normalité)
from statsmodels.stats.stattools import jarque_bera
jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(data)
print(f"Jarque-Bera: JB={jb_stat:.4f}, p={jb_pvalue:.4f}")

# 6. Test de Kolmogorov-Smirnov
stat, p_value = stats.kstest(data, 'norm')
print(f"KS: {stat:.4f}, p={p_value:.4f}")

# 7. Test de Levene (égalité des variances)
stat, p_value = stats.levene(group1, group2)
print(f"Levene: {stat:.4f}, p={p_value:.4f}")


[OK] 13. TESTS DE PROPORTIONS

from statsmodels.stats.proportion import proportions_ztest, proportion_confint

# Test de proportion
# H0: proportion = p0
count = 45  # succès
nobs = 100  # observations
p0 = 0.5    # proportion hypothétique

stat, p_value = proportions_ztest(count, nobs, p0)
print(f"Test Z: {stat:.4f}, p={p_value:.4f}")

# Intervalle de confiance pour proportion
ci = proportion_confint(count, nobs, alpha=0.05, method='wilson')
print(f"IC 95%: [{ci[0]:.4f}, {ci[1]:.4f}]")

# Comparaison de deux proportions
count1, nobs1 = 45, 100
count2, nobs2 = 30, 100

stat, p_value = proportions_ztest([count1, count2], [nobs1, nobs2])
print(f"Test Z (2 proportions): {stat:.4f}, p={p_value:.4f}")


[OK] 14. PUISSANCE STATISTIQUE ET TAILLE D'ÉCHANTILLON

from statsmodels.stats.power import TTestIndPower, TTestPower

# Calculer la puissance d'un test
analysis = TTestIndPower()

# Puissance pour un test t à deux échantillons
power = analysis.solve_power(
    effect_size=0.5,  # Cohen's d
    nobs1=50,         # taille échantillon 1
    alpha=0.05,       # seuil de significativité
    ratio=1.0         # ratio des tailles d'échantillon
)
print(f"Puissance: {power:.4f}")

# Taille d'échantillon nécessaire
sample_size = analysis.solve_power(
    effect_size=0.5,
    power=0.80,      # puissance désirée
    alpha=0.05,
    ratio=1.0
)
print(f"Taille échantillon nécessaire: {sample_size:.0f}")

# Effet détectable minimal
effect_size = analysis.solve_power(
    nobs1=50,
    power=0.80,
    alpha=0.05,
    ratio=1.0
)
print(f"Taille d'effet détectable: {effect_size:.4f}")


[OK] 15. RÉGRESSION ROBUSTE

from statsmodels.robust.robust_linear_model import RLM

# Générer données avec outliers
np.random.seed(42)
X = np.random.randn(100)
y = 2 + 3 * X + np.random.randn(100) * 0.5

# Ajouter des outliers
y[95:] = y[95:] + 10

X_with_const = sm.add_constant(X)

# Régression OLS classique (sensible aux outliers)
model_ols = sm.OLS(y, X_with_const)
results_ols = model_ols.fit()

# Régression robuste (résistante aux outliers)
model_rlm = RLM(y, X_with_const, M=sm.robust.norms.HuberT())
results_rlm = model_rlm.fit()

print("OLS coefficients:", results_ols.params)
print("Robust coefficients:", results_rlm.params)

# Poids des observations (faibles pour outliers)
weights = results_rlm.weights
print(f"Poids min: {weights.min():.4f}, max: {weights.max():.4f}")


[OK] 16. RÉGRESSION QUANTILE

from statsmodels.regression.quantile_regression import QuantReg

# Régression sur différents quantiles
np.random.seed(42)
X = np.random.randn(100)
y = 2 + 3 * X + np.random.randn(100) * (1 + 0.5 * X)

X_with_const = sm.add_constant(X)

# Médiane (quantile 0.5)
model_q50 = QuantReg(y, X_with_const)
results_q50 = model_q50.fit(q=0.5)

# Quantiles 0.25 et 0.75
results_q25 = model_q50.fit(q=0.25)
results_q75 = model_q50.fit(q=0.75)

print("Quantile 0.25:", results_q25.params)
print("Quantile 0.50:", results_q50.params)
print("Quantile 0.75:", results_q75.params)

# Visualisation
plt.scatter(X, y, alpha=0.5)
plt.plot(X, results_q25.predict(X_with_const), 'r-', label='Q25')
plt.plot(X, results_q50.predict(X_with_const), 'g-', label='Q50')
plt.plot(X, results_q75.predict(X_with_const), 'b-', label='Q75')
plt.legend()
plt.show()


[OK] 17. MODÈLES MIXTES (EFFETS ALÉATOIRES)

from statsmodels.regression.mixed_linear_model import MixedLM

# Données avec structure hiérarchique
np.random.seed(42)
groups = np.repeat(range(20), 10)  # 20 groupes de 10 observations
X = np.random.randn(200)
group_effects = np.random.randn(20)[groups]  # Effet aléatoire par groupe
y = 2 + 3*X + group_effects + np.random.randn(200) * 0.5

data = pd.DataFrame({
    'y': y,
    'x': X,
    'group': groups
})

# Modèle mixte avec intercept aléatoire
model = MixedLM.from_formula('y ~ x', data=data, groups=data['group'])
results = model.fit()
print(results.summary())

# Effets aléatoires estimés
random_effects = results.random_effects
print(random_effects)


[OK] 18. RÉGRESSION PONDÉRÉE (WLS)

from statsmodels.regression.linear_model import WLS

# Générer données hétéroscédastiques
np.random.seed(42)
X = np.random.randn(100)
# Variance augmente avec X
y = 2 + 3*X + np.random.randn(100) * (1 + np.abs(X))

X_with_const = sm.add_constant(X)

# OLS classique
model_ols = sm.OLS(y, X_with_const)
results_ols = model_ols.fit()

# WLS avec poids inversement proportionnels à la variance
weights = 1 / (1 + np.abs(X))**2
model_wls = WLS(y, X_with_const, weights=weights)
results_wls = model_wls.fit()

print("OLS params:", results_ols.params)
print("WLS params:", results_wls.params)


[OK] 19. ANALYSE DE SURVIE (KAPLAN-MEIER)

from statsmodels.duration.survfunc import SurvfuncRight

# Données de survie
durations = np.array([5, 6, 6, 2, 4, 7, 10, 9, 8, 3])
# 1 = événement observé, 0 = censuré
event_observed = np.array([1, 0, 1, 1, 1, 0, 1, 1, 0, 1])

# Fonction de survie de Kaplan-Meier
surv = SurvfuncRight(durations, event_observed)

# Probabilités de survie
print("Temps:", surv.surv_times)
print("Survie:", surv.surv_prob)

# Intervalle de confiance
print("IC inf:", surv.surv_prob_lb)
print("IC sup:", surv.surv_prob_ub)

# Visualisation
plt.step(surv.surv_times, surv.surv_prob, where='post')
plt.fill_between(surv.surv_times, surv.surv_prob_lb, surv.surv_prob_ub,
                 alpha=0.3, step='post')
plt.xlabel('Temps')
plt.ylabel('Probabilité de survie')
plt.title('Courbe de Kaplan-Meier')
plt.show()


[OK] 20. MODÈLE DE COX (HAZARD PROPORTIONNEL)

from statsmodels.duration.hazard_regression import PHReg

# Données pour modèle de Cox
data = pd.DataFrame({
    'time': [5, 6, 6, 2, 4, 7, 10, 9, 8, 3],
    'event': [1, 0, 1, 1, 1, 0, 1, 1, 0, 1],
    'age': [45, 50, 55, 60, 48, 52, 58, 43, 61, 49],
    'treatment': [0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
})

# Modèle de Cox
model = PHReg.from_formula('time ~ age + treatment', data=data,
                           status='event')
results = model.fit()
print(results.summary())

# Hazard ratios
hazard_ratios = np.exp(results.params)
print("Hazard ratios:", hazard_ratios)


[OK] 21. STATISTIQUES DESCRIPTIVES

import statsmodels.api as sm

# Données
data = np.random.randn(100)

# Statistiques descriptives complètes
desc = sm.stats.DescrStatsW(data)

print(f"Moyenne: {desc.mean:.4f}")
print(f"Écart-type: {desc.std:.4f}")
print(f"Variance: {desc.var:.4f}")
print(f"Médiane: {np.median(data):.4f}")
print(f"Min: {data.min():.4f}")
print(f"Max: {data.max():.4f}")

# Intervalle de confiance pour la moyenne
ci = desc.tconfint_mean(alpha=0.05)
print(f"IC 95% moyenne: [{ci[0]:.4f}, {ci[1]:.4f}]")

# Asymétrie et aplatissement
from scipy.stats import skew, kurtosis
print(f"Skewness: {skew(data):.4f}")
print(f"Kurtosis: {kurtosis(data):.4f}")


[OK] 22. CORRÉLATIONS ET TESTS

from scipy.stats import pearsonr, spearmanr, kendalltau

# Deux variables
x = np.random.randn(100)
y = 2*x + np.random.randn(100)

# Corrélation de Pearson (linéaire)
r, p_value = pearsonr(x, y)
print(f"Pearson r: {r:.4f}, p={p_value:.4f}")

# Corrélation de Spearman (rang)
rho, p_value = spearmanr(x, y)
print(f"Spearman rho: {rho:.4f}, p={p_value:.4f}")

# Tau de Kendall
tau, p_value = kendalltau(x, y)
print(f"Kendall tau: {tau:.4f}, p={p_value:.4f}")

# Corrélation partielle (contrôlant pour z)
from statsmodels.stats.stattools import durbin_watson

z = np.random.randn(100)
# Régression de x sur z
model_xz = sm.OLS(x, sm.add_constant(z))
resid_x = model_xz.fit().resid

# Régression de y sur z
model_yz = sm.OLS(y, sm.add_constant(z))
resid_y = model_yz.fit().resid

# Corrélation des résidus = corrélation partielle
partial_r, p_value = pearsonr(resid_x, resid_y)
print(f"Corrélation partielle: {partial_r:.4f}")


[OK] 23. EXEMPLE COMPLET: ANALYSE DE RÉGRESSION

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns

# Générer dataset
np.random.seed(42)
n = 200

data = pd.DataFrame({
    'age': np.random.randint(20, 70, n),
    'experience': np.random.randint(0, 30, n),
    'education': np.random.choice(['HS', 'Bachelor', 'Master', 'PhD'], n),
    'gender': np.random.choice(['M', 'F'], n)
})

# Variable dépendante (salaire)
data['salary'] = (
    30000 +
    500 * data['age'] +
    1000 * data['experience'] +
    np.where(data['education'] == 'Bachelor', 10000, 0) +
    np.where(data['education'] == 'Master', 20000, 0) +
    np.where(data['education'] == 'PhD', 30000, 0) +
    np.random.randn(n) * 5000
)

# 1. Exploration des données
print(data.describe())
print(data.head())

# 2. Visualisation des relations
sns.pairplot(data, hue='education')
plt.show()

# 3. Matrice de corrélation
numeric_cols = ['age', 'experience', 'salary']
corr_matrix = data[numeric_cols].corr()
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm')
plt.show()

# 4. Modèle de régression multiple
model = smf.ols('salary ~ age + experience + C(education) + C(gender)', 
                data=data)
results = model.fit()

print("\n" + "="*80)
print("RÉSULTATS DE LA RÉGRESSION")
print("="*80)
print(results.summary())

# 5. Diagnostics
# Résidus vs fitted values
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Plot 1: Résidus vs valeurs prédites
axes[0, 0].scatter(results.fittedvalues, results.resid, alpha=0.5)
axes[0, 0].axhline(y=0, color='r', linestyle='--')
axes[0, 0].set_xlabel('Valeurs prédites')
axes[0, 0].set_ylabel('Résidus')
axes[0, 0].set_title('Résidus vs Valeurs prédites')

# Plot 2: Q-Q plot (normalité des résidus)
from scipy import stats
stats.probplot(results.resid, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')

# Plot 3: Distribution des résidus
axes[1, 0].hist(results.resid, bins=30, edgecolor='black')
axes[1, 0].set_xlabel('Résidus')
axes[1, 0].set_ylabel('Fréquence')
axes[1, 0].set_title('Distribution des résidus')

# Plot 4: Résidus standardisés
standardized_resid = results.resid / results.resid.std()
axes[1, 1].scatter(range(len(standardized_resid)), standardized_resid, alpha=0.5)
axes[1, 1].axhline(y=0, color='r', linestyle='--')
axes[1, 1].axhline(y=2, color='r', linestyle=':')
axes[1, 1].axhline(y=-2, color='r', linestyle=':')
axes[1, 1].set_xlabel('Index')
axes[1, 1].set_ylabel('Résidus standardisés')
axes[1, 1].set_title('Résidus standardisés')

plt.tight_layout()
plt.show()

# 6. Tests de diagnostic
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson

# Test d'hétéroscédasticité
lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(results.resid, 
                                                     results.model.exog)
print(f"\nTest Breusch-Pagan (hétéroscédasticité)")
print(f"  LM Statistic: {lm:.4f}")
print(f"  p-value: {lm_pvalue:.4f}")
print(f"  Interprétation: {'Hétéroscédasticité détectée' if lm_pvalue < 0.05 else 'Homoscédasticité OK'}")

# Test de normalité
jb_stat, jb_pvalue, skew, kurt = sm.stats.stattools.jarque_bera(results.resid)
print(f"\nTest Jarque-Bera (normalité)")
print(f"  JB Statistic: {jb_stat:.4f}")
print(f"  p-value: {jb_pvalue:.4f}")
print(f"  Interprétation: {'Non-normalité détectée' if jb_pvalue < 0.05 else 'Normalité OK'}")

# Test d'autocorrélation
dw = durbin_watson(results.resid)
print(f"\nDurbin-Watson (autocorrélation)")
print(f"  Statistic: {dw:.4f}")
print(f"  Interprétation: {'Pas d\'autocorrélation' if 1.5 < dw < 2.5 else 'Autocorrélation possible'}")

# 7. VIF (multicolinéarité)
from statsmodels.stats.outliers_influence import variance_inflation_factor

# Créer matrice X sans la constante et variables catégorielles encodées
X_numeric = data[['age', 'experience']].values
vif_data = pd.DataFrame()
vif_data["Variable"] = ['age', 'experience']
vif_data["VIF"] = [variance_inflation_factor(X_numeric, i) for i in range(2)]
print(f"\n{vif_data}")
print("Interprétation: VIF > 5 indique multicolinéarité")

# 8. Influence des observations
influence = results.get_influence()
cooks_d = influence.cooks_distance[0]
leverage = influence.hat_matrix_diag

# Identifier observations influentes
influential_threshold = 4 / len(data)
influential_obs = np.where(cooks_d > influential_threshold)[0]
print(f"\nObservations influentes (Cook's D > {influential_threshold:.4f}):")
print(f"  Indices: {influential_obs[:10]}")  # Afficher les 10 premières

# 9. Prédictions
# Créer nouvelles observations
new_data = pd.DataFrame({
    'age': [30, 45, 50],
    'experience': [5, 15, 20],
    'education': ['Bachelor', 'Master', 'PhD'],
    'gender': ['M', 'F', 'M']
})

predictions = results.predict(new_data)
pred_ci = results.get_prediction(new_data).conf_int()

print(f"\nPrédictions pour nouvelles observations:")
for i, row in new_data.iterrows():
    print(f"  {row['age']} ans, {row['experience']} ans exp., {row['education']}: "
          f"${predictions.iloc[i]:,.0f} (IC: ${pred_ci.iloc[i, 0]:,.0f} - ${pred_ci.iloc[i, 1]:,.0f})")


[OK] 24. EXEMPLE COMPLET: SÉRIES TEMPORELLES

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose

# 1. Générer série temporelle avec tendance et saisonnalité
np.random.seed(42)
n = 365 * 3  # 3 ans de données quotidiennes
dates = pd.date_range('2020-01-01', periods=n, freq='D')

# Composantes
trend = 0.5 * np.arange(n)
seasonal = 20 * np.sin(2 * np.pi * np.arange(n) / 365)
noise = np.random.randn(n) * 5
ts = pd.Series(trend + seasonal + noise + 100, index=dates, name='value')

print("Série temporelle créée:")
print(ts.head())
print(f"Période: {ts.index.min()} à {ts.index.max()}")

# 2. Visualisation
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(ts)
axes[0].set_title('Série temporelle originale')
axes[0].set_xlabel('Date')
axes[0].set_ylabel('Valeur')
axes[0].grid(True, alpha=0.3)

# Moyenne mobile
ts_ma = ts.rolling(window=30).mean()
axes[1].plot(ts, alpha=0.5, label='Original')
axes[1].plot(ts_ma, color='red', label='Moyenne mobile (30j)')
axes[1].set_title('Avec moyenne mobile')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 3. Test de stationnarité
print("\n" + "="*80)
print("TEST DE STATIONNARITÉ")
print("="*80)

result = adfuller(ts)
print(f"ADF Statistic: {result[0]:.4f}")
print(f"p-value: {result[1]:.4f}")
print(f"Valeurs critiques:")
for key, value in result[4].items():
    print(f"  {key}: {value:.4f}")

if result[1] < 0.05:
    print("Conclusion: Série STATIONNAIRE")
else:
    print("Conclusion: Série NON STATIONNAIRE - différenciation nécessaire")

# 4. Décomposition
print("\n" + "="*80)
print("DÉCOMPOSITION DE LA SÉRIE")
print("="*80)

decomposition = seasonal_decompose(ts, model='additive', period=365)

fig, axes = plt.subplots(4, 1, figsize=(12, 10))
decomposition.observed.plot(ax=axes[0], title='Original')
decomposition.trend.plot(ax=axes[1], title='Tendance')
decomposition.seasonal.plot(ax=axes[2], title='Saisonnalité')
decomposition.resid.plot(ax=axes[3], title='Résidus')

for ax in axes:
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 5. ACF et PACF pour déterminer paramètres ARIMA
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

plot_acf(ts, lags=40, ax=axes[0])
axes[0].set_title('Autocorrélation (ACF)')

plot_pacf(ts, lags=40, ax=axes[1])
axes[1].set_title('Autocorrélation partielle (PACF)')

plt.tight_layout()
plt.show()

# 6. Différenciation pour rendre stationnaire
ts_diff = ts.diff().dropna()

# Test sur série différenciée
result_diff = adfuller(ts_diff)
print(f"\nAprès différenciation:")
print(f"  ADF p-value: {result_diff[1]:.4f}")
print(f"  Stationnaire: {result_diff[1] < 0.05}")

# 7. Split train/test
train_size = int(len(ts) * 0.8)
train, test = ts[:train_size], ts[train_size:]

print(f"\nTrain: {len(train)} observations ({train.index.min()} à {train.index.max()})")
print(f"Test:  {len(test)} observations ({test.index.min()} à {test.index.max()})")

# 8. Modèle ARIMA
print("\n" + "="*80)
print("MODÈLE ARIMA")
print("="*80)

model_arima = ARIMA(train, order=(1, 1, 1))
results_arima = model_arima.fit()
print(results_arima.summary())

# Prédictions sur test set
forecast_arima = results_arima.forecast(steps=len(test))

# Métriques
from sklearn.metrics import mean_squared_error, mean_absolute_error

mse = mean_squared_error(test, forecast_arima)
mae = mean_absolute_error(test, forecast_arima)
rmse = np.sqrt(mse)

print(f"\nPerformance ARIMA:")
print(f"  RMSE: {rmse:.2f}")
print(f"  MAE:  {mae:.2f}")

# 9. Modèle SARIMA (avec saisonnalité)
print("\n" + "="*80)
print("MODÈLE SARIMA")
print("="*80)

model_sarima = SARIMAX(train, 
                       order=(1, 1, 1),
                       seasonal_order=(1, 1, 1, 365))
results_sarima = model_sarima.fit(disp=False)
print(results_sarima.summary())

# Prédictions SARIMA
forecast_sarima = results_sarima.forecast(steps=len(test))

mse_sarima = mean_squared_error(test, forecast_sarima)
mae_sarima = mean_absolute_error(test, forecast_sarima)
rmse_sarima = np.sqrt(mse_sarima)

print(f"\nPerformance SARIMA:")
print(f"  RMSE: {rmse_sarima:.2f}")
print(f"  MAE:  {mae_sarima:.2f}")

# 10. Visualisation des prédictions
fig, ax = plt.subplots(figsize=(14, 6))

ax.plot(train.index, train, label='Train', color='blue')
ax.plot(test.index, test, label='Test', color='green')
ax.plot(test.index, forecast_arima, label='ARIMA', color='red', linestyle='--')
ax.plot(test.index, forecast_sarima, label='SARIMA', color='orange', linestyle='--')

ax.set_title('Prédictions vs Valeurs réelles')
ax.set_xlabel('Date')
ax.set_ylabel('Valeur')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 11. Diagnostics du modèle
results_sarima.plot_diagnostics(figsize=(14, 8))
plt.show()

# 12. Prédictions futures
future_steps = 60
forecast_future = results_sarima.forecast(steps=future_steps)
forecast_obj = results_sarima.get_forecast(steps=future_steps)
forecast_ci = forecast_obj.conf_int()

# Visualisation des prédictions futures
fig, ax = plt.subplots(figsize=(14, 6))

# Données historiques (derniers 180 jours)
ax.plot(ts[-180:], label='Historique', color='blue')

# Prédictions futures avec IC
future_dates = pd.date_range(ts.index[-1] + pd.Timedelta(days=1), 
                             periods=future_steps, freq='D')
ax.plot(future_dates, forecast_future, label='Prévisions', color='red')
ax.fill_between(future_dates,
                forecast_ci.iloc[:, 0],
                forecast_ci.iloc[:, 1],
                alpha=0.3, color='red', label='IC 95%')

ax.set_title('Prévisions futures (60 jours)')
ax.set_xlabel('Date')
ax.set_ylabel('Valeur')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nPremières prédictions futures:")
for i in range(min(10, len(forecast_future))):
    date = future_dates[i]
    pred = forecast_future.iloc[i]
    ci_low = forecast_ci.iloc[i, 0]
    ci_high = forecast_ci.iloc[i, 1]
    print(f"  {date.date()}: {pred:.2f} (IC: {ci_low:.2f} - {ci_high:.2f})")


[OK] 25. SÉLECTION DE MODÈLES

from statsmodels.stats.anova import anova_lm

# Comparer plusieurs modèles
np.random.seed(42)
n = 100
data = pd.DataFrame({
    'y': np.random.randn(n) + 10,
    'x1': np.random.randn(n),
    'x2': np.random.randn(n),
    'x3': np.random.randn(n)
})

# Modèles de complexité croissante
model1 = smf.ols('y ~ x1', data=data).fit()
model2 = smf.ols('y ~ x1 + x2', data=data).fit()
model3 = smf.ols('y ~ x1 + x2 + x3', data=data).fit()

# Comparaison avec ANOVA
anova_results = anova_lm(model1, model2, model3)
print("\nComparaison des modèles (ANOVA):")
print(anova_results)

# Critères d'information
print("\nCritères d'information:")
print(f"Modèle 1 - AIC: {model1.aic:.2f}, BIC: {model1.bic:.2f}")
print(f"Modèle 2 - AIC: {model2.aic:.2f}, BIC: {model2.bic:.2f}")
print(f"Modèle 3 - AIC: {model3.aic:.2f}, BIC: {model3.bic:.2f}")
print("\nNote: Plus faible AIC/BIC = meilleur modèle")

# R² et R² ajusté
print("\nR² et R² ajusté:")
print(f"Modèle 1 - R²: {model1.rsquared:.4f}, R² adj: {model1.rsquared_adj:.4f}")
print(f"Modèle 2 - R²: {model2.rsquared:.4f}, R² adj: {model2.rsquared_adj:.4f}")
print(f"Modèle 3 - R²: {model3.rsquared:.4f}, R² adj: {model3.rsquared_adj:.4f}")


[OK] 26. VALIDATION CROISÉE

from sklearn.model_selection import KFold

# K-Fold validation pour régression
np.random.seed(42)
n = 200
X = np.random.randn(n, 2)
y = 2 + 3*X[:, 0] - 1.5*X[:, 1] + np.random.randn(n)

X_with_const = sm.add_constant(X)

# 5-fold cross-validation
kf = KFold(n_splits=5, shuffle=True, random_state=42)
r2_scores = []
rmse_scores = []

for train_idx, test_idx in kf.split(X):
    # Train
    X_train = X_with_const[train_idx]
    y_train = y[train_idx]
    
    # Test
    X_test = X_with_const[test_idx]
    y_test = y[test_idx]
    
    # Modèle
    model = sm.OLS(y_train, X_train)
    results = model.fit()
    
    # Prédictions
    y_pred = results.predict(X_test)
    
    # Métriques
    r2 = 1 - np.sum((y_test - y_pred)**2) / np.sum((y_test - y_test.mean())**2)
    rmse = np.sqrt(np.mean((y_test - y_pred)**2))
    
    r2_scores.append(r2)
    rmse_scores.append(rmse)

print(f"\n5-Fold Cross-Validation:")
print(f"  R² moyen: {np.mean(r2_scores):.4f} (±{np.std(r2_scores):.4f})")
print(f"  RMSE moyen: {np.mean(rmse_scores):.4f} (±{np.std(rmse_scores):.4f})")


[OK] 27. BONNES PRATIQUES

"""
[OK] 1. TOUJOURS VÉRIFIER LES HYPOTHÈSES
    - Linéarité (scatter plots)
    - Indépendance (Durbin-Watson)
    - Homoscédasticité (Breusch-Pagan)
    - Normalité des résidus (Jarque-Bera, Q-Q plot)

[OK] 2. DIAGNOSTICS DE RÉGRESSION
    - Examiner les résidus
    - Identifier les outliers (Cook's distance)
    - VIF pour multicolinéarité (VIF > 5 = problème)
    - Leverage points

[OK] 3. SÉLECTION DE VARIABLES
    - Ne pas inclure trop de variables (overfitting)
    - Utiliser R² ajusté, AIC, BIC
    - Backward/forward selection
    - Tests de significativité (p-values)

[OK] 4. VALIDATION
    - Split train/test (80/20 ou 70/30)
    - Cross-validation pour petits datasets
    - Out-of-sample predictions

[OK] 5. INTERPRÉTATION
    - Toujours regarder summary()
    - Intervalles de confiance des coefficients
    - Taille d'effet vs significativité statistique
    - Contexte métier

[OK] 6. SÉRIES TEMPORELLES
    - Vérifier stationnarité (ADF test)
    - Différencier si nécessaire
    - ACF/PACF pour choisir paramètres
    - Diagnostics (Ljung-Box)

[OK] 7. REPORTING
    - Inclure summary complet
    - Graphiques de diagnostics
    - Métriques de performance
    - Interprétation en langage simple
"""


[OK] 28. RÉSUMÉ DES CONCEPTS CLÉS

"""
STATSMODELS = Modélisation statistique en Python

RÉGRESSION:
    OLS: Moindres carrés ordinaires
    WLS: Moindres carrés pondérés
    RLM: Régression robuste
    Logit: Régression logistique
    GLM: Modèles linéaires généralisés
    
SÉRIES TEMPORELLES:
    ARIMA(p,d,q): Autorégressif intégré moyenne mobile
    SARIMA: ARIMA avec saisonnalité
    Décomposition: Tendance + Saisonnalité + Résidu
    Tests: ADF (stationnarité)
    
TESTS STATISTIQUES:
    t-test: Comparaison de moyennes
    Chi2: Tests d'indépendance
    ANOVA: Analyse de variance
    Normalité: Shapiro-Wilk, Jarque-Bera
    Corrélations: Pearson, Spearman, Kendall
    
DIAGNOSTICS:
    Résidus: Plots, tests d'autocorrélation
    Hétéroscédasticité: Breusch-Pagan
    Multicolinéarité: VIF
    Influence: Cook's distance, leverage
    
SÉLECTION:
    AIC/BIC: Critères d'information
    R² ajusté: Ajusté pour nombre de variables
    ANOVA: Comparaison de modèles
    Cross-validation: Performance out-of-sample
    
FORMULES: