Migliori pratiche per il debugging degli errori nella regressione logistica con Python

Le migliori strategie per il debug degli errori nella regressione logistica con Python

Ottimizzazione delle prestazioni utilizzando dati non strutturati del mondo reale

È stato scritto molto sulle basi della Regressione Logistica (LR) – la sua versatilità, le prestazioni collaudate nel tempo, anche la matematica sottostante. Ma sapere come implementare con successo l’LR e risolvere gli errori inevitabili è molto più difficile. Queste informazioni vivono nelle profondità dei siti QA, negli articoli accademici o semplicemente vengono con l’esperienza.

La realtà è che non ogni output sarà pulito come l’iconico dataset dell’iris, che classifica rapidamente i tipi di fiori. Su dataset di grandi dimensioni (i dataset che probabilmente usi per il lavoro), i modelli di LR avranno probabilmente problemi. Coefficienti incredibilmente alti. Errori standard NaN. Manca la convergenza.

A cosa possiamo attribuire la colpa? Dati errati? Il tuo modello?

Per organizzare un po’ lo scenario, ho condotto alcune ricerche e ho compilato una lista degli errori comuni dell’LR, delle cause e delle possibili soluzioni.

• Esposizione del potere del filtro di Kalman
• Comprensione delle Variabili Strumentali
• Il problema della colorazione dei grafi soluzioni esatte ed euristiche

La tabella sopra non è esaustiva, ma è tutta in un unico posto. Di seguito utilizzerò dati fabbricati (che sono sbilanciati, sparsi, categorici… ma tipici) per ricreare artificialmente questi errori e poi correggerli nuovamente.

Ma prima.

Sfondo Rapido

Mentre inizi a pensare all’LR, ci sono alcune cose uniche di questi modelli che dovresti ricordare:

1. L’LR è tecnicamente un modello di probabilità, non un classificatore
2. L’LR richiede un confine di decisione lineare; assume linearità tra le caratteristiche e la variabile target, che puoi determinare con visualizzazioni, tabelle incrociate o analisi ConvexHull da scipy.spatial
3. Non puoi avere valori mancanti o outlier significativi
4. All’aumentare del numero di caratteristiche, è più probabile che si verifichi la multicollinearità e l’eccessivo adattamento (risolvibili con VIF e regolarizzazione, rispettivamente)
5. I pacchetti LR più popolari in Python provengono da sklearn e statsmodels

Ora entriamo in alcuni problemi comuni.

Problema 1: Avviso di convergenza, ottimizzazione fallita

Utilizzando il pacchetto python statsmodels, potresti eseguire un’LR di base e ricevere l’avviso “ConvergenceWarning: Massima verosimiglianza, ottimizzazione fallita nella convergenza.”

import statsmodels.formula.api as smfmodel = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \    C(educ,Treatment('Some High School')) + \    C(emp,Treatment('Not Employed')) + \    C(gender,Treatment('Male')) + \    C(hispanic,Treatment('N')) + \    C(dept,Treatment('Medical Oncology')) + \    C(sign_status,Treatment('Activated')) + \    C(mode,Treatment('Internet')) + \    C(inf,Treatment('Missing'))", data=df_a)    results = model.fit(method='bfgs')print(results.summary())

I coefficienti e gli errori standard possono anche popolarsi, come al solito, inclusi quelli della LogisticRegression nel pacchetto sklearn.linear_model.

Dati i risultati, potresti essere tentato di fermarti qui. Ma non farlo. Vuoi assicurarti che il tuo modello converga per produrre la migliore (più bassa) funzione di costo e aderenza del modello [1].

Puoi risolvere questo problema in statsmodels o sklearn modificando il solver/metodo o aumentando il parametro maxiter. In sklearn, il parametro di regolazione C può anche essere abbassato per applicare una maggiore regolarizzazione L2 e può essere testato iterativamente lungo una scala logaritmica (100, 10, 1, 0.1, 0.01, ecc.) [2].

Nel mio caso specifico, ho aumentato maxiter=300000 e il modello ha convergato. Il motivo per cui ha funzionato è che ho fornito al modello (ad esempio, solvers) più tentativi di convergenza e di trovare i migliori parametri [3]. E quei parametri, come i coefficienti e i valori p, diventano effettivamente più accurati.

Problema 2: Aggiunto una caratteristica, ma gli output del LR non si sono aggiornati

Questo è facile da trascurare, ma facile da diagnosticare. Se stai lavorando con molte caratteristiche o non l’hai notato durante la pulizia dei dati, potresti includere accidentalmente una caratteristica categorica nel tuo modello di regressione logistica che è quasi costante o ha solo un livello… cose brutte. L’inclusione di tale caratteristica renderà i coefficienti e gli errori standard senza alcun avviso.

Ecco il nostro modello e l’output senza la cattiva caratteristica.

model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \    C(educ,Treatment('Some High School')) + \    C(emp,Treatment('Not Employed')) + \    C(gender,Treatment('Maschile')) + \    C(hispanic,Treatment('N')) + \    C(dept,Treatment('Medical Oncology')) + \    C(sign_status,Treatment('Attivato')) + \    C(mode,Treatment('Internet')) + \    C(inf,Treatment('Mancante'))", data=df_a)    results = model.fit(method='bfgs',maxiter=30000)print(results.summary())

Ora, aggiungiamo una caratteristica con un solo livello, con la categoria di riferimento impostata su ‘Mancante’.

model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \    C(educ,Treatment('Some High School')) + \    C(emp,Treatment('Not Employed')) + \    C(gender,Treatment('Maschile')) + \    C(hispanic,Treatment('N')) + \    C(dept,Treatment('Medical Oncology')) + \    C(sign_status,Treatment('Attivato')) + \    C(mode,Treatment('Internet')) + \    C(inf,Treatment('Mancante')) + \    C(type,Treatment('Mancante'))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)print(results.summary())

E qui vediamo l’impatto di quella caratteristica superflua… che non ha avuto alcun impatto. Il valore R-quadrato rimane lo stesso, evidenziato sopra e sotto.

La regressione logistica non può fare stime significative su una caratteristica con un livello o valori costanti e potrebbe scartarla dal modello. Questo non influisce sul modello stesso, ma la pratica migliore è includere solo le caratteristiche nel modello che hanno un impatto positivo.

Una semplice value_counts() della nostra caratteristica di tipo rivela che ha un solo livello, indicando che si desidera eliminarla dal modello.

Problema 3: “Inverting Hessian failed” (e errori standard NaN)

Questo è un grande errore e si verifica spesso. Creiamo questo nuovo errore includendo altre quattro caratteristiche nel nostro modello LR: VA, VB, VC e VI. Sono tutte categoriche, ognuna con 3 livelli (0, 1 e “Mancante” come riferimento).

modello = smf.logit("target_soddisfazione ~ C(età,Trattamento('0-13')) + \    C(istruzione,Trattamento('Qualche scuola superiore')) + \    C(impiego,Trattamento('Non impiegato')) + \    C(genere,Trattamento('Maschio')) + \    C(ispanico,Trattamento('N')) + \    C(dipartimento,Trattamento('Oncologia medica')) + \    C(stato_firma,Trattamento('Attivato')) + \    C(modalità,Trattamento('Internet')) + \    C(info,Trattamento('Mancante')) + \    C(tipo,Trattamento('Mancante')) + \    C(VA,Trattamento('Mancante')) + \    C(VB,Trattamento('Mancante')) + \    C(VC,Trattamento('Mancante')) + \    C(VI,Trattamento('Mancante'))", data=df_a) risultati = modello.fit(method='bfgs', maxiter=300000)print(risultati.summary())

Ecco il nostro nuovo errore:

Esplorando un po’ gli output, vediamo che i coefficienti sono visualizzati ma non gli errori standard o i valori p. Sklearn fornisce anche i coefficienti.

Ma perché i valori NaN? E perché hai bisogno di invertire il Hessiano, comunque?

Molti algoritmi di ottimizzazione utilizzano l’inverso del Hessiano (o una stima di esso) per massimizzare la funzione di verosimiglianza. Quindi quando l’inversione fallisce, significa che l’Hessiano (tecnicamente, una derivata seconda del logaritmo di verosimiglianza) non è definito positivo [4]. Visivamente, alcune caratteristiche hanno curvature molto elevate mentre altre hanno curvature più basse, e il risultato è che il modello non può produrre errori standard per i parametri.

Più precisamente, quando un Hessiano non è invertibile, nessuna modifica computazionale può renderlo invertibile perché semplicemente non esiste [5].

Cosa sono le cause di questo problema? Ci sono tre più comuni:

1. Ci sono più variabili caratteristiche rispetto alle osservazioni
2. La/e caratteristica/e categorica/e ha/hanno frequenza molto bassa
3. Esiste una multicollinearità tra le caratteristiche

Dato che il nostro dataset ha molte righe (~25.000) rispetto alle caratteristiche (14) possiamo tranquillamente ignorare la prima possibilità (anche se esistono soluzioni per questo preciso problema [6]). La terza possibilità potrebbe essere in gioco e possiamo verificarla con il fattore di inflazione della varianza (VIF). La seconda possibilità è un po’ più facile da diagnosticare, quindi cominciamo da lì.

Nota che le caratteristiche VA, VB, VC e VI sono principalmente composte da 1 e 0.

In particolare, la caratteristica VI ha una frequenza molto bassa (relativamente) per la categoria ‘Mancante’, effettivamente dell’ordine dello 0,03% (9/24.874).

Supponiamo di confermare con il contesto aziendale che possiamo unire 1 e “Mancante” insieme. Oppure, perlomeno, qualsiasi potenziale conseguenza di refactoring dei dati in questo modo sarebbe molto inferiore rispetto all’accettazione di un modello con errori noti (e senza errori standard da considerare).

Quindi abbiamo creato VI_alt, che ha 2 livelli, e 0 può fungere da riferimento.

Sostituendo VI_alt con VI nel modello,

model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \    C(educ,Treatment('Some High School')) + \    C(emp,Treatment('Not Employed')) + \    C(gender,Treatment('Male')) + \    C(hispanic,Treatment('N')) + \    C(dept,Treatment('Medical Oncology')) + \    C(sign_status,Treatment('Activated')) + \    C(mode,Treatment('Internet')) + \    C(inf,Treatment('Missing')) + \    C(type,Treatment('Missing')) + \    C(VA,Treatment('Missing')) + \    C(VB,Treatment('Missing')) + \    C(VC,Treatment('Missing')) + \    C(VI_alt,Treatment(0))", data=df_a) risultati = model.fit(method='bfgs', maxiter=300000)print(risultati.summary())

C’è un leggero miglioramento nel modello complessivo perché converge senza errori, i coefficienti vengono visualizzati e gli errori standard sono ora presenti. Ancora una volta, è comunque un modello poco adattato, ma è ora un modello funzionante, che è il nostro obiettivo qui.

La terza e ultima possibilità per cui l’inversione della Hessiana fallisce è rappresentata dalla multicollinearità. Ora, la multicollinearità è un argomento scottante nell’apprendimento automatico, penso perché (1) affligge molti modelli popolari come LR, regressione lineare, KNN e Naive Bayes, (2) l’euristica VIF per rimuovere le feature si rivela essere più un’arte che una scienza e (3) in definitiva, gli addetti ai lavori non sono d’accordo sul fatto di rimuovere o meno tali feature in primo luogo se ciò comporta l’introduzione di un bias di selezione o la perdita di informazioni chiave del modello [5-9].

Non mi addentrerò in questa questione. È complicata.

Ma mostrerò come calcolare i VIF e forse possiamo trarre le nostre conclusioni.

Innanzitutto, i VIF chiedono “Quanto bene una delle feature è spiegata congiuntamente da tutte le altre feature?” La stima del VIF di una feature si “gonfia” quando c’è una dipendenza lineare con altre feature.

Dati tutti i fattori nel nostro dataset, calcoliamo i VIF di seguito.

from statsmodels.stats.outliers_influence import variance_inflation_factorvariables = results.model.exog vif = [variance_inflation_factor(variables, i) for i in range(variables.shape[1])]vifs = pd.DataFrame({'variabili':results.model.exog_names,'vif':[ '%.2f' % elem for elem in vif ]})vifs.sort_index(ascending=False).head(14)

Da questa lista parziale, si può notare che le feature VA, VB, VD mostrano tutte VIF verso l’infinito, che supera di gran lunga la soglia “rule of thumb” di 5. Ma dobbiamo fare attenzione a euristiche come queste, poiché le soglie VIF hanno due accortezze:

1. Il valore di 5 è relativo ad altre feature – ad esempio, se la grande maggioranza dei VIF delle feature è inferiore a 7 e un subset più piccolo dei VIF delle feature è superiore a 7, allora 7 potrebbe essere una soglia più ragionevole.
2. Le feature categoriche con un numero minore di casi nella categoria di riferimento rispetto agli altri livelli produrranno VIF elevati anche se tali feature non sono correlate con altre variabili nel modello [8].

Nel nostro caso, le feature VA, VB, VC sono tutte altamente collineari. Le tabelle incrociate confermano questo, e una matrice di correlazione di Pearson lo farebbe anche se fossero variabili continue.

Il consenso generale per risolvere questo problema è quello di eliminare sistematicamente una feature alla volta, valutare tutti i VIF, annotare possibili miglioramenti e continuare fino a quando tutti i VIF scendono al di sotto della soglia selezionata. Bisogna prestare attenzione alla perdita di potenziale potere esplicativo delle feature, sia in relazione al contesto aziendale che alla variabile target. Test statistici confermativi come il chi-quadro e la visualizzazione possono aiutare a decidere quale feature eliminare tra due scelte possibili.

Eliminiamo VB e impostiamo i cambiamenti dei VIF:

model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \    C(educ,Treatment('Some High School')) + \    C(emp,Treatment('Not Employed')) + \    C(gender,Treatment('Male')) + \    C(hispanic,Treatment('N')) + \    C(dept,Treatment('Medical Oncology')) + \    C(sign_status,Treatment('Activated')) + \    C(mode,Treatment('Internet')) + \    C(inf,Treatment('Missing')) + \    C(type,Treatment('Missing')) + \    C(VA,Treatment('Missing')) + \    C(VC,Treatment('Missing')) + \    C(VI_alt,Treatment(0))", data=df_a) risultati = model.fit(method='bfgs', maxiter=300000)print(risultati.summary())

VA e VC hanno ancora dei VIF all’infinito. Ancora peggio, il modello sta ancora generando errori standard NaN (non mostrati). Lasciamo cadere VC.

model = smf.logit("soddisfazione_obiettivo ~ C(eta, Trattamento('0-13')) + \    C(educazione, Trattamento('Diploma di scuola media')) + \    C(occupazione, Trattamento('Non occupato')) + \    C(sesso, Trattamento('Maschio')) + \    C(ispanico, Trattamento('N')) + \    C(reparto, Trattamento('Oncologia medica')) + \    C(stato_firma, Trattamento('Attivata')) + \    C(_, Trattamento('Internet')) + \    C(inf, Trattamento('Mancante')) + \    C(tipo, Trattamento('Mancante')) + \    C(VA, Trattamento('Mancante')) + \    C(VI_alt, Trattamento(0))", data=df_a) risultati = model.fit(method='bfgs', maxiter=300000)print(risultati.summary())

Infine, il modello produce errori standard, quindi anche se i VIF della variabile VA sono ancora > 5, non è un problema a causa della seconda avvertenza sopra, la categoria di riferimento che ha un piccolo numero di casi rispetto agli altri livelli.

Extra Credit: Diciamo che sai assolutamente che VA, VB e VC sono fondamentali per spiegare la variabile target e che hai bisogno di averle nel modello. Date le nuove caratteristiche, assumi che l’ottimizzazione stia lavorando su uno spazio complesso e che “la conversione in Hessian ha fallito” possa essere aggirata scegliendo nuovi risolutori e punti di partenza (start_params in sklearn). La formazione di un nuovo modello che non assuma limiti lineari sarebbe anche un’opzione.

Problema 4: Errore di “quasi-separazione” (possibilmente completo) e coefficienti/e errori standard improcedibili

Potremmo entusiasmarci nel vedere coefficienti elevati e punteggi di precisione elevati. Ma spesso, queste stime sono incredibilmente grandi e sono causate da un altro problema comune: la separazione perfetta.

La separazione perfetta (o quasi perfetta) si verifica quando una o più caratteristiche sono fortemente associate alla variabile target. In effetti, possono essere quasi identiche ad essa.

Posso produrre questo errore prendendo la variabile target soddisfazione_obiettivo e creando una nuova caratteristica chiamata nuova_soddisfazione_obiettivo che è del 95% simile ad essa.

df_a['nuova_soddisfazione_obiettivo'] = pd.concat([pd.DataFrame(np.where(df_a.soddisfazione_obiettivo.iloc[:1000]==1,0,df_a.soddisfazione_obiettivo[:1000]),columns=["soddisfazione_obiettivo"]),pd.DataFrame(df_a.soddisfazione_obiettivo.iloc[1000:],columns=['soddisfazione_obiettivo'])],axis=0)

E inseriamo nuova_soddisfazione_obiettivo nel modello.

model = smf.logit("soddisfazione_obiettivo ~ C(eta, Trattamento('0-13')) + \    C(educazione, Trattamento('Diploma di scuola media')) + \    C(occupazione, Trattamento('Non occupato')) + \    C(sesso, Trattamento('Maschio')) + \    C(ispanico, Trattamento('N')) + \    C(reparto, Trattamento('Oncologia medica')) + \    C(stato_firma, Trattamento('Attivata')) + \    C(_, Trattamento('Internet')) + \    C(inf, Trattamento('Mancante')) + \    C(tipo, Trattamento('Mancante')) + \    C(VA, Trattamento('Mancante')) + \    C(VI_alt, Trattamento(0)) + \    C(nuova_soddisfazione_obiettivo, Trattamento(0))", data=df_a) risultati = model.fit(method='lbfgs', maxiter=1000000)print(risultati.summary())

I coefficienti e gli errori standard vengono visualizzati, ma otteniamo questa avvertenza di quasi-separazione:

La caratteristica ha un coefficiente ridicolmente alto e un rapporto di probabilità di circa 70.000.000, il che sembra implausibile.

Sotto il cofano, ciò accade perché le caratteristiche che separano “troppo bene” creano una pendenza ingrandita, quindi un coefficiente e un errore standard elevati. È possibile che il modello LR potrebbe anche non convergere [10].

Quei due punti rossi, i casi classificati in modo errato rimossi, avrebbero effettivamente impedito una separazione perfetta, aiutando il modello LR a convergere e le stime dell’errore standard a essere più realistiche.

La cosa da ricordare sulla separazione perfetta in sklearn è che può generare un modello che sembra avere una precisione quasi perfetta, nel nostro caso 98%, ma in realtà ha una caratteristica new_target_satisfaction che è quasi duplicata del target target_satisfaction.

categorical_features = ['educ','emp','gender','hispanic','dept','sign_status','mode','inf','type','VA','VI_alt','new_target_satisfaction']df1_y = df_a['target_satisfaction']df1_x = df_a[['educ','emp','gender','hispanic','dept','sign_status','mode','inf','type','VA','VI_alt','new_target_satisfaction']]# crea un pipeline per tutte le caratteristiche categorichetransformer categorici = Pipeline(steps=[    ('ohe', OneHotEncoder(handle_unknown='ignore'))])# crea la trasformazione delle colonne globalicol_transformer = ColumnTransformer(transformers=[    ('ohe', OneHotEncoder(handle_unknown='ignore'), categorical_features)],   # ('num', numeric_transformer, numeric_features)],                                    remainder='passthrough')lr = Pipeline(steps = [('preprocessor', col_transformer),                            ('classifier', LogisticRegression(solver='lbfgs',max_iter=200000))])#['liblinear', 'newton-cg', 'lbfgs', 'sag', 'saga']X_train, X_test, y_train, y_test = train_test_split(df1_x, df1_y, test_size=0.2, random_state=101)lr_model = lr.fit(X_train, y_train)y_pred = lr_model.predict(X_test)# per sfruttare il ripristino della soglia# THRESHOLD = 0.5# y_pred = np.where(lr_model.predict_proba(X_test)[:,1] > THRESHOLD, 1, 0)print(classification_report(y_test, y_pred)) 

La soluzione più comune sarebbe semplicemente eliminare la caratteristica. Ma ci sono un numero crescente di soluzioni alternative:

1. Applicare la correzione di Firth, che massimizza una funzione di verosimiglianza “penalizzata”. Attualmente, statsmodels non dispone di questa funzionalità ma R sì [11]
2. La regressione penalizzata può anche funzionare, in particolare testando combinazioni di risolutori, penalizzazioni e tassi di apprendimento [2]. Ho mantenuto new_target_satisfaction nel modello e ho provato varie combinazioni, ma ha fatto poca differenza in questo caso specifico
3. Alcuni professionisti sostituiscono manualmente alcune osservazioni selezionate casualmente nella caratteristica problematica per renderla meno perfettamente separata dal target, come aggiungere di nuovo quei cerchi rossi nell’immagine sopra [8, 10]. Eseguire una crosstab su quella caratteristica con il target può aiutarti a capire quale percentuale di casi scambiare. Mentre lo fai, potresti chiederti Perché? Perché sto rifattorizzando una caratteristica in questo modo solo perché il modello LR può accettarla? Per aiutarti a dormire meglio, alcune ricerche sostengono che la separazione perfetta sia solo sintomatica del modello, non dei nostri dati [8, 11]
4. Infine, alcuni professionisti contradittori non vedono nulla di sbagliato in coefficienti così elevati [8]. Un rapporto di probabilità molto alto suggerisce semplicemente che è fortemente suggestivo di una correlazione e predice quasi perfettamente. Fai attenzione alle conclusioni e lascialo così. La base di questo argomento è che i coefficienti elevati sono una conseguenza sfortunata del test di Wald e del rapporto di verosimiglianza che richiedono una valutazione delle informazioni con un’ipotesi alternativa

Conclusione

La regressione logistica è sicuramente versatile e potente, purché si riescano a superare le sfide dei dataset del mondo reale. Spero che questa panoramica fornisca una buona base per possibili soluzioni. Quale consiglio hai trovato più interessante? Ne hai altri da suggerire?

Grazie per la lettura. Condividi le tue opinioni nella sezione Commenti e fammi sapere quali altri problemi hai incontrato con la regressione logistica. Sarà un piacere connettersi e confrontarsi su LinkedIn.

Dai un’occhiata ad alcuni dei miei altri articoli:

Usare le reti bayesiane per prevedere il volume dei servizi ausiliari negli ospedali

Perché bilanciare le classi è sopravvalutato

Feature Engineering con i codici CPT

Citazioni

1. Allison, Paul. Convergence Failures in Logistic Regression. University of Pennsylvania, Philadelphia, PA. https://www.people.vcu.edu/~dbandyop/BIOS625/Convergence_Logistic.pdf
2. Geron, Aurelien. Hands-On Machine Learning with Scikit-Learn, Keras & Tensorflow. Seconda Edizione. Pubblicato da: O’Reilly, 2019.
3. Documentazione di scikit-learn, sklearn.linear_model.LogisticRegression. https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
4. Discussione su Google Groups. “MLE error: Warning: Inverting hessian failed: Maybe I cant use ‘matrix’ containers?” 2014. https://groups.google.com/g/pystatsmodels/c/aELcgNpg5f8
5. Gill, Jeff & King, Gary. What to Do When Your Hessian Is Not Invertible. Alternatives to Model Respecification in Nonlinear Estimation. SOCIOLOGICAL METHODS & RESEARCH, Vol. 33, №1, agosto 2004 54–87 DOI: 10.1177/0049124103262681. https://gking.harvard.edu/files/help.pdf
6. “Selezione delle variabili per un problema di classificazione binaria.” StackExchange, CrossValidated. Online at: https://stats.stackexchange.com/questions/64271/variable-selection-for-a-binary-classification-problem/64637#64637
7. “Nell’apprendimento supervisionato, perché è sbagliato avere variabili correlate.” StackExchange, Data Science. Online at: https://datascience.stackexchange.com/questions/24452/in-supervised-learning-why-is-it-bad-to-have-correlated-features
8. “Come gestire la separazione perfetta nella regressione logistica”. StackExchange, CrossValidated. Online at: https://stats.stackexchange.com/questions/11109/how-to-deal-with-perfect-separation-in-logistic-regression
9. Allison, Paul. “When Can You Safely Ignore Multicollinearity”. Statistical Horizons, 10 settembre 2012. Online at: https://statisticalhorizons.com/multicollinearity/
10. Cornell Statistical Consulting Unit. Separation and Convergence Issues in Logistic Regression. Statnews #82. Pubblicato: febbraio 2012. Online at: https://cscu.cornell.edu/wp-content/uploads/82_lgsbias.pdf
11. Logistf: Firth’s Bias-Reduced Logistic Regression. Documentazione di R. Pubblicato il 18 agosto 2023. Online at: https://cran.r-project.org/web/packages/logistf/index.html

Nota: Tutte le immagini sono dell’autore, salvo diversa indicazione.