6.6 Applicazioni

La creazione della tabella dati può essere fatta come indicato di seguito, da una tabella creata con il comando matrix ():

>tab < – as.ping(matrix(c(815,115,208,327), nrow = 2, byrow = TRUE,

dimnames = list(EST = c(“+”,”−”),

CAD = c(“+”,”-”))))

> scheda

Si noti che la tabella è stata riorganizzata in modo da corrispondere alle solite notazioni, con gli eventi “positivi” sulla prima riga (in genere, l’esposizione) e la prima colonna (in genere, la malattia) della tabella.

Sarebbe quindi possibile calcolare tutte le quantità richieste utilizzando pochissime operazioni aritmetiche con R. Tuttavia, il pacchetto epiR contiene tutti gli strumenti necessari per rispondere alle domande epidemiologiche nelle tabelle 2 × 2. Quindi, il comando epi.tests () fornisce valori di prevalenza, sensibilità / specificità, nonché valori predittivi positivi e negativi:

>library(epiR)

> epi.test (tab)

Per quanto riguarda l’area sotto la curva, è possibile utilizzare anche un comando esterno, roc.da.tabella(), nel pacchetto epicalc:

> libreria(epicalc)

> roc.da.tabella (tab, grafico = FALSO)

Importazione dei dati hdis.dat è raggiunto come di seguito:

> bp< – read.tabella (“hdis.dat”, header = TRUE)

> str(bp)

Come può essere verificato, non esiste un’etichetta associata ai livelli delle variabili di interesse (bpress per la pressione sanguigna, chol per il tasso di colesterolo). Per generare e associare etichette, è possibile inserire i seguenti comandi:

> blab < – c(“<117”,”117-126”,”127-136”,”137-146”,

“147-156”,”157-166”,”167-186”,”>186”)

> clab < – c(“<200”,”200-209”,”210-219”,”220-244”,

“245-259”,”260-284”,”>284”)

> bp$bpress < – factor(bp$bpress, labels = blab)

> bp$chol < – factor(bp cho chol, labels = clab)

L’ultimo comando converte le variabili originali in variabili qualitative e associa alle loro modalità le etichette definite da blab e clab. Per verificare che il database sia ora nella forma desiderata, sono utili i seguenti comandi:

> str(bp)

> summary(bp)

La tabella delle frequenze relative esposte nell’istruzione può ora essere riprodotta:

>round(xtabs(hdis/total bpress + chol, data = bp), 2)

Poiché ci concentreremo solo sulla relazione tra pressione sanguigna e infarto miocardico, è necessario aggregare i dati sui livelli di colesterolo. In altre parole, è necessario accumulare i valori per ciascun livello di pressione sanguigna, indipendentemente dai livelli di colesterolo. L’opportunità dovrebbe anche essere presa per rinominare i livelli della variabile bpress usando i centri degli intervalli di classe:

> blab2< – c(111.5,121.5,131.5,141.5,151.5,161.5,176.5,191.5)

> bp$bpress < rep(blab2, ogni = 7)

> dfrm < – aggregato(bp,

elenco(bpress = bp), somma)

e ‘ l’ultimo comando, di aggregazione(), che rende possibile l’aggregazione dei dati: è la somma di tutti i valori (somma) della variabile chol che non sono inclusi nell’elenco di variabili che devono essere mantenute per l’analisi. Allo stesso tempo, i risultati vengono memorizzati in un nuovo database denominato dfrm. Si noti che in questo caso non si fa uso della notazione della formula ma di quella della sintassi alternativa per aggregate () (variabile di risposta seguita da una lista composta da una o più variabili esplicative). Una panoramica dei dati aggregati può essere ottenuta dando head() come input:

> head(dfrm, 5)

Quando è disponibile una proporzione p, il suo valore su una scala le cui unità sono logits è dato dal log delle relazioni(p/(1 − p)). Da ciò i seguenti comandi possono essere derivati per convertire le proporzioni, calcolate come hdis / totale in unità logit:

> logit < – function(x) log(x/(1-x))

> dfrm$prop < dfrm$clive/dfrm$totale

> dfrm$logit < – logit(dfrm$clive/dfrm$totale)

Si deve osservare che una piccola funzione è stata definita per la conversione di valori x, che qui viene considerato un proporzioni, in un suo equivalente log(x/(1 − x)). Allo stesso modo, si potrebbe scrivere:

>log((dfrm hd hdis/dfrm total total)/(1-dfrm hd hdis/dfrm total total))

I valori calcolati sono stati aggiunti direttamente al frame dati con il nome prop e logit.

Il modello di regressione logistica è scritta nel modo seguente:

> glm(rbind(clive,totale-clive) bpress,dati = dfrm,family = binomial)

La formulazione utilizzata, rbind(clive, totale-clive) bpress, tiene conto del fatto che siamo l’accesso a dati raggruppati e non le singole risposte. Il comando glm () con l’opzione family = binomial corrisponde ad una regressione logistica, che, senza entrare molto nel dettaglio, utilizza per impostazione predefinita la funzione logit come collegamento canonico.

I coefficienti di regressione può essere visualizzato facendo uso di riepilogo():

> summary(glm(rbind(clive, totale-clive) bpress, dati = dfrm,

family = binomial))

Il risultato precedente comprende le informazioni essenziali per rispondere alla domanda della significatività statistica di associazione tra pressione arteriosa e il rischio di attacco di cuore: il coefficiente di regressione (sulla scala log-odds) è pari a 0,024 ed è significativo alla soglia usuale del 5% (vedi colonna Pr (>|z|)).

La probabilità di avere un attacco di cuore in base ai diversi livelli di pressione sanguigna in esame si ottiene come segue:

> glm.res < – glm(cbind(hdis, total-hdis) bpress, data = dfrm,

famiglia = binomio)

> predict(glm.res)

Va notato che i risultati intermedi generati da R sono stati memorizzati con il nome glm.res prima di utilizzare il comando predict (). Le previsioni generate da R sono espresse in forma logit ed è possibile confrontare i logit osservati e previsti:

> cbind(dfrm, logit.predit = predire(glm . res))

Praticamente tutti gli elementi sono disponibili per rappresentare le proporzioni di infarto miocardico osservate e previste graficamente in base al livello di pressione sanguigna. Mancano le previsioni del modello espresse sotto forma di proporzioni e non sulla scala logodds. Può quindi essere auspicabile tracciare la curva di previsione, vale a dire la probabilità di avere un infarto in base alla pressione sanguigna, senza limitare quest’ultima agli otto valori osservati per la variabile bpress. Di seguito una possibile soluzione:

> dfrm prop prop.predit < – predict (glm.la funzione(x)

> f < – function(x)

1 / (1 + exp (- (coef (glm.res) + coef(glm.res)*x)))

> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,

xlab = “Blood pressure”,ylab = “Infarction likelihood”,

panel = function(x, y, …) {

panel.xyplot(x, y, col = “gray30”, pch = 19, …)

panel.curve(f, lty = 3, col = “gray70”)

panel.points(x, dfrm$prop.predit, col = “gray70”, …)

})

Poiché i dati grezzi (individuali) non sono disponibili, è necessario lavorare direttamente con la tabella delle frequenze fornita nell’istruzione del problema:

> alcohol < – matrix(c(666,104,109,96), nr = 2,

dimnames = list(c(“Control”,”Case”),

c(“<> = 80”)))

> alcohol

Regarding the odds ratio, the command oddsratio() from the package vcd will be used:

>library(vcd)

>oddsratio(alcohol, log = FALSE)

L’opzione log = FALSE assicura che il risultato restituito corrisponda correttamente a un odds ratio e non al log del odds ratio. Per ottenere un intervallo di confidenza asintotico verrà impiegato quanto segue:

>confint(oddsratio(alcohol, log = FALSE))

Allo stesso modo, summary(oddsratio(alcohol)) potrebbe essere usato per eseguire un test di ipotesi sul log odds ratio (H0 : log(θ) = 0).

Il comando prop.test () può essere utilizzato per verificare l’ipotesi che la percentuale di persone con un consumo giornaliero sia ≥ 80 g ed è identica nei casi e nei controlli, indicando i valori osservati dalla tabulazione incrociata riportata nell’istruzione del problema:

> prop.test (c (96,109), c(200,775), correct = FALSE)

Questo test è esattamente equivalente a un test Z per testare la differenza tra due proporzioni stimate dai dati (se non viene utilizzata alcuna correzione di continuità).

Per quanto riguarda il modello di regressione logistica, la tabella di contingenza deve essere trasformata in una tabella di dati in cui le due variabili qualitative (malattia ed esposizione) sono rappresentate correttamente, vale a dire con un frame di dati. Questo può essere fatto usando il comando melt () dal pacchetto reshape2 che rende possibile trasformare i dati in forma tabellare in formato “lungo”. In questa fase, i livelli della variabile di malattia deve essere registrato con 0 (controllo) e 1 (casi), anche se questo non è davvero necessario e il livello 0 dovrebbe essere considerata come la categoria di riferimento (che facilita l’interpretazione dei risultati):

> library(reshape2)

> alcol.df < – melt(alcohol)

> nomi(alcohol.df)< – c(“malattia”, “esposizione”, “n”)

> livelli(alcol.df disease malattia)< – c(1,0)

> alcol.df disease disease < – relevel (alcohol.df disease disease,”0″)

Il modello di regressione logistica viene quindi scritto nel modo seguente:

> glm.res < – glm(esposizione alla malattia, dati = alcol.df,

famiglia = binomio, pesi = n)

> sommario (glm.res)

Il risultato di interesse è la riga associata con l’esposizione > = 80, in quanto fornisce informazioni circa il valore del coefficiente di regressione associato alla variabile di esposizione, come stimato da R, con il suo errore standard, così come il valore della statistica test. Qui, il coefficiente di regressione viene interpretato come il registro del rapporto di probabilità. Si noti che si otterrebbe esattamente gli stessi risultati scambiando il ruolo delle variabili nella formulazione precedente, malattia da esposizione.

Il rapporto di probabilità calcolato sopra può essere trovato dal coefficiente di regressione associato al fattore di interesse (esposizione), così come il suo intervallo di confidenza del 95%, nel modo seguente:

> exp(coef(glm.res))

> exp (confint (glm.res))

La seconda soluzione consiste nel considerare il numero di casi e il numero totale di individui:

>alcohol2< – dati.telaio(expos = c(“<> = 80”), cassa = c(104,96),

totale = c(104 + 666, 96 + 109))

> summary(glm(rbind(caso, totale-caso) fiere, dati = alcohol2,

family = binomial))

4) i dati arrivano da uno studio che cerca di stabilire la prognosi validità di creatina chinasi concentrazione nel corpo, relativa alla prevenzione dell’insorgenza di infarto del miocardio .

Poiché il nome delle variabili non è nel file di dati, sarà necessario assegnarle immediatamente dopo l’importazione dei dati:

>sck < – read.tabella (“sck.dat”, header = FALSE)

> nomi(sck) < – c(“ck”, “pres”, “abs”)

> summary(sck)

Il numero totale di soggetti corrisponde alla somma delle frequenze per le due variabili pres e abs, che è:

> sum(sck)

o, equivalentemente: sum (sck pre pres) + sum (sck abs abs) (ma non sck sc pres + sck abs abs).

Per calcolare le frequenze relative di queste due variabili, è necessario conoscere le quantità totali per variabile. Questi possono essere ottenuti utilizzando il comando applica e operativo da colonne:

> ni < applica(sck, 2, somma)

sulla Base di questo, è sufficiente dividere ogni valore delle variabili pres e abs per i valori calcolati in precedenza. I valori ottenuti verranno memorizzati in due nuove variabili nella stessa tabella di dati:

> sck pre pres.prop < – sck$pres/ni

> sck$abs.prop < – sck$abs/ni

è possibile verificare che i calcoli sono corretti: la somma dei valori per ogni variabile deve ora essere pari a 1:

> applica(sck, 2, somma)

Quindi le proporzioni ottenuti sono semplicemente rappresentato in un unico grafico, considerando i valori della variabile ck come ascisse:

> xyplot(pres.prop + abs.prop ck, dati = sck, tipo = c (”b”,”g”),

auto.key = TRUE, ylab = “Frequency”)

L’istruzione type = c (”b”,” g”) significa che vogliamo che i punti vengano visualizzati collegati da linee (”b “= ” o ” + “l’) e una griglia (”g”).

Il modello di regressione per i dati raggruppati è scritto come:

> glm.res < – glm(cbind (pres,abs) ck, data = sck, family = binomiale)

> sommario (glm.res)

Le previsioni, espresse in forma di probabilità e non sulla scala log-odds, sono ottenute mediante il comando predict specificando l’opzione type = “response” come segue:

> glm.pred < – predict (glm.res, type = “response”)

> nomi (glm.pred) < – sck c ck

Considerando che le probabilità ≥ 0,5 designano individui malati, si ottiene così la seguente distribuzione:

> glm.pred

Si conclude che le persone saranno considerate malate, secondo questo modello, per valori ck di 80 o più. Questo è facilmente verificabile con un grafico in cui la probabilità stimabile sulla base di valori della variabile di ck sono rappresentati:

> sck$malati < – sck$pres/(sck$pres + sck$abs)

> xyplot(glm.in questo caso, la funzione è quella di utilizzare il pannello di controllo.xyplot (panel)

pannello.xyplot (sck c ck, sck sick sick, pch = 19, col = “grey”)

})

In una prima fase, è necessario “decomprimere” i dati raggruppati e creare una tabella con due colonne: la prima che rappresenta la variabile ck e la seconda che rappresenta la presenza o l’assenza di malattia. Useremo i conteggi per sottogruppo precedentemente calcolati e disponibili nella variabile ni:

> sck.espandere < – dati.la struttura(ck = c(rep(sck c ck, sck pre pres),

rep(sck c ck, sck abs abs)),

sick = c(rep(1,ni), rep(0,ni)))

> tabella(sck.per maggiori informazioni clicca qui.expand, tapply (sick, ck, sum))

Gli ultimi due comandi hanno lo scopo di garantire che vengano rilevati gli stessi conteggi e che la distribuzione delle persone malate in base al valore di ck sia corretta. Si può anche verificare che si ottengano gli stessi risultati per quanto riguarda la regressione logistica e allo stesso tempo si possano aggiungere i valori previsti nella tabella dati precedente (si potrebbe seguire la stessa procedura di prima e replicare le previsioni, ma questo è più semplice in questo modo):

> glm.res2 < – glm(ck malato, dati = sck.espandi, famiglia = binomio)

> sck.expand prediction prediction < – ifelse (predire (glm.res2,

type = “response”) > = 0.5,

1, 0)

> con(sck.expand, table (sick, prediction))

L’ultimo comando consente di visualizzare una matrice di confusione in cui vengono incrociate la diagnostica reale e quelle previste dal modello di regressione. I tassi di classificazione corretti possono essere confrontati quando variiamo la soglia di riferimento:

> classif.scheda < – con(sck.expand, table (sick, prediction))

> sum (diag (classif.tab))/sum(classif.tab)

Per visualizzare la curva ROC, utilizzeremo il pacchetto ROCR:

>library(ROCR)

> pred < – prediction(predict(glm.res2, tipo = “risposta”),

sck.expand$sick)

> perf < – performance(pred, “tpr”,”fpr”)

> plot(perf, ylab = “Sensitivity”, xlab = “1-Specificity”)

> grid()

> abline(0, 1, lty = 2)

The value of the area under the curve is obtained as follows:

> performance(pred, “auc”)@”y.values”