6.6 Aplicații

crearea tabelului de date se poate face după cum este indicat în continuare, dintr – un tabel creat cu matricea de comandă ():

> tab< – as.tabel (matrice (c (815,115,208,327), nrow = 2, byrow = adevărat,

dimnames = listă (EST = c(„+”,”−”),

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

> tab

trebuie remarcat faptul că tabelul a fost reorganizat astfel încât să corespundă notațiilor obișnuite, cu evenimente „pozitive” pe primul rând (de obicei, expunerea) și prima coloană (de obicei, boala) a tabelului.

atunci ar fi posibil să se calculeze toate cantitățile necesare folosind foarte puține operații aritmetice cu R. Cu toate acestea, pachetul epiR conține toate instrumentele necesare pentru a răspunde la întrebările epidemiologice din cele 2 tabele de la 2 la sută. Astfel, comanda epi.teste () oferă valori de prevalență, sensibilitate / specificitate, precum și valori predictive pozitive și negative:

> bibliotecă(epiR)

> epi.teste (tab)

în ceea ce privește zona de sub curbă, se poate folosi și o comandă externă, roc.din.tabelul(), în pachetul epicalc:

> biblioteca(epicalc)

> Roc.din.tabel (tab, grafic = fals)

importul hdis de date.dat se realizează după cum urmează:

> bp< – read.tabel („hdis.dat”, header = TRUE)

> str(bp)

după cum se poate verifica, nu există nicio etichetă asociată cu nivelurile variabilelor de interes (bpress pentru tensiunea arterială, chol pentru rata colesterolului). Pentru a genera și asocia etichete, pot fi inserate următoarele comenzi:

> 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$chol, etichete = clab)

Ultima comandă convertește variabilele originale în variabile calitative și asociază la modalitățile lor etichetele definite de blab și clab. Pentru a verifica dacă baza de date este acum în forma dorită, următoarele comenzi sunt utile:

> str(bp)

> rezumat(bp)

tabelul frecvențelor relative expuse în declarație poate fi reprodus acum:

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

deoarece ne vom concentra doar pe relația dintre tensiunea arterială și infarctul miocardic, este necesară agregarea datelor privind nivelurile de colesterol. Cu alte cuvinte, este necesar să se acumuleze valorile pentru fiecare nivel al tensiunii arteriale, indiferent de nivelul colesterolului. De asemenea,trebuie luată posibilitatea de a redenumi nivelurile variabilei bpress folosind centrele intervalelor de clasă:

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

>BP$bpress < – rep(blab2, fiecare = 7)

>dfrm < – agregat(bp,

list(bpress = BP), sum)

este ultima comandă, aggregate(), care face posibilă agregarea datelor: însumează toate valorile (suma) variabilei chol care nu sunt incluse în lista variabilelor care trebuie păstrate pentru analiză. În același timp, rezultatele sunt stocate într-o nouă bază de date numită dfrm. Rețineți că în acest caz nu folosim notația formulei, ci cea a sintaxei alternative pentru agregat() (variabilă de răspuns urmată de o listă formată din una sau mai multe variabile explicative). O prezentare generală a datelor agregate poate fi obținută dând head () ca intrare:

> head(dfrm, 5)

atunci când o proporție p este disponibilă, valoarea sa pe o scară ale cărei unități sunt logits este dată de Jurnalul de relații(p/(1 − p)). Din acestea pot fi derivate următoarele comenzi pentru a converti proporțiile, calculate ca hdis / total în unități logit:

> logit < – funcție(x) log(x/(1-x))

> dfrm$prop < – dfrm$hdis/dfrm$total

> dfrm$logit < – logit(dfrm$hdis/dfrm$total)

trebuie observat că o funcție mică a fost definită pentru conversia valorilor X, care aici se presupune a fi proporții, în jurnalul său echivalent(x/(1 − x)). De asemenea, s-ar putea scrie:

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

Valorile calculate au fost adăugate direct în cadrul de date sub numele prop și logit.

modelul de regresie logistică este scris în felul următor:

> glm(cbind(hdis,total-hdis) bpress,data = dfrm,family = binomial)

formularea utilizată, cbind(hdis, total-hdis) bpress, ține cont de faptul că accesăm date grupate și nu răspunsuri individuale. Comanda glm () cu opțiunea family = binomial corespunde unei regresii logistice, care, fără a intra prea mult în detalii, folosește în mod implicit funcția logit ca legătură canonică.

coeficienții de regresie pot fi afișați folosind rezumatul():

> rezumat(glm(cbind(hdis, total-hdis) bpress, data = dfrm,

family = binomial))

rezultatul anterior include informațiile esențiale pentru a răspunde la întrebarea semnificației statistice a asocierii dintre tensiunea arterială și probabilitatea atacului de cord: coeficientul de regresie (pe scara log-odds) este egal cu 0,024 și este semnificativ la pragul obișnuit de 5% (a se vedea coloana Pr(>|z|)).

probabilitatea de a avea un atac de cord în funcție de diferitele niveluri ale tensiunii arteriale luate în considerare se obține după cum urmează:

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

familie = binom)

> prezice(glm.res)

trebuie remarcat faptul că rezultatele intermediare generate de R au fost stocate cu numele glm.res înainte de a utiliza comanda prezice(). Predicțiile generate de R sunt exprimate în formă logit și este posibil să se compare logitele observate și prezise:

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

practic toate elementele sunt disponibile pentru a reprezenta proporțiile infarctului miocardic observate și prezise grafic în funcție de nivelul tensiunii arteriale. Predicțiile modelului exprimate sub formă de proporții și nu pe scara logodds lipsesc. Prin urmare, poate fi de dorit să se traseze curba de predicție, adică probabilitatea de a avea un atac de cord în funcție de tensiunea arterială, fără a se limita la cele opt valori observate pentru variabila bpress. Aici urmează o posibilă soluție:

> dfrm$prop.predit < – prezice(glm.res, type = „response”)

> 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”,…)

})

