Analisi della sopravvivenza: curve di Kaplan-Meier e modello di Cox

Le funzioni di sopravvivenza possono essere studiate con due metodi: il metodo attuariale e il metodo di Kaplan-Meier. Quest'ultimo viene utilizzato quando:
- si conosce il tempo esatto in cui avviene l'evento considerato (e non l'intervallo, come nel metodo atturiale)
- i persi al follow-up (truncated o censored, in inglese) partecipano alla determinazione della funzione di sopravvivenza fino al loro ritiro, perchè considerati a rischio fino a quel momento
Così la curva di Kaplan-Meier (e quindi la probabilità di non-evento) cambia ogni volta che accade l'evento nel gruppo considerato.

Supponiamo di voler analizzare la curva di sopravvivenza di 10 soggetti, con il metodo di Kaplan-Meier. Supponiamo che questo gruppo sia affetto da una certa malattia, e siamo quindi interessati a stimare la probabilità che un soggetto ha di sopravvivere fino ad un certo tempo.
I dati vengono raccolti annualmente, e sono i seguenti (per assurdo ammettiamo che gli eventi trocato o deceduto avvengano esattamente alla fine di ogni anno; in realtà si dovrebbe considerare il momento esatto, come si fa nell'esempio successivo):

  Tempo   Numero.Morti   Numero.Troncati
     1            0               1
     3            1               0
     4            0               1
     5            2               0
     6            0               1
     7            2               1
     8            1               0

Una tabella di questo genere viene letta nel modo seguente: alla fine del primo anno non è morto nessun paziente, ed è stato perso al follow-up 1 paziente; al termine del terzo anno è morto 1 paziente; e così via.

In R il pacchetto per l'analisi delle funzione di sopravvivenza è la library survival. Il formato con cui sono stati presentati i dati, non è però quello ottimale per l'analisi in R. La funzione che tra breve useremo richiede in input questi stessi dati, presentati in questo modo:

   time status
     1      0
     3      1
     4      0
     5      1
     5      1
     6      0
     7      1
     7      1
     7      0
     8      1

Per ogni paziente in osservazione, viene definito uno status:
0 per i troncati
1 per i deceduti
Questa tabella viene letta nel modo seguente: al termine del primo anno un paziente è stato perso al follow-up; se in un certo tempo, più di un paziente viene perso o è deceduto, si scrive una riga per ognuno di essi. Quindi osservando il settimo anno, si scrivono 3 righe, due relative ai 2 pazienti deceduti, e una relativa al paziente troncato
In pratica ogni volta che accade un certo evento (paziente troncato o deceduto), si scrive una riga per paziente, specificando il tempo in cui è avvenuto e cosa è avvenuto.

Creiamo quindi un dataframe in R con i seguenti dati:

time <- c(1,3,4,5,5,6,7,7,7,8)
status <- c(0,1,0,1,1,0,1,1,0,1)
popol <- data.frame(time, status)

A questo punto possiamo generare la curva di Kaplan-Meier:

library(survival)
fit <- survfit(Surv(time, status) ~ 1, data=popol)
plot(fit)
text( fit$time, 0, format(fit$n.risk), cex = 0.7 )

Osserviamo il grafico. Ciò che interessa a noi è la linea continua, che indica la probabilità di sopravvivenza in funzione del tempo; le linee discontinue fanno riferimento all'intervallo di confidenza della probabilità di sopravvivenza; se non vogliamo l'intervallo di confidenza, scriviamo plot(fit, conf.int=F). I trattini che possiamo vedere, in corrispondenza dei tempi 1,4,6 indicano i soggetti che sono stati censurati (persi). I numeri posizionati col comando text sull'asse X indicano il numero di soggetti in osservazione ai vari tempi.
Con la funzione summary(fit) otteniamo le probabilità in corrispondenza dei punti dove la curva cambia:

 fit <- survfit(Surv(time, status) ~ 1, data=popol, se.fit=TRUE, conf.int=.95, type="kaplan-meier")
 summary(fit)

Call: survfit(formula = Surv(time, status) ~ 1, data = popol, se.fit = TRUE, 
    conf.int = 0.95)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    3      9       1    0.889   0.105        0.706        1.000
    5      7       2    0.635   0.169        0.377        1.000
    7      4       2    0.317   0.180        0.105        0.964
    8      1       1    0.000     NaN           NA           NA

La colonna survival della tabella indica la probabilità di sopravvivenza (l'asse Y del grafico di prima). Ad esempio, deduciamo da questa che la probabilità (cumulativa, e non puntuale) di sopravvivere al termine del 5° anno è del 63.5%. Essendo il gruppo A costituito da 10 soggetti, che al termine dell'ottavo anno sono tutti scomparsi (persi oppure deceduti), risulta ovvio che la probabilità di sopravvivere al termine dell'ottavo anno è pari a zero.

---

Come ulteriore esercizio delle funzioni viste fin qui, utilizziamo i dati estratti dalla tabella pubblicata a pagina 2 di questo pdf.
Riassiumiamo rapidamente i comandi per generare la Kaplan-Meier survival curve, e per ottenere l'ultima colonna (Cumulative proportion surviving).

time <- c(0.909, 1.112, 1.322, 1.328, 1.536, 2.713, 2.741, 2.743, 3.524, 4.079)
status <- c(1,1,0,1,1,1,0,1,0,0)
table1 <- data.frame(time, status)

fit <- survfit(Surv(time, status) ~ 1, data=table1)

summary(fit)

Call: survfit(formula = Surv(time, status) ~ 1, data = table1)

  time n.risk n.event survival std.err lower 95% CI upper 95% CI
 0.909     10       1    0.900  0.0949        0.732        1.000
 1.112      9       1    0.800  0.1265        0.587        1.000
 1.328      7       1    0.686  0.1515        0.445        1.000
 1.536      6       1    0.571  0.1638        0.326        1.000
 2.713      5       1    0.457  0.1662        0.224        0.932
 2.743      3       1    0.305  0.1666        0.104        0.890


plot(fit, col="red", lwd=2, conf.int=F)

---

Vediamo adesso come confrontare due curve di sopravvivenza di Kaplan-Meier utilizzando il Log-rank test (o test di Mantel-Cox).
Si dispone di due gruppi (A e B), e dei seguenti dati:

Gruppo A
  Tempo Numero.Morti Numero.Troncati
     1            0               1
     3            1               0
     4            0               1
     5            2               0
     6            0               1
     7            2               1
     8            1               0

Gruppo B
  Tempo Numero.Morti Numero.Troncati
     2            2               0
     3            0               1
     4            1               0
     6            0               2
     7            1               0
    10            1               0

Come abbiamo visto prima, questa forma non è quella idonea per le funzioni di R. Leggiamo le due tabelle in questo modo:

Gruppo A
  1 soggetto troncato al tempo 1
  1 soggetto morto al tempo 3
  1 soggetto troncato al tempo 4
  2 soggetti morti al tempo 5
  1 soggetto troncato al tempo 6
  2 soggetti morti e 1 soggetto troncato al tempo 7
  1 soggetto morto al tempo 8

Gruppo B
  2 soggetti morti al tempo 2
  1 soggetto troncato al tempo 3
  1 soggetto morto al tempo 4
  2 soggetti troncati al tempo 6
  1 soggetto morto al tempo 7
  1 soggetto morto al tempo 10

Traduciamo quindi le tabelle nel seguente formato:

Gruppo A
   time status
     1      0
     3      1
     4      0
     5      1
     5      1
     6      0
     7      1
     7      1
     7      0
     8      1

Gruppo B
  time status
    2      1
    2      1
    3      0
    4      1
    6      0
    6      0
    7      1
   10      1

Adesso creiamo un dataframe contenente questi valori, più la variabile che specifica il gruppo di appartenenza (che trasformiamo in un fattore):

time <- c(1,3,4,5,5,6,7,7,7,8,2,2,3,4,6,6,7,10)
status <- c(0,1,0,1,1,0,1,1,0,1,1,1,0,1,0,0,1,1)
groups <- c(rep("GroupA",10), rep("GroupB",8))
pop <- data.frame(time, status, groups)
pop$groups <- factor(pop$groups)

Adesso disegnano su un unico grafico le due curve di Kaplan-Meier. Dapprima fittiamo il modello: anzichè specificare ~1, il modello deve essere in funzione del gruppo di appartenenza, e quindi abbiamo ~groups:

 fit <- survfit(Surv(time, status) ~ groups, data=pop)
 summary(fit)

Call: survfit(formula = Surv(time, status) ~ groups, data = pop)

                groups=GroupA 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    3      9       1    0.889   0.105        0.706        1.000
    5      7       2    0.635   0.169        0.377        1.000
    7      4       2    0.317   0.180        0.105        0.964
    8      1       1    0.000     NaN           NA           NA

                groups=GroupB 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    2      8       2     0.75   0.153       0.5027            1
    4      5       1     0.60   0.182       0.3315            1
    7      2       1     0.30   0.231       0.0664            1
   10      1       1     0.00     NaN           NA           NA

 plot(fit, col=2:3, lwd=2)

Ora confrontiamo le due curve con il Log-rank test:

survdiff(Surv(time, status) ~ groups, data=pop)

Call:
survdiff(formula = Surv(time, status) ~ groups, data = pop)

               N Observed Expected (O-E)^2/E (O-E)^2/V
groups=GroupA 10        6     6.05  0.000364   0.00110
groups=GroupB  8        5     4.95  0.000445   0.00110

 Chisq= 0  on 1 degrees of freedom, p= 0.974 

Commentiamo l'output: la colonna N indica il numero di soggetti analizzati per ciascun gruppo; la colonna observed è il numero di morti osservate (O); la colonna expected è il numero di morti attese per ciascun gruppo (E). Infine vengono riportati i valori per il calcolo del chi-quadro della statistica test.
Il chi-quadro calcolato è pari alla somma dei valori della colonna (O-E)^2/E; esso ci viene dato anche alla fine, e ha valore 0 (approssimato, in realtà è 0.000809). Il relativo p-value (0.974) è abbondantemente superiore a 0.05, quindi accettiamo l'ipotesi nulla: l'andamento della sopravvivenza nei due gruppi è significativamente uguale.

Per calcolare il relative risk del gruppoA rispetto al gruppoB, si applica questa formula:

$$RR=\frac{O_1/E_1}{O_2/E_2}$$

Calcoliamo in R utilizzando i valori della tabella ottenuta:

RR <- (6/6.05) / (5/4.95)
RR
[1] 0.9818182

Avendo ottenuto un Log-rank test non significativo, è abbastanza prevedibile ottenere anche un relative-risk molto vicino ad 1 (il rischio di morte del gruppo A è pressochè uguale a quello del gruppo B).

---

Vediamo ora un altro esempio di confronto di due curve KM. Utilizzo qui i dati presenti in questo link. Vediamo subito che uno dei modi di presentare i dati è il seguente:

Group 1: 143, 165, 188, 188, 190, 192, 206, 208, 212, 216, 220, 227, 230, 235, 246, 265, 303, 216*, 244*
Group 2: 142, 157, 163, 198, 205, 232, 232, 232, 233, 233, 233, 233, 239, 240, 261, 280, 280, 295, 295, 323, 204*, 344*

Per ogni individuo in analisi, viene fornito il tempo di sopravvivenza prima dell'evento, mentre vengono asteriscati i soggetti persi al follow-up (censored). Traduciamo questi dati in R:

Time1 <- c(143, 165, 188, 188, 190, 192, 206, 208, 212, 216, 220, 227, 230, 235, 246, 265, 303, 216, 244)
Time2 <- c(142, 157, 163, 198, 205, 232, 232, 232, 233, 233, 233, 233, 239, 240, 261, 280, 280, 295, 295, 323, 204, 344)

Status1 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0)
Status2 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0)

