Regresní modelování v analýze přežití Tomáš Pavlík pavlik@iba.muni.cz Bi7491 Regresní modelování 9. května 2012 Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 1 / 121 Motivace Motivace V klasické regresi nás zajímá vliv vysvětlujících proměnných (sledovaných faktorů) na hodnoty sledované náhodné veličiny. V logistické regresi nás zajímá vliv vysvětlujících proměnných na výskyt nebo nevýskyt nějakého jevu. Tyto metody jsou však nevhodné ve chvíli, kdy 1 studujeme čas do výskytu nějakého jevu 2 může se stát, že u všech subjektů nemáme o výskytu tohoto jevu kompletní informace. Analýza přežití (survival analysis) – soubor metod pro popis a modelování doby do výskytu sledovaného jevu v čase. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 2 / 121 Motivace Motivace - příklad Jméno pacienta Zemřel Délka sledování Pan Silný Ne 12,0 měsíců Pan Slabý Ano 6,8 měsíců Pan Moucha Ano 4,8 měsíců Pan Komár Ano 9,8 měsíců Pan Skála Ne 10,8 měsíců Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 3 / 121 Motivace Příklady Čas od diagnózy do úmrtí onkologického pacienta Čas od zahájení léčby do progrese onemocnění Čas od propuštění z nemocnice do opakované hospitalizace Čas od infekce HIV do propuknutí AIDS Čas od narození do prvního požití alkoholu/drog Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 4 / 121 Motivace Modelová data - projekt Alert Projekt Alert - Klinický registr pro standardizovaný sběr diagnostických, prognostických a léčebných dat u pacientů s akutní leukémií. charakteristika pacientů s AML, ALL a APL v ČR a SR sběr a vyhodnocování základních epidemiologických a klinických dat sběr a vyhodnocování molekulárně-biologických a genetických dat modelování přežití Definice souboru pacientů: Dospělí pacienti s de novo AML i sekundární AML, diagnostikovaní v ČR (5 hematologických center) v období 1999–2009 Kurativně léčení (léčba se snahou o odstranění nemoci) N = 1091 Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 5 / 121 Literatura Literatura Marubini & Valsecchi: Analysing Survival Data from Clinical Trials and Observational Studies Collett: Modelling Survival Data in Medical Research Hosmer & Lemeshow: Applied Survival Analysis Klein & Moeschberger: Survival Analysis: Techniques for Censored and Truncated Data Therneau & Grambsch: Modeling Survival Data: Extending the Cox Model Andersen, Borgan, Gill & Keiding: Statistical Models Based on Counting Processes Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 6 / 121 Úvod a definice pojmů Úvod a definice pojmů Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 7 / 121 Úvod a definice pojmů Definice pojmů Čas do výskytu sledované události (time to event), jinak také čas přežití (survival time, failure time, event time), je nezáporná náhodná veličina. Budeme ji značit T. Lze ji jednoznačně popsat následujícími charakteristikami: Distribuční funkcí: F(t) = P(T ≤ t), t ≥ 0 Hustotou pravděpodobnosti: f (t) = F′ (t), t ≥ 0 Funkcí přežití: S(t) = P(T > t) = 1 − F(t), t ≥ 0 Náhodná veličina T může být jak spojitá, definovaná na intervalu (0, ∞), tak diskrétní, nabývající nejvýše spočetně mnoha hodnot, např. a1, a2, . . . , an. My se zde budeme zabývat pouze spojitým případem. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 8 / 121 Úvod a definice pojmů Klíčové prvky v analýze přežití Pro definici času přežití jako náhodné veličiny potřebujeme stanovit následující: Jednoznačný počátek sledování (např. datum diagnózy, narození, zahájení léčby) Časovou škálu (např. reálný čas v měsících nebo letech) Koncovou událost (např. úmrtí, progrese onemocnění, rehospitalizace, otěhotnění) Sledovaná událost by měla být snadno pozorovatelná či měřitelná (úmrtí × progrese onemocnění). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 9 / 121 Úvod a definice pojmů Úlohy v analýze přežití 1 Popis – bodové a intervalové odhady pravděpodobnosti přežití v čase t 2 Srovnání – hodnocení rozdílů v přežívání dvou a více skupin (pacientů, osob, zvířat, věcí, bakterií, apod.) 3 Modelování – hodnocení vlivu vysvětlujících proměnných na pozorovaný čas přežití Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 10 / 121 Úvod a definice pojmů Specifikum dat přežití - Cenzorování Definovaná událost se nemusí v průběhu sledování vyskytnout u všech subjektů (pacientů) → Nevíme tedy, kdy a jestli vůbec daná událost nastala. Víme pouze, že nenastala před ukončením sledování. Cenzorování je ztrátou určité informace, cenzorované časy přežití však nelze z analýzy vyřadit. Typy cenzorování: 1 Cenzorování zprava 2 Cenzorování zleva 3 Intervalové cenzorování Bez ohledu na typ cenzorování, hlavním předpokladem analýzy přežití je, že cenzorování je tzv. neinformativní (non-informative). To znamená, že výskyt cenzorování nijak nesouvisí s výskytem sledovaných událostí. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 11 / 121 Úvod a definice pojmů Typy cenzorování 1 Cenzorování zprava – pozorujeme minimum z hodnot skutečného a cenzorovaného času přežití. Klinické studie Observační studie 2 Cenzorování zleva – pozorujeme maximum z hodnot skutečného a cenzorovaného času přežití. Věk prvního požití drog, kouření 3 Intervalové cenzorování – víme pouze, že skutečný čas přežití leží v intervalu (D, H). Sledování pacientů v dlouhodobějších intervalech (screening karcinomu prostaty, pacienti s CML) Nejčastějším typem cenzorování je v biolgii a medicíně cenzorování zprava. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 12 / 121 Úvod a definice pojmů Cenzorování zprava - příklad Délka sledování a výskyt události X X X X zahájení sledování ukončení sledování = cenzorování X = výskyt události Obrázek 1. Sledování času přežití v přítomnosti cenzorování. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 13 / 121 Úvod a definice pojmů Cenzorování - vliv procenta cenzorovaných hodnot na odhad funkce přežití 0 20 40 60 80 0.00.20.40.60.81.0 Mesíce Podílzijícíchpacientu 0 20 40 60 80 0.00.20.40.60.81.0 Mesíce Podílzijícíchpacientu Obrázek 2. Křivka přežití bez cenzorovaných hodnot a s 50% cenzorovaných hodnot. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 14 / 121 Úvod a definice pojmů Délka sledování Design (plán) studie má vždy vliv na její analýzu. Informace o výskytu sledované události v čase je primárně zachycena v kompletních (skutečných) časech přežití. Délka sledování (studie) proto musí být dostatečná, abychom pozorovali dostatek událostí a tedy i informace o sledovaném procesu. Souvisí s četností výskytu události (např. mortalitou u rakoviny). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 15 / 121 Úvod a definice pojmů Čas přežití v přítomnosti cenzorování Jak jsou data o přežití reprezentována v přítomnosti cenzorování? Ti je zaznamenaný čas přežití pro pacienta i Označme 𝛿i indikátor pozorované události. Pak 𝛿i = {︃ 1 když pozorujeme skutečný čas do události 0 když je čas Ti cenzorován Data přežití n subjektů pak budeme reprezentovat dvojicemi (ti , 𝛿i ), i = 1, . . . , n. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 16 / 121 Úvod a definice pojmů Hustota pravděpodobnosti (density function) Hustota udává pravděpodobnost výskytu sledované události v čase t. T jako spojitá náhodná veličina: f (t) = lim u→0 1 u P(t ≤ T ≤ t + u) T jako diskrétní náhodná veličina: f (t) = P(T = t) = {︃ fj když t = aj , j = 1, . . . , n 0 když t ̸= aj , j = 1, . . . , n Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 17 / 121 Úvod a definice pojmů Funkce přežití (survival function) Hlavní charakteristikou je tzv. funkce přežití. Ta udává pravděpodobnost, že jedinec přežije (vzhledem k výskytu sledované události) déle než do času t. T jako spojitá náhodná veličina: S(t) = P(T > t) = ∫︁ ∞ t f (x)dx = 1 − F(t) T jako diskrétní náhodná veličina: S(t) = ∑︁ aj ≥t f (aj ) = ∑︁ aj ≥t fj Nabývá hodnot mezi 1 a 0 (respektive 100% a 0%), kdy hodnotu 1 má v počátečním čase a hodnotou 0 při výskytu poslední události Je to nerostoucí funkce Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 18 / 121 Úvod a definice pojmů Funkce přežití - příklad 0 20 40 60 80 100 120 0.00.20.40.60.81.0 Mesíce Podílzijícíchpacientu Obrázek 3. Odhad funkce přežití pacientů s AML léčených kurativně (N = 1091). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 19 / 121 Úvod a definice pojmů Riziková funkce (hazard function) Druhou důležitou charakteristikou je tzv. riziková funkce. Ta vyjadřuje intenzitu výskytu sledované události v čase t za podmínky, že subjekt přežil do času t. Riziková funkce je klíčová pro modelování přežití - řada modelů je definována a interpretována pomocí rizikové funkce. T jako spojitá náhodná veličina: h(t) = lim u→0 P{t < T ≤ t + u|T > t} u = lim u→0 P{t < T ≤ t + u}/P{T > t} u = lim u→0 [F(t + u) − F(t)]/u S(t) = dF(t)/dt S(t) = f (t) S(t) = − d ln S(t) dt Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 20 / 121 Úvod a definice pojmů Riziková funkce (hazard function) T jako diskrétní náhodná veličina: h(aj ) = hj = P(T = aj | T ≥ aj ) = P(T = aj ) P(T ≥ aj ) = f (aj ) S(aj ) = f (aj ) ∑︀ k:ak ≥aj f (ak ) Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 21 / 121 Úvod a definice pojmů Riziková funkce - příklad 0 20 40 60 80 100 120 0.00.20.40.60.81.0 Mesíce Riziko Obrázek 4. Modelové příklady rizikové funkce, respektive vývoje rizika výskytu sledované události v čase. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 22 / 121 Úvod a definice pojmů Kumulativní riziková funkce (cumulative hazard function) Kumulativní riziková funkce odpovídá celkovému riziku výskytu sledované události od začátku sledování až do času t. T jako spojitá náhodná veličina: H(t) = ∫︁ t 0 h(x)dx T jako diskrétní náhodná veličina: H(t) = ∑︁ j:aj 1, riziková funkce je monotónně rostoucí ve všech podskupinách definovaných vysvětlujícími proměnnými. Vliv proměnných je multiplikativní. S nárůstem xk o 1 se zvýší riziko exp(𝛽k )-krát. Pro dva subjekty s vektory x1 a x2 platí HR = h(t, x2) h(t, x1) = exp ((x2 − x1)′ 𝛽) Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 71 / 121 Coxův model proporcionálních rizik Coxův model proporcionálních rizik Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 72 / 121 Coxův model proporcionálních rizik Coxův model proporcionálních rizik Coxův model proporcionálních rizik je dán vztahem h(t, xi ) = h0(t) exp(xi1 𝛽1 + xi2 𝛽2 + . . . + xip 𝛽p), Předpoklady vysvětlující proměnné mají aditivní vliv na škále log h(t, xi ). vliv vysvětlujících proměnných je stejný pro ∀t. Na rozdíl od parametrických modelů nespecifikujeme h0(t). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 73 / 121 Coxův model proporcionálních rizik Poměr rizik (Hazard ratio, HR) Poměr rizik pro subjekty s vektory vysvětlujících proměnných x1 a x2 lze vyjádřit takto: HR = h(t, xi ) h(t, xj ) = h0(t) exp(x′ i 𝛽) h0(t) exp(x′ j 𝛽) = exp((xi − xj )′ 𝛽), Bodový odhad pro poměr rizik je ^HR = h(t, xi ) h(t, xj ) = h0(t) exp(x′ i ^𝛽) h0(t) exp(x′ j ^𝛽) = exp((xi − xj )′ ^𝛽), kde ^𝛽 je maximálně věrohodný odhad 𝛽. 100(1 − 𝛼)% interval spolehlivosti lze získat takto exp((x1 − x2)′ ^𝛽 ± z1−𝛼/2 ^SE((x1 − x2)′ ^𝛽)). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 74 / 121 Coxův model proporcionálních rizik Základní riziková funkce Základní riziková funkce, h0(t), vyjadřuje riziko pro subjekty, pro které platí x1 = 0, . . . , xp = 0. Jedná se o riziko odpovídající referenční skupině subjektů: h(t, xi = 0) = h0(t) exp(0 * 𝛽1 + 0 * 𝛽2 + . . . + 0 * 𝛽p) = h0(t). Pro odhad vlivu vysvětlujících proměnných, HR, základní rizikovou funkci nepotřebujeme. Potřebujeme ji ale pro odhad funkce přežití a rizikové funkce pro každé x. Můžeme použít odhad kumulativní rizikové funkce, H0(t), dle Breslowa: ^H0(t) = ∑︁ ti ≤t ^h0(ti ) = ∑︁ ti ≤t di ∑︀ j∈Ri exp(x′ j ^𝛽) , kde di je počet úmrtí v čase ti a Ri je příslušný počet pacientů v riziku sledované události. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 75 / 121 Coxův model proporcionálních rizik Odhad regresních koeficientů Pro odhad regresních koeficientů navrhl Cox metodu parciální věrohodnosti (partial likelihood), kdy je místo věrohodnostní funkce maximalizována tzv. parciální věrohodnostní funkce. Ta má pro vektor regresních koeficientů 𝛽 tvar: L(𝛽) = k∏︁ i=1 exp(x′ i 𝛽) ∑︀ j∈Ri exp(x′ j 𝛽) . Logaritmus parciální věrohodnostní funkce pak má tvar: log L(𝛽) = k∑︁ i=1 ⎧ ⎨ ⎩ x′ i 𝛽 − log ⎡ ⎣ ∑︁ j∈Ri exp(x′ j 𝛽) ⎤ ⎦ ⎫ ⎬ ⎭ = k∑︁ i=1 li . Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 76 / 121 Sestavení modelu Sestavení modelu Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 77 / 121 Sestavení modelu Cíl studie ovlivňuje výběr vysvětlujících proměnných Před modelováním je třeba si ujasnit, proč vlastně modelujeme a co od modelu požadujeme. Modely umí zároveň zohlednit více proměnných – jejich výběr ale vždy leží na analytikovi! Scénáře mohou být následující: 1 Cílem hodnocení je jediná proměnná 2 Cílem je prediktivní model 3 Cílem je identifikace potenciálních prediktorů Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 78 / 121 Sestavení modelu Cílem hodnocení je jediná proměnná Jde nám o kvantifikaci a statistickou významnost vlivu jedné vysvětlující proměnné, např. podané léčby (předpokládáme totiž její vliv na přežití). Typické je modelování výsledků klinických studií. Ostatní vysvětlující proměnné vystupují v modelu v roli adjustačních faktorů - chceme odstranit jejich vliv. V modelu by měly kvůli adjustaci zůstat i proměnné s nevýznamným vlivem - i ty mohou hrát roli v identifikaci vlivu sledované proměnné. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 79 / 121 Sestavení modelu Význam regresního modelování Regresní model umožňuje současně uvažovat vliv více proměnných a vzájemně tak adjustovat jejich vlivy: Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 80 / 121 Sestavení modelu Cílem je prediktivní model Z množiny zaznamenávaných proměnných chceme vybrat sadu s významným vlivem na přežití a schopností jeho predikce. U všech předpokládáme relevanci vzhledem k přežití, ale ve výsledném modelu se vyskytují pouze ty “nejsilnější”. Klíčové téma je složitost modelu (počet vysvětlujících proměnných a jejich interakcí). Statistická významnost proměnné nemusí zaručovat její přínos pro predikci. Je třeba rozlišovat regresní model a prediktivní model. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 81 / 121 Sestavení modelu Cílem je identifikace potenciálních prediktorů Částečně exploratorní práce – z množiny zaznamenávaných proměnných chceme vybrat ty významně asociované s přežitím, např. expresní profily zkoumaných genů (mohou jich být až tisíce). Jde nám o redukci počtu vysvětlujících proměnných a identifikaci těch nejvýznamnějších. Je třeba dávat pozor na falešně pozitivní výsledky. Problém n ≪ p. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 82 / 121 Sestavení modelu Výběr vysvětlujících proměnných I “The data analyst knows more than computer” Klíčové téma, roli zde hraje i úplnost dat - které proměnné si vůbec můžeme dovolit do modelování zahrnout. Je-li vysvětlujících proměnných mnoho, je třeba zapojit vícerozměrné metody identifikovat shluky proměnných, které jsou korelované. Následně vybrat pouze reprezentativní zástupce. Cílem použití vícerozměrných metod je odstranit redundantní informaci. Tento postup nelze automatizovat, závisí na znalosti konkrétního problému. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 83 / 121 Sestavení modelu Výběr vysvětlujících proměnných II “The data analyst knows more than computer” Nejoblíbenější jsou tzv. stepwise procedury: Backward elimination - vysvětlující proměnné z modelu postupně ubíráme. Forward selection - vysvětlující proměnné do modelu postupně přidáváme. Nelze je používat slepě, analytik musí vždy nad modelováním přemýšlet a pracovat s literaturou/odborníkem. Zvlášť obtížné je modelování interakcí. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 84 / 121 Sestavení modelu Rozvaha nad velikostí vzorku Velikost vzorku je zvlášť v analýze přežití extrémně důležitá ⇒ jde nám hlavně o dostatečný počet sledovaných událostí (problém cenzorování). Právě čas do sledované události představuje v analýze přežití informaci! Velikost vzorku by měla být vždy plánována před zahájením experimentu (pomocí software). Peduzzi a kol. (1995) navrhli na základě simulací alespoň 10 událostí na 1 vysvětlující proměnnou zahrnutou do modelu. Co když nemáme alespoň 10 událostí na 1 vysvětlující proměnnou? Bodové i intervalové odhady regresních koeficientů budou nevěrohodné. Model může selhat ve výběru významných proměnných. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 85 / 121 Sestavení modelu Kategorizace spojitých proměnných Kategorizace spojitých proměnných by měla být prováděna co nejméně, protože se jedná o ztrátu informace, která může zvýšit variabilitu a vést ke zkreslení výsledků. Někdy je však důležitá z hlediska interpretace, např. věkové kategorie (0 – 50 let, 51 – 60 let, 61 – 70 let, nad 70 let). Kategorizace může být použita pro řešení nelinearity vlivu vysvětlující proměnné. Doporučení pro kategorizaci spojitých proměnných: 1 Kategorie by měly mít logiku, případně odrážet biologickou podstatu. 2 Je třeba dávat pozor na počty subjektů v kategoriích, nejlépe je mít kategorie vyvážené. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 86 / 121 Regresní diagnostika Regresní diagnostika Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 87 / 121 Regresní diagnostika Co je regresní diagnostika? Nástroje regresní diagnostiky slouží pro hodnocení vhodnosti modelu (goodness of fit) a splnění předpokladů modelu (model assumptions). Cílem je zjistit, zda data nejsou v rozporu s předpoklady modelu a zda náš model adekvátně vystihuje přežití sledovaného souboru. Hodnocení vhodnosti modelu je těsně spjato s výběrem modelu a výběrem vysvětlujících proměnných. Nástroje regresní diagnostiky: 1 Vizualizace neparametrických odhadů charakteristik přežití 2 Výpočet reziduí modelu a jejich vizualizace (grafy reziduí) 3 Akaikeho informační kritérium (AIC) 4 Statistické testy Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 88 / 121 Regresní diagnostika Hodnocení vhodnosti modelu “Každý model je špatný, ale některý může být i užitečný.” Velmi náročný úkol, protože objektivně není dáno, co je vhodný model a co už ne. Pro srovnání dvou modelů lze použít tzv. Akaikeho informační kritérium. Pro hodnocení celkového fitu lze použít test dle Parzena & Lipsitze (1999) založený na podobném principu jako Pearsonův chí-kvadrát test pro kontingenční tabulky. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 89 / 121 Regresní diagnostika Akaikeho informační kritérium (AIC) AIC (Akaike, 1974) je statistika zahrnující věrohodnost modelu a jeho složitost. Slouží k posouzení schoposti různých modelů fitovat pozorovaná data. AIC = −2ℓ(𝛽, (t1, 𝛿1), . . . , (tn, 𝛿n)) + 2(c + a) kde ℓ(𝛽, (t1, 𝛿1), . . . , (tn, 𝛿n)) je logaritmus věrohodnostní funkce modelu, c je počet vysvětlujících proměnných v modelu a a je počet parametrů uvažovaného rozdělení pravděpodobnosti. Nižší hodnoty AIC indikují lepší model. Nelze srovnávat AIC parametrických modelů a Coxova modelu. Proč? Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 90 / 121 Regresní diagnostika Rezidua modelu Obecně lze rezidua definovat jako rozdíl mezi pozorovanou a predikovanou hodnotou sledované veličiny. Velké hodnoty reziduí, případně jejich “systematické chování”, jsou indikátorem špatného modelu. Kvůli cenzorování není jednoduché hodnoty reziduí interpretovat – přínosem je grafická vizualizace a vyhlazení trendu (např. jádrovým vyhlazováním). Obecně lze říci, že by grafy reziduí neměly vykazovat žádný trend (rezidua by měla tovřit rovnoměrný horizontální pás). Vybrané typy reziduí v analýze přežití: Martingale rezidua Deviance rezidua Skórová rezidua Schoenfeldova rezidua Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 91 / 121 Regresní diagnostika Martingale rezidua Martingale reziduum představuje rozdíl mezi pozorovaným a předpokládaným počtem událostí (dle daného modelu) u subjektu i. Předpokládaný počet událostí (^Ei ) je reprezentován kumulativním rizikem do času ti . ^rM i = 𝛿i − ^Ei = 𝛿i − ^H0(^𝛽, ti ) exp(x′ i ^𝛽) Martingale reziduum je jedno číslo charakterizující shodu pozorování s předpokládaným rizikem. Reziduum deviance je také jedno číslo. Skórové reziduum i Schoenfeldovo reziduum je vektor hodnot - každý subjekt má jednu hodnotu pro každou z vysvětlujících proměnných. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 92 / 121 Regresní diagnostika Význam jednotlivých reziduí Rezidua Použití v regresní diagnostice Martingale Identifikace nelineárního vlivu proměnné Martingale Identifikace nesprávně vyloučené proměnné Deviance Identifikace odlehlých pozorování/jedinců Skórová Identifikace vlivných pozorování/jedinců Schoenfeldova Testování předpokladu proporcionality rizik Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 93 / 121 Regresní diagnostika Martingale rezidua - příklad 1 Martingaleresidual Age at diagnosis 30 40 50 60 70 80 −5 −4 −3 −2 −1 0 1 Zdroj: Bradburn et al. (2003) Survival Analysis Part III: Multivariate data analysis – choosing a model and assessing its adequacy and fit. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 94 / 121 Regresní diagnostika Martingale rezidua - příklad 2 Zdroj: Therneau et al. (1990) Martingale-Based Residuals for Survival Models. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 95 / 121 Regresní diagnostika Skórová rezidua - příklad 0 50 100 150 200 -0.06-0.020.020.06 Patient ID Sokalscore(Highrisk) 0 50 100 150 200 -0.6-0.4-0.20.00.2 Patient ID Imatinibdosage x Zdroj: Analýza dat registru CAMELIA. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 96 / 121 Regresní diagnostika Ověření předpokladu proporcionality rizik 1 Grafické ověření - nejjednodušší forma ověření, K-M křivky by se měly “rovnoměrně vzdalovat”, v žádném případě se nesmí křížit. Ještě vhodnější je použití logaritmu kumulativní rizikové funkce. 2 Test pomocí časově závislé proměnné 3 Test založený na škálovaných Schoenfeldových reziduích Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 97 / 121 Regresní diagnostika Grafické ověření proporcionality rizik Máme k dispozici jednoduchou pomůcku pro ověření jednorozměrné proporcionality, neoť funkce přežití vzhledem k proměnné k musí v případě Coxova modelu splňovat Sk (t) = exp(−H0(t) exp(xk 𝛽k )), k = 1, . . . , p, a tedy platí log[− log(Sk (t))] = log H0(t) + xk 𝛽k , Pokud je předpoklad proporcionality splněn, křivky log[− log(Sk (t))], pro jednotlivé hodnoty proměnné k budou přibližně paralelní. Grafické ověření však nezohledňuje více faktorů. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 98 / 121 Accelerated Failure Time model Accelerated Failure Time model Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 99 / 121 Accelerated Failure Time model Accelerated Failure Time (AFT) model 1 Čas přežití itého subjektu, Ti , je nezáporný → můžeme modelovat jeho logaritmus. 2 AFT model je pak definován jako: log Ti = x′ i 𝛽 + 𝜖i nebo log Ti = x′ i 𝛽 + 𝜇 + 𝜎𝜖i , kde 𝜖i je reziduální člen s daným rozdělením pravděpodobnosti. Ti = exp(x′ i 𝛽) exp(𝜖i ) = exp(x′ i 𝛽)T0i Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 100 / 121 Accelerated Failure Time model AFT model - alternativní zápis Jinou formou zápisu AFT modelu je zápis pomocí funkce přežití: Si (t) = S0(t/ exp(x′ i 𝛽)), kde S0(t) představuje základní funkci přežití odpovídající referenční skupině s vektorem vysvětlujících proměnných xi = 0. Regresní AFT model je vhodnou alternativou pro model proporcionálního rizika tehdy, když je předpoklad proporcionality rizik porušen. Odhad regresních koeficientů je opět založen na metodě maximální věrohodnosti. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 101 / 121 Relativní přežití Relativní přežití Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 102 / 121 Relativní přežití Vyjádření pravděpodobnosti přežití pacientů Nejčastěji jako tzv. celkové přežití (overall survival), někdy označováno také jako pozorované přežití (observed survival). Celkové přežití odráží celkovou mortalitu pacientů bez ohledu na přesnou příčinu úmrtí. Chceme-li kvantifikovat mortalitu spojenou pouze se sledovanou diagnózou, musíme zohlednit pouze vybrané příčiny úmrtí. Jak? 1 Výpočet přežití specifického dle diagnózy (cause-specific survival) 2 Výpočet relativního přežití (relative survival) Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 103 / 121 Relativní přežití Přežití specifické pro sledovanou diagnózu Výpočet je jednoduchý, jediný rozdíl proti výpočtu celkového přežití je v tom, že časy pacientů, kteří zemřeli z jiné příčiny než z příčiny základního onemocnění, jsou cenzorovány. Problém s výpočtem přežití specifického pro sledovanou diagnózu však nastává při hodnocení populačních dat, která svojí kvalitou záznamů nemohou monitorovaným klinickým studiím konkurovat. Problém s kódováním přesné příčiny úmrtí u onkologických pacientů nemusí být v administrativě záznamu, ale spíše v jeho nejednoznačnosti, protože ne vždy je klinicky zřejmé, jestli pacient zemřel v souvislosti se sledovaným onemocněním nebo ne. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 104 / 121 Relativní přežití Relativní přežití Relativní přežití představuje poměr celkového a tzv. očekávaného přežití (expected survival). Očekávané přežití odráží mortalitu v obecné populaci, která odpovídá sledované skupině pacientů věkem, pohlavím a obdobím diagnózy. Relativní přežití ≈ vážený ekvivalent celkového přežití, váhou je přežití obecné populace. Relativní přežití ≈ celkové přežití korigované na mortalitu spojenou s dalšími chorobami, na něž může pacient zemřít. Relativní přežití ≈ odhad pravděpodobnosti přežití, který odpovídá pouze zátěži představované sledovanou diagnózou. Relativní přežití představuje míru tzv. excess mortality, která je asociována s danou chorobou bez ohledu na to, zda přímo či nepřímo. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 105 / 121 Relativní přežití Relativní přežití Označme h* (t) a S* (t) očekávanou rizikovou funkci a očekávanou funkci přežití obecné populace. Obdobně označme h(t) a S(t) pozorovanou rizikovou funkci a funkci celkového přežití. Relativní rizikovou funkci, označme ji hR (t), pak vypočteme jako hR (t) = h(t) − h* (t), a relativní funkci přežití, SR (t), jako SR (t) = S(t) S*(t) . Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 106 / 121 Relativní přežití Relativní přežití - příklad Zdroj: Data NOR. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 107 / 121 Relativní přežití Metody odhadu očekávaného přežití Principem výpočtu očekávaného přežití je výpočet přežití z populačních mortalitních tabulek odpovídajícího srovnatelné skupině (vzhledem k věku, pohlaví a roku diagnózy) z obecné populace, u které předpokládáme, že prakticky není zasažena sledovaným onemocněním. Metody výpočtu očekávaného přežití: Edererova metoda I Edererova metoda II Hakulinenova metoda Všechny metody vycházejí ze stejných datových podkladů, tedy úmrtnostních tabulek pro danou populaci (stát). Mezinárodní zdroj: www.mortality.org. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 108 / 121 Relativní přežití Rozdíl v metodách odhadu očekávaného přežití Edererova metoda I: srovnatelná skupina z obecné populace je uvažována v riziku bez omezení, nebereme tedy ohled na mortalitu a cenzorování v hodnocené kohortě. Edererova metoda II: srovnatelná skupina z obecné populace je uvažována v riziku pouze do úmrtí nebo cenzorování odpovídajícího jedince v hodnocené kohortě. Hakulinenova metoda: zařazení do skupiny osob v riziku bere ohled na cenzorování, ale v případě úmrtí je srovnatelná skupina z obecné populace v riziku až do ukončení sledování. Metody jsou sice výpočetně odlišné, nicméně pokud odhadujeme pouze krátkodobé přežití (např. 5leté relativní přežití), jejich výsledky jsou velmi podobné. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 109 / 121 Relativní přežití Interval spolehlivosti pro relativní přežití Bodový odhad relativního přežití v čase t je nutné doplnit intervalovým odhadem. Pointou výpočtu je, že očekávané přežití bereme za konstantní (jeho variabilitu zanedbáváme). var (︀ SR (t) )︀ = var (︂ S(t) S*(t) )︂ = var(S(t)) S*(t)2 = SE(S(t))2 S*(t)2 . Můžeme si to dovolit? Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 110 / 121 Relativní přežití Intervalově specifické relativní přežití Relativní přežití nejčastěji počítáme pomocí metody úmrtnostních tabulek. Rozdělíme-li délku sledování do l intervalů, lze kumulativní formu relativního přežití vyjádřit jako SR (t) = l∏︁ i=1 pR i , kde pR i je tzv. intervalově specifické relativní přežití. Hodnoty intervalově specifického relativního přežití souvisí s pojmem statistické vyléčení. Jak? Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 111 / 121 Relativní přežití Statistické vyléčení Ve chvíli, kdy se odhady intervalově specifického relativního přežití dosáhnou hodnoty 1, lze říci, že se mortalita sledovaných pacientů v daném intervalu dostala na úroveň mortality populační. V tomto případě pak mluvíme o tzv. statistickém vyléčení. Nelze zaměňovat pojmy klinické a statistické vyléčení. Pojem klinické vyléčení chápeme na úrovní jedince jako vymizení všech klinických projevů nemoci. Pojem statistické vyléčení chápeme na úrovni skupiny pacientů jako srovnání mortality s populační úrovní. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 112 / 121 Relativní přežití Statistické vyléčení - příklad Cancer of the stomach r 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 Annual follow-up interval 1 2 3 4 5 6 7 8 9 10 Cancer of the breast r 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 Annual follow-up interval 1 2 3 4 5 6 7 8 9 10 Zdroj: Dickman PW. (2004) Population-based cancer survival analysis. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 113 / 121 Relativní přežití Statistické vyléčení - příklad Cancer of the stomach rsr 0.0 0.2 0.4 0.6 0.8 1.0 Annual follow-up interval 0 1 2 3 4 5 6 7 8 9 10 Cancer of the breast rsr 0.0 0.2 0.4 0.6 0.8 1.0 Annual follow-up interval 0 1 2 3 4 5 6 7 8 9 10 Zdroj: Dickman PW. (2004) Population-based cancer survival analysis. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 114 / 121 Modely s podílem vyléčených pacientů Modely s podílem vyléčených pacientů Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 115 / 121 Modely s podílem vyléčených pacientů Modely s podílem vyléčených pacientů Pointou modelů s podílem vyléčených pacientů (cure fraction models) je rozdělení sledovaných pacientů (na základě rizika úmrtí) do dvou skupin: Nevyléčení pacienti, jejichž úmrtí lze přičítat sledovanému onemocnění. Vyléčení pacienti, u nichž k úmrtí ve sledovaném období nedošlo (nebo případně zemřeli z jiných příčin, než bylo sledované onemocnění). Jsou-li v dané kohortě pacientů přítomni tzv. statisticky vyléčení pacienti, křivka funkce přežití, která jde zpravidla k nule, se narovná k pomyslné nenulové asymptotě. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 116 / 121 Modely s podílem vyléčených pacientů Modelování celkového a relativního přežití Modely s podílem vyléčených pacientů lze použít jak pro modelování celkového přežití (využití v klinických analýzách), tak i relativního přežití (využití v populačních analýzách). Nejpoužívanějšími modely jsou tzv. smíšené modely, které jsou definovány jako Model pro celkové přežití: S(t) = c + (1 − c)SU(t), kde c je tzv. podíl statisticky vyléčených pacientů a SU(t) je funkce přežití pro tzv. nevyléčené pacienty (uncured, bound to die). Model pro relativní přežití: S(t) = S* (t)SR (t) = S* (t) (c + (1 − c)SU(t)) , kde S* (t) je očekávané přežití odpovídající populace. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 117 / 121 Modely s podílem vyléčených pacientů Smíšené a nesmíšené modely Kromě smíšeného modelu s podílem vyléčených pacientů byly definovány i tzv. smíšené modely s podílem vyléčených pacientů, které jsou definovány jako Model pro celkové přežití: S(t) = cF(t) , kde c je tzv. podíl statisticky vyléčených pacientů a F(t) je vhodně zvolená distribuční funkce. Model pro relativní přežití: S(t) = S* (t)cF(t) , kde S* (t) je očekávané přežití odpovídající populace. Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 118 / 121 Modely s podílem vyléčených pacientů Používaná rozdělení pravděpodobnosti Lze definovat i semi-parametrické modely s podílem vyléčených pacientů, ale většina těchto modelů je parametrických s některým z následujících rozdělení pravděpodobnosti: Weibullovo rozdělení - 2 parametry Log-normální rozdělení - 2 parametry Gamma rozdělení - 3 parametry Směs dvou Weibullových rozdělení Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 119 / 121 Modely s podílem vyléčených pacientů Co vlastně u těchto modelů modelujeme? Naším cílem je modelování charakteristik týkajících se 1 vyléčených pacientů - modelujeme c, tedy podíl statisticky vyléčených pacientů. 2 nevyléčených pacientů - modelujeme parametry rozdělení pravděpodobnosti (např. 𝜆 a 𝛾 v případě Weibullova rozdělení). Odhad regresních koeficientů je opět založen na metodě maximální věrohodnosti. Podle toho, co všechno může záviset na vysvětlujících proměnných, roste náročnost na počet pozorování (sledovaných událostí). Tomáš Pavlík (MU, Brno) Analýza přežití 9. května 2012 120 / 121