mercoledì 19 agosto 2009

Regressione logistica ordinale

In vari post precedenti (tra cui questi: [1^][2^][3^]), abbiamo discusso su come eseguire un modello di regressione logistica qualosa la variabile dipendente (Y) sia di tipo dicotomico. Se la variabile dipendente può invece assumente più di 2 valori (ossia la variabile dipendente è policotomica), si ricorre alla regressione logistica ordinale, in inglese ordered logistic regression. Esistono numerosi approcci a questo problema, e i principali modelli che sono stati formulati sono: cumulative logit o proportional odds model (Walker and Duncan 1967; McCullagh 1980), il continuation ratio model (Feinberg 1980), il constrained and unconstrained partial proportional odds models (Peterson and Harrell 1990), ed il adjacent-category logistic model (Agresti 1984). In questo articolo vengono esaminati alcuni di questi: LINK.
Eseguire in R una ordered logistic regression è alquanto difficoltoso; mancano difatti pacchetti dedicati esclusivamente all'analisi di dati categoriali (nominali). Inoltre, a quanto ne so, l'unico modello che è stato implementato è il proportional-odds model, ed è questo che utilizzerò negli esempi che ora vedremo.

Utilizzerò il dataset disponibile in questo pdf, ottima dispensa che ho trovato in rete sulla regressione logistica.
La seguente tabella di contingenza presenta i dati raccolti da un sondaggio universitario sottoposto ad alcuni studenti. L'obiettivo che ci proponiamo è quello di realizzare un modello che spieghi la variabile dipendente Y ordinale "giudizio complessivo), in base a "chiarezza" (variabile indipendente X1) e "interesse" (variabile indipendente X2).
Disponiamo quindi di:
- variabile dipendente policotomica
- 2 variabili indipendenti policotomiche



Per prima cosa dobbiamo cercare un modo per importare i dati di questa tabella in R. Io utilizzo il seguente schema: creo un file.txt sostituendo ai descrittori delle variabili ("decisamente no"..."decisamente si"), i numeri, da 1 a 4; si creano 3 colonne (X1, X2, Y) più una quarta colonna contenente le frequenze osservate. Il contenuto del file è il seguente:


chiarezza interesse giudizio frequenza
1 1 1 10
1 1 2 2
1 1 3 3
1 1 4 0
1 2 1 6
1 2 2 2
1 2 3 0
1 2 4 0
1 3 1 4
1 3 2 10
1 3 3 1
1 3 4 0
1 4 1 2
1 4 2 1
1 4 3 3
1 4 4 1
2 1 1 10
2 1 2 6
2 1 3 3
2 1 4 2
2 2 1 3
2 2 2 6
2 2 3 4
2 2 4 1
2 3 1 4
2 3 2 32
2 3 3 31
2 3 4 1
2 4 1 6
2 4 2 18
2 4 3 20
2 4 4 6
3 1 1 10
3 1 2 10
3 1 3 21
3 1 4 2
3 2 1 1
3 2 2 27
3 2 3 73
3 2 4 4
3 3 1 1
3 3 2 41
3 3 3 375
3 3 4 47
3 4 1 3
3 4 2 12
3 4 3 195
3 4 4 138
4 1 1 2
4 1 2 2
4 1 3 19
4 1 4 10
4 2 1 1
4 2 2 12
4 2 3 48
4 2 4 26
4 3 1 0
4 3 2 14
4 3 3 288
4 3 4 240
4 4 1 1
4 4 2 7
4 4 3 182
4 4 4 1224


Ogni elemento della riga deve essere spaziato (potremmo anche utilizzare la virgola). Salvato il file dataset.txt, importiamo i dati in R con il seguente comando:


dati <- read.table("/datasetR/dataset.txt", header=T)


All'interno di dati si trova ora la tabella, ripresa dal file txt di cui specifichiamo l'indirizzo (attenzione ad usare lo slash /, anzichè il backslash \, altrimenti R comunica un errore). Specificando header=TRUE, informiamo R che la prima riga del file da leggere è costituita dal nome delle colonne.
In che formato sono le colonne della tabella? Verifichiamo con il seguente comando:


sapply(dati, class)

chiarezza interesse giudizio frequenza
integer" "integer" "integer" "integer"


Le 4 colonne sono interpretate come numeriche. Eppure solo la colonna frequenza è di tipo numerico; le altre devono essere trasformate in fattori. Per far questo utilizziamo la funcione factor(); possiamo anche specificare il nome di ciascun elemento dei fattori (non è obbligatorio). Procediamo in questo modo e osserviamo la trasformazione, e verifichiamo con la funzione sapply():


dati$giudizio <- factor(dati$giudizio, labels=c("fantastico","buono","discreto","pessimo"))
dati$chiarezza <- factor(dati$chiarezza, labels=c("ottima","buona","discreta","pessima"))
dati$interesse <- factor(dati$interesse, labels=c("moltissimo","molto","poco","pochissimo"))

