6.6 Aplikace

vytvoření tabulky dat může být provedeno, jak je dále uvedeno, z tabulky vytvořené pomocí příkazu matrix() :

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

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

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

> tab

Je třeba poznamenat, že tabulka byl reorganizován tak, aby odpovídala obvyklé notace, s „pozitivní“ události na první řádek (typicky, expozice) a prvním sloupci (typicky, nemoc) tabulky.

bylo by pak možné vypočítat všechny požadované veličiny pomocí velmi málo aritmetických operací s R. Balíček epiR však obsahuje všechny nástroje potřebné k odpovědi na epidemiologické otázky v tabulkách 2 × 2. To znamená, že příkaz epi.testy() poskytuje prevalence, citlivost/specificita hodnoty, stejně jako pozitivní a negativní prediktivní hodnoty:

> library(epiR)

> epi.testy (tab)

Pokud jde o oblast pod křivkou, lze také použít externí příkaz, roc.z.tabulka (), v balíčku epicalc:

> knihovna (epicalc)

> roc.z.tabulka(karta, graf = FALSE)

Import dat hdis.dat je dosaženo jako dále:

> bp < – číst.tabulka („hdis.dat“, header = TRUE)

> str(bp)

Jak to může být ověřena, je bez štítku spojené s úrovní proměnných zájmu (bpress pro krevní tlak, chol pro cholesterolu). Pro generování a přidružení štítků lze vložit následující příkazy:

> 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 < – faktor(bp$chol, popisky = clab)

poslední příkaz převádí původní proměnné do kvalitativních proměnných a společníků k jejich způsoby, štítky, definované žvanit a clab. K ověření, že databáze je nyní v požadované formě, následující příkazy jsou užitečné:

> str(bp)

> summary(bp)

tabulka relativní četnosti vystaveny v prohlášení, které nyní mohou být reprodukovány:

> round(xtabs(jednotce/celkem bpress + chol, data = bp), 2)

Protože my se zaměříme na vztah mezi krevním tlakem a infarktem myokardu, je nutné, aby souhrnné údaje o cholesterolu. Jinými slovy, je nutné akumulovat hodnoty pro každou hladinu krevního tlaku bez ohledu na hladinu cholesterolu. Příležitost by také měla být přijata, abyste přejmenovat úrovně proměnné bpress pomocí centrech třídy intervaly:

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

> bp$bpress < – rep(blab2, každý = 7)

> dfrm < – agregátní(bp,

seznam(bpress = bp), sum)

To je poslední příkaz, aggregate(), která umožňuje agregovat data: sumy všech hodnot (sum) proměnná chol, které nejsou zahrnuty do seznamu proměnných, které mají být zachovány pro analýzu. Současně jsou výsledky uloženy v nové databázi s názvem dfrm. Všimněte si, že v tomto případě nepoužíváme zápis vzorce, ale alternativní syntaxi pro agregát () (proměnná odezvy následovaná seznamem sestávajícím z jedné nebo více vysvětlujících proměnných). Přehled souhrnných údajů lze získat tím, že hlavy() jako vstupní:

> hlava(dfrm, 5)

Když se podíl p je k dispozici, její hodnota na stupnici, jejíž jednotky jsou logits je dána vztahem log(p/(1 − p)). Z toho lze odvodit následující příkazy pro převod proporcí, vypočtených jako hdis / total v jednotkách logit:

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

> dfrm$prop < – dfrm$jednotce/dfrm$celkem

> dfrm$logit < – logit(dfrm$jednotce/dfrm$celkem)

Je třeba poznamenat, že malé funkce byla definována pro převod hodnot x, které je zde předpokládá, že v poměru, v jeho ekvivalent log(x/(1 − x)). Stejně tak by se dalo psát:

> log((dfrm$jednotce/dfrm$celkem)/(1-dfrm$jednotce/dfrm$celkem))

vypočtené hodnoty byly přidány přímo do dat snímku pod názvem prop a logit.

V modelu logistické regrese je psán následujícím způsobem:

> glm(cbind(jednotce,celková-jednotce) bpress,data = dfrm,rodina = binomické)

složení používán, cbind(jednotce, celková-jednotce) bpress, bere v úvahu skutečnost, že jsme přístup k seskupené data a ne jednotlivé odpovědi. Příkaz glm() s možností rodiny = binomické odpovídá logistické regrese, který, aniž by se moc do detailů, používá ve výchozím nastavení na logit funkce jako kanonický link.

regresní koeficienty mohou být zobrazeny s využitím summary():

