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.