Un altro modo per valutare l'associazione è quello di sviluppare un modello di regressione, in questo caso un modello di regressione logistica multipla, poiché Y, X1 e X2 sono bernoulliane (1 o 0). Abbiamo già visto, in quest'altro post, come creare un modello di regressione logistica semplice. Cerchiamo di applicare gli stessi principi anche ai dati che abbiamo già analizzato. In virtù dell'analisi precedente, ci aspettiamo che la variabile dipendente Y dipenda con basso grado dalle due variabili indipendenti X1 (fumo) e X2 (alcol), ossia ci aspettiamo di ottenere dei bassi valori beta1 e beta2. Cerchiamone la conferma.
La tabella sottostante riassume le due tabelle già considerate:
(la tabella è orribile, ma ho qualche problemino con le tabelle in LaTeX !!).
Per prima cosa creiamo un dataframe con tutti i dati, nel modo seguente:
fumo <- c(rep("NO", 2), rep("SI", 2))
alcol <- c(rep("NO", 1),rep("SI", 1))
tumore <- c(rep("NO", 4),rep("SI", 4))
freq <- c(79,16,115,5,76,22,98,8)
data <- data.frame(fumo, alcol, tumore, freq)
> data
fumo alcol tumore freq
1 NO NO NO 79
2 NO SI NO 16
3 SI NO NO 115
4 SI SI NO 5
5 NO NO SI 76
6 NO SI SI 22
7 SI NO SI 98
8 SI SI SI 8
A questo punto possiamo costruire il nostro modello di regressione logistica multipla, con le stesse regole viste nel post sulla regressione logistica semplice. Vi è solo da aggiungere che quando si dispone di due variabili indipendenti, nella formula vengono entrambe dichiarate con un 'più', come potete vedere qui sotto:
fit <- glm(tumore ~ fumo + alcol, weights = freq, data = data, family = binomial)
> summary(fit)
Call:
glm(formula = tumore ~ fumo + alcol, family = binomial, data = data,
weights = freq)
Deviance Residuals:
1 2 3 4 5 6 7 8
-10.263 -5.372 -11.955 -2.909 10.464 4.786 12.289 2.995
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.05363 0.15620 -0.343 0.731
fumoSI -0.09542 0.20109 -0.475 0.635
alcolSI 0.43477 0.30988 1.403 0.161
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 580.57 on 7 degrees of freedom
Residual deviance: 577.95 on 5 degrees of freedom
AIC: 583.95
Number of Fisher Scoring iterations: 4
La formula del modello di regressione logistica multipla che stiamo cercando è:
Quindi la probabilità che Y assuma valore 1 (presenza di tumore), in funzione delle due variabili dipendenti X1 e X2, è data da quella formula. I 3 coefficienti sono stati stimati, e sono pari a:
Quindi la formula è:
Essendo noti questi parametri, possiamo calcolare diversi valori. Innanzitutto le varia probabilità.
La probabilità di avere tumore in assenza di fumo e alcol è:
La probabilità di avere tumore in presenza di fumo, ma non di alcol, è:
La probabilità di avere tumore in presenza di alcol, ma non di fumo, è:
La probabilità di avere tumore in presenza di alcol e fumo è:
Inoltre possiamo calcolare gli odds, con le seguenti formule:
Infine possiamo stimare gli odds ratio delle due variabili:
Dopo questo piccolo riassunto di formule, un paio di considerazioni. Innanzitutto il modello qui ricreato non è statisticamente significativo, essendo i p-value di tutt'e 3 i coefficienti maggiore di 0.05. Ho proceduto nei calcoli, perchè serva come esempio. In secondo luogo, vediamo che i 3 coefficienti beta hanno valori molto piccoli (e piccoli sono i valori degli OR), proprio come ci aspettavamo dall'analisi col chi-quadro di Mantel-Haenszel. Infine ricordo che tutti i valori riportati in tabella sono assolutamente inventati così come le conclusioni (paradossali: il fumo qui appare un fattore di protezione contro il tumore, cosa assolutamente non vera!).
Per calcolare tutte le formule presentate in R, ricordo come richiamare i valori dei coefficienti (così da poter essere utilizzati nelle formule):
coef(fit)[1], oppure fit$coefficient[1]coef(fit)[2], oppure fit$coefficient[2]coef(fit)[3], oppure fit$coefficient[3]Se abbiamo modo di ritenere che le due variabili indipendenti X1 e X2 interagiscano tra di loro, è possibile effettuare un modello di regressione logistica multipla con interazioni. In R questo si realizza semplicemente sostituendo al 'più', il segno moltiplicativo nella formula del modello. Vediamo subito cosa otteniamo:
fit <- glm(tumore ~ fumo * alcol, weights = freq, data = data, family = binomial)
> summary(fit)
Call:
glm(formula = tumore ~ fumo * alcol, family = binomial, data = data,
weights = freq)
Deviance Residuals:
1 2 3 4 5 6 7 8
-10.319 -5.261 -11.906 -3.091 10.408 4.904 12.335 2.787
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03871 0.16067 -0.241 0.810
fumoSI -0.12125 0.21146 -0.573 0.566
alcolSI 0.35717 0.36575 0.977 0.329
fumoSI:alcolSI 0.27280 0.69114 0.395 0.693
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 580.57 on 7 degrees of freedom
Residual deviance: 577.79 on 4 degrees of freedom
AIC: 585.79
Number of Fisher Scoring iterations: 4
In questo caso abbiamo ottenuto 4 coefficienti beta, che sono quelli relativi alla seguente formula:
In parole semplici, se questo modello fosse veritiero, significherebbe che la presenza simulatanea di fumo E alcol aumenta la probabilità di contrarre tumore molto più della semplice somma (le due condizioni aggravano vicendevolmente la salute dell'individuo).
Le formule e le considerazioni fatte precedentemente sono valide anche qui.
