mercoledì 12 agosto 2009

Operazioni con l'odd-ratio in R

Il software R non mette a disposizione molte funzioni per lavorare con le tabelle di contingenza 2x2, e quindi con gli odd-ratio. Tuttavia vedremo in questo post come calcolare l'OR, come stimare l'intervallo di confidenza dell'OR con i metodi di Wolff e Miettinen, come verificare la significatività dell'OR, e come confrontare tra loro due odd-ratios.

Consideriamo la seguente tabelle 2x2:



Immettiamo questi dati in R:


dati <- c(235, 103, 525, 635)
tabella <- matrix(dati, nrow=2)
colnames(tabella) <- c("malSI", "malNO")
rownames(tabella) <- c("InquiSI", "inquiNO")
table <- as.table(tabella)
table

malSI malNO
InquiSI 235 525
inquiNO 103 635


Abbiamo creato così la tabella a doppia entrata. La formula per ottenere l'OR è questa:



In R possiamo calcolarlo, partendo dalla tabella, in questo modo:


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


In R è quindi abbastanza facile richiamare i valori contenuti in una tabella, semplicemente specificando il nome della tabella, e le coordinate (numero riga, numero colonna) della cella da utilizzare.




Verifichiamo adesso l'ipotesi della significatività dell'OR appena calcolato. Si può procedere con uno Z-test, oppure in termini di chi-quadro corretto.

Cominciamo a vedere la verifica di ipotesi dell'OR con lo Z-test.

L'ipotesi nulla è che OR = 1, il che significa, in termini di logaritmo naturale, che H0: ln(OR) = 0.
La statistica test è:



L'errore standard di ln(OR) è:



Quindi la statistica test in R diventa:


esOR <- sqrt(sum( 1/dati ))
z <- log(OR) / esOR
z

[1] 7.685698


Poichè Z è maggiore del valore critico 1.96, rifiuto l'ipotesi nulla: l'odd ratio è significativo.

Adesso verifichiamo l'OR in termini di chi-quadro corretto. L'ipotesi nulla è H0: OR=1, e la statistica test è:



E(a) è il valore atteso della cella a, calcolato con la formula:
var(a) è la varianza della cella a, calcolata con questa formula:
Portiamo le formule in R e calcoliamo il chi-quadro corretto:


E.a <- (table[1,1]+table[1,2])*(table[1,1]+table[2,1]) / sum(dati)
E.a
[1] 171.4820

var.a <- (table[1,1]+table[1,2])*(table[1,1]+table[2,1])*
(table[2,1]+table[2,2])*(table[1,2]+table[2,2]) /
(sum(dati)*sum(dati)*(sum(dati)-1))
var.a
[1] 65.4635

chi.square <- ((abs(table[1,1]-E.a)-0.5)^2) / var.a
chi.square
[1] 60.6639


Anche con il test del chi-quadro corretto, viene rifiutata l'ipotesi nulla: l'OR è significativo.




Calcoliamo ora l'intervallo di confidenza dell'odd-ratio con il metodo di Wolff. Per far questo ci occorre calcolare l'errore standard del logaritmo naturale dell'OR. Semplicemente questo valore è pari alla somma dei reciproci delle celle, sotto radice quadrata. Quindi in R diventa:


esOR <- sqrt(sum( 1/dati ))


L'intervallo di confidenza con il metodo di Wolff (alpha = 0.05), è:


ciMin <- log(OR) - (1.96 * esOR)
> ciMin
[1] 0.7562176

ciMax <- log(OR) + (1.96 * esOR)
> ciMax
[1] 1.273949

ORmin <- exp(ciMin)
ORmax <- exp(ciMax)

> ORmin
[1] 2.130204
> ORmax
[1] 3.574942


Quindi l'odd-ratio risulta essere compreso tra 2.13 e 3.57.

Possiamo calcolare inoltre l'intervallo di confidenza dell'OR con il metodo di Miettinen:



Il chi quadro che qui compare è il chi-quadro corretto già calcolato; Z è il calore critico per alpha stabilito. In R i due limiti dell'intervallo di confidenza sono:


espMin <- 1 - 1.96/(sqrt(chi.square))
espMax <- 1 + 1.96/(sqrt(chi.square))

mORmin <- OR^espMin
mORmax <- OR^espMax

> mORmin
[1] 2.137509

> mORmax
[1] 3.562724


Con il metodo di Miettinen, l'OR risulta compreso tra 2.14 e 3.56, valori molti vicini all'intervallo di confidenza calcolato con il metodo di Wolff.




Possiamo calcolare l'OR anche partendo dalla regressione logistica semplice, OR è pari all'esponenziale del regressore beta1:


table

malSI malNO
InquiSI 235 525
inquiNO 103 635

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

Var1 Var2 Freq
1 InquiSI malSI 235
2 inquiNO malSI 103
3 InquiSI malNO 525
4 inquiNO malNO 635

fit <- glm(Var2 ~ Var1, weights = Freq, data = dft, family = binomial(logit))
OR <- exp(fit$coefficient[2])
OR

Var1inquiNO
2.759593


Per calcolare l'intervallo di confidenza con il metodo di Wolff, si utilizza la funzione confint():


ci <- confint(fit)

2.5 % 97.5 %
(Intercept) 0.6514110 0.9592309
Var1inquiNO 0.7588576 1.2770034


Calcoliamo l'esponenziale dei valori ottenuti in questa tabella, e otteniamo l'intervallo di confidenza:


> exp(ci[2,1])
[1] 2.135835

> exp(ci[2,2])
[1] 3.585878


L'intervallo di confidenza è (2.13, 3.58). Possiamo specificare anche il metodo con cui calcolare l'intervallo di confidenza: confint(fit, method="miettinen"), oppure confint(fit, method="wolff").




Il test Z utilizzato prima per verificare la significatività dell'OR può anche essere utilizzato per confrontare tra loro gli OR di due tabelle. Gli OR vengono calcolati indipendentemente; l'ipotesi nulla è OR1 = OR2, e la statistica test di confronto è la seguente:



La varianza dell'odd ratio è pari alla somma dei reciproci delle celle.
Supponiamo di avere le seguenti due generiche tabelle, di cui si vogliono confrontare gli odd-ratios:



Immettiamo i dati in R, calcoliamo gli OR, quindi le loro varianze, e infine il valore Z:


dati1 <- c(47, 21, 62, 77)
tabella1 <- matrix(dati1, nrow=2)
colnames(tabella1) <- c("SI", "NO")
rownames(tabella1) <- c("A", "B")
table1 <- as.table(tabella1)

dati2 <- c(54, 21, 95, 61)
tabella2 <- matrix(dati2, nrow=2)
colnames(tabella2) <- c("SI", "NO")
rownames(tabella2) <- c("A", "B")
table2 <- as.table(tabella2)

dft1 <- as.data.frame(table1)
dft2 <- as.data.frame(table2)

fit1 <- glm(Var2 ~ Var1, weights = Freq, data = dft1, family = binomial(logit))
OR1 <- exp(fit1$coefficient[2])

fit2 <- glm(Var2 ~ Var1, weights = Freq, data = dft2, family = binomial(logit))
OR2 <- exp(fit2$coefficient[2])

varOR1 <- sum(1/dati1)
varOR2 <- sum(1/dati2)

#ricordo che ln(OR) = coefficiente beta1 della regressione logistica
num <- fit1$coefficient[2] - fit2$coefficient[2]
num <- OR1-OR2
den <- sqrt(varOR1 + varOR2)
Z.OR <- num / den

pnorm(-abs(Z.OR))
[1] 0.004917582


Il p-value è 0.005, inferiore all'alpha prefissato di 0.05. Quindi rifiutiamo l'ipotesi nulla: gli OR delle due tabelle sono significativamente diversi.

Nessun commento:

Posta un commento