Modelli di regressione multinivello e il paradosso di Simpson

Modelli multinivello e paradosso di Simpson.

Evitare conclusioni false con gli strumenti appropriati

L'analisi dei dati influenza le nostre conclusioni, ma dovremmo utilizzare gli strumenti appropriati per intraprendere il percorso giusto. Foto di Brendan Church su Unsplash

L’analisi dei dati è – come indica il nome di questa professione – una parte integrante del lavoro di un Data Scientist, che va dalle statistiche descrittive e dai modelli di regressione semplici fino ad approcci sofisticati di machine learning. Tuttavia, questi approcci devono essere gestiti con attenzione e la selezione dell’approccio corretto è tutto tranne che banale. I dati complessi spesso includono strutture nascoste che, se non considerate adeguatamente, possono portare a fallacie che possono finire in conclusioni non valide.

In questo articolo, voglio dare un esempio del paradosso di Simpson e mostrare come un’analisi semplice ma miope possa portare a conclusioni false che sembrano essere supportate dai dati, anche se non sono altro che una cattiva interpretazione. In questo modo, dimostro l’utilizzo di modelli di regressione multilivello come un modo appropriato di analisi dei dati per dati gerarchici, ovvero nidificati.

Il problema

Iniziamo! Immaginiamo di aver pubblicato un’app per smartphone e vogliamo conoscere meglio i nostri utenti e quanto sono soddisfatti. Pertanto, effettuiamo un piccolo sondaggio e chiediamo a alcuni dei nostri utenti di dirci quanto sono soddisfatti della nostra app su una scala da 1 (molto insoddisfatto) a 4 (molto soddisfatto). Inoltre, abbiamo misurato quanto tempo hanno trascorso nella nostra app nell’ultima settimana e, per ottenere un campione ricco, abbiamo chiesto agli utenti di diversi paesi. Quindi i nostri dati possono apparire così (utilizzo dati generati per questo articolo):

   Soddisfazione  Tempo_trascorso  Paese0      2.140440    1.585295        01      2.053545    0.636235        02      1.589258    1.468033        13      1.853545    0.968651        24      1.449286    0.967104        2.      .            .              ..      .            .              ..      .            .              .

Siamo interessati alla relazione tra il tempo trascorso nella nostra app e la soddisfazione riportata. Per essere precisi, vogliamo sapere se trascorrere più tempo nella nostra app è associato a una maggiore o minore soddisfazione e vogliamo quantificare questa associazione, ovvero vogliamo fare una dichiarazione del tipo “trascorrere un’ora in più nella nostra app è associato a essere x volte più/meno soddisfatti”. Quando osserviamo i dati, potremmo avere già una prima intuizione, ovvero che trascorrere più tempo nell’app è correlato a una minore soddisfazione:

Regressione lineare

Facciamo una regressione lineare per vedere se abbiamo ragione. Con la regressione lineare, cerchiamo di prevedere la soddisfazione data il tempo trascorso come una funzione lineare della forma soddisfazione = intercetta + coefficiente_di_regressione * tempo_trascorso. Possiamo farlo facilmente con il pacchetto statsmodels utilizzando la funzione OLS (Ordinary Least Squares).

import statsmodels.api as smresult = sm.OLS(df["Soddisfazione"], sm.add_constant(df["Tempo_trascorso"])).fit()print(result.params)

Il metodo add_constant è solo un dettaglio tecnico che usiamo per dire al modello che vogliamo avere un’intercetta nella nostra equazione (che è richiesta, fintanto che i nostri dati non sono standardizzati). Il risultato.params ci fornisce due valori, ovvero l’intercetta (const) e il coefficiente_di_regressione per la variabile Tempo_trascorso.

const         3.229412Tempo_trascorso   -0.655470

Quindi, il nostro modello ci dice che la soddisfazione può essere prevista come 3,229 – 0,655 * tempo_trascorso. In altre parole, trascorrere un’ora in più nell’app porta a una diminuzione di 0,655 punti nella soddisfazione (a causa del segno meno). Tuttavia, non si parte da zero, ma la soddisfazione media di una persona solo dalla loro prima impressione (cioè tempo_trascorso = 0) è 3,229. Possiamo rendere ciò visibile anche come una linea con intercetta 3,229 e pendenza -0,665:

Certo, questa previsione non è perfetta, ma almeno ci fornisce una tendenza. Ok, quindi il caso è chiaro, giusto? Trascorrere più tempo nell’app porta a una diminuzione della soddisfazione, e possiamo anche quantificare questa diminuzione. Potremmo trarre le nostre conclusioni da ciò e pensare a come migliorare l’applicazione (vogliamo che i nostri utenti siano più soddisfatti man mano che utilizzano l’app, ovviamente), o fare un sondaggio più dettagliato per scoprire perché gli utenti non sono soddisfatti.

Bene, non così veloce!

Raggruppamento per paese

Ricordate, abbiamo raccolto dati da utenti di diversi paesi? Cosa succede se guardiamo i dati separati per paese? Nel grafico seguente vediamo gli stessi dati di prima, ma ora abbiamo evidenziato ogni paese con un colore diverso.

Ci sono due osservazioni che possiamo fare da quel grafico. Primo, sembra che i paesi differiscano nella loro soddisfazione e nel tempo trascorso nell’app. I soggetti del paese blu trascorrono più tempo nell’app ma sono meno soddisfatti rispetto ai soggetti degli altri paesi, in media. Ancora di più, quando guardiamo i tre paesi separatamente, potremmo pensare che l’associazione tra il tempo trascorso nell’app e la soddisfazione sia effettivamente positiva. Non contraddice questa analisi precedente?

Paradosso di Simpson

Beh, in realtà non è chiamato così a causa di quei Simpson... Foto di Stefan Grage su Unsplash

L’effetto che abbiamo appena visto si chiama paradosso di Simpson. Si verifica quando una correlazione nei dati è diversa tra i gruppi rispetto all’interno dei gruppi. Sebbene sia molto controintuitivo, questo può accadere effettivamente (come abbiamo appena visto), e la ragione è rappresentata da variabili di confondimento. Spieghiamo questo con il nostro esempio di cui sopra. Osservando ogni paese isolatamente, vediamo una tendenza positiva: più tempo trascorso nell’app è associato a una maggiore soddisfazione. Tuttavia, come abbiamo già visto, i paesi differiscono nella loro soddisfazione media e nel tempo trascorso nell’app. Nel paese blu, la soddisfazione media è più bassa ma il tempo medio trascorso nell’app è più alto rispetto al paese arancione o verde; una tendenza che si oppone alla tendenza all’interno dei paesi. Tuttavia, potrebbe esserci un’altra variabile che causa questo. Ad esempio, si potrebbe immaginare che nel paese blu le persone siano più spesso annoiate, il che porta a una minore soddisfazione in generale (e quindi a un umore meno positivo nei confronti della nostra app) ma più tempo da trascorrere nell’app. Ovviamente, questa è solo una possibile spiegazione e ce ne possono essere molte altre. Tuttavia, la spiegazione corretta non è così importante al momento. Per noi è importante capire che ci sono differenze sistematiche tra i paesi.

Allora, perché non abbiamo scoperto questo nella nostra analisi precedente? Abbiamo commesso un errore nell’eseguire la regressione lineare? Beh, sì, perché era sbagliato fare una regressione lineare perché una delle assunzioni fondamentali della regressione lineare è stata violata: una regressione lineare assume che tutti i punti dati siano campionati in modo indipendente e dalla stessa distribuzione. Questo non è il caso nel nostro esempio! Ovviamente, le distribuzioni del tempo trascorso nell’app e della soddisfazione sono diverse nei diversi paesi. Ora, se le assunzioni della regressione lineare vengono violate, la regressione lineare non è lo strumento giusto per l’analisi dei dati.

Modelli gerarchici

