Évaluation de l'Hypothèse des Odds Proportionnels dans la Régression Logistique Ordinale : Méthodes et Implémentation en Python

Exploration du Modèle des Odds Proportionnels pour la Régression Logistique Ordinale Le modèle des odds proportionnels, également connu sous le nom de modèle des probabilités cumulatives, a été introduit pour la première fois par McCullagh en 1980. Ce modèle généralise la régression logistique binaire aux situations où la variable dépendante est ordinales, c'est-à-dire composée de valeurs catégorielles ordonnées. Le modèle des odds proportionnels repose sur plusieurs hypothèses, notamment l'indépendance des observations, la linéarité des log-odds, l'absence de multicollinéarité parmi les prédicteurs et l'hypothèse des odds proportionnels. Cette dernière hypothèse stipule que les coefficients de régression sont constants pour toutes les seuils de la variable dépendante ordinate. Pour garantir la validité et l'interprétabilité du modèle, il est crucial de vérifier cette hypothèse. L'article est organisé en quatre sections principales : 1. Introduction au modèle des odds proportionnels 2. Évaluation de l'hypothèse des odds proportionnels via le test du rapport de vraisemblance 3. Évaluation de l'hypothèse des odds proportionnels via les ajustements séparés 4. Exemples d'implémentation en Python 1. Introduction au Modèle des Odds Proportionnels La structure des données suppose que nous avons N observations indépendantes, chacune représentée par un vecteur de p variables explicatives ( \mathbf{X}i = (X{i1}, X_{i2}, \ldots, X_{ip}) ) et une variable dépendante ou réponse Y qui prend des valeurs ordinale de 1 à K. Le modèle des odds proportionnels modélise les probabilités de distribution cumulées de la variable réponse Y, définies par ( \gamma_j = P(Y \leq j | \mathbf{X}_i) ) pour ( j = 1, 2, \ldots, K − 1 ), comme fonctions des variables explicatives ( \mathbf{X}_i ). La formule du modèle est : [ \text{logit}(\gamma_j) = \theta_j - \mathbf{X}_i^\top \beta \quad \text{(1)} ] où ( \theta_j ) sont les intercepts pour chaque catégorie j, et ( \beta ) est le vecteur des coefficients de régression, identique pour toutes les catégories. On observe une tendance monotone dans les coefficients ( \theta_j ) à travers les catégories de la variable réponse Y. Ce modèle peut également être derivé en supposant l'existence d'une variable latente continue ( Y^ ), qui suit un modèle de régression linéaire avec une moyenne conditionnelle ( \eta = \mathbf{X}_i^\top \beta ). La variable observée Y est dérivée de ( Y^ ) à travers des seuils ( \theta_j ) comme suit : [ Y^* = \mathbf{X}_i^\top \beta + \epsilon \quad \text{(2)} ] où ( \epsilon ) est un terme d'erreur, généralement supposé suivre une distribution logistique standard dans le modèle des odds proportionnels. 2. Évaluation de l'Hypothèse des Odds Proportionnels : Le Test du Rapport de Vraisemblance Pour évaluer l'hypothèse des odds proportionnels, Brant (1990) propose l'utilisation du test du rapport de vraisemblance. Cela commence par ajuster un modèle moins restrictif où les coefficients de régression sont autorisés à varier selon les catégories. Ce modèle est défini par : [ \text{logit}(\gamma_j) = \theta_j - \mathbf{X}_i^\top \beta_j \quad \text{pour} \ j = 1, \ldots, K - 1 \quad \text{(4)} ] où ( \beta_j ) est le vecteur des coefficients de régression pour chaque catégorie j. Ensuite, on effectue un test du rapport de vraisemblance pour comparer ce modèle non-proportionnel (ou saturé) avec le modèle proportionnel (ou réduit) où l'hypothèse des odds proportionnels tient. Le test statistique du rapport de vraisemblance est : [ \lambda = -2 \left( \ln \mathcal{L}_0(\hat{\theta}_0) - \ln \mathcal{L}_1(\hat{\theta}) \right) ] où ( \mathcal{L}_0 ) et ( \mathcal{L}_1 ) sont les fonctions de vraisemblance des modèles contraint et non-contraint, respectivement, et ( \hat{\theta}_0 ) et ( \hat{\theta} ) sont les estimations maximales de vraisemblance (EMV) sous ces modèles. Le nombre total de paramètres dans le modèle non-proportionnel est ( (K - 1)(p + 1) ), tandis que dans le modèle proportionnel, il est ( (K - 1) + p ). Ainsi, le test du rapport de vraisemblance suit une distribution khi-deux avec ( (K - 2)p ) degrés de liberté. En appliquant ce test au dataset "Qualité des vins rouges," on obtient une statistique du rapport de vraisemblance de ( 53.207 ) et une p-value de ( 1.066 \times 10^{-9} ). Cette p-value indique que l'hypothèse des odds proportionnels est violée au niveau de significativité de 5%. 3. Évaluation de l'Hypothèse des Odds Proportionnels : L'Approche des Ajustements Séparés Une autre méthode pour évaluer l'hypothèse des odds proportionnels est l'approche des ajustements séparés, basée sur la distance de Mahalanobis. On ajuste une série de ( K - 1 ) modèles de régression logistique binaire, chacun pour un seuil j de la variable réponse Y. Pour chaque seuil j, on définit une variable binaire ( Z_j ) qui vaut 1 si l'observation dépasse le seuil j, et 0 sinon. Les modèles sont définis par : [ \text{logit}(\pi_j) = \theta_j - \mathbf{X}_i^\top \beta_j \quad \text{(10)} ] où ( \pi_j = P(Z_j = 1 | \mathbf{X}_i) = 1 - \gamma_j ). L'hypothèse que l'on teste est : [ H_0 : \beta_1 = \beta_2 = \cdots = \beta_{K-1} \quad \text{(11)} ] Pour évaluer cette hypothèse, on construit une matrice ( \mathbf{D} ) qui capture les différences entre les coefficients ( \beta_j ). La statistique de Wald ( X^2 ) est ensuite calculée comme suit : [ X^2 = (\mathbf{D} \beta)^T (\mathbf{D} \hat{\mathbf{V}} \mathbf{D}^T)^{-1} (\mathbf{D} \beta) ] où ( \beta ) est le vecteur des coefficients estimés et ( \hat{\mathbf{V}} ) est la matrice de variance-covariance globale. En appliquant cette méthode au dataset "Qualité des vins rouges," on trouve une statistique de Wald de ( 41.880 ) et une p-value de ( 1.232 \times 10^{-7} ). Ces résultats confirment également la violation de l'hypothèse des odds proportionnels au niveau de significativité de 5%. 4. Exemples d'Implémentation en Python Les exemples suivants illustrent comment implementer les méthodes d'évaluation de l'hypothèse des odds proportionnels en Python, en utilisant le package statsmodels pour les régressions et les tests statistiques. Charger les données et prétraitement initial : ```python import pandas as pd data = pd.read_csv("winequality-red.csv", sep=";") data['quality'] = data['quality'].replace({3: 4, 8: 7}) print("Nombre d'observations:", data.shape[0]) ``` Gérer les valeurs aberrantes des variables explicatives : ```python def remove_outliers_iqr(df, column): Q1 = df[column].quantile(0.25) Q3 = df[column].quantile(0.75) IQR = Q3 - Q1 lower_bound = Q1 - 1.5 * IQR upper_bound = Q3 + 1.5 * IQR return df[(df[column] >= lower_bound) & (df[column] <= upper_bound)] explanatory_vars = ['volatile acidity', 'free sulfur dioxide', 'total sulfur dioxide'] for col in explanatory_vars: data = remove_outliers_iqr(data, col) ``` Créer des diagrammes en boîte pour chaque variable par groupe de qualité : ```python import matplotlib.pyplot as plt import seaborn as sns for i, var in enumerate(explanatory_vars): plt.subplot(1, 3, i + 1) sns.boxplot(x='quality', y=var, data=data) plt.title(f'Diagramme en boîte de {var} par qualité') plt.xlabel('Qualité') plt.ylabel(var) plt.tight_layout() plt.show() ``` Ajuster le modèle ordinal pour le test des odds proportionnels : ```python from statsmodels.miscmodels.ordinal_model import OrderedModel from sklearn.preprocessing import StandardScaler data[explanatory_vars] = StandardScaler().fit_transform(data[explanatory_vars]) def fit_ordered_logistic_regression(data, response_var, explanatory_vars): model = OrderedModel(data[response_var], data[explanatory_vars], distr='logit') result = model.fit(method='bfgs', disp=False) return result response_var = 'quality' result = fit_ordered_logistic_regression(data, response_var, explanatory_vars) print(result.summary()) log_reduced = result.llf print(f"Log-vraisemblance du modèle réduit: {log_reduced}") ``` Calculer le modèle multinomial complet : ```python import statsmodels.api as sm data_sm = sm.add_constant(data[explanatory_vars]) model_full = sm.MNLogit(data[response_var], data_sm) result_full = model_full.fit(method='bfgs', disp=False) print(result_full.summary()) log_full = result_full.llf print(f"Log-vraisemblance du modèle complet: {log_full}") Calculer la statistique du rapport de vraisemblance LR_statistic = 2 * (log_full - log_reduced) print(f"Statistique du rapport de vraisemblance: {LR_statistic}") df1 = (data['quality'].nunique() - 1) * len(explanatory_vars) p_value = 1 - chi2.cdf(LR_statistic, df1) print(f"P-value: {p_value}") if p_value < 0.05: print("Rejet de l'hypothèse nulle : L'hypothèse des odds proportionnels est violée.") else: print("Échec de rejet de l'hypothèse nulle : L'hypothèse des odds proportionnels est satisfaite.") ``` Ajuster les modèles logistiques binaires pour vérifier l'hypothèse des odds proportionnels : ```python def fit_binary_models(data, explanatory_vars, y): qualities = np.sort(np.unique(y)) thresholds = qualities[:-1] p = len(explanatory_vars) n = len(y) K_1 = len(thresholds) binary_models = [] beta_hat = np.full((K_1, p + 1), np.nan) var_hat = [] for j, t in enumerate(thresholds): z_j = (y > t).astype(int) model = sm.Logit(z_j, sm.add_constant(data[explanatory_vars])) res = model.fit(disp=0) binary_models.append(res) beta_hat[j, :] = res.params.values var_hat.append(res.cov_params().values) return binary_models, beta_hat, sm.add_constant(data[explanatory_vars]), var_hat, thresholds binary_models, beta_hat, X_with_const, var_hat, thresholds = fit_binary_models(data, explanatory_vars, data[response_var]) print("Coefficients estimés (beta_hat):") print(beta_hat) ``` Calculer les probabilités ajustées pour construire la matrice diagonale Wjl : ```python def compute_pi_hat(binary_models, X_with_const): n = X_with_const.shape[0] K_1 = len(binary_models) pi_hat = np.full((n, K_1), np.nan) for m, model in enumerate(binary_models): pi_hat[:, m] = model.predict(X_with_const) return pi_hat pi_hat = compute_pi_hat(binary_models, X_with_const) print("Forme de pi_hat:", pi_hat.shape) print("Aperçu de pi_hat:\n", pi_hat[:5, :]) ``` Assembler la matrice de variance-covariance globale V̂(β̃) : ```python def assemble_varBeta(pi_hat, X_with_const): X = X_with_const.values n, p1 = X.shape p = p1 - 1 K_1 = pi_hat.shape[1] varBeta = np.zeros((K_1 * p, K_1 * p)) for j in range(K_1): pi_j = pi_hat[:, j] Wj = np.diag(pi_j * (1 - pi_j)) XtWjX_inv = np.linalg.pinv(X.T @ Wj @ X) var_block_diag = XtWjX_inv[1:, 1:] row_start = j * p row_end = (j + 1) * p varBeta[row_start:row_end, row_start:row_end] = var_block_diag for l in range(j + 1, K_1): pi_l = pi_hat[:, l] Wml = np.diag(pi_l - pi_j * pi_l) Wl = np.diag(pi_l * (1 - pi_l)) XtWlX_inv = np.linalg.pinv(X.T @ Wl @ X) block_cov = (XtWjX_inv @ (X.T @ Wml @ X) @ XtWlX_inv)[1:, 1:] col_start = l * p col_end = (l + 1) * p varBeta[row_start:row_end, col_start:col_end] = block_cov varBeta[col_start:col_end, row_start:row_end] = block_cov.T return varBeta varBeta = assemble_varBeta(pi_hat, X_with_const) print("Forme de varBeta:", varBeta.shape) print("Aperçu de varBeta:\n", varBeta[:5, :5]) ``` Construire la matrice de contraste (D) : ```python def construct_D(K_1, p): D = np.zeros(((K_1 - 1) * p, K_1 * p)) I = np.eye(p) for i in range(K_1 - 1): for j in range(K_1): temp = I if j == i else (-I if j == i + 1 else np.zeros((p, p))) row_start = i * p row_end = (i + 1) * p col_start = j * p col_end = (j + 1) * p D[row_start:row_end, col_start:col_end] += temp return D D = construct_D(len(thresholds), len(explanatory_vars)) print("Forme de D:", D.shape) print("Aperçu de D:\n", D[:5, :5]) ``` Calculer la statistique de Wald pour tester l'hypothèse des odds proportionnels : ```python def wald_statistic(D, betaStar, varBeta): Db = D @ betaStar V = D @ varBeta @ D.T inv_V = np.linalg.inv(V) X2 = float(Db.T @ inv_V @ Db) return X2 K_1 = len(binary_models) p = len(explanatory_vars) D = construct_D(K_1, p) betaStar = beta_hat[:, 1:].flatten() X2 = wald_statistic(D, betaStar, varBeta) ddl = (K_1 - 1) * p pval = 1 - chi2.cdf(X2, ddl) print(f"Statistique de Wald X² = {X2:.4f}") print(f"Degrés de liberté = {ddl}") print(f"P-value = {pval:.4g}") if pval < 0.05: print("Rejet de l'hypothèse nulle : L'hypothèse des odds proportionnels est violée.") else: print("Échec de rejet de l'hypothèse nulle : L'hypothèse des odds proportionnels est satisfaite.") ``` Conclusion Cet article visait à illustrer comment tester l'hypothèse des odds proportionnels dans le cadre de la régression logistique ordinale et à encourager les lecteurs à consulter l'article de Brant (1990) pour obtenir une compréhension plus approfondie. Brant aborde également la validation générale du modèle de régression logistique ordrale, comme le test de l'ajustement de la variable latente ( Y^* ) à une distribution logistique ou l'évaluation de l'adéquation d'une fonction de lien alternative. Dans cet article, nous avons effectué une évaluation globale de l'hypothèse des odds proportionnels sans investiguer les coefficients spécifiques qui pourraient être responsables des violations. Brant propose une analyse plus fine, ce qui motive notre suggestion de lecture intégrale de son travail. Notes sur les Données Les données de cet article sont issues du dataset "Wine Quality" et sont licensed sous la licence Creative Commons Attribution 4.0 International (CC BY 4.0). elles peuvent être utilisées, distribuées et adaptées, même à des fins commerciales, à condition d'accorder le crédit approprié. Références Brant, R. (1990). "Assessing Proportionality in the Proportional Odds Model for Ordinal Logistic Regression." Biometrics, 1171–78. McCullagh, P. (1980). "Regression Models for Ordinal Data." Journal of the Royal Statistical Society: Series B (Methodological), 42(2), 109–27. Wasserman, L. (2013). All of Statistics: A Concise Course in Statistical Inference. Springer Science & Business Media. Cortez, P., Cerdeira, A., Almeida, F., Matos, T., & Reis, J. (2009). Wine Quality [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C56S3T. Évaluation de l'Événement par les Professionnels de l'Industrie Les professionnels de la statistique et de la science des données soulignent l'importance de tester rigoureusement l'hypothèse des odds proportionnels pour assurer la validité des modèles de régression logistique ordinale. Les méthodes proposées par Brant (1990) offrent un moyen robuste et reconnu pour évaluer cette hypothèse. L'implémentation pratique de ces tests en Python, comme démontrée dans cet article, permet aux professionnels de la data science d'intégrer ces vérifications dans leurs analyses, améliorant ainsi la fiabilité des modèles et l'interprétation des résultats. Profil de l'Entreprise L'entreprise mentionnée est spécialisée dans l'analyse de données pour l'industrie vinicole. Elle utilise des techniques avancées de régression logistique ordinale pour prédire la qualité des vins basée sur divers facteurs chimiques et sensoriels. La validation statistique de ces modèles est cruciale pour fournir des insights precis et pertinents aux vignerons et aux amateurs de vin.