Gru1 <- rep("Group1", length(Status1))
Gru2 <- rep("Group2", length(Status2))

Time <- c(Time1,Time2)
Status <- c(Status1,Status2)
Group <- factor(c(Gru1,Gru2))

Osserviamo il formato della funzione Surv:

> Surv(Time1, Status1)
 [1] 143  165  188  188  190  192  206  208  212  216  220  227  230  235  246  265 
[17] 303  216+ 244+

> Surv(Time2, Status2)
 [1] 142  157  163  198  205  232  232  232  233  233  233  233  239  240  261  280 
[17] 280  295  295  323  204+ 344+

Come possiamo vedere, la funzione Surv trasforma i dati nel formato descritto prima. Adesso fittiamo il modello:

KM.fit <- survfit(Surv(Time,Status) ~ Group)
summary(KM.fit)

Call: survfit(formula = Surv(Time, Status) ~ Group)

                Group=Group1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  143     19       1    0.947  0.0512       0.8521        1.000
  165     18       1    0.895  0.0704       0.7669        1.000
  188     17       2    0.789  0.0935       0.6259        0.996
  190     15       1    0.737  0.1010       0.5632        0.964
  192     14       1    0.684  0.1066       0.5041        0.929
  206     13       1    0.632  0.1107       0.4480        0.890
  208     12       1    0.579  0.1133       0.3946        0.850
  212     11       1    0.526  0.1145       0.3435        0.806
  216     10       1    0.474  0.1145       0.2949        0.761
  220      8       1    0.414  0.1145       0.2412        0.712
  227      7       1    0.355  0.1124       0.1911        0.661
  230      6       1    0.296  0.1082       0.1447        0.606
  235      5       1    0.237  0.1015       0.1023        0.548
  246      3       1    0.158  0.0934       0.0495        0.504
  265      2       1    0.079  0.0728       0.0130        0.481
  303      1       1    0.000     NaN           NA           NA

                Group=Group2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  142     22       1   0.9545  0.0444      0.87136        1.000
  157     21       1   0.9091  0.0613      0.79656        1.000
  163     20       1   0.8636  0.0732      0.73151        1.000
  198     19       1   0.8182  0.0822      0.67189        0.996
  205     17       1   0.7701  0.0904      0.61180        0.969
  232     16       3   0.6257  0.1051      0.45020        0.870
  233     13       4   0.4332  0.1082      0.26548        0.707
  239      9       1   0.3850  0.1063      0.22408        0.662
  240      8       1   0.3369  0.1034      0.18465        0.615
  261      7       1   0.2888  0.0992      0.14731        0.566
  280      6       2   0.1925  0.0864      0.07991        0.464
  295      4       2   0.0963  0.0647      0.02580        0.359
  323      2       1   0.0481  0.0469      0.00712        0.326

Se siamo interessati a conoscere la mediana del tempo di sopravvivenza per ciascun gruppo, e il suo intervallo di confidenza, scriviamo:

KM.fit

Call: survfit(formula = Surv(Time, Status) ~ Group)

             records n.max n.start events median 0.95LCL 0.95UCL
Group=Group1      19    19      19     17    216     206     265
Group=Group2      22    22      22     20    233     232     280

Per il gruppo1 la mediana è 216, con intervallo di confidenza (206,265), mentre per il gruppo2 è 233, con CI = (232,280). L'algoritmo per calcolare l'intervallo di confidenza della mediana può essere scelto con il comando conf.type="" (nell'help della funzione survfit sono descritti i vari algoritmi che possiamo scegliere).
In R non ho trovato funzioni per calcolare la media; ma del resto la media nell'analisi delle curve di sopravvivenza di Kaplan-Meier non viene accettata da buona parte degli statistici.

Adesso poniamo in grafico le curve di Kaplan-Meier:

plot(KM.fit, conf.int=F, col=2:3, lwd=2, xlab="time", ylab="Survival probability")
mtext("K-M survival curve for Group1 and Group2", 3, line=2, cex=1.5)
abline(h=0)

text(300, 1, expression(bold("Group 1 - median = 216")), col=2, cex=.8, lwd=2)
text(300, 0.95, expression(bold("Group 2 - median = 233")), col=3, cex=.8, lwd=2)

Adesso eseguiamo il Log-rank test sui due gruppi:

survdiff(Surv(Time, Status) ~ Group)

Call:
survdiff(formula = Surv(Time, Status) ~ Group)

              N Observed Expected (O-E)^2/E (O-E)^2/V
Group=Group1 19       17     12.2     1.885      3.12
Group=Group2 22       20     24.8     0.928      3.12

 Chisq= 3.1  on 1 degrees of freedom, p= 0.0771

Abbiamo ottenuto p-value > 0.05, quindi l'andamento della sopravvivenza nei due gruppi è significativamente (95%) uguale. Il rischio relativo del gruppo1 rispetto al gruppo2 è:

RR <- (17/12.2) / (20/24.8)
RR
[1] 1.727869

Il gruppo1 ha un rischio di morte 1.73 volte più elevato del gruppo2.
Se vogliamo esaminare la differenza delle probabilità di sopravvivenza calcolata per ogni tempo esaminato, e il suo intervallo in confidenza, utilizziamo la funzione survdiffplot disponibile nel pacchetto Design:

library(Design)
survdiffplot(KM.fit)

Per ottenere il grafico dell'hazard in funzione del tempo, dobbiamo utilizzare le funzioni coxph() e quindi basehaz() e plottare il risultato:

temp <- coxph(Surv(Time, Status)~strata(Group))
plot(basehaz(temp))

---

Proportional hazard model: analisi della sopravvivenza col modello di Cox

Il Log-rank test che abbiamo visto precedentemente permette di confrontare l'andamento della sopravvivenza in più gruppi, che potrebbero anche essere suddivisi in base ad un certo predittore (ad esempio si potrebbero confrontare le curve di sopravvivenza nei maschi e nelle femmine). Esso però non ci fornisce alcuna idea sul peso che assume la regola con la quale dividiamo i gruppi (ossia non ci dice, ad esempio, quanto influisce la variabile "sesso" sulla sopravvivenza).
Il modello di Cox è il modello di regressione sviluppato proprio per poter analizzare il peso di certe variabili sulla sopravvivenza. Esso altro non è che una regressione multipla, di cui segue gli stessi principi.
In presenza di un certo numero di fattori di rischio (probabili), chiamate covariate, è possibile stimare il rischio di morte (hazard) ad un certo tempo t, in presenza delle covariate z. La funzione di rischio (che rappresenta la probabilità di morire al tempo t per un individuo con vettore di covariate z, sopravvissuto fino al tempo t) viene generalmente fattorizzata in due parti:
$$\h(t,z)=\lambda_0(t)\cdot H(z,\beta)$$

Qui $\lambda_0(t)$ è la funzione di rischio di base (baseline hazard) per un individuo con vettore di covariate z=0 che dipende solo da t ed è lasciata non specificata; H(z,\beta) è la funzione hazard, con beta i regressori (il peso di ciascuna covariate z).
Nel modello più semplice, la funzione hazard H è quella esponenziale:
$$H(z,\beta)=exp(\beta_1x_1+\beta_2x_2+...)$$

Il modello di Cox è definito semi-parametrico. Ma cominciamo subito con un esempio.

Carichiamo il pacchetto survival, e utilizziamo il dataset ovarian, ivi disponibile:

library(survival)
Carico il pacchetto richiesto: splines

ovarian
   futime fustat     age resid.ds rx ecog.ps
1      59      1 72.3315        2  1       1
2     115      1 74.4932        2  1       1
3     156      1 66.4658        2  1       2
4     421      0 53.3644        2  2       1
5     431      1 50.3397        2  1       1
6     448      0 56.4301        1  1       2
7     464      1 56.9370        2  2       2
8     475      1 59.8548        2  2       2
9     477      0 64.1753        2  1       1
10    563      1 55.1781        1  2       2
11    638      1 56.7562        1  1       2
12    744      0 50.1096        1  2       1
13    769      0 59.6301        2  2       2
14    770      0 57.0521        2  2       1
15    803      0 39.2712        1  1       1
16    855      0 43.1233        1  1       2
17   1040      0 38.8932        2  1       2
18   1106      0 44.6000        1  1       1
19   1129      0 53.9068        1  2       1
20   1206      0 44.2055        2  2       1
21   1227      0 59.5890        1  2       2
22    268      1 74.5041        2  1       2
23    329      1 43.1370        2  1       1
24    353      1 63.2192        1  2       2
25    365      1 64.4247        2  2       1
26    377      0 58.3096        1  2       1

Le colonne futime e fustat sono rispettivamente il tempo di analisi dei 26 soggetti, e lo status (morto e censored); rx rappresenta i due gruppi. Vogliamo stimare il peso dei regressori age (età) nei due gruppi rx. Utilizziamo la funzione coxph():

cph1 <- coxph(Surv(futime,fustat) ~ rx+age, data=ovarian)
summary(cph1)

Call:
coxph(formula = Surv(futime, fustat) ~ rx + age, data = ovarian)

  n= 26 

        coef exp(coef) se(coef)      z Pr(>|z|)   
rx  -0.80397   0.44755  0.63205 -1.272  0.20337   
age  0.14733   1.15873  0.04615  3.193  0.00141 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

    exp(coef) exp(-coef) lower .95 upper .95
rx     0.4475      2.234    0.1297     1.545
age    1.1587      0.863    1.0585     1.268

Rsquare= 0.457   (max possible= 0.932 )
Likelihood ratio test= 15.89  on 2 df,   p=0.0003551
Wald test            = 13.47  on 2 df,   p=0.00119
Score (logrank) test = 18.56  on 2 df,   p=9.34e-05

Scriviamo subito l'equazione del nostro modello:
$$h(t|rx,age)=\lambda_0(t)exp(\beta_1\cdot rx + \beta_2\cdot age)$$
$$h(t|rx,age)=h(t)exp(-0.80397\cdot rx + 0.14733\cdot age)$$

Qual è il rischio relativo di un soggetto appartenente al gruppo rx=1 di 50 anni, rispetto ad un soggetto appartenente al gruppo rx=2 di 60 anni?
$$hr(rx=1,age=50:rx=2,age=60)=exp(-0.804(1-2)+0.147(50-60))=0.514$$

Possiamo mettere in grafico il modello, con la funzione plot(survfit(cph1)). Vengono consideato tutte le covariate prese in considerazione dal modello.

Nell'output vengono anche forniti i valori z e p-value per verificare se i coefficienti beta delle covariate sono significativi. La colonna zeta è il test di Wald per ogni regressore, ottenuta dividendo il valore beta per lo standard error del coefficiente stesso.
In questo caso è significativo solamente il regressore relativo alla covariata age.
Il likelihood ratio test, il Wald test e la statistica score chi-square che osserviamo alla fine rappresentano i test statistici per verificare se tutti i regressori considerati insieme sono differenti da zero (serve quindi per validare il modello fittato). In questo caso p-value < 0.05 in tutti i test, e quindi l'ipotesi nulla viene rifiutata.

---

Analizziamo un altro esempio, utilizzando il dataset qui disponibile. Poichè nel txt non è specificato il nome delle colonne, lo assegnamo noi con la funzione names():

data <- read.table("http://www.sph.emory.edu/~cdckms/CoxPH/Anderson_i.txt", F, sep=",")
names(data) <- c("pred1", "pred2", "pred3", "time", "status")
data[1:10, 1:5]

   pred1 pred2 pred3 time status
1      0  1.45     0   35      0
2      0  1.47     0   34      0
3      0  2.20     0   32      0
4      0  2.53     0   32      0
5      0  1.78     0   25      0
6      0  2.57     0   23      1
7      0  2.32     0   22      1
8      0  2.01     0   20      0
9      0  2.05     0   19      0
10     0  2.16     0   17      0

Fittiamo il modello di regressione di Cox, conoscendo che esiste un'interazione tra pred2 e pred3:

coxph1 <- coxph(Surv(time,status)~pred1+pred2*pred3, data=data)
summary(coxph1)

Call:
coxph(formula = Surv(time, status) ~ pred1 + pred2 * pred3, data = data)

  n= 42 

                  coef  exp(coef)   se(coef)      z Pr(>|z|)    
pred1          8.01642 3030.32167    3.06634  2.614  0.00894 ** 
pred2          1.84990    6.35920    0.44718  4.137 3.52e-05 ***
pred3         -4.12837    0.01611    1.88621 -2.189  0.02862 *  
pred2:pred3    0.60680    1.83455    0.29746  2.040  0.04136 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

            exp(coef) exp(-coef) lower .95 upper .95
pred1       3.030e+03    0.00033 7.4372579 1.235e+06
pred2       6.359e+00    0.15725 2.6470659 1.528e+01
pred3       1.611e-02   62.07646 0.0003995 6.496e-01
pred2:pred3 1.835e+00    0.54509 1.0240708 3.286e+00

Rsquare= 0.704   (max possible= 0.988 )
Likelihood ratio test= 51.15  on 4 df,   p=2.072e-10
Wald test            = 35.56  on 4 df,   p=3.57e-07
Score (logrank) test = 71.79  on 4 df,   p=9.548e-15

Il p-value di tutti i regressori è minore di 0.05, così come quelli relativi ai test sul modello completo.

Una delle assunzioni del modello è il "proportional hazards" (i coefficienti del modello non dipendono dal tempo), il quale può essere esaminato con la funzione cox.zph(), per effettuare il test di rischio proporzionale:

 cox.zph(coxph1)

                rho chisq     p
pred1        0.1303 0.587 0.444
pred2        0.2117 1.089 0.297
pred3       -0.1165 0.544 0.461
pred2:pred3  0.0936 0.394 0.530
GLOBAL           NA 1.716 0.788

Essendo il p-value > 0.05 per tutte le covariate e per il modello globale, l'assunzione è stata verificata.
Possiamo anche verificare graficamente i residui (di Schoenfeld) plottando la funzione:

par(mfrow=c(2,2))
plot(cox.zph(coxph1))

I punti rispettano una certa linearità (come richiesto dall'assunzione); se i punti fossero stati maggiormente dispersi, non avremmo potuto verificare l'assuzione di proporzionalità, e quindi si dovrebbe procedere con trasformazioni della variabile.

Un'altra verifica che andrebbe affettuata è quella della non-linearità delle covariate (salvo quelle dicotomiche). In questo pdf nella sezione 5.3 sono presenti i test da effettuare. -- Continua a leggere! --

Calcolo della tabella ANOVA partendo da dati aggregati

Può capitare di dover eseguire un test per verificare l'uguaglianze delle medie tra più gruppi, di cui però non disponiamo dei dati grezzi (il valore di ogni singolo dato), ma solamente della media, della deviazione standard e della numerosità campionaria di ciascun gruppo. Questo capita frequentemente nel momento in cui vogliamo fare ulteriori analisi sui dati pubblicati in articoli scientifici, che ci forniscono i dati descrittivi dei gruppi, e non i gruppi in sè.
Si possono ricavare facilmente le formule per ottenere la tabella ANOVA partendo dai dati descrittivi (in questo pdf sono presenti le formule per eseguire la one-way ANOVA, e le ANOVA bi- e tridimensionale partendo dai dati descrittivi). Supponiamo di avere i seguenti 4 gruppi:

a <- c(7.4, 6.6, 6.7, 6.1, 6.5, 7.2)
b <- c(7.1, 7.3, 6.8, 6.9, 7)
c <- c(6.8, 6.3, 6.4, 6.7, 6.5, 6.8)
d <- c(6.4, 6.9, 7.6, 6.8, 7.3)

Da questi dati calcoliamo: numerosità campionaria, media e deviazione standard, e raccogliamo tutti i dati in una tabella di riepilogo:

> mean(a)
[1] 6.75
> mean(b)
[1] 7.02
> mean(c)
[1] 6.583333
> mean(d)
[1] 7
> 
> length(a)
[1] 6
> length(b)
[1] 5
> length(c)
[1] 6
> length(d)
[1] 5
> 
> sd(a)
[1] 0.4764452
> sd(b)
[1] 0.1923538
> sd(c)
[1] 0.2136976
> sd(d)
[1] 0.4636809

Adesso dimentichiamo i dati campionari, e facciamo una ANOVA parametrica sui dati di quest'ultima tabella. Le formule sono:

$$GrandMean=\overline{\overline{X}}=\sum \overline{x_i}\cdot n_i$$

$$SSB=\sum (\overline{x_i}-\overline{\overline{X}})^2\cdot n_i$$

$$SSE=\sum s_i^2\cdot (n-1)$$

La GrandMean è la media di tutti i gruppi; la SSB è la devianza tra gruppi, e la SSE è la devianza entro gruppi (o d'errore). Le altre formule sono quelle consuete, che già conoscevamo:

$$df_{SSB}=k-1$$

$$df_{SSE}=N-k$$

$$MSB=\frac{SSB}{df_{SSB}}$$

$$MSE=\frac{SSE}{df_{SSE}}$$

$$F=\frac{MSB}{MSE}$$

I df sono i gradi di liberta; MSB e MSE sono la varianza tra gruppi e d'errore; F è il valore della statistica test.

Traduciamo questi calcoli in R, e procediamo:

nA <- 6; nB <- 5; nC <- 6; nD <- 5
meanA <- 6.75; meanB <- 7.02; meanC <- 6.58; meanD <- 7 
sdA <- .476; sdB <- .192; sdC <- .214; sdD <- .464

GrandMean <- (meanA * nA + meanB * nB + meanC * nC + meanD * nD) / (nA + nB + nC + nD)
> GrandMean
[1] 6.821818

SSB <- (((meanA - GrandMean)^2)*nA) + (((meanB - GrandMean)^2)*nB) + (((meanC - GrandMean)^2)*nC) + (((meanD - GrandMean)^2)*nD)
> SSB
[1] 0.7369273

SSE <- (sdA^2 * (nA-1)) + (sdB^2 * (nB-1)) + (sdC^2 * (nC-1)) + (sdD^2 * (nD-1))
> SSE
[1] 2.3705

dfSSB <- 4 - 1
dfSSE <- (nA+nB+nC+nD) - 4

MSB <- SSB / dfSSB
> MSB
[1] 0.2456424
 
MSE <- SSE / dfSSE
> MSE
[1] 0.1316944
 
F <- MSB / MSE
> F
[1] 1.865245

Possiamo calcolare il p-value di una distribuzione F, per prendere la nostra decisione statistica, in questo modo:

p.value <- 1 - pf(F, dfSSB, dfSSE)
> p.value
[1] 0.1716690

Raccogliamo tutti i dati nella classica tabella ANOVA:


$$\begin{tabular}{c|ccccc}\hline Source of variation&df&SS&MS&F-statistic&p-value\\ \hline Between&3&.737&.246&1.86&.172\\Error&18&2.37&.132&&\\ \end{tabular}$$

Essendo p-value > 0.05, accettiamo l'ipotesi nulla: i 4 gruppi hanno medie uguali.

Adesso vediamo i valori che avremmo ottenuto, se avessimo utilizzato i dati campionari:

a <- c(7.4, 6.6, 6.7, 6.1, 6.5, 7.2)
b <- c(7.1, 7.3, 6.8, 6.9, 7)
c <- c(6.8, 6.3, 6.4, 6.7, 6.5, 6.8)
d <- c(6.4, 6.9, 7.6, 6.8, 7.3)

dati <- c(a,b,c,d)
gruppi = factor(rep(letters[1:4], times=c(length(a),length(b),length(c),length(d))))

anova(lm(dati ~ gruppi))

Analysis of Variance Table

Response: dati
          Df  Sum Sq Mean Sq F value Pr(>F)
gruppi     3 0.72730 0.24243  1.8402 0.1760
Residuals 18 2.37133 0.13174 

E questa è la tabella ottenuta con la funzione xtable, disponibile nel pacchetto xtable (converte le tabelle di R in codice LaTeX):


$$\begin{tabular}{lrrrrr} \hline & Df & Sum Sq & Mean Sq & F value & Pr($>$F) \\ \hline gruppi & 3 & 0.73 & 0.24 & 1.84 & .1760 \\ Residuals & 18 & 2.37 & 0.13 & & \\ \hline\end{tabular}$$

Confrontando i valori, essi risultano molto simili, approssimazioni a parte. -- Continua a leggere! --

ANOVA post-hoc tests #2

In questo post continuiamo a vedere come utilizzare R per i test post-hoc di una ANOVA significativa. Raccomando la lettura del primo post in riguardo allo stesso argomento, dove si è parlato del test LSD di Fisher, e del test HSD di Tukey. Qui utilizzerò gli stessi dati, per esaminare insieme i seguenti test:

• test t di Bonferroni


• test S di Scheffé

---

Test t di Bonferroni

Il test per i confronti multipli simultanei proposto dallo statistico italiano Bonferroni, si basa sul confronto tra il valore t di Bonferroni, e una quantità che viene calcolata, di volta in volta, per ogni coppia di gruppi che si esamina.
Il t di Bonferroni si ricerca sulle tavole t di Student, in corrispondenza dei gradi di libertà della varianza d'errore (N-k), e del valore alpha ottenuto dividendo l'alpha di scelta (solitamente 0.05) per il doppio del numero di confronti da effettuare. Nel nostro caso abbiamo 4 gruppi; il numero di confronti da effettuare è quindi 6, che moltiplicato per 2 fa 12. In R tale t si calcola in questo modo. Supponiamo di avere una lista, con i k gruppi; il doppio del numero di confronti da effettuare è semplicemente k*(k-1). Quindi

a <- c(13.47, 10.21, 15.10, 14.65, 9.03, 15.14, 6.09, 3.43, 5.95, 10.72, 10.01, 8.17)
b <- c(4.02, 14.03, 3.09, 10.25, 6.25, 1.00, 7.50, 10.18, 6.07, 2.03, 4.17, 7.28)
c <- c(10.56, 14.61, 11.88, 11.43, 9.73, 12.37, 3.82, 13.04, 13.28, 15.28, 12.90, 10.98)
d <- c(7.74, 2.77, 10.29, 4.03, 10.23, 7.67, 7.93, 6.75, 5.60, 8.19, 2.02, 7.88)
 
dati <- list(g1=a, g2=b, g3=c, g4=d)

confronti2 <- length(dati) * (length(dati) - 1)
alpha.bonferr <- 0.05 / confronti2

abs(qt(alpha.bonferr, 44))
[1] 2.762815

Questo valore va confrontato con il risultato della seguente formula. Siano r ed s due generici gruppi in esame:

$$t_{bonferroni}=\frac{|\overline{x_r}-\overline{x_s}|}{\sqrt{S_e^2(\frac{1}{n_r}+\frac{1}{n_s})}}$$

con $S_e^2$ = varianza entro gruppi.
Se il t di Bonferroni calcolato è maggiore di quello tabulato (nel nostro caso 2.76), allora la differenza è significativa (le medie sono significativamente differenti); altrimenti sono uguali.
Procedendo manualmente, quindi, si crea una matrice triangolare con i t di Bonferroni calcolati, e si individuano quelli maggiori (in valore assoluto) del t-tabulato, ossia le differenze significative (qui le indico con un asterisco). Con i nostri dati avremmo:

$$\begin{tabular}{c|c|c|c}&b&c&d\\\hline a&2.79*&10.86&2.479\\\hline b&&3.88*&0.317\\\hline c&&&3.56* \end{tabular}$$

In R non è disponibile una funzione per eseguire direttamente il test di Bonferroni. Tuttavia sono riuscito a trovare in rete una bozza della funzione, che ho modificato per svolgere i calcoli visti finora. Ecco il codice, quindi, per eseguire il test di Bonferroni in R:

bonferroni.test = function(x, y)
{
# x is a list of vectors
# y is the alpha 
 num <- 0
 denom <- 0
 grandmean.num <- 0
 sb2.num <- 0
 nn <- numeric(length(x))
 ss2 <- numeric(length(x))
 for(i in 1:length(x)) {
  nn[i] <- length(x[[i]])
  ss2[i] <- (sum((x[[i]] - mean(x[[i]]))^2))/(nn[i] - 1
   )
  num <- ((nn[i] - 1) * ss2[i]) + num
  denom <- nn[i] + denom
  grandmean.num <- (nn[i] * mean(x[[i]])) + 
   grandmean.num
 }
 sw2 <- num/(denom - length(x))
 grandmean <- grandmean.num/denom
 for(j in 1:length(x)) {
  sb2.num <- (nn[j] * ((mean(x[[j]]) - grandmean)^2)) + 
   sb2.num
 }
 sb2 <- sb2.num/(length(x) - 1)
 FF <- sb2/sw2
 p.val <- 1 - pf(FF, length(x) - 1, denom - length(x))
 effort <- numeric(length(x) * (length(x) - 1))
 effort1 <- character(length(x) * (length(x) - 1))
 tt <- matrix(effort, ncol = length(x))
 pp <- matrix(effort, ncol = length(x))
 qq <- matrix(effort1, ncol = length(x))
 alpha.star <- y/((length(x) * (length(x) - 1)))
 deg.free <- denom - length(x)
 for(k in 1:(length(x) - 1)) {
  for(l in 2:length(x)) {
   if(k >= l) {
    tt[k, l] <- 0
   }
   else {
    tt[k, l] <- abs((mean(x[[k]]) - mean(x[[l
      ]]))/sqrt(sw2 * ((1/nn[k]) + (1/nn[
      l]))))
    alpha <- 1 - (alpha.star)
    t.tab <- qt(alpha, deg.free)
    pp[k, k] <- 0
    if(tt[k, l] == 0) {
      qq[k, l] <- c("Zero")
    }
    else if(tt[k, l] < t.tab) {
      qq[k, l] <- c("Not")
    }
    else if(tt[k, l] > t.tab) {
      qq[k, l] <- c("Signif")
    }
   }
  }
 }
 results <- list(grandmean, sw2, sb2, FF, p.val, tt, 
  alpha.star, qq)
 results <- list(Grand.mean = grandmean, Within.var = sw2, 
  Between.var = sb2, F.value = FF, p.value = p.val, 
  bonferroni.t = tt, alpha.star = alpha.star, 
  bonferroni.sig = qq)
 results
}

L'input è costituito dalla lista dei gruppi (x), e dal valore di alpha da cui vogliamo partire (qui userò alpha = 0.05). Applichiamo subito la funzione, e studiamo l'output:

a <- c(13.47, 10.21, 15.10, 14.65, 9.03, 15.14, 6.09, 3.43, 5.95, 10.72, 10.01, 8.17)
b <- c(4.02, 14.03, 3.09, 10.25, 6.25, 1.00, 7.50, 10.18, 6.07, 2.03, 4.17, 7.28)
c <- c(10.56, 14.61, 11.88, 11.43, 9.73, 12.37, 3.82, 13.04, 13.28, 15.28, 12.90, 10.98)
d <- c(7.74, 2.77, 10.29, 4.03, 10.23, 7.67, 7.93, 6.75, 5.60, 8.19, 2.02, 7.88)
 
dati <- list(g1=a, g2=b, g3=c, g4=d)

bonferroni.test(dati, 0.05)

$Grand.mean
[1] 8.725417

$Within.var
[1] 11.32240

$Between.var
[1] 81.22261

$F.value
[1] 7.173622

$p.value
[1] 0.0005046445

$bonferroni.t
     [,1]     [,2]     [,3]      [,4]
[1,]    0 2.796572 1.086478 2.4793041
[2,]    0 0.000000 3.883050 0.3172684
[3,]    0 0.000000 0.000000 3.5657816

$alpha.star
[1] 0.004166667

$bonferroni.sig
     [,1] [,2]     [,3]     [,4]    
[1,] ""   "Signif" "Not"    "Not"   
[2,] ""   ""       "Signif" "Not"   
[3,] ""   ""       ""       "Signif"

Per comodità la funzione fornisce in output anche le varianze (tra gruppi ed entro gruppi), il risultato della statistica F di una ANOVA parametrica, e il suo corrispettivo p-value (per controllare se è utile eseguire i confronti multipli oppure no).
Quindi viene presentata la matrice triangolare con i valori t di Bonferroni calcolati come sopra. E' in pratica la stessa tabella che abbiamo scritto manualmente. Quindi viene riportato il valore alpha (0.05 / 12), ed infine una matrice triangolare contenente la significatività del test.
Risultano quindi significative (ossia le medie sono differenti) i confronti tra i gruppi 1-2, 2-3, 3-4. Osservate che il test HSD di Tukey (applicato in questo post) ha dato lo stesso risultato, ossia ha individuato come significative le stesse coppie.

---

Test di Scheffé per i confronti multipli

Il test implementato da Scheffé è molto simile al testo di Bonferroni. Per ogni coppia di dati si calcola la seguente quantità:

$$S=\frac{|\overline{x_r}-\overline{x_s}|}{\sqrt{S_e^2\cdot(\frac{1}{n_r}+\frac{1}{n_s})}}$$

con $S_e^2$ = varianza entro gruppi.
Questo valore (che è lo stesso della procedura di Bonferroni) viene confrontato con il valore critico:

$$S_\alpha=\sqrt{(k-1)F_{(alpha, k-1, N-k)}}$$

Se la quantità calcolata supera il valore critico, allora la differenza è significativa.
Il test è quindi molto simile al test di Bonferroni, eccetto il valore critico di confronto. Consideriamo il consueto dataset, e calcoliamo il valore critico manualmente:

a <- c(13.47, 10.21, 15.10, 14.65, 9.03, 15.14, 6.09, 3.43, 5.95, 10.72, 10.01, 8.17)
b <- c(4.02, 14.03, 3.09, 10.25, 6.25, 1.00, 7.50, 10.18, 6.07, 2.03, 4.17, 7.28)
c <- c(10.56, 14.61, 11.88, 11.43, 9.73, 12.37, 3.82, 13.04, 13.28, 15.28, 12.90, 10.98)
d <- c(7.74, 2.77, 10.29, 4.03, 10.23, 7.67, 7.93, 6.75, 5.60, 8.19, 2.02, 7.88)

dati <- list(g1=a, g2=b, g3=c, g4=d)

N <- length(a) + length(b) + length(c) + length(d)
alpha <- 1 - .05
F.cv <- sqrt( (length(dati)-1) * qf(alpha, (length(dati)-1), (N-length(dati))) )

> F.cv
[1] 2.906785

Proviamo ora a calcolare la quantità S per la coppia di gruppi a-b:

num <- abs( mean(a) - mean(b) )

dati.tot <- c(a,b,c,d)
gruppi = factor(rep(letters[1:4], each = 12))
within.var <-  anova(lm(dati.tot~gruppi))[2,3]

den <- sqrt( within.var*(1/length(a)+1/length(b)) )

S <- num/den

> S
[1] 2.796572

Per la coppia a-b, S < S-critico, quindi la media di questi due gruppi è significativamente uguale.
Possiamo procedere in modo analogo per le altre coppie. Oppure utilizziamo la funzione sottostante, che è sostanzialmente uguale a quella del test di Bonferroni visto poco sopra, salvo le modifiche del test di Scheffé. Questo è il codice per eseguire il test di Scheffé in R:

scheffe.test = function(x, y)
{
# x is a list of vectors
# y is the alpha 
 num <- 0
 denom <- 0
 grandmean.num <- 0
 sb2.num <- 0
 nn <- numeric(length(x))
 ss2 <- numeric(length(x))
 for(i in 1:length(x)) {
  nn[i] <- length(x[[i]])
  ss2[i] <- (sum((x[[i]] - mean(x[[i]]))^2))/(nn[i] - 1
   )
  num <- ((nn[i] - 1) * ss2[i]) + num
  denom <- nn[i] + denom
  grandmean.num <- (nn[i] * mean(x[[i]])) + 
   grandmean.num
 }
 sw2 <- num/(denom - length(x))
 grandmean <- grandmean.num/denom
 for(j in 1:length(x)) {
  sb2.num <- (nn[j] * ((mean(x[[j]]) - grandmean)^2)) + 
   sb2.num
 }
 sb2 <- sb2.num/(length(x) - 1)
 FF <- sb2/sw2
 p.val <- 1 - pf(FF, length(x) - 1, denom - length(x))
 effort <- numeric(length(x) * (length(x) - 1))
 effort1 <- character(length(x) * (length(x) - 1))
 tt <- matrix(effort, ncol = length(x))
 pp <- matrix(effort, ncol = length(x))
 qq <- matrix(effort1, ncol = length(x))
 alpha.star <- y/((length(x) * (length(x) - 1)))
 deg.free <- denom - length(x)
 for(k in 1:(length(x) - 1)) {
  for(l in 2:length(x)) {
   if(k >= l) {
    tt[k, l] <- 0
   }
   else {
    tt[k, l] <- abs((mean(x[[k]]) - mean(x[[l
      ]]))/sqrt(sw2 * ((1/nn[k]) + (1/nn[
      l]))))
    alpha <- 1 - y
    f.tab <- sqrt( (length(x)-1) * qf(alpha, (length(x)-1), (denom-length(x))))
    pp[k, k] <- 0
    if(tt[k, l] == 0) {
      qq[k, l] <- c("Zero")
    }
    else if(tt[k, l] < f.tab) {
      qq[k, l] <- c("Not")
    }
    else if(tt[k, l] > f.tab) {
      qq[k, l] <- c("Signif")
    }
   }
  }
 }
 results <- list(grandmean, sw2, sb2, FF, p.val, tt, 
  alpha.star, qq)
 results <- list(Within.var = sw2, 
  Between.var = sb2, F.value = FF, p.value = p.val, 
  Scheffé.S = tt, alpha = y, S.critical.value = f.tab,
  Scheffé.sig = qq)
 results
}

Adesso studiamo l'output della funzione:

a <- c(13.47, 10.21, 15.10, 14.65, 9.03, 15.14, 6.09, 3.43, 5.95, 10.72, 10.01, 8.17)
b <- c(4.02, 14.03, 3.09, 10.25, 6.25, 1.00, 7.50, 10.18, 6.07, 2.03, 4.17, 7.28)
c <- c(10.56, 14.61, 11.88, 11.43, 9.73, 12.37, 3.82, 13.04, 13.28, 15.28, 12.90, 10.98)
d <- c(7.74, 2.77, 10.29, 4.03, 10.23, 7.67, 7.93, 6.75, 5.60, 8.19, 2.02, 7.88)

dati <- list(g1=a, g2=b, g3=c, g4=d)

scheffe.test(dati, .05)

$Within.var
[1] 11.32240

$Between.var
[1] 81.22261

$F.value
[1] 7.173622

$p.value
[1] 0.0005046445

$Scheffé.S
     [,1]     [,2]     [,3]      [,4]
[1,]    0 2.796572 1.086478 2.4793041
[2,]    0 0.000000 3.883050 0.3172684
[3,]    0 0.000000 0.000000 3.5657816

$alpha
[1] 0.05

$F.critical.value
[1] 2.906785

$Scheffé.sig
     [,1] [,2]  [,3]     [,4]    
[1,] ""   "Not" "Not"    "Not"   
[2,] ""   ""    "Signif" "Not"   
[3,] ""   ""    ""       "Signif"

Nell'ordine ci vengono forniti:
- la varianza entro gruppi
- la varianza tra gruppi
- il valore della statistica F dell'ANOVA e il suo p-value
- la quantità S calcolata per ogni coppia
- il valore alpha prescelto
- il valore critico di confronto
- la significativà (una coppia è significativa quando le medie sono differenti)

In base al test di Scheffé sono significative le coppie 2-3 e 3-4. -- Continua a leggere! --

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. -- Continua a leggere! --

Regressione lineare multivariata

Accolgo volentieri l'invito di Fabio, e mi accingo a cominciare alcuni post sulla statistica multivariata.
Nella regressione lineare semplice, abbiamo immaginato che una certa variabile Y dipendesse dall'andamento di un'altra variabile (X), in maniera lineare con andamento crescente o decrescente. Abbiamo quindi visto come realizzare e disegnare la retta che pone in relazione le due variabili, e come valutare la bontà del modello.
Nella realtà (scientifica, economica, psicometrica, etc.), quasi mai un evento dipende solamente dall'andamento di un certo fattore. Tutti gli eventi (anche i più comuni) sono influenzati da numerosissimi elementi. Risulta pertanto molto più utile formulare un modello che tenga conto di tutte queste influenze (o, come vedremo in seguito, delle influenze maggiori). Ciò si ottiene con lo studio della regressione lineare multipla (o multivariata). In generale si indica con Y la variabile dipendente, e con X seguito da un numero in pedice le variabili indipendenti che si suppone abbiano un effetto. Le X vengono chiamate predittori e la formula generale del modello che cerchiamo è:

$$Y=\beta_0+\beta_1X_1+\beta_2X_2+...+\beta_nX_n$$

$\beta_0$ è l'intercetta, mentre $\beta_1$, $\beta_2$ eccetera, sono i regressori, i quali rappresentano il coefficiente angolare della retta che otterremmo al variare del predittore corrispondente, qualora tutti gli altri predittori fossero costanti.
La rappresentazione grafica dipende dal numero di predittori che si vogliono considerare. Con un solo predittore, si ottiene (come abbiamo visto) una retta; con 2 predittori si ottiene un piano nello spazio a 3 dimensioni; con 3 o più predittori non è possibile la rappresentazione grafica, in quanto occorrerebbe rappresentare uno spazio a più di 3 dimensioni.

Cominciamo a vedere i passaggi da effettuare quando si vuole realizzare un modello di regressione multipla. Utilizzerò qui il dataset disponibile su questo sito. In esso viene raccolto il consumo di gelati pro-capite, considerando il periodo dell'anno, il prezzo del gelato, i guadagni delle famiglie e la temperatura media ambientale. Potremmo ad esempio chiederci: il consumo di gelato aumenta o diminuisce quando la temperatura aumenta? Sul consumo di gelati influisce maggiormente la temperatura o il suo prezzo?
Importiamo i dati. Sul sito è disponibile un file in formato .dat. Possiamo aprirlo con il blocco note per visualizzarlo, e renderci conto di come è strutturato. In R è facile importare questo formato di dati, anche senza la necessità di salvare il file sul nostro computer, utilizzando questo comando:

dati <- read.table("http://www.sci.usq.edu.au/staff/dunn/Datasets/Books/Hand/Hand-R/icecream-R.dat", header = TRUE, row.names = 1)

Adesso in dati è contenuta una tabella, con i valori riportati nel file specificato. Perchè R possa accedere al file .dat, non ci deve essere alcuna limitazione di connessione; se non riesce ad accedere ad internet, conviene controllare il proprio firewall (che spesso blocca la connessione a internet da parte di alcuni programmi). Per comodità trasformiamo la tabella in un dataframe:

datiDF <- data.frame(dati)

Disponiamo quindi di 4 variabili e 30 osservazioni. Per prima cosa cerchiamo di indagare il consumo di gelato in base alla temperatura e al prezzo (ignoriamo le entrate). Le colonne che ci interessano, per ora, sono datiDF$Consumption, datiDF$Price, datiDF$Temp. Cerchiamo un modello di regressione multipla a 3 predittori. La tecnica per stimare i regressori è detta ordinary least square OLS, che ripercorre gli stessi principi della tecnica dei minimi quadrati della regressione semplice.
In R stimiamo il modello in questo modo:

model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
summary(model)

Call:
lm(formula = datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08226 -0.02051  0.00184  0.02272  0.10076 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.59655    0.25831   2.309   0.0288 *  
datiDF$Price -1.40176    0.92509  -1.515   0.1413    
datiDF$Temp   0.00303    0.00047   6.448 6.56e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.04132 on 27 degrees of freedom
Multiple R-squared: 0.6328,     Adjusted R-squared: 0.6056 
F-statistic: 23.27 on 2 and 27 DF,  p-value: 1.336e-06 

Il model contiene la formula che vogliamo indagare. Dapprima si specifica la variabile dipendente (il consumo), quindi si specificano le variabili indipendenti, tra loro unite da un più.
Con la funzione summary() otteniamo una serie di informazioni. Innanzitutto il valore dell'intercetta e dei regressori; in questo caso:
$\beta_0=0.59655$
$\beta_1=-1.40176$
$\beta_2=0.00303$
Pertanto il modello di regressione multipla stimato è:
$Y=0.59655-1.40176\cdot X_1+0.00303\cdot X_2$

Osserviamo il segno del primo regressore, che descrive il predittore "prezzo": esso è negativo, il che significa che per ogni aumento unitario del prezzo del gelato, il consumo procapite diminuisce di circa una unità e mezzo (1.40); viceversa la temperatura ha un effetto positivo (seppure molto scarso) sul consumo di gelato (all'aumentare della temperatura, aumenta anche il consumo).

Possiamo rappresentare graficamente il modello creato. Installiamo il pacchetto scatterplot3d, e utilizziamo l'unica funzione in esso compreso, specificando le variabili nell'ordine x, y, z (nel nostro caso meglio specificare (X1, X2, Y)):

library(scatterplot3d)
sc <- scatterplot3d(datiDF$Price, datiDF$Temp, datiDF$Consumption, pch=16)

Possiamo anche intersecare il piano che abbiamo trovato, in questo modo:

model <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
sc$plane(model, col="red")

Analizziamo adesso il modello. Il coefficiente di determinazione (Multiple R-squared) consente di valutare la bontà del modello, ossia la proporzione di variabilità della Y spiegata dalle variabili esplicative considerate. In questo caso R-squared è pari a 0.63 (il 63% dei dati viene spiegato dalle variabili esplicative). Il valore R-squared adjusted è una correzione del coefficiente di determinazione che tiene conto del numero di predittori utilizzati (secondo alcuni statistici è più preciso quest'ultimo valore per definire la bontà del modello o goodness-of-fit).
Inoltre è opportuno effettuare una verifica dell'intercetta e dei regressori, ossia un t-test con ipotesi nulla H0: beta = 0. Ossia si deve verificare se il valore campionario possa essere esteso all'intera popolazione. Nel nostro caso i p-value associati alle variabili sono:
p-value $\beta_0$ = 0.0288 (significativo)
p-value $\beta_1$ = 0.1413 (non significativo - dovremmo ignorarla)
p-value $\beta_2$ = 6.56e-07 (significativo)

La statistica F ci indica se il modello è da scartare nella sua interessa, oppure se può essere ritenuto valido. Essendo F-statistic maggiore dell'F-tabulato, rifiutiamo l'ipotesi nulla che il modello sia da scartare nella sua interezza.

Possiamo calcolare in questo modo gli intervalli di confidenza dell'intercetta e dei 2 regressori:

confint(model)
                    2.5 %      97.5 %
(Intercept)   0.066549244 1.126551076
datiDF$Price -3.299896373 0.496375437
datiDF$Temp   0.002066035 0.003994569

Se vogliamo invece calcolare l'intervallo di confidenza di un valore predetto (ossia se vogliamo prevedere il valore che assumerà Y, dati i predittori), utilizziamo la funzione predict():

model <- lm(Consumption ~ Price + Temp, data=datiDF)
predict(model, data.frame(Price=0.280, Temp=30), interval="confidence")

        fit       lwr       upr
1 0.2949663 0.2700105 0.3199220

---

Proviamo a complicare un pò le cose, e supponiamo di considerare nel nostro modello anche la variabile "entrate". All'aumentare dei guadagni, aumenta anche il consumo di gelato?

model <- lm(Consumption ~ Price + Temp + Income, data=datiDF)
summary(model)

Call:
lm(formula = Consumption ~ Price + Temp + Income, data = datiDF)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.065302 -0.011873  0.002737  0.015953  0.078986 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.1973151  0.2702162   0.730  0.47179    
Price       -1.0444140  0.8343573  -1.252  0.22180    
Temp         0.0034584  0.0004455   7.762  3.1e-08 ***
Income       0.0033078  0.0011714   2.824  0.00899 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.03683 on 26 degrees of freedom
Multiple R-squared: 0.719,      Adjusted R-squared: 0.6866 
F-statistic: 22.17 on 3 and 26 DF,  p-value: 2.451e-07 

Il regressore beta3 associato al predittore entrate risulta essere significativo, e di segno positivo. Quindi all'aumentare delle entrate, aumenta anche il consumo di gelati.
Il modello a 3 predittori viene complessivamente accettato (la statistica F risulta significativa), e si assiste ad un miglioramento della goodness-of-fit (è aumentato il valore dell'R quadro aggiustato).

Una domanda che potremmo porci è la seguente: tra le tre variabili considerate, qual è quella che maggiormente influisce sulla variabile consumo? Per far questo possiamo confrontare i valori dei regressori beta, ma solo dopo averli standardizzati (altrimenti non è possibile il confronto). Nel pacchetto QuantPsyc è presenta la funzione lm.beta() che fa al caso nostro:

library(QuantPsyc)
lm.beta(model)

     Price       Temp     Income 
-0.1324351  0.8632558  0.3140085 

Quelli qui riportati sono i coefficienti beta standardizzati. Osserviamo che la variabile esplicativa che maggiormente influisce sul consumo di gelati è la temperatura.
Se non volete scaricare il pacchetto per utilizzare solo questa funzione, è sufficiente sapere che il calcolo di coefficienti beta standardizzati è pari al prodotto del valore beta per il rapporto tra la deviazione standard della X considerata, e la deviazione standard della Y:

coef(model)[2]*(sd(datiDF$Price))/(sd(datiDF$Consumption))
     Price 
-0.1324351 

coef(model)[3]*(sd(datiDF$Temp))/(sd(datiDF$Consumption))
     Temp 
0.8632558 

coef(model)[4]*(sd(datiDF$Income))/(sd(datiDF$Consumption))
   Income 
0.3140085 

In maniera analoga (ma anche un pò più rischiosa), si potrebbe valutare la dipendenza anche osservando la correlazione tra le variabili. Se calcoliamo il coefficiente di correlazione (con la funzione cor())tra il consumo e le variabili prezzo, entrate, temperatura, osserviamo che il maggior valore è relativo alla variabile temperatura (è maggiore la correlazione tra queste due variabili).

Spesso si finisce col cercare una quantità di possibili variabili esplicative, da rendere poco agevole lo sviluppo di un modello di regressione multipla, che teoricamente dovrebbe contenere il minor numero di predittori possibile. La scelta non si basa solo sulla valutazione "ad occhio" del peso di ciascun predittore (come abbiamo fatto prima), ma viene utilizzato il test F parziale, che in pratica va a testare se c'è un incremento dell'R-quadro nel momento in cui si aggiunge un'altra variabile esplicativa al modello.
In R per eseguire un test F-parziale, si confronta la tabella dell'analisi della varianza dei due modelli. Supponiamo di voler confrontare le bontà dei modelli (Price + Temp) e (Price + Temp + Income):

model1 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)
model2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income)
anova(model1, model2)

Analysis of Variance Table

Model 1: datiDF$Consumption ~ datiDF$Price + datiDF$Temp
Model 2: datiDF$Consumption ~ datiDF$Price + datiDF$Temp + datiDF$Income
  Res.Df      RSS Df Sum of Sq      F   Pr(>F)   
1     27 0.046090                                
2     26 0.035273  1  0.010817 7.9734 0.008989 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Se rifiutiamo l'ipotesi H0, vuol dire che l'inserimento della variabile esplicativa ha apportato un miglioramento al modello. In questo caso F = 7.97, con p-value < 0.05. Rifiutiamo quindi l'ipotesi nulla: l'aggiunta del predittore entrate migliora il modello.

Possiamo infine calcolare i coefficiente di determinazioni parziali. Fissata una delle variabili esplicative, a quanto ammonta la bontà del modello considerando l'altra variabile?
Riprendiamo a considerare solo 2 predittori (Temp e Price). Costruiamo prima i due modelli lineari semplici, quindi il modello a 2 predittori. Tenendo presente le formule per calcolare il coefficiente di determinazione parziale, abbiamo:

lm.x1 <- lm(datiDF$Consumption ~ datiDF$Price)
lm.x2 <- lm(datiDF$Consumption ~ datiDF$Temp)
lm.x1.x2 <- lm(datiDF$Consumption ~ datiDF$Price + datiDF$Temp)

ssr.x1.x2 <- sum(anova(lm.x1.x2)$"Sum Sq"[1:2])
ssr.x1    <- anova(lm.x1)$"Sum Sq"[1]
ssr.x2    <- anova(lm.x2)$"Sum Sq"[1]
ssr.x1.x2 - ssr.x1
sse.x1 <- anova(lm.x1)$"Sum Sq"[2]
sse.x2 <- anova(lm.x2)$"Sum Sq"[2]

(ssr.x1.x2 - ssr.x2) / sse.x2
[1] 0.07837315

(ssr.x1.x2 - ssr.x1) / sse.x1 
[1] 0.6062858

Il coefficiente di determinazione parziale corrispondente alla variabile X1, quando X2 è tenuta costante, ci dice che per una data temperatura, il 7.83% della variabilità del consumo di gelati può essere spiegata dalla variabilità del prezzo nei negozi. Il coefficiente di determinazione parziale corrispondente alla variabile X2, quando X1 è tenuta costante, ci dice che per un dato prezzo, il 60.63% della variabilità del consumo di gelati può essere spiegata dalla variabilità della temperatura (sono valori in accordo con il peso definito con il calcolo dei regressori standardizzati).


Queste viste finora solo alcune delle analisi preliminari che dobbiamo effettuare quando si cerca un modello di regressione lineare multipla. Nei post successivi vedremo altri esempi, e andremo ad effettuare altre analisi, come l'analisi dei residui, la scelta dei predittori e il confronto tra più modelli. -- Continua a leggere! --