Cosa possiamo fare ora per analizzare i nostri dati in modo più appropriato? Fortunatamente, sono disponibili modelli statistici che estendono l’idea della regressione lineare ai dati gerarchici. Parliamo di dati gerarchici se i punti dati che abbiamo campionato sono nidificati in una struttura gerarchica, come nel nostro caso, in cui le persone che abbiamo interrogato sono nidificate nei paesi. Tali modelli statistici sono chiamati modelli lineari gerarchici, modelli multilivello o modelli lineari a effetti misti. Tali modelli tengono conto delle strutture di gruppo introducendo cosiddetti effetti fissi ed effetti casuali. In un semplice esempio, in cui vogliamo prevedere una variabile data un’altra variabile singola (come vogliamo prevedere la soddisfazione data il tempo trascorso nell’app), gli effetti fissi consistono in un’intercetta e una pendenza per tutti i gruppi insieme. Finora è esattamente lo stesso della regressione lineare.

Ora gli effetti casuali possono introdurre una deviazione da quell’intercetta per ogni gruppo separatamente. Ad esempio, l’intercetta per il paese blu potrebbe essere leggermente più bassa e l’intercetta per il paese verde potrebbe essere leggermente più alta rispetto all’intercetta fissa. Ciò spiegherebbe le differenze nei livelli medi di soddisfazione dei paesi.

Inoltre, gli effetti casuali possono introdurre una deviazione della pendenza per ogni gruppo. Ad esempio, nel gruppo arancione, la pendenza potrebbe essere più alta rispetto alla pendenza fissa (cioè l’associazione tra soddisfazione e tempo trascorso è più forte), e nel paese verde potrebbe essere più bassa.

Modelli gerarchici in azione

Vediamo questo in azione per capire cosa succede realmente. Facciamo un’analisi nuova, ma questa volta utilizziamo la funzione mixedlm di statsmodels. Specifichiamo che vogliamo predire la soddisfazione data il tempo_trascorso (e non viceversa) mediante la formula “Soddisfazione ~ Tempo_trascorso” e indichiamo che la colonna “Paese” del nostro dataframe determina i diversi gruppi. Inoltre, il parametro re_formula=”Tempo_trascorso” indica al modello che vogliamo avere una pendenza separata per ogni gruppo. Senza di ciò, gli effetti casuali considererebbero solo un’intercetta specifica del gruppo, ma non una pendenza specifica del gruppo.

import statsmodels.formula.api as smfresult = smf.mixedlm("Soddisfazione ~ Tempo_trascorso", data=df, groups=df["Paese"], re_formula="Tempo_trascorso").fit()print(result.fe_params)print(result.random_effects)

Se stampiamo i fixed_effects (fe_params) e i random_effects, otteniamo valori come questi:

Effetti fissi  Intercetta     2.286638  Tempo_trascorso    0.497657Effetti casuali  {0: Gruppo -0.958805, Tempo_trascorso -0.018178,   1: Gruppo 0.155233,  Tempo_trascorso 0.274222,   2: Gruppo 0.803572,  Tempo_trascorso -0.256044}

Allora, cosa significa? Per gli effetti fissi, abbiamo un valore per l’intercetta e un valore per la nostra variabile tempo_trascorso. Per gli effetti casuali, invece, abbiamo due valori per paese (0,1,2): uno per l’intercetta (Gruppo) e uno per la pendenza della nostra variabile (Tempo_trascorso). Come abbiamo visto prima, gli effetti casuali descrivono la deviazione rispetto agli effetti medi per ogni gruppo. Per i nostri tre gruppi, possiamo costruire tre diverse equazioni lineari aggiungendo gli effetti casuali agli effetti fissi sia per l’intercetta che per la pendenza.

soddisfazione_0 = (2.286638 - 0.958805) + (0.497657 - 0.018178) * tempo_trascorso = 1.327833 + 0.479479 * tempo_trascorsosoddisfazione_1 = (2.286638 + 0.155233) + (0.497657 + 0.274222) * tempo_trascorso = 2.441871 + 0.771879 * tempo_trascorsosoddisfazione_2 = (2.286638 + 0.803572) + (0.497657 - 0.256044) * tempo_trascorso = 3.090210 + 0.241613 * tempo_trascorso

Vediamo che l’intercetta casuale per il gruppo 0 è negativa (-0.958) e l’intercetta casuale per il gruppo 2 è positiva (0.803), quindi il gruppo 0 si trova al di sotto dell’intercetta fissa e il gruppo 2 al di sopra. Di conseguenza, il gruppo 0 ha l’intercetta più bassa nella sua funzione lineare (1.327) e il gruppo 2 ha la più alta (3.090). In altre parole, nel paese 0, la soddisfazione inizia a un livello inferiore rispetto al paese 2.