deoarece datele brute (individuale) nu sunt disponibile, este necesar să lucrați direct cu tabelul de frecvență furnizat în instrucțiunea problemei:

> 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(alcool, log = FALSE)

opțiunea log = FALSE asigură că rezultatul returnat corespunde corect unui raport de cote și nu jurnalului raportului de cote. Pentru a obține un interval de încredere asimptotic se vor folosi următoarele:

> confint(oddsratio(alcool, log = fals))

în mod similar, summary(oddsratio(alcool)) ar putea fi utilizat pentru a efectua un test de ipoteză privind raportul cotelor log (H0 : log(inkt) = 0).

comanda prop.testul () poate fi folosit pentru a testa ipoteza că proporția persoanelor cu un consum zilnic este de 80 g și este identică atât în cazuri, cât și în controale, indicând valorile observate din tabularea încrucișată dată în enunțul problemei:

> prop.test(c (96.109), c(200.775), corect = fals)

Acest test este exact echivalent cu un test Z pentru a testa diferența dintre două proporții estimate din date (dacă nu se folosește nicio corecție de continuitate).

în ceea ce privește modelul de regresie logistică, tabelul de contingență trebuie transformat într-un tabel de date în care cele două variabile calitative (boala și expunerea) sunt reprezentate corespunzător, adică cu un cadru de date. Acest lucru se poate face folosind comanda melt() din pachetul reshape2 care face posibilă transformarea datelor în formă tabelară în format „lung”. În acest stadiu, nivelurile bolii variabile ar trebui recodificate cu 0 (control) și 1 (cazuri), deși acest lucru nu este cu adevărat necesar, iar nivelul 0 ar trebui considerat drept categoria de referință (care facilitează interpretarea rezultatelor):

> bibliotecă(reshape2)

> alcool.df < – topiți(alcool)

> nume(alcool.df) < – c(„boală”, „expunere”, „n”)

> niveluri(alcool.DF $ boala) < – c(1,0)

> alcool.DF$boala < – relevel(alcool.DF$boala, „0”)

modelul de regresie logistică este apoi scris în felul următor:

> glm.res < – glm(expunere la boală, date = alcool.DF,

familie = binom, greutăți = n)

> rezumat(glm.res)

rezultatul interesului este rândul asociat expunerii> = 80 deoarece oferă informații despre valoarea coeficientului de regresie asociat variabilei de expunere estimată de R, cu eroarea sa standard, precum și valoarea statisticii testului. Aici, coeficientul de regresie este interpretat ca Jurnalul raportului de cote. Rețineți că s-ar obține exact aceleași rezultate prin schimbarea rolului variabilelor în formularea anterioară, boala de expunere.

raportul de cote calculat mai sus poate fi găsit din coeficientul de regresie asociat factorului de interes (expunere), precum și din intervalul său de încredere de 95%, în felul următor:

> exp(coef(glm.res))

