martedì 21 luglio 2009

Regressione logistica semplice su variabili qualitative dicotomiche

In questo post vedremo brevemente come realizzare un modello di regressione logistica qualora si disponga di variabili categoriali, o qualitative, organizzate in tabelle di contingenza a doppia entrata. In questo modello la variabile dipendente Y è una variabile bernoulliana, che può assumere i valori 0 oppure 1. La probabilità che assuma valore 1 dipende dai regressori, ossia dalle variabili indipendenti che si prendono in considerazione; anch'esse sono variabili dicotomiche.

In generale, la probabilità che sia Y = 1 in funzione dei regressori è:



Il nostro obiettivo è quello di stimare i valori dei parametri beta.

Cominciamo ad esaminare un modello di regressione logistica semplice (un solo regressore). Sia la variabile dipendente che quella indipendente può assumere valori 0 o 1.

Si consideri il seguente esempio. La tabella sottostante riporta i risultati di uno studio sul reflusso gastroesofageo. Si vuole valutare quanto la presenza di una fattore di stress possa influenzare l'insorgenza di questa patologia.



Innanzitutto importiamo i valori in R. Dobbiamo creare una tabella a doppia entrata; procediamo in questo modo:


reflusso <- matrix(c(251,131,4,33), nrow=2)
colnames(reflusso) <- c("reflNO", "reflSI")
rownames(reflusso) <- c("stressNO", "stressSI")
table <- as.table(reflusso)


Osserviamo la tabella richiamando la variabile:


table

reflNO reflSI
stressNO 251 4
stressSI 131 33


Ora aggiustiamo i dati per effettuare la regressione logistica. Dobbiamo creare un data frame:


dft <- as.data.frame(table)
dft

Var1 Var2 Freq
1 stressNO reflNO 251
2 stressSI reflNO 131
3 stressNO reflSI 4
4 stressSI reflSI 33


Possiamo ore procedere fittando i dati, e quindi effetuando la regressione logistica:


fit <- glm(Var2 ~ Var1, weights = Freq, data = dft, family = binomial(logit))
summary(fit)

Call:
glm(formula = Var2 ~ Var1, family = binomial(logit), data = dft,
weights = Freq)

Deviance Residuals:
1 2 3 4
-2.817 -7.672 5.765 10.287

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.1392 0.5040 -8.213 < 2e-16 ***
Var1stressSI 2.7605 0.5403 5.109 3.23e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 250.23 on 3 degrees of freedom
Residual deviance: 205.86 on 2 degrees of freedom
AIC: 209.86

Number of Fisher Scoring iterations: 6



Questo che vediamo è l'output della funzione. Anzitutto commentiamo il codice per effettuare la regressione. La regressione logistica viene chiamata imponendo la famiglia: family = binomial(logit). Il codice Var2 ~ Var1 significa che vogliamo creare un modello che ci spieghi la variabile Var2 (presenza o assenza del reflusso) in funzione della variabile Var1 (presenza o assenza di eventi stressogeni). In pratica Var2 è la variabile indipendente Y e Var1 è la variabile dipendente X (il regressore). Fornita la formula da analizzare, si specifica il peso di ciascuna variabile, dati contenuti nella colonna Freq del dataframe dft (quindi si scrive weights = Freq e data = dft per specificare la posizione dove sono contenuti i valori). Con summary(fit) si ottengono una serie di dati utili.

I valori dei parametri $\beta_0$ e $\beta_1$ sono rispettivamente i valori (intercept) e Var1stress1. Possiamo quindi scrivere il nostro modello empirico:



La variabile indipendente x può assumere valore zero oppure uno. Se assume valore 0 (quindi in assenza di eventi stressanti), allora la probabilità di avere il reflusso è pari a:



Se sono presenti eventi stressanti (x=1), la probabilità di avere reflusso è:



Le formula per gli odds sono:



Possiamo infine calcolare l'odd ratio OR:



Una persona che ha vissuto un evento stressante ha una propensione a sviluppare il reflusso gastroesofageo 15.807 volte più grande della persona che non ha subito eventi stressanti.

Gli odds, le probabilità, e l'OR, possono essere rapidamente calcolati in R ricordando che:

