6.6 zastosowania

wytworzenie tabeli danych może być wykonane w sposób wskazany poniżej, z tabeli utworzonej przy użyciu macierzy poleceń() :

> zakładka < – as.table (matrix (c (815,115,208,327), nrow = 2, byrow = TRUE,

dimnames = list(EST = C(„+”,”−”),

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

> tab

należy zauważyć, że tabela została zreorganizowana tak, aby odpowiadała zwykłym notacjom, z „dodatnimi” zdarzeniami w pierwszym wierszu (zazwyczaj ekspozycja) i pierwszej kolumnie (zazwyczaj choroba) tabeli.

możliwe byłoby wówczas obliczenie wszystkich wymaganych wielkości przy użyciu bardzo niewielu operacji arytmetycznych z R. Pakiet epiR zawiera jednak wszystkie narzędzia potrzebne do odpowiedzi na pytania epidemiologiczne w tabelach 2 × 2. Tak więc, Komenda epi.tests() podaje częstość występowania, czułość / specyficzność, a także pozytywne i negatywne wartości predykcyjne:

> biblioteka(epiR)

> epi.tests (tab)

Jeśli chodzi o obszar pod krzywą, można również użyć polecenia zewnętrznego, roc.od.table(), w pakiecie epicalc:

> biblioteka(epicalc)

> roc.od.tabela (tab, graph = FALSE)

Importowanie danych.dat jest osiągany w następujący sposób:

> bp < – Czytaj.stół („hdis.dat”, header = TRUE)

> str(bp)

Jak można zweryfikować, nie ma etykiety związanej z poziomami zmiennych zainteresowania (bpress dla ciśnienia krwi, chol dla poziomu cholesterolu). Aby wygenerować i powiązać etykiety, można wstawić następujące polecenia:

> 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, labels = clab)

ostatnie polecenie konwertuje pierwotne zmienne NA zmienne jakościowe i przyporządkowuje do nich etykiety zdefiniowane przez blab i clab. Aby sprawdzić, czy baza danych jest teraz w pożądanej formie, pomocne są następujące polecenia:

> str(bp)

> podsumowanie(bp)

tabela względnych częstotliwości ujawnionych w instrukcji może być teraz odtworzona:

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

ponieważ skupimy się wyłącznie na związku między ciśnieniem krwi a zawałem mięśnia sercowego, konieczne jest agregowanie danych dotyczących poziomu cholesterolu. Innymi słowy, konieczne jest gromadzenie wartości dla każdego poziomu ciśnienia krwi, niezależnie od poziomu cholesterolu. Należy również skorzystać z okazji, aby zmienić nazwy poziomów zmiennej bpress za pomocą środków przedziałów klasowych:

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

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

>dfrm < – agregat(bp,

list(bpress = BP), Sum)

jest to ostatnie polecenie, aggregate(), które umożliwia agregację danych: sumuje wszystkie wartości (sumy) zmiennej chol, które nie są zawarte w liście zmiennych, które mają być zachowane do analizy. Jednocześnie wyniki są przechowywane w nowej bazie danych o nazwie dfrm. Zauważ, że w tym przypadku nie używamy notacji formuły, ale alternatywnej składni dla aggregate () (zmienna odpowiedzi, po której następuje lista składająca się z jednej lub więcej zmiennych objaśniających). Przegląd zagregowanych danych można uzyskać, podając head () jako dane wejściowe:

> head(dfrm, 5)

gdy dostępna jest proporcja p, jej wartość na skali, której jednostki są logitami, jest podawana przez dziennik relacji(p/(1 − p)). Z tego można uzyskać następujące polecenia do konwersji proporcji, obliczonych jak HDI / total w jednostkach logit:

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

> dfrm$prop < – dfrm$HDIS/dfrm$Total

> dfrm$logit < – logit(dfrm$HDIS/dfrm$Total)

należy zauważyć, że mała funkcja została zdefiniowana do konwersji wartości x, które tutaj przyjmuje się proporcje, w równoważnym logu(x/(1 − x)). Równie dobrze można by napisać:

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

obliczone wartości zostały dodane bezpośrednio do ramki danych pod nazwami prop i logit.

model regresji logistycznej zapisywany jest w następujący sposób:

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

stosowany preparat, cbind(hdis, total-hdis) bpress, uwzględnia fakt, że uzyskujemy dostęp do danych zgrupowanych, a nie indywidualnych odpowiedzi. Polecenie glm () z opcją family = binomial odpowiada regresji logistycznej, która, nie wchodząc w szczegóły, domyślnie używa funkcji logit jako łącza kanonicznego.

współczynniki regresji mogą być wyświetlane za pomocą summary ():

> summary(glm(cbind (HDI, total-HDI) bpress, data = dfrm,

family = dwumian))

poprzedni wynik zawiera podstawowe informacje, aby odpowiedzieć na pytanie o znaczenie statystyczne związku między ciśnieniem krwi a prawdopodobieństwem zawału serca: współczynnik regresji (w skali log-odds) wynosi 0,024 i jest znaczący przy zwykłym progu 5% (patrz kolumna Pr (>|z|)).

prawdopodobieństwo wystąpienia zawału serca w zależności od różnych rozważanych poziomów ciśnienia krwi jest następujące:

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

family = dwumian)

> predict(glm.res)

należy zauważyć, że pośrednie wyniki generowane przez R zostały zapisane pod nazwą glm.res przed użyciem polecenia predict (). Przewidywania generowane przez R są wyrażone w postaci logitu i można porównać zaobserwowane i przewidywane logity:

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

praktycznie wszystkie elementy są dostępne, aby przedstawić proporcje obserwowanego i przewidywanego zawału mięśnia sercowego w sposób graficzny w zależności od poziomu ciśnienia krwi. Brakuje przewidywań modelu wyrażonych w postaci proporcji, a nie w skali logoddsa. Może zatem być pożądane narysowanie krzywej przewidywania, to znaczy prawdopodobieństwa wystąpienia zawału serca w zależności od ciśnienia krwi, bez ograniczania tego ostatniego do ośmiu wartości obserwowanych dla zmiennej bpress. Oto możliwe rozwiązanie:

> dfrm$prop.predit < – predict(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”, …)

})

ponieważ (pojedyncze) nieprzetworzone dane nie są dostępne, należy pracować bezpośrednio z tabelą częstotliwości podaną w instrukcji problemu:

> 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:

> biblioteka(VCD)

> oddsratio(alcohol, log = FALSE)

opcja log = FALSE zapewnia, że zwracany wynik odpowiada poprawnie ilorazowi szans, a nie logowi ilorazu szans. Aby uzyskać asymptotyczny przedział ufności, stosuje się następujące metody:

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

podobnie, podsumowanie(oddsratio(alkohol)) może być użyte do przeprowadzenia testu hipotezy na log odds ratio (H0 : log(θ) = 0).

właściwość polecenia.test () może być użyty do przetestowania hipotezy, że odsetek osób z dziennym spożyciem wynosi ≥ 80 g i jest identyczny zarówno w przypadkach, jak i w kontrolach, wskazując wartości obserwowane z tabeli krzyżowej podanej w oświadczeniu problemowym:

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

Ten test jest dokładnie równoważny testowi Z, aby sprawdzić różnicę między dwoma proporcjami oszacowanymi na podstawie danych (jeśli nie stosuje się korekcji ciągłości).

w odniesieniu do modelu regresji logistycznej tabela awaryjności musi zostać przekształcona w tabelę danych, w której dwie zmienne jakościowe (choroba i ekspozycja) są właściwie reprezentowane, to znaczy z ramką danych. Można to zrobić za pomocą polecenia melt () z pakietu reshape2, które umożliwia przekształcenie danych w formie tabelarycznej do formatu „long”. Na tym etapie poziomy zmiennej choroby należy przekodować za pomocą 0 (kontrola) i 1 (przypadki), chociaż nie jest to naprawdę konieczne, a poziom 0 należy uznać za kategorię odniesienia (co ułatwia interpretację wyników):

> library(reshape2)

> alcohol.df < – stop(alkohol)

> nazwy(alkohol.df) < – C(„choroba”, „ekspozycja”, „n”)

> poziomy(alkohol.df$) < – c(1,0)

> alkohol.df$choroba < – (alkohol.df$, „0”)

model regresji logistycznej zapisywany jest w następujący sposób:

> glm.res < – glm (disease exposure, data = alkohol.df,

family = dwumian, weights = n)

> podsumowanie(glm.res)

wynikiem zainteresowania jest wiersz związany z ekspozycją > = 80, ponieważ dostarcza informacji o wartości współczynnika regresji związanego ze zmienną ekspozycji oszacowaną przez R, z jej błędem standardowym, a także o wartości statystyki badania. Tutaj współczynnik regresji jest interpretowany jako log ilorazu szans. Należy zauważyć, że można uzyskać dokładnie takie same wyniki poprzez zamianę roli zmiennych w poprzednim sformułowaniu, choroba ekspozycji.

iloraz szans obliczony powyżej można znaleźć na podstawie współczynnika regresji związanego z czynnikiem zainteresowania (ekspozycją), a także jego 95% przedziału ufności, w następujący sposób:

> exp(coef(glm.res))

> exp(confint(glm.res))

drugie rozwiązanie polega na uwzględnieniu liczby przypadków i ogólnej liczby osób:

>< – dane.frame(expos = C(„<> = 80”), case = C(104,96),

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

> podsumowanie (glm (cbind (case, total-case) expos, data = alcohol2,

family = dwumian))

4) Dane te pochodzą z badania mającego na celu ustalenie prognostycznej ważności stężenia kinazy kreatynowej w organizmie w zapobieganiu wystąpieniu zawału mięśnia sercowego .

ponieważ nazwy zmiennych nie ma w pliku danych, konieczne będzie ich przypisanie natychmiast po zaimportowaniu danych:

> sck < – read.stół („sck.dat”, header = FALSE)

> nazwy(SCK) < – c(„CK”, „pres”, „abs”)

> podsumowanie(sck)

suma liczba badanych odpowiada sumie częstotliwości dla dwóch zmiennych pres i ABS, czyli:

> Sum(SCK)

lub, równoważnie: sum (sck$pres) + sum (SCK$abs) (ale nie sck$pres + SCK$abs).

aby obliczyć względne częstotliwości tych dwóch zmiennych, konieczne jest poznanie całkowitych wielkości według zmiennej. Można je uzyskać za pomocą polecenia apply i operować kolumnami:

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

na tej podstawie wystarczy podzielić każdą wartość zmiennych pres i abs przez wartości wcześniej obliczone. Uzyskane wartości zostaną zapisane w dwóch nowych zmiennych w tej samej tabeli danych:

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

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

można sprawdzić, czy obliczenia są poprawne: suma wartości dla każdej zmiennej musi być teraz równa 1:

> apply(sck, 2, sum)

wtedy uzyskane proporcje są po prostu reprezentowane na jednym wykresie, biorąc pod uwagę wartości wartości zmiennej.zmienna CK jako abscissae:

> xyplot(Pres.prop + abs.prop ck, data = sck, type = C („b”, „g”),

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

Instrukcja type = c(„b”, „g”) oznacza, że chcemy, aby punkty były wyświetlane połączone liniami („b” = „o” + „l”) oraz siatką („g”).

model regresji dla grupowanych danych zapisany jest jako:

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

> podsumowanie(glm.res)

przewidywania, wyrażone w formie prawdopodobieństwa, a nie w skali log-odds, są uzyskiwane za pomocą polecenia predict poprzez określenie opcji type = „response” w następujący sposób:

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

> nazwy(glm.pred) < – sck$CK

biorąc pod uwagę, że prawdopodobieństwo ≥ 0,5 oznacza osoby chore, otrzymuje się następujący rozkład:

> glm.pred

wnioskuje się, że ludzie będą uważani za chorych, według tego modelu, dla wartości ck 80 lub więcej. Można to łatwo zweryfikować za pomocą wykresu, na którym reprezentowane są prawdopodobieństwa przewidywane na podstawie wartości zmiennej ck:

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

> div > xyplot(GLM.pred ~ sck$CK, type = „l”,

ylab = „prawdopodobieństwo”, xlab = „ck”,

panel = function(…) {

panel.xyplot(…)

panel.xyplot (sck$CK, sck$sick, pch = 19, col = „grey”)

})

w pierwszym kroku konieczne jest „dekompresja” zgrupowanych danych i utworzenie tabeli z dwiema kolumnami: pierwsza reprezentująca zmienną ck i druga reprezentująca obecność lub brak choroby. Użyjemy zliczeń według podgrup wcześniej obliczonych i dostępnych w zmiennej ni:

> sck.rozwiń < – dane.frame (CK = c(rep(sck$CK, sck$pres),

rep(sck$CK, sck$abs)),

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

> table(SCK.expand$sick)

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

dwa ostatnie polecenia mają na celu zapewnienie, że te same liczby są napotkane i że rozkład chorych osób zgodnie z wartością ck jest prawidłowy. Można również zweryfikować, że te same wyniki są uzyskiwane w odniesieniu do regresji logistycznej i jednocześnie można dodać wartości przewidywane w poprzedniej tabeli danych (można zastosować tę samą procedurę, jak poprzednio, i można powielać prognozy, ale jest to prostsze w ten sposób):

> glm.res2 < – glm(chory ck, data = sck.expand, family=dwumian)

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

type = „response”) > = 0.5,

1, 0)

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

ostatnie polecenie umożliwia wyświetlenie macierzy pomieszania, w której krzyżują się rzeczywiste diagnostyki i przewidywane przez model regresji. Prawidłowe wskaźniki klasyfikacji można porównać, gdy zmieniamy próg odniesienia:

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

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

aby wyświetlić krzywą ROC, użyjemy pakietu ROCR:

> library(ROCR)

>prediction(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”