Vediamo anche che le pendenze differiscono tra i gruppi. Nel gruppo 1, la pendenza è la più alta con 0.771, mentre per il gruppo 2 è solo 0.241. Ciò significa che l’associazione tra soddisfazione e tempo trascorso nell’app è molto più alta nel paese 1 rispetto al paese 2. In altre parole, nel paese 1 un aumento di un’ora di tempo trascorso nell’app porta a un incremento di 0.771 punti in soddisfazione (in media), mentre per il paese 2 aggiunge solo altri 0.241 punti in soddisfazione. Inoltre, tutte le pendenze sono positive, come ci aspettavamo dal grafico precedente, ma ciò contrasta con la pendenza negativa della regressione lineare che abbiamo fatto all’inizio.

Possiamo tracciare una linea di regressione per ogni paese ora:

Ora vediamo chiaramente il trend positivo in ogni paese e gli intercetti diversi (cioè le posizioni in cui le linee sarebbero al tempo_spent=0).

Conclusioni

Nell’esempio sopra, abbiamo visto come un’analisi miope possa facilmente portarci a conclusioni false. Ignorando la struttura annidata dei dati, cioè gli utenti provenienti da paesi diversi, avremmo potuto facilmente fermarci dopo la regressione lineare e avremmo concluso che più tempo trascorso nell’app era associato a una minore soddisfazione. Solo comprendendo che i nostri dati non soddisfano le assunzioni principali della regressione lineare, poiché i punti dati non sono estratti dalla stessa distribuzione, siamo stati motivati a fare ulteriori analisi che hanno rivelato che, in effetti, è il contrario: più tempo trascorso nell’app è effettivamente associato a una maggiore soddisfazione.

Quindi, traiamo alcune conclusioni da questo esempio:

  • Prima di utilizzare un metodo statistico per l’analisi, le sue assunzioni dovrebbero essere validate con i dati.
  • Per i dati annidati, l’assunzione che tutti i punti dati siano estratti dalla stessa distribuzione potrebbe non sempre essere vera.
  • Può accadere che un trend nei dati nel complesso sia diverso da un trend all’interno di singoli gruppi che compongono tali dati. Questo è chiamato paradosso di Simpson.
  • I modelli lineari multilivello sono un modo per gestire strutture dati annidate e evitare conclusioni false dal paradosso di Simpson.

Approfondimenti

Abbiamo utilizzato la seguente implementazione di modelli gerarchici in statsmodels:

  • https://www.statsmodels.org/stable/mixed_linear.html

Ho utilizzato il seguente libro di testo di statistica (disponibile solo in tedesco, sfortunatamente).

  • Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden.

Informazioni di base sui modelli multilivello possono essere trovate anche qui:

  • Snijders, T. A. B.; Bosker, R. J. (2011). Multilevel Analysis: an Introduction to Basic and Advanced Multilevel Modeling (2nd ed.). London: Sage. ISBN 9781446254332.

Se desideri riprodurre i risultati, ecco come sono stati generati i dati:

import numpy as npimport pandas as pdgroup_1_x = np.random.uniform(0.5, 1.8, 25)group_1_y = (1 + 0.3 * group_1_x) + np.random.rand(len(group_1_x))#start_2, end_2, step_2 = 0.3, 1.3, 0.04group_2_x = np.random.uniform(0.3, 1.3, 22)group_2_y = (2 + 0.7*group_2_x) + np.random.rand(len(group_2_x))#start_3, end_3, step_3 = 0, 1, 0.04group_3_x = np.random.uniform(0, 1, 32)group_3_y = (2.5 + 0.3*group_3_x) + np.random.rand(len(group_3_x))all_x = np.concatenate([group_1_x, group_2_x, group_3_x])all_y = np.concatenate([group_1_y, group_2_y, group_3_y])df = pd.DataFrame({"Soddisfazione": all_y, "Tempo_trascorso":all_x, "Paese":[0]*len(group_1_x) + [1]*len(group_2_x) + [2]*len(group_3_x)})

Ti è piaciuto questo articolo? Seguimi per essere notificato dei miei futuri post.