> summary(glm(cbind(jednotce, celková-jednotce) bpress, data = dfrm,

rodina = binomické))

předchozí výsledek obsahuje základní informace, odpovědět na otázku, statistická významnost asociace mezi krevní tlak a infarkt pravděpodobnost: regresní koeficient (na stupnici log-odds) se rovná 0,024 a je významný při obvyklém Prahu 5% (viz sloupec Pr (>|z|)).

pravděpodobnost infarktu podle různých hladin uvažovaného krevního tlaku se získá následujícím způsobem:

> glm.res < – glm(cbind(jednotce, celková-jednotce) bpress, data = dfrm,

rodina = binomické)

> predict(glm.res)

je třeba poznamenat, že mezilehlé výsledky generované R byly uloženy s názvem glm.res před použitím příkazu predict(). Předpovědi generované R jsou vyjádřeny v logit podobě a je možné porovnat pozorované a předpokládané logits:

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

Prakticky všechny prvky jsou k dispozici, představují rozměrů infarktu myokardu pozorované a předpokládané graficky podle úrovně krevního tlaku. Předpovědi modelu vyjádřené ve formě proporcí a nikoli na stupnici logodds chybí. Pravděpodobnost infarktu podle krevního tlaku, aniž by byla omezena na osm hodnot pozorovaných pro proměnnou bpress. Zde následuje možné řešení:

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

> f < – funkce(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“, …)

})

Od (individuální), raw data nejsou k dispozici, je nutné pracovat přímo s frekvenční tabulky uvedené v prohlášení problém:

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

možnost log = FALSE zajišťuje, že výsledek se vrátil odpovídá správně na odds ratio a ne log odds ratio. Získat asymptotický interval spolehlivosti následující bude zaměstnán:

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

Podobně, shrnutí(oddsratio(alkohol)) by mohly být použity k provedení test hypotézy o log odds ratio (H0 : log(θ) = 0).

příkaz prop.test() může být použita k testování hypotézy, že podíl osob s denní spotřebou je ≥ 80 g a je shodná v případech, stejně jako u kontrol, s uvedením hodnot, pozorované z cross-tabulation uvedeny v prohlášení problém:

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

Tento test je přesně stejné, jako Z-test pro testování rozdílu mezi dvěma poměry odhadnout z dat (pokud ne kontinuity korekce je zaměstnán).

S ohledem na modelu logistické regrese, kontingenční tabulce musí být transformována do datové tabulky, v nichž dvě kvalitativní proměnné (nemoc a expozice) jsou zastoupeny správně, to znamená, že se datový rámec. To lze provést pomocí příkazu melt () z balíčku reshape2, který umožňuje transformovat data v tabulkové formě do“ dlouhého “ formátu. V této fázi, úrovně proměnné onemocnění by měla být překódovat s 0 (kontrola) a 1 (případy), i když to není opravdu nutné a úroveň 0 by měla být považována za referenční kategorii (což usnadňuje interpretaci výsledků):

> library(reshape2)

> alkoholu.df < – melt(alcohol)

> names (alcohol.df) < – c(„nemoc“, „expozice“, „n“)

> úrovně(alkoholu.df$disease) < – c (1,0)

> alkohol.df$disease < – relevel(alcohol.df$nemoc, „0“)

logistické regresní model je pak zapsán následujícím způsobem:

> glm.res < – glm(expozice nemoci, data = alkohol.df,

rodina = binomické, závaží = n)

> summary(glm.res)

výsledkem zájmu je na řadě spojené s expozicí > = 80, protože poskytuje informace o hodnotě regresního koeficientu spojené s expozicí variabilní podle odhadu R, s její standardní chyba, stejně jako hodnota testovací statistika. Zde je regresní koeficient interpretován jako protokol poměru šancí. Všimněte si, že jeden by získal přesně stejné výsledky výměnou role proměnných v předchozí formulaci, expoziční nemoc.

šance poměr vypočítaný výše, lze nalézt z regresní koeficient asociovaný s faktorem zájmu (expozice), stejně jako jeho 95% interval spolehlivosti, v následujícím způsobem:

> exp(coef(glm.res))

> exp (confint (glm. res))

druhé řešení se skládá z ohledem na počet případů a celkový počet jedinců:

> alcohol2 < – data.rám(expos = c(„<> = 80“), case = c(104,96),

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

> summary(glm(cbind(pouzdro, celková-případ), veletrhů, data = alcohol2,

rodina = binomické))

4) Tyto údaje pochází ze studie snaží zjistit prognostickou platnost kreatin kinázy koncentrace v těle na prevenci vzniku infarktu myokardu .

Protože jméno proměnné není v datovém souboru, bude nutné přiřadit jim okamžitě poté, co údaje jsou importovány:

> sck < – číst.tabulka („sck.dat“, header = FALSE)

> jména(sck) < – c(„ck“, „pres“, „abs“)

> summary(sck)

celkový počet předmětů odpovídá součtu četností dvou proměnných pres a abs, to je:

> sum(sck)

nebo, ekvivalentně: součet (sck$pres) + součet (sck$abs) (ale ne sck$pres + sck$abs).

pro výpočet relativních frekvencí těchto dvou proměnných je nutné znát celkové množství proměnnou. Tyto mohou být získány pomocí příkazu použít a provozní podle sloupce:

> ni < apply(sck, 2, sum)

na tomto Základě, to stačí rozdělit každou hodnotu proměnné pres a abs hodnoty vypočítané dříve. Získané hodnoty budou uloženy ve dvou nových proměnných ve stejné datové tabulce:

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

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

je možné, ověřte, zda výpočty jsou správné: součet hodnot pro každou proměnnou musí nyní být rovna 1.

> použít(sck, 2, sum)

rozměry získané jsou jednoduše zastoupeny v jediném grafu, s ohledem na hodnoty proměnné ck jako abscissae:

> xyplot(pres.prop + abs.prop ck, data = sck, type = c („b“,“g“),

auto.key = TRUE, ylab = „Četnost“)

instrukce, typ = c(„b“, „g“) znamená, že chceme body, které mají být zobrazeny propojené linky („b“ = „o“ + „l‘), stejně jako mřížka („g“).

regresní model pro seskupená data je zapsán jako:

> glm.res < – glm(cbind(pres,abs, airbag řidiče) ck, data = sck, rodina = binomické)

> summary(glm.res)

předpovědi, vyjádřené ve formě pravděpodobnosti a ne na log-odds měřítku, jsou získány pomocí příkazu předpovídat pomocí parametru type = „response“ takto:

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

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

vzhledem k tomu, že pravděpodobnosti ≥ 0.5 určit jednotlivce, kteří jsou nemocní, následující rozdělení je tedy získat:

> glm.před

dospělo se k závěru, že lidé budou podle tohoto modelu považováni za nemocné pro hodnoty CK 80 nebo více. To je snadno ověřit pomocí grafu, kde pravděpodobnosti předpovědět na základě hodnoty proměnné ck jsou zastoupeny:

> sck$nemocný < – sck$pres/(sck$pres + sck$abs)

> xyplot(glm.pred ~ sck$ck, type = „l“,

ylab = „Pravděpodobnost“, xlab = „ck“,

panel = funkce ( … ) {

panel.xyplot ( … )

panel.xyplot(sck$ck, sck$nemocný, pch = 19, col = „grey“)

})

V prvním kroku, je nutné, aby „rozbalit“ skupinových dat a vytvořit tabulku se dvěma sloupci: první představuje proměnnou ck a druhá představuje přítomnost nebo nepřítomnost choroby. Použijeme počty podle podskupiny dříve vypočtené a dostupné v proměnné ni:

> sck.rozbalte < – data.rám(ck = c(rep(sck$ck, sck$pres),

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

nemocný = c(rep(1,ni), rep(0,ni)))

> table(sck.expand$sick)

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

poslední dva příkazy mají zajistit, aby se vyskytly stejné počty a aby rozdělení nemocných podle hodnoty ck bylo správné. To může být také ověřeno, že stejné výsledky jsou získány, pokud jde o logistické regrese a zároveň hodnoty předpokládané v předchozí tabulka dat může být přidáno (stejný postup jako předtím mohl být dodržovány, a předpovědi, mohla by být zopakována, ale je to jednodušší tímto způsobem):

> glm.res2 < – glm(sick ck, data = sck.expand, family = binomial)

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

type = „response“) > = 0.5,

1, 0)

> s(sck.rozbalit, tabulka (nemocný, predikce))

poslední příkaz umožňuje zobrazit matici zmatku, ve které jsou překročeny skutečné diagnostiky a diagnostiky předpovídané regresním modelem. Správné míry klasifikace lze porovnat, když měníme referenční prahovou hodnotu:

> classif.tab < – s(sck.expand, table (sick, prediction))

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

K zobrazení ROC křivky, použijeme balíček ROCR:

> library(ROCR)

> pred < – predikce(predict(glm.res2, type = „response“),

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”