fit$coefficient[1] = (intercept)
fit$coefficient[2] =

Inoltre:

summary(fit)$coefficient[1,2] = standard error del coefficiente $\beta_0$
summary(fit)$coefficient[2,2] = standard error del coefficiente $\beta_1$

E quindi abbiamo:


pi0 <- exp(fit$coefficient[1]) / (1 + exp(fit$coefficient[1]))
pi1 <- exp(fit$coefficient[1] + fit$coefficient[2]) / (1 + exp(fit$coefficient[1]+fit$coefficient[2]))

odds0 <- pi0 / (1 - pi0)
odds1 <- pi1 / (1 - pi1)

OR <- odds1 / odds0

#stesso risultato con:
OR <- exp(fit$coefficient[2])

#l'intervallo di confidenza per l'odd-ratio è:
ORmin <- exp( fit$coefficient[2] - qnorm(.975) * summary(fit)$coefficient[2,2] )

ORmax <- exp( fit$coefficient[2] + qnorm(.975) * summary(fit)$coefficient[2,2] )


La più nota formula per l'odd-ratio è la seguente (stesso risultato):



che in R calcoliamo così:


OR <- (table[1,1]*table[2,2]) / (table[1,2]*table[2,1])


L'acronimo AIC sta per Akaike's information criterion. Questo parametro non ci fornisce alcun dato sul modello appena creato. Esso può però essere utile nel confrontare questo modello con altri eventualmente presi in considerazione (il modello con AIC più basso è il migliore).

La regressione logistica multipla viene trattata in questo post

5 commenti:

  1. Ciao! questo blog è molto chiaro e utile, ti ringrazio. Vorrei chiederti se sai la differenza tra queste due funzioni:
    glm(Var2 ~ Var1,family = binomial(link =logit))
    e
    glm(Var2 ~ Var1,family = binomial(link =probit))

    ti ringrazio!
    Alessandra

    RispondiElimina
  2. Ciao Alessandra,
    ti ricambio i complimenti, in quanto mi hai fatto proprio una bella domanda :)
    Tra logit e probit le differenze teoriche sono molte, quelle pratiche un pò meno. Matematicamente parlando la differenza sta nella distribuzione degli errori, logistica (nel caso della logit) e normale (nel caso della probit). Dal punto di vista pratico questo si ripercuote sulla pendenza della curva sigmoidea, e quindi avremo diversi parametri stimati della stessa funzione.

    Come puoi immaginare, questo può determinare grandi cambiamenti, qualora si voglia predire un certo parametro. Anche se, nella pratica, se provi a plottare con gli stessi dati una logit e una probit, le due sigmoidi saranno ""quasi"" sovrapponibili (dipende dai casi).

    Per quanto riguarda la scelta, ti posso dire che in linea di massima per una multivariata è preferibile la logit, mentre per una bivariata è meglio una probit. Ovviamente ci sono le eccezioni (e di ciò te ne accorgi validando il modello e le assunzioni).
    Se stai esaminando una serie di dati molto particolare, ti conviene comunque andare a vedere i modelli già utilizzati in altri articoli scientifici, in quanto ogni distribuzione può avere delle sue peculiarità che giustificano l'uso dell'uno o dell'altro metodo.

    Ti segnalo questo pdf, per approfondire un poco le differenze tra logit e probit :)

    http://www.iasri.res.in/ebook/EBADAT/6-Other%20Useful%20Techniques/5-Logit%20and%20Probit%20Analysis%20Lecture.pdf

    RispondiElimina
  3. Ciao,
    non mi è chiaro se c'è un errore nel calcolo di pi1.
    Se ho capito bene quando calcoliamo pi1 dobbiamo moltiplicare fit$coefficient[2] * Var1, corretto?

    Grazie.
    Massimo

    RispondiElimina
  4. Nel caso pi1 Var1=1.
    Ci sono arrivato dopo aver postato il commento.
    Scusa.

    Ciao
    Massimo

    RispondiElimina
  5. Ciao
    ho una variabile che varia da -1 a +1. Posso fare una regressione logistica per vedere la sua variazione nel tempo?
    Se sì, si può utilizzare un approccio di maximum likelihood?
    grazie
    Daniele

    RispondiElimina