Problema 1: avviso di convergenza, ottimizzazione non riuscita
Utilizzando il pitone modelli statali pacchetto, è possibile eseguire un LR di base e ricevere l’avviso “ConvergenceWarning: l’ottimizzazione della massima probabilità non è riuscita a convergere”.
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 essere compilati, come al solito, compresi quelli di Regressione logistica nel sklearn.linear_model pacchetto.
Visti 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 l’adattamento del modello (1).
Puoi risolverlo in modelli statali O imparato modificando il risolutore/metodo o aumentando il valore maxiter parametro. In imparatoil parametro di ottimizzazione 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 particolare, sono aumentato a maxiter=300000 e il modello converge. Il motivo per cui ha funzionato è perché ho fornito al modello (ad esempio, i solutori) più tentativi di convergenza e di ricerca dei parametri migliori (3). E questi parametri, come i coefficienti e i valori p, diventano effettivamente più accurati.
Problema 2: aggiunta una funzionalità, ma le uscite LR non sono state aggiornate
Questo è facile da perdere, ma facile da diagnosticare. Se stai lavorando con molte funzionalità o non l’hai notato durante la pulizia dei dati, potresti includere accidentalmente una funzionalità categoriale nel tuo modello LR che è quasi costante o ha un solo livello… cose brutte. L’inclusione di tale funzionalità restituirà coefficienti ed errori standard senza alcun avviso.
Ecco il nostro modello e l’output senza la caratteristica negativa.
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'))", data=df_a)results = model.fit(method='bfgs',maxiter=30000)
print(results.summary())
Ora aggiungiamo una funzionalità di un livello tipocon 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('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'))", 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ò effettuare stime significative su una caratteristica con un livello o valori costanti e potrebbe eliminarla dal modello. Ciò non è determinante per il modello stesso, ma la pratica migliore è comunque quella di includere nel modello solo le funzionalità che hanno un impatto positivo.
Un semplice valore_conta() della nostra caratteristica tipo rivela che ha un livello, indicando che desideri eliminare la funzionalità dal tuo modello.
Problema 3: “Inversione dell’Hessian non riuscita” (ed errori standard NaN)
Questo è un grosso problema e succede spesso. Creiamo questo nuovo errore includendo altre quattro funzionalità nel nostro modello LR: VA, VB, VC e VI. Sono tutti categorici, ciascuno con 3 livelli (0, 1 e “Mancante” come riferimento).
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,Treatment('Missing'))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Ed ecco il nostro nuovo errore:
Esplorando un po’ gli output, vediamo i coefficienti renderizzati ma nessuno degli errori standard o dei valori p. Imparare fornisce anche i coefficienti.
Ma perché i valori NaN? E comunque perché hai bisogno di invertire l’Assia?
Molti algoritmi di ottimizzazione utilizzano l’inverso dell’Assia (o una sua stima) per massimizzare la funzione di verosimiglianza. Quindi, quando l’inversione fallisce, significa che l’Assia (tecnicamente, una derivata seconda della logverosimiglianza) non è un definito positivo (4). Visivamente, alcune caratteristiche hanno curvature enormi mentre altre hanno curvature più piccole 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 inverso perché semplicemente non esiste (5).
Quali sono le cause di ciò? Ce ne sono tre più comuni:
- Ci sono più variabili di funzionalità che osservazioni
- Le caratteristiche categoriali hanno livelli di frequenza molto bassi
- Esiste multicollinearità tra le caratteristiche
Poiché il nostro set di dati ha molte righe (~25.000) relative alle caratteristiche (14), possiamo tranquillamente ignorare la prima possibilità (sebbene esistano soluzioni per questo esatto problema (6)). Potrebbe essere in gioco la terza possibilità e possiamo verificarla con il fattore di inflazione della varianza (VIF). La seconda possibilità è un po’ più semplice da diagnosticare, quindi iniziamo da lì.
Si noti che le caratteristiche VA, VB, VC e VI sono costituite principalmente da 1 e 0.
In particolare, la caratteristica VI ha una frequenza (relativa) molto bassa per la categoria “Mancante”, in realtà dell’ordine dello 0,03% (9/24.874).
Diciamo che confermiamo con il nostro contesto aziendale che possiamo comprimere insieme 1 e “Mancanti”. O, per lo meno, qualsiasi potenziale conseguenza del refactoring dei dati in questo modo sarebbe molto inferiore all’accettazione di un modello con errori noti (e senza errori standard di cui parlare).
Quindi abbiamo creato VI_alt, che ha 2 livelli e 0 può servire come riferimento.
Scambiare VI_alt per 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) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
C’è un leggero miglioramento nel modello complessivo perché converge senza errori, vengono ora visualizzati i coefficienti di rendering e gli errori standard. Ancora una volta, è ancora un modello poco adattato, ma ora funziona modello, che è il nostro obiettivo qui.
La nostra terza e ultima possibilità per l’inversione del fallimento dell’Assia è la multi-collinearità. Ora, la multi-collinearità è un argomento caldo nell’apprendimento automatico, penso perché (1) affligge molti modelli popolari come LR, regressione lineare, KNN e Naive Bayes (2) l’euristica VIF per la rimozione delle funzionalità risulta essere più un’arte che scienza e (3) in definitiva, i professionisti non sono d’accordo sull’opportunità o meno di rimuovere tali funzionalità in primo luogo se ciò avviene a scapito dell’iniezione di bias di selezione o della perdita di informazioni chiave sul modello (5–9).
Non scenderò in quella tana del coniglio. È profondo.
Ma mostrerò come calcolare i VIF e forse potremo trarre le nostre conclusioni.
Innanzitutto, il VIF chiede davvero “Quanto bene una delle funzionalità è spiegata congiuntamente da tutte le altre mie funzionalità?” La stima VIF di una caratteristica verrà “gonfiata” quando esiste una dipendenza lineare con altre caratteristiche.
Date tutte le funzionalità del nostro set di dati, 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({'variables':results.model.exog_names,'vif':( '%.2f' % elem for elem in vif )})
vifs.sort_index(ascending=False).head(14)
Da questo elenco parziale, le caratteristiche VA, VB, VD mostrano tutte VIF verso l’infinito, che supera di gran lunga la soglia della “regola pratica” di 5. Ma dobbiamo fare attenzione alle euristiche come questa, poiché le soglie VIF hanno due avvertenze:
- Il limite 5 è relativo ad altre funzionalità: ad esempio, se la grande maggioranza dei VIF delle funzionalità scende al di sotto di 7 e un sottoinsieme più piccolo di VIF delle funzionalità è superiore a 7, allora 7 potrebbe essere una soglia più ragionevole.
- Le caratteristiche categoriche con un numero minore di casi nella categoria di riferimento rispetto ad altri livelli produrranno VIF elevati anche se tali caratteristiche non sono correlate con altre variabili nel modello (8).
Nel nostro caso, le caratteristiche VA, VB, VC sono tutte altamente colineari. Le tabelle incrociate lo confermano, e anche la matrice di correlazione di Pearson lo farebbe se fossero variabili continue.
Il consenso generale per risolvere questo problema è quello di eliminare sistematicamente una funzionalità alla volta, rivedere tutti i VIF, annotare i possibili miglioramenti e continuare finché tutti i VIF non scendono al di sotto della soglia selezionata. Occorre prestare particolare attenzione a perdere ogni possibile potere esplicativo delle caratteristiche, relative sia al contesto aziendale che alla variabile target. Test statistici di conferma come il chi quadrato e la visualizzazione possono aiutare a decidere quale caratteristica eliminare tra due possibili scelte.
Lasciamo cadere VB e notiamo le modifiche 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) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
VA e VC hanno ancora VIF all’infinito. Ancor peggio, il modello sta ancora eseguendo il rendering degli errori standard NaN (non mostrati). Lasciamo perdere VC.
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(VI_alt,Treatment(0))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Infine, il modello produce errori standard, quindi anche se i VIF della caratteristica VA sono ancora > 5, non è un problema a causa del secondo avvertimento precedente, poiché la categoria di riferimento ha un numero limitato di casi rispetto agli altri livelli.
Credito extra: supponiamo che tu sappia assolutamente che VA, VB e VC sono fondamentali per spiegare la variabile target e ne hai bisogno nel modello. Date le funzionalità aggiuntive, presupponiamo che l’ottimizzazione stia funzionando su uno spazio complesso e che il problema “Inversione dell’Hessian non riuscita” possa essere aggirato scegliendo nuovi solutori e punti di partenza (start_params in sklearn). Anche la formazione di un nuovo modello che non presuppone confini lineari sarebbe un’opzione.
Problema 4: Errore di “quasi-separazione possibilmente completa” (e coefficienti e/o errori standard irragionevolmente grandi)
Potremmo essere entusiasti di vedere coefficienti elevati e punteggi di alta precisione. Ma spesso queste stime sono inverosimili e sono causate da un altro problema comune: la perfetta separazione.
La separazione perfetta (o quasi perfetta) si verifica quando una o più caratteristiche sono fortemente associate alla variabile target. In effetti, potrebbero essere quasi identici ad esso.
Posso produrre questo errore prendendo la variabile target target_soddisfazione e creando una nuova funzionalità da esso chiamata nuovo_obiettivo_soddisfazione è simile al 95% ad esso.
df_a('new_target_satisfaction') =
pd.concat((pd.DataFrame(np.where(df_a.target_satisfaction.iloc(:1000)==1,0,df_a.target_satisfaction(:1000)),columns=("target_satisfaction")),
pd.DataFrame(df_a.target_satisfaction.iloc(1000:),columns=('target_satisfaction'))),axis=0)
E metti nuovo_obiettivo_soddisfazione 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(VI_alt,Treatment(0)) + \
C(new_target_satisfaction,Treatment(0))", data=df_a) results = model.fit(method='lbfgs', maxiter=1000000)
print(results.summary())
Rendering di coefficienti ed errori standard, ma riceviamo questo avviso di quasi-separazione:
La funzione ha un coefficiente ridicolmente alto e un rapporto di probabilità di circa 70.000.000, il che non è plausibile.
Dietro le quinte, ciò accade perché le caratteristiche che si separano “troppo bene” creano una pendenza gonfiata, quindi un coefficiente e un errore standard elevati. È possibile che anche il modello LR non converga (10).
Quei due punti rossi, eliminati i casi classificati erroneamente, avrebbero effettivamente impedito una perfetta separazione, aiutando il modello LR a convergere e le stime dell’errore standard a essere più realistiche.
La cosa da ricordare sulla perfetta separazione in imparato è che può produrre un modello che sembra avere una precisione quasi perfetta, nel nostro caso il 98%, ma in realtà ha una caratteristica nuovo_obiettivo_soddisfazione questo è un quasi duplicato dell’obiettivo target_soddisfazione.
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'))
# create a pipeline for all categorical features
categorical_transformer = Pipeline(steps=(
('ohe', OneHotEncoder(handle_unknown='ignore'))))
# create the overall column transformer
col_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)
# to leverage threshold resetting
# 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 funzionalità. Ma ci sono un numero crescente di soluzioni alternative:
- Applicare la correzione di Firth, che massimizza una funzione di verosimiglianza “penalizzata”. Attualmente, modelli statali non ha questa funzionalità disponibile ma R sì (11)
- Anche la regressione penalizzata può funzionare, in particolare testando combinazioni di risolutori, penalità e tassi di apprendimento (2). Ho continuato nuovo_obiettivo_soddisfazione nel modello e ho provato varie combinazioni, ma in questo caso particolare ha fatto poca differenza
- Alcuni praticanti ne scambieranno manualmente alcuni osservazioni selezionate casualmente nella caratteristica problematica per renderla meno perfettamente separata dal bersaglio, come aggiungere nuovamente quei cerchi rossi nell’immagine sopra (8, 10). L’esecuzione di una tabella incrociata su tale funzionalità con la destinazione può aiutare a determinare quale percentuale di casi scambiare. Mentre lo fai potresti chiederti Perché? Perché sto eseguendo il refactoring di una funzionalità come questa solo in modo che il modello LR possa accettarla? Per aiutarti a dormire meglio, alcune ricerche sostengono che la separazione perfetta è solo sintomatica del modello, non dei nostri dati (8, 11)
- Infine, alcuni praticanti contrarian in realtà non vedono nulla di sbagliato in coefficienti così alti (8). Un rapporto di probabilità molto elevato suggerisce semplicemente che è fortemente indicativo di un’associazione e che la previsione è quasi perfetta. Avvertite i risultati e lasciate le cose come stanno. La base di questa argomentazione è che i coefficienti elevati sono una sfortunata conseguenza del test di Wald e del rapporto di verosimiglianza che richiede una valutazione delle informazioni con un’ipotesi alternativa
Conclusione
La regressione logistica è sicuramente versatile e potente se riesci a superare le sfide derivanti dai set di dati del mondo reale. Spero che questa panoramica fornisca una buona base per le possibili soluzioni. Quale suggerimento hai trovato più interessante? Quali sono gli altri?
Fonte: towardsdatascience.com