> dati
chiarezza interesse giudizio frequenza
1 ottima moltissimo fantastico 10
2 ottima moltissimo buono 2
3 ottima moltissimo discreto 3
4 ottima moltissimo pessimo 0
5 ottima molto fantastico 6
6 ottima molto buono 2
7 ottima molto discreto 0
8 ottima molto pessimo 0
9 ottima poco fantastico 4
10 ottima poco buono 10
11 ottima poco discreto 1
12 ottima poco pessimo 0
13 ottima pochissimo fantastico 2
14 ottima pochissimo buono 1
15 ottima pochissimo discreto 3
16 ottima pochissimo pessimo 1
17 buona moltissimo fantastico 10
18 buona moltissimo buono 6
19 buona moltissimo discreto 3
20 buona moltissimo pessimo 2
21 buona molto fantastico 3
22 buona molto buono 6
23 buona molto discreto 4
24 buona molto pessimo 1
25 buona poco fantastico 4
26 buona poco buono 32
27 buona poco discreto 31
28 buona poco pessimo 1
29 buona pochissimo fantastico 6
30 buona pochissimo buono 18
31 buona pochissimo discreto 20
32 buona pochissimo pessimo 6
33 discreta moltissimo fantastico 10
34 discreta moltissimo buono 10
35 discreta moltissimo discreto 21
36 discreta moltissimo pessimo 2
37 discreta molto fantastico 1
38 discreta molto buono 27
39 discreta molto discreto 73
40 discreta molto pessimo 4
41 discreta poco fantastico 1
42 discreta poco buono 41
43 discreta poco discreto 375
44 discreta poco pessimo 47
45 discreta pochissimo fantastico 3
46 discreta pochissimo buono 12
47 discreta pochissimo discreto 195
48 discreta pochissimo pessimo 138
49 pessima moltissimo fantastico 2
50 pessima moltissimo buono 2
51 pessima moltissimo discreto 19
52 pessima moltissimo pessimo 10
53 pessima molto fantastico 1
54 pessima molto buono 12
55 pessima molto discreto 48
56 pessima molto pessimo 26
57 pessima poco fantastico 0
58 pessima poco buono 14
59 pessima poco discreto 288
60 pessima poco pessimo 240
61 pessima pochissimo fantastico 1
62 pessima pochissimo buono 7
63 pessima pochissimo discreto 182
64 pessima pochissimo pessimo 1224


> sapply(dati, class)
chiarezza interesse giudizio frequenza
"factor" "factor" "factor" "integer"


Ora il formato dei dati è quello opportuno per proseguire. In realtà, se il nostro file txt contenesse direttamente il nome, e non il numero, per ogni categoria di ogni variabile, R avrebbe interpretato direttamente come un fattore le varie colonne.

Per eseguire una regressione (logistica) ordinale con il proportional-odds model, installiamo e carichiamo la libreria MASS, e utilizziamo la funzione polr():


library(MASS)
model <- polr(giudizio ~ chiarezza + interesse, data=dati, weights=frequenza, Hess=T)
summary(model)

Call:
polr(formula = giudizio ~ chiarezza + interesse, data = dati,
weights = frequenza, Hess = T)

Coefficients:
Value Std. Error t value
chiarezzabuona 1.1720 0.3509 3.340
chiarezzadiscreta 3.3901 0.3340 10.151
chiarezzapessima 5.3735 0.3404 15.784
interessemolto 0.5239 0.2546 2.058
interessepoco 1.4395 0.2240 6.426
interessepochissimo 3.1973 0.2296 13.924

Intercepts:
Value Std. Error t value
fantastico|buono 0.873 0.345 2.528
buono|discreto 3.041 0.369 8.229
discreto|pessimo 6.862 0.390 17.576

Residual Deviance: 4277.66
AIC: 4295.66


Osserviamo l'input di polr(): viene prima dichiarata la formula da analizzare (il giudizio in funzione di chiarezza e interesse), quindi si specifica la sorgente dei dati (data=dati), le frequenza (weights), e con la dicitura Hess=T abilitiamo la funzione sotto-utilizzata summary() sul modello (se non specifichiamo Hess=TRUE, ci viene restituito errore).

In output invece vengono forniti i valori dei parametri beta [coefficients] (relativi alle due variabili indipendenti) e dei parametri tau [intercepts] (rispettivi delle 4 modalità della variabile dipendente). Per quanto riguarda l'interpretazione dei coefficienti, rimando alla già citata dispensa sulla regressione logistica.

Concentrandoci invece sul programma R, sono disponibili altre funzioni per la ordered logit:

- Nella library Design è disponibile la funzione lrm().
- Nella library VGAM è disponibile la funzione vglm().
- Nella library Zelig è disponibile la funzione zelig().
- Nella library nnet è disponibile la funzione multinom().

Per ognuna di queste funzioni è disponibile l'help relativo (basta caricare la library e quindi digitare help(nome_funzione)). Tutte dovrebbero eseguire delle regressioni logistiche ordinali, ma gli output di tutte le funzioni, sugli stessi dati, sono completamente differenti. L'eccessiva confusione che si è creata (numerosi modelli + diverse modalità di input dati + diversi output), rendono l'argomento alquanto nebuloso.

4 commenti:

  1. Ciao,
    complimenti per gli argomenti trattati, preciso e chiaro.
    Posso permetterti di chiederti una cosa?
    Come tratteresti il problema della multicollinearità in una regressione logitica?

    RispondiElimina
  2. Sono molto contento di avere la possibilità di trovare informazioni utili visitare il vostro blog. Continuate a fare il lavoro impressionante.

    RispondiElimina
  3. ciao, potresti postare il pdf in quanto il link non funziona?
    grazie mille

    RispondiElimina
  4. Ciao, purtroppo il link per la dispensa non funziona più

    RispondiElimina