> exp(confint(glm.a doua soluție constă în luarea în considerare a numărului de cazuri și a numărului total de persoane:

> alcohol2 < – date.cadru (expos = c („<> = 80”), caz = c (104,96),

total = C(104 + 666, 96 + 109))

> rezumat(glm(cbind(caz, caz total) expos, date = alcool2,

familie = binom))

4) aceste date provin dintr-un studiu care urmărește stabilirea valabilității prognostice a concentrației creatinkinazei în organism cu privire la prevenirea apariției infarctului miocardic .

deoarece numele variabilelor nu se află în fișierul de date, va fi necesar să le atribuiți imediat după importul datelor:

> sck< – read.tabel („sck.dat”, header = FALSE)

> nume(sck) < – c(„ck”, „pres”, „abs”)

> rezumat(sck)

totalul numărul de subiecți corespunde sumei frecvențelor pentru cele două variabile pres și ABS, adică:

> sumă(SCK)

sau, echivalent: sumă ( sck$pres) + sumă(sck$abs) (dar nu sck$pres + sck$abs).

pentru a calcula frecvențele relative ale acestor două variabile, este necesar să cunoaștem cantitățile totale pe variabile. Acestea pot fi obținute folosind comanda apply și funcționând pe coloane:

> ni< – apply(sck, 2, sum)

pe baza acestui lucru, este suficient să împărțiți fiecare valoare a variabilelor pres și abs la valorile calculate anterior. Valorile obținute vor fi stocate în două variabile noi în același tabel de date:

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

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

este posibil să se verifice dacă calculele sunt corecte: suma valorilor pentru fiecare variabilă trebuie să fie acum egală cu 1:

> apply(sck, 2, sum)

atunci proporțiile obținute sunt pur și simplu reprezentate într-o singură diagramă, având în vedere valorile variabila CK ca abscisă:

> xyplot(pres.prop + abs.prop ck, date = sck, tip = c („b”,”g”),

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

instrucțiunea type = C(„b”, „g”) înseamnă că dorim ca punctele să fie afișate conectate prin linii („b” = „o” + „l’), precum și o grilă („g”).

modelul de regresie pentru datele grupate este scris ca:

> glm.res < – glm(cbind(pres,abs) ck, date = sck, familie = binom)

> rezumat(glm.res)

predicțiile, exprimate sub formă de probabilități și nu pe scara log-odds, sunt obținute prin intermediul comenzii prezic prin specificarea opțiunii type = „response” după cum urmează:

> glm.pred < – prezice(glm.res, type= „response”)

> nume(glm.pred) < – sck$ck

având în vedere că probabilitățile 0.5 0.5 desemnează persoanele bolnave, se obține astfel următoarea distribuție:

> glm.pred

se concluzionează că oamenii vor fi considerați bolnavi, conform acestui model, pentru valori ck de 80 sau mai mult. Acest lucru este ușor de verificat cu o diagramă în care probabilitățile prezise pe baza valorilor variabilei ck sunt reprezentate:

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

> xyplot(GLM.pred ~ sck$ck, tip = „l”,

ylab = „probabilitate”, xlab = „ck”,

panou = funcție(…) {

panou.xyplot ( … )

panou.xyplot (sck$ck, sck$sick, pch = 19, col = „grey”)

})

într-o primă etapă, este necesar să „decomprimați” datele grupate și să creați un tabel cu două coloane: prima reprezentând variabila ck și a doua reprezentând prezența sau absența bolii. Vom folosi numărările pe subgrupuri calculate anterior și disponibile în variabila ni:

> sck.expand < – date.frame (ck = c(rep(sck$ck, SCK$pres),

rep(sck$ck, sck$abs)),

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

> tabel (sck.expand$bolnav)

> cu(sck.expand, tapply (sick, ck, sum))

ultimele două comenzi sunt destinate să se asigure că se întâlnesc aceleași numărări și că distribuția bolnavilor în funcție de valoarea ck este corectă. De asemenea, se poate verifica dacă se obțin aceleași rezultate în ceea ce privește regresia logistică și în același timp se pot adăuga valorile prezise în tabelul de date anterior (aceeași procedură ca și anterior ar putea fi urmată și predicțiile ar putea fi reproduse, dar acest lucru este mai simplu în acest fel):

> glm.res2 < – glm(CK bolnav, date = sck.extinde, familie = binom)

> sck.expand $ predicție < – ifelse(prezice(glm.res2,

type = „response”)> = 0.5,

1, 0)

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

Ultima comandă face posibilă afișarea unei matrice de confuzie în care sunt încrucișate diagnosticul real și cele prezise de modelul de regresie. Ratele corecte de clasificare pot fi comparate atunci când variem pragul de referință:

> classif.tab < – cu(sck.expand, tabel(bolnav, predicție))

> sum(diag(classif.tab)) / sumă (classif.pentru a afișa curba ROC, vom folosi pachetul ROCR:

> bibliotecă(ROCR)

> pred< – predicție(prezice(glm.res2, type = „răspuns”),

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”