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