W »In mir jmytraf Jury statystyczne
1'iml/N/y !«*<! przedstawia jak przeprowadzić taki test metodą bootstrap.
b 0 definiujemy funkcję do wyznaczania statystyki T dla replikacji
bootstapowych, drugim argumentem tej funkcji są indeksy wylosowanych elementów w danej replikacji
> NtAtystykaT <- function(x,i) (mean(x[i]> - raean(x)}/sd(x[i])
> 0 r. a implementujmy funkcje wyznaczającą p-wartość dla dwustronnego test u
bootstrapowego opartego na statystyce T b boostrapowyTestT <- function(wektor, mu0*0) {.
♦ 0 wyznaczamy 999 replikacji z próby i liczymy dla nich wartości statystyki testowej
♦ wartosciStatystykiT = wynikBoot <- boot(wektor, statystyka!*, R=999, stype="i")
♦ 0 wyznaczamy dystrybuantę dla rozkładu statystyki T na replikacjąch .
♦ dystrybuantaStatystykiT ■ ecdf(wartosciStatystykiT$t)
♦ T ■ mean(wektor - muO)/sd(wektor)
♦ 0 wyznaczamy kwantyl odpowiadający wartości statystyki T w rozkładzie statystyki testowej dla replikacji bootstrpowych
♦ kwantyl » dystrybuantaStatystykiT(T)
+ 0 wyznaczamy p-wartość. dla testu dwustronnego
$ l-2*abs(kwantyl-0.5)
+ }
> 0 testujemy czy Średnia VEGF jest istotnie różna od 3000 (nie jest), ,
wynikiem funkcji boostrapowyTestT O jest ocena p-wartości
> boostrapowyTestT(daneO$VEGF, 3000)
[1] 0.4204204
> 0 testujemy czy średnia VEGF jest istotnie różna od Ą000 (jest) 3
> boostrapowyTestT(daneO$VEGF, 4000)
[1] 0.02002002
Mtt' Im I • I*m« N»me ntuidnr iliitwil t mmii w 4*lnM i upy Iliom AimI miinlnrlng llfelM llintiuli not as «t||nli», *.Im miel ftluo
Im* i*v<i|ile<l.
hllllti K llarrcll (iiiiHWTring w i|«nwUon on • Ht^Morlantion of ••lit Iiiiuhim vni ittblos Itl Mirvival modolling) fnrtune(132)
W tyra rozdziale przedstawione będą skrótowo funkcje pozwalające na wykonanie podstawowej analizy przeżycia z użyciem pakietu R. Funkc je do przeprowadzania takiej analizy dostępne są w pakiecie survival. Zostaną ono przedstawione na przykładzie zbioru danych o pacjeutkach oddziału onkologii (zbiór danych daneO, patrz załącznik).
Zmienną, którą będziemy tutaj analizować, jest przedział czasu od operacji do wznowienia się choroby. Ponieważ w przypadku łagodnych postaci nowotworów piersi; u znacznego procentu kobiet nie dochodzi do wznowienia się choroby, dlatego też pomiary podzielić będzie można na dwa zbiory. Pierwszy, w którym będą pacjentki, u których po zakończeniu leczenia doszło do wznowienia się choroby (w skrócie:; doszło do wznowy) i drugi, w którym będą pacjentki, które są pod obserwacją przez określony okres czasu i nie doszło u nich do wznowienia się choroby. Dla drugiej grupy pacjentek nie wiemy kiedy dojdzie do wznowy i czy w ogóle do niej dojdzie; Mamy jedynie informacje, że w określonym okresie czasu nic złego się nie działo. O takich obserwacjach będziemy mówić, że są cenzurowane, a więc mamy dla nich jedynie część informacji.
Zanim przystąpimy do jakichkolwiek czynności związanych z analizą przeżycia musimy stworzyć obiekty klasy Surv, który będzie przechowywał informacje o czasach przeżycia (czasach do wznowienia się choroby) i o cenzurowaniu. Do tego celu służy funkcja Surv(survival), która jako pierwszy argument przyjmuje wektor czasów. a jako drugi argument przyjmuje wektor wartości logicznych lub identyfikatorów {0,1} określający, które pomiary są cenzurowane.
> tt wczytujemy zbiór danych
> daneO - read.table("http://wuw.bieeek.pl/R/dane/dane0.csv", sep-";",
header»T)
> # funkcją Sury określamy, które obserwacje są cenzurowane a które nie
> czasy = withCdaneO, Surv(0kres.bez.wznowy, Niepowodzenia=="wznowa"))
> czasy[1:10]
. [1] 22+ 53+ 38+ 26+ 19+ 36 33+ 38+ 38 37+
Dla obiektu klasy Surv dostępna jest przeciążona funkcja printO. Odpowiada ona za to, że wyświetlając obiekty tej klasy obserwacje cenzurowane zaznaczane są plusem. Obiekt klasy Surv (w tym przypadku jest to obiekt czasy) możemy wykorzystać w dalszych analizach.
Estymator Kapłana Meyera oceniający prawdopodobieństwo przeżycia można wyznaczyć funkcją survfit(survival). Prawdopodobieństwo przeżycia czasu I określone jest następującym wzorem
S(t) = Pr(T > t),
gdzie T to rozkład czasu życia. Innymi słowy funkcja przeżycia określa jakie jest prawdopodobieństwo życia dłużej niż do czasu 1.
Pierwszym argumentem funkcji suryfitO jest formula, w której zmienną objaśnianą powinna być zmienna klasy Surv będąca wynikiem funkcji Surv(survival) (więcej o składni formuł przeczytać można w podrozdziale 2.1.8). Prawa strona formuły jest nieobowiązkowa może być pominięta lub może wskazywać na zmienną : jakościową. Jeżeli zostanie wskazana zmienna jakościowa, to krzywa przeżycia będzie wyznaczona dla każdego poziomu danej zmiennej. Wynikiem działania funkcji suryfitO jest klasy survf it i dla takich obiektów dostępne są przeciążone funkcjo plot() (rysująca krzywą przeżycia) oraz summaryO (wypisującą wartości krzywej przeżycia we wskazanych punktach określonych argumentem times. Jeżeli ten argument nie zostanie wskazany wypisywane będą wartości krzywej przeżycia w punktach odpowiadającym wszystkim niecenzurowanym obserwacjom).
Funkcja suryfitO domyślnie wyznacza estymator Kaplana-Meyera, ale dostępny jest w niej również estymator Fleminga-Harringtona (wybrać go można zmieniając. wartość argumentu type="f leming-harrington"). Określając argument weights możemy również wskazać wektor wag dla poszczególnych olrserwacji. Poniżej przedstawiamy przykład użycia funkcji surfyitO oraz prezentacji jej wyników.
> # wyznaczamy estymator Kapłana Meyera
> krzywaKM - suryfit(czasy)