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