Proportional Odds Modell zur ordinalen logistischen Regression wird vorgestellt und getestet.

Das proportionale Odds-Modell für ordinale logistische Regression wurde erstmals von McCullagh (1980) eingeführt. Dieses Modell erweitert die binäre logistische Regression auf Situationen, in denen die abhängige Variable ordinal ist, was bedeutet, dass sie geordnete kategoriale Werte enthält. Das proportionale Odds-Modell basiert auf mehreren Annahmen, darunter die Unabhängigkeit der Beobachtungen, die Linearität der Log-Odds, die Abwesenheit von Multikollinearität unter den Prädiktoren und die proportionale Odds-Annahme. Diese letzte Annahme besagt, dass die Regressionskoeffizienten über alle Schwellenwerte der ordinalen abhängigen Variable konstant sind. Es ist entscheidend, dass die proportionale Odds-Annahme erfüllt ist, um die Gültigkeit und Interpretierbarkeit des Modells zu gewährleisten. In dieser Arbeit konzentrieren wir uns auf zwei Ansätze, die Brant (1990) in seinem Artikel "Assessing Proportionality in the Proportional Odds Model for Ordinal Logistic Regression" vorgeschlagen hat, um die proportionale Odds-Annahme zu bewerten. Wir demonstrieren auch, wie diese Techniken in Python implementiert werden können, indem wir sie auf reale Daten anwenden. Ob Sie aus einem Hintergrund in Data Science, maschinellem Lernen oder Statistik kommen, soll dieser Artikel Ihnen helfen, das Modellfit in der ordinalen logistischen Regression zu verstehen. Die Arbeit ist in vier Hauptabschnitte gegliedert: Einführung in das proportionale Odds-Modell und seine Annahmen. Bewertung der proportionalen Odds-Annahme mithilfe des Likelihood-Verhältnis-Tests. Bewertung der proportionalen Odds-Annahme durch Separatfits. Beispiele zur Implementierung dieser Bewertungsmethoden in Python mit Daten. 1. Einführung in das proportionale Odds-Modell Wir nehmen an, wir haben N unabhängige Beobachtungen. Jede Beobachtung wird durch einen Vektor von p erklärenden Variablen ( X_i = (X_{i1}, X_{i2}, \ldots, X_{ip}) ) und eine abhängige oder Antwortvariable ( Y ) repräsentiert, die ordinale Werte von 1 bis K annimmt. Das proportionale Odds-Modell modelliert die kumulativen Verteilungswahrscheinlichkeiten der Antwortvariable ( Y ), definiert als ( \gamma_j = P(Y \leq j | X_i) ) für ( j = 1, 2, \ldots, K - 1 ), als Funktionen der erklärenden Variablen ( X_i ). Das Modell wird wie folgt formuliert: [ \text{logit}(\gamma_j) = \log\left(\frac{\gamma_j}{1 - \gamma_j}\right) = \theta_j - \beta^T X_i \quad (1) ] Dabei sind ( \theta_j ) die Achsenabschnitte für jede Kategorie ( j ) und erfüllen die Bedingung ( \theta_1 < \theta_2 < \cdots < \theta_{K-1} ), und ( \beta ) ist der Vektor der Regressionskoeffizienten, der für alle Kategorien gleich ist. Man beobachtet einen monotonen Trend in den Koeffizienten ( \theta_j ) über die Kategorien der Antwortvariable ( Y ). Dieses Modell wird auch als gruppiertes kontinuierliches Modell bezeichnet, da es durch die Annahme eines kontinuierlichen latenteen Variables ( Y^* ) abgeleitet werden kann. Dieses latente Variable folgt einem linearen Regressionsmodell mit bedingtem Mittelwert ( \eta = \beta^T X_i ) und ist durch Schwellenwerte ( \theta_j ) mit der Antwortvariable ( Y ) in Verbindung gebracht: [ Y^* = \beta^T X_i + \epsilon \quad (2) ] Dabei ist ( \epsilon ) ein Fehlerterm (zufälliger Rausch), der in der Regel einer standardlogistischen Verteilung folgt. Das latente Variable ( Y^* ) ist nicht direkt beobachtbar und in Intervalle durch die Schwellenwerte ( \theta_1, \theta_2, \ldots, \theta_{K-1} ) unterteilt, was die beobachtete ordinale Variable ( Y ) erzeugt. 2. Bewertung der proportionalen Odds-Annahme: Der Likelihood-Verhältnis-Test Um die proportionale Odds-Annahme in einem ordinalen logistischen Regressionsmodell zu evaluieren, schlägt Brant (1990) den Likelihood-Verhältnis-Test vor. Dieser Ansatz beginnt damit, ein weniger restriktives Modell zu fitten, in dem die Regressionskoeffizienten über die Kategorien variieren dürfen. Dieses Modell ist wie folgt definiert: [ \text{logit}(\gamma_j) = \theta_j - \beta_j^T X_i \quad \text{für } j = 1, \ldots, K-1 \quad (4) ] Dabei sind ( \beta_j ) die Vektoren der Regressionskoeffizienten für jede Kategorie ( j ). Die Koeffizienten ( \beta_j ) dürfen hierbei über die Kategorien variieren, was bedeutet, dass die proportionale Odds-Annahme nicht erfüllt ist. Wir führen dann den konventionellen Likelihood-Verhältnis-Test durch, um die Hypothese zu evaluieren: [ H_0 : \beta_j = \beta \quad \text{für alle } j = 1, 2, \ldots, K-1 \quad (5) ] Dieser Test vergleicht das unbeschränkte (nicht-proportionale oder gesättigte) Modell mit dem beschränkten (proportionalen Odds oder reduzierten) Modell. Der Likelihood-Verhältnis-Teststatistik ( \lambda ) folgt einer Chi-Quadrat-Verteilung mit Freiheitsgraden, die der Differenz der Anzahl der Parameter zwischen dem vollständigen und dem beschränkten Modell entsprechen. Angenommen, der vollständige Parameterspace ist: [ \Theta = (\theta_1, \theta_2, \ldots, \theta_q, \ldots, \theta_p) ] und der restriktive Parameterspace unter der Nullhypothese ist: [ \Theta_0 = (\theta_1, \theta_2, \ldots, \theta_q) ] Dann folgt die Likelihood-Verhältnis-Teststatistik ( \lambda ) einer Chi-Quadrat-Verteilung mit ( p - q ) Freiheitsgraden. Hierbei stellt ( p ) die Gesamtanzahl der Parameter im vollständigen (unbeschränkten oder "gesättigten") Modell dar, während ( K - 1 ) der Anzahl der Parameter im reduzierten (restriktiven) Modell entspricht. Für das ordinales logistische Regressionsmodell mit der proportionalen Odds-Annahme vergleichen wir zwei Modelle: Unbeschränktes Modell (nicht-proportionale Odds): Dieses Modell erlaubt es, dass jeder Ausgangsschwelle eigene Regressionskoeffizienten hat, was bedeutet, dass wir nicht davon ausgehen, dass die Regressionskoeffizienten über alle Schwellen gleich sind. Das Modell ist definiert als: [ \text{logit}(\gamma_j) = \theta_j - \beta_j^T X_i \quad \text{für } j = 1, \ldots, K-1 \quad (7) ] Die Gesamtanzahl der Parameter im unbeschränkten Modell beträgt: [ (K-1) \text{ Schwellen} + (K-1) \times p \text{ Steigungen} = (K-1)(p+1) ] Proportionales Odds-Modell: Dieses Modell nimmt an, dass ein einziger Satz von Regressionskoeffizienten für alle Schwellen gilt: [ \text{logit}(\gamma_j) = \theta_j - \beta^T X_i \quad \text{für } j = 1, \ldots, K-1 \quad (8) ] Die Gesamtanzahl der Parameter im proportionalen Odds-Modell beträgt: [ (K-1) \text{ Schwellen} + p \text{ Steigungen} = (K-1) + p ] Die Likelihood-Verhältnis-Teststatistik folgt einer Chi-Quadrat-Verteilung mit Freiheitsgraden: [ \text{df} = [(K-1) \times (p+1)] - [(K-1) + p] = (K-2) \times p ] Dieser Test bietet eine formale Möglichkeit, die proportionale Odds-Annahme für die gegebenen Daten zu bewerten. Bei einem Signifikanzniveau von 1%, 5% oder einem anderen konventionellen Schwellenwert wird die proportionale Odds-Annahme verworfen, wenn der Teststatistik ( \lambda ) den kritischen Wert der Chi-Quadrat-Verteilung mit ( (K-2) \times p ) Freiheitsgraden überschreitet. 3. Bewertung der proportionalen Odds-Annahme: Der Separate-Fits-Ansatz Um diesen Abschnitt zu verstehen, ist es wichtig, das Konzept und die Eigenschaften der Mahalanobis-Distanz zu kennen. Diese Distanz wird verwendet, um die Unterschiede zwischen zwei Vektoren in einem multivariaten Raum zu quantifizieren, der die gleiche Verteilung teilt. Seien ( x = (x_1, x_2, \ldots, x_p)^T ) und ( y = (y_1, y_2, \ldots, y_p)^T ). Die Mahalanobis-Distanz zwischen ihnen ist definiert als: [ D_M(x, y) = \sqrt{(x - y)^T \Sigma^{-1} (x - y)} ] wobei ( \Sigma ) die Kovarianzmatrix der Verteilung ist. Das quadrierte Mahalanobis-Distanz ( D_M(X, \mu)^2 ) folgt einer Chi-Quadrat-Verteilung mit ( p ) Freiheitsgraden, wenn ( X \sim N(\mu, \Sigma) ) ein p-dimensionaler normalverteilter Zufallsvektor mit Mittelwert ( \mu ) und Kovarianzmatrix ( \Sigma ) ist. Die natürliche Methode zur Evaluierung der proportionalen Odds-Annahme besteht darin, ( K - 1 ) binäre logistische Regressionsmodelle (wobei ( K ) die Anzahl der Kategorien der Antwortvariable ist) zu fitten und dann die statistischen Eigenschaften der geschätzten Parameter zu nutzen, um eine Teststatistik für die proportionale Odds-Hypothese zu konstruieren. Das Verfahren ist wie folgt: Separate binäre logistische Regressionsmodelle für jede Schwelle ( j = 1, 2, \ldots, K-1 ) der ordinalen Antwortvariable ( Y ) konstruieren. Für jede Schwelle ( j ) definieren wir eine binäre Variable ( Z_j ), die den Wert 1 hat, wenn die Beobachtung die Schwelle ( j ) überschreitet, und 0, wenn sie dies nicht tut: [ Z_j = \begin{cases} 1 & \text{wenn } Y > j \ 0 & \text{sonst} \end{cases} ] Die Wahrscheinlichkeit ( \pi_j = P(Z_j = 1 | X) = 1 - \gamma_j ) erfüllt das binäre logistische Regressionsmodell: [ \text{logit}(\pi_j) = \theta_j - \beta^T X_i \quad \text{für } j = 1, \ldots, K-1 \quad (10) ] Die Hypothese, dass die Regressionskoeffizienten ( \beta_j ) über alle ( K - 1 ) Modelle gleich sind, testen. Dies ist äquivalent zum Testen der Hypothese: [ H_0 : \beta_1 = \beta_2 = \cdots = \beta_{K-1} \quad (11) ] Seien ( \hat{\beta}j ) die Maximum-Likelihood-Schätzer der Regressionskoeffizienten für jedes binäre Modell und ( \hat{\beta} = (\hat{\beta}_1^T, \hat{\beta}_2^T, \ldots, \hat{\beta}{K-1}^T)^T ) der globale Vektor der Schätzer. Dieser Vektor ist asymptotisch normalverteilt, sodass ( \hat{\beta}_j \approx \beta ) mit Varianz-Kovarianz-Matrix ( \hat{\Sigma}(\hat{\beta}_j) ). Die allgemeinen Terme dieser Matrix, ( \text{cov}(\hat{\beta}_j, \hat{\beta}_k) ), müssen bestimmt werden und sind gegeben durch: [ \text{cov}(\hat{\beta}j, \hat{\beta}_k) = (X^T W{jj} X)^{-1} X^T W_{jk} X (X^T W_{kk} X)^{-1} ] wobei ( W_{jk} ) eine ( n \times n )-diagonale Matrix mit typischer Eintrag ( \pi_k - \pi_j \pi_k ) für ( j \leq k ) ist, und ( X ) eine ( n \times p )-Matrix der erklärenden Variablen, die links mit einer Spalte von Einsen erweitert ist. Eine Matrix ( \Delta ) konstruieren, die die Unterschiede zwischen den Koeffizienten ( \hat{\beta}_j ) erfasst. Jeder Vektor ( \hat{\beta}_j ) hat die Dimension ( p ). Die Matrix ( \Delta ) ist definiert als: [ \Delta = \begin{pmatrix} I & -I & 0 & \cdots & 0 \ 0 & I & -I & \cdots & 0 \ \vdots & \vdots & \vdots & \ddots & \vdots \ 0 & 0 & 0 & \cdots & I & -I \end{pmatrix} ] wobei ( I ) die ( p \times p )-Einheitsmatrix ist. Das Produkt ( \Delta \hat{\beta} ) ergibt einen Vektor der Unterschiede zwischen den Koeffizienten ( \hat{\beta}_j ). Die Wald-Teststatistik ( X^2 ) definieren, um die proportionale Odds-Annahme zu testen. Diese Statistik kann als Mahalanobis-Distanz zwischen dem Vektor ( \Delta \hat{\beta} ) und dem Nullvektor interpretiert werden. Die Wald-Teststatistik ist definiert als: [ X^2 = (\Delta \hat{\beta})^T \hat{V}(\Delta \hat{\beta}) ] und folgt asymptotisch einer Chi-Quadrat-Verteilung mit ( (K-2) \times p ) Freiheitsgraden unter der Nullhypothese. Die Herausforderung besteht darin, die Varianz-Kovarianz-Matrix ( \hat{V}(\hat{\beta}) ) zu bestimmen. In seinem Artikel bereitet Brant einen expliziten Schätzer für diese Matrix vor, der auf den Maximum-Likelihood-Schätzern ( \hat{\beta}_j ) aus jedem binären Modell basiert. Beispiel: Anwendung der beiden proportionalen Odds-Tests Die Daten für dieses Beispiel stammen aus dem "Wine Quality"-Datensatz, der Informationen über Rotweinproben und deren Qualitätseinschätzungen enthält. Der Datensatz umfasst 1.599 Beobachtungen und 12 Variablen. Die Zielvariable "Qualität" ist ordinal und ursprünglich von 3 bis 8 skaliert. Um genügend Beobachtungen in jeder Gruppe zu gewährleisten, kombinieren wir die Kategorien 3 und 4 (zu Kategorie 4) sowie die Kategorien 7 und 8 (zu Kategorie 7), sodass die Antwortvariable vier Ebenen hat. Wir behandeln Ausreißer in den erklärenden Variablen mit der Interquartilsdistanz-Methode (IQR) und wählen drei Prädiktoren—flache Säure, freies Schwefeldioxid und totales Schwefeldioxid—für unser ordinales logistisches Regressionsmodell, die wir standardisieren, sodass sie den Mittelwert 0 und die Standardabweichung 1 haben. Die Tabellen 1 und 2 unten präsentieren die Ergebnisse der drei binären logistischen Regressionsmodelle und des proportionalen Odds-Modells, jeweils. Es gibt einige Diskrepanzen in diesen Tabellen, insbesondere in den Koeffizienten für "flache Säure". Zum Beispiel beträgt der Unterschied im Koeffizienten für "flache Säure" zwischen dem ersten und zweiten binären Modell -0.280, während der Unterschied zwischen dem zweiten und dritten Modell 0.361 ist. Diese Unterschiede deuten darauf hin, dass die proportionale Odds-Annahme möglicherweise nicht erfüllt ist (es wird auch empfohlen, die Standardfehler der Unterschiede zu berechnen, um mehr Sicherheit zu erlangen). Um die Gesamtsignifikanz der proportionalen Odds-Annahme zu bewerten, führen wir den Likelihood-Verhältnis-Test durch, der einen Teststatistik von ( \text{LR} = 53.207 ) und einen p-Wert von ( 1.066 \times 10^{-9} ) bei Vergleich mit der Chi-Quadrat-Verteilung mit 6 Freiheitsgraden ergibt. Dies deutet darauf hin, dass die proportionale Odds-Annahme bei einem Signifikanzniveau von 5% verletzt ist, was darauf schließen lässt, dass das Modell für die Daten möglicherweise nicht angemessen ist. Wir verwenden auch den Separate-Fits-Ansatz, um diese Annahme weiter zu untersuchen. Die Wald-Teststatistik wird berechnet als ( X^2 = 41.880 ), mit einem p-Wert von ( 1.232 \times 10^{-7} ), ebenfalls basierend auf der Chi-Quadrat-Verteilung mit 6 Freiheitsgraden. Dies bestätigt, dass die proportionale Odds-Annahme bei einem Signifikanzniveau von 5% verletzt ist. Schlussfolgerung Dieses Paper hatte zwei Hauptziele: erstens, zu zeigen, wie man die proportionale Odds-Annahme im Kontext der ordinalen logistischen Regression testet, und zweitens, Leser dazu zu ermutigen, den Artikel von Brant (1990) für ein tieferes Verständnis des Themas zu lesen. Brants Arbeit reicht weiter als die Evaluation der proportionalen Odds-Annahme. Er diskutiert auch Methoden zur Bewertung der Gesamtheadäquanz des ordinalen logistischen Regressionsmodells. Zum Beispiel diskutiert er, wie man testen kann, ob die latente Variable ( Y^* ) tatsächlich einer logistischen Verteilung folgt oder ob eine alternative Link-Funktion angemessener sein könnte. In diesem Artikel haben wir uns auf eine globale Bewertung der proportionalen Odds-Annahme konzentriert, ohne detailliert zu untersuchen, welche spezifischen Koeffizienten für Verletzungen verantwortlich sein könnten. Brant behandelt jedoch auch diese feinere Analyse, weshalb wir Leser stark dazu ermutigen, seinen Artikel in voller Länge zu lesen. Bildrechte Alle Bilder und Visualisierungen in diesem Artikel wurden vom Autor unter Verwendung von Python (pandas, matplotlib, seaborn, plotly) und Excel erstellt, es sei denn, anders angegeben. Referenzen Brant, Rollin. 1990. "Assessing Proportionality in the Proportional Odds Model for Ordinal Logistic Regression." Biometrics, 1171–78. McCullagh, Peter. 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. Daten und Lizenz Der in diesem Artikel verwendete Datensatz steht unter der Creative Commons Attribution 4.0 International (CC BY 4.0) Lizenz. Er erlaubt die Nutzung, Vervielfältigung und Adaption, auch zu kommerziellen Zwecken, vorausgesetzt, angemessene Anerkennung wird gewährt. Der ursprüngliche Datensatz wurde veröffentlicht von [Name des Autors / Instituts] und ist hier verfügbar. Python-Code Importieren der Daten und Anzahl der Beobachtungen ```python import pandas as pd data = pd.read_csv("winequality-red.csv", sep=";") data.head() Repartition der Zielvariable "Qualität" data['quality'].value_counts(normalize=False).sort_index() Kombinieren der Kategorien 3 und 4, sowie 7 und 8 data['quality'] = data['quality'].replace({3: 4, 8: 7}) data['quality'].value_counts(normalize=False).sort_index() print("Anzahl der Beobachtungen:", data.shape[0]) ``` Ausreißerbehandlung der erklärenden Variablen (ohne die Zielvariable "Qualität") mit der IQR-Methode ```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)] for col in data.columns: if col != 'quality': data = remove_outliers_iqr(data, col) ``` Erstellen von Boxplots für jede Variable pro Gruppe der Qualität ```python var_names_without_quality = [col for col in data.columns if col != 'quality'] import matplotlib.pyplot as plt import seaborn as sns plt.figure(figsize=(15, 10)) for i, var in enumerate(var_names_without_quality): plt.subplot(3, 4, i + 1) sns.boxplot(x='quality', y=var, data=data) plt.title(f'Boxplot von {var} nach Qualität') plt.xlabel('Qualität') plt.ylabel(var) plt.tight_layout() plt.show() ``` Implementierung der ordinalen Regression für den Proportionalitäts-Test ```python from statsmodels.miscmodels.ordinal_model import OrderedModel from sklearn.preprocessing import StandardScaler explanatory_vars = ['flache Säure', 'freies Schwefeldioxid', 'totales Schwefeldioxid'] 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()) Berechnung der Log-Likelihood des reduzierten Modells log_reduced = result.llf print(f"Log-Likelihood des reduzierten Modells: {log_reduced}") ``` Berechnung des vollständigen multinomialen Modells ```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()) Berechnung der Log-Likelihood des vollständigen Modells log_full = result_full.llf print(f"Log-Likelihood des vollständigen Modells: {log_full}") Berechnung der Likelihood-Verhältnis-Statistik LR_statistic = 2 * (log_full - log_reduced) print(f"Likelihood-Verhältnis-Statistik: {LR_statistic}") Berechnung der Freiheitsgrade df1 = (num_of_thresholds - 1) * len(explanatory_vars) df2 = result_full.df_model - OrderedModel( data[response_var], data[explanatory_vars], distr='logit' ).fit().df_model print(f"Freiheitsgrade: {df1}") print(f"Freiheitsgrade für das vollständige Modell: {df2}") Berechnung des p-Werts from scipy.stats import chi2 print("Die LR-Statistik:", LR_statistic) p_value = chi2.sf(LR_statistic, df1) print(f"p-Wert: {p_value}") if p_value < 0.05: print("Nullhypothese ablehnen: Die proportionale Odds-Annahme wird verletzt.") else: print("Nullhypothese nicht ablehnen: Die proportionale Odds-Annahme gilt.") ``` (K - 1) binäre Logit-Schätzungen zur Überprüfung der proportionalen Odds-Annahme ```python import numpy as np import statsmodels.api as sm import pandas as pd def fit_binary_models(data, explanatory_vars, y): qualities = np.sort(np.unique(y)) # Sortierte einzigartige Kategorien von y thresholds = qualities[:-1] # Schwellenwerte zur Definition der binären Modelle (K - 1) p = len(explanatory_vars) n = len(y) K_1 = len(thresholds) binary_models = [] beta_hat = np.full((K_1, p+1), np.nan) # Speicherung der geschätzten Koeffizienten p_values_beta_hat = np.full((K_1, p+1), np.nan) var_hat = [] z_mat = pd.DataFrame(index=np.arange(n)) # Speicherung der binären Zielvariablen X_with_const = sm.add_constant(data[explanatory_vars]) # Designmatrix mit Achsenabschnitt for j, t in enumerate(thresholds): z_j = (y > t).astype(int) # Binäre Zielvariable: 1, wenn y > t, sonst 0 z_mat[f'z>{t}'] = z_j model = sm.Logit(z_j, X_with_const) res = model.fit(disp=0) binary_models.append(res) beta_hat[j, :] = res.params.values # Speicherung der Koeffizienten (mit Achsenabschnitt) p_values_beta_hat[j, :] = res.pvalues.values var_hat.append(res.cov_params().values) # Speicherung der vollständigen (p+1) × (p+1) Kovarianzmatrix return binary_models, beta_hat, X_with_const, var_hat, z_mat, thresholds Anwendung der Funktion auf die Daten binary_models, beta_hat, X_with_const, var_hat, z_mat, thresholds = fit_binary_models( data, explanatory_vars, data[response_var] ) Anzeige der geschätzten Koeffizienten print("Geschätzte Koeffizienten (beta_hat):") print(beta_hat) Anzeige der Designmatrix (Prädiktoren mit Achsenabschnitt) print("Designmatrix X (mit Achsenabschnitt):") print(X_with_const) Anzeige der Schwellenwerte print("Schwellenwerte:") print(thresholds) Anzeige der ersten Zeilen der binären Zielvariablenmatrix print("z_mat (erstellte binäre Zielvariablen):") print(z_mat.head()) ``` Berechnung der gefitteten Wahrscheinlichkeiten (π̂) für die Konstruktion der Wjl n × n Diagonalmatrix ```python def compute_pi_hat(binary_models, X_with_const): n = X_with_const.shape[0] # Anzahl der Beobachtungen K_1 = len(binary_models) # Anzahl der binären Modelle (K - 1) pi_hat = np.full((n, K_1), np.nan) # Initialisierung der Vorhersagematrix for m, model in enumerate(binary_models): pi_hat[:, m] = model.predict(X_with_const) return pi_hat Annahme: Sie haben bereits - binary_models: Liste der gefitteten Modelle von der vorherigen Funktion - X_with_const: Designmatrix mit Achsenabschnitt (n, p+1) pi_hat = compute_pi_hat(binary_models, X_with_const) Anzeige der Form und einer Vorschau der Vorhersagematrix print("Form von pi_hat:", pi_hat.shape) # Erwartete Form: (n, K - 1) print("Vorschau von pi_hat:\n", pi_hat[:5, :]) # Erste 5 Zeilen ``` Berechnung der Gesamtgeschätzten Kovarianzmatrix V̂(β̃) in zwei Schritten ```python def assemble_varBeta(pi_hat, X_with_const): X = X_with_const.values if hasattr(X_with_const, 'values') else np.asarray(X_with_const) n, p1 = X.shape # p1 = p + 1 (mit Achsenabschnitt) 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)) # Diagonale Gewichtsmatrix für Modell j XtWjX_inv = np.linalg.pinv(X.T @ Wj @ X) var_block_diag = XtWjX_inv[1:, 1:] # Entfernen des Achsenabschnitts 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 Berechnung der Gesamtgeschätzten Kovarianzmatrix varBeta = assemble_varBeta(pi_hat, X_with_const) Anzeige der Form und einer Vorschau print("Form von varBeta:", varBeta.shape) # Erwartet: ((K - 1) * p, (K - 1) * p) print("Vorschau von varBeta:\n", varBeta[:5, :5]) # Obere linke Ecke ``` Füllen der Diagonalblöcke der Gesamtgeschätzten Kovarianzmatrix (ohne Achsenabschnitt) ```python def fill_varBeta_diagonal(varBeta, var_hat): K_1 = len(var_hat) # Anzahl der binären Modelle p = var_hat[0].shape[0] - 1 # Anzahl der Prädiktoren (ohne Achsenabschnitt) for m in range(K_1): block = var_hat[m][1:, 1:] # Entfernen des Achsenabschnitts row_start = m * p row_end = (m + 1) * p varBeta[row_start:row_end, row_start:row_end] = block return varBeta Flachlegen der Steigungsparameter (ohne Achsenabschnitt) in betaStar betaStar = beta_hat[:, 1:].flatten() Füllen der Diagonalblöcke der Gesamtgeschätzten Kovarianzmatrix varBeta = fill_varBeta_diagonal(varBeta, var_hat) Diagnostische Ausgaben print("Form von varBeta nach Diagonalblockfüllung:", varBeta.shape) # Erwartet: ((K - 1) * p, (K - 1) * p) print("Vorschau von varBeta nach Diagonalblockfüllung:\n", varBeta[:5, :5]) # Obere linke Ecke ``` Konstruktion der (K - 2)p × (K - 1)p Kontrastmatrix ```python def construct_D(K_1, p): D = np.zeros(((K_1 - 1) * p, K_1 * p)) # Initialisierung der D-Matrix I = np.eye(p) # Einheitsmatrix für Blockeinfügung for i in range(K_1 - 1): # i = 0 bis (K - 2) for j in range(K_1): # j = 0 bis (K - 1) if j == i: temp = I # Positiver Einheitsblock elif j == i + 1: temp = -I # Negativer Einheitsblock else: temp = np.zeros((p, p)) # Nullblock sonst row_start = i * p row_end = (i + 1) * p col_start = j * p col_end = (j + 1) * p D[row_start:row_end