ניתוח הישרדות

**ניתוח הישרדות** (באנגלית: Survival analysis) הוא ענף בסטטיסטיקה, העוסק באופן מסורתי בניתוח במשך הזמן הצפוי לעבור עד שיתרחש אירוע מסוים, בטרמינולוגיה סטטיסטית - זמן-למאורע (Time-to-event). אירוע יכול להיות, לדוגמה, מוות של חולה או קלקול במערכת מכנית כגון מכונה כלשהי. הטיפול הייחודי באמדן של זמן-למאורע נובע מהקושי לאסוף מדגמים בהם מועד התרחשות האירוע ידוע עבור כל התצפיות. עבור התצפיות שלא נאסף עבורם מועד ההתרחשות, מודלי הישרדות מתחשבים בזמן שידוע שעד אליו האירוע לא התרחש. בעזרת ניתוח הישרדות ניתן לענות על שאלות כגון: איזה אחוז מאוכלוסיית החולים יישארו בחיים לאחר פרק זמן מסוים? איזה גורמים משפיעים על הסיכון למות? האם טיפול תרופתי מסוים מאריך את תוחלת החיים? כמה אנשים מתים מנזקי העישון?

באופן מסורתי, ניתוח הישרדות עוסק במידול מאורע בודד. עם זאת, עם פיתוח הענף פותחו מבחנים סטטיסטים העוסקים במידול מצבים בהם האירוע הוא חוזרני (*Recurring event*), כלומר יכול להתרחש יותר מפעם אחת. מצבו הקליני של החולה יכול להחמיר מספר פעמים במשך חייו, ומכונה יכולה להתקלקל, לחזור לפעילות לאחר תיקון ואז להתקלקל שוב.

אם לכל תצפית ידוע הזמן-למאורע, ניתן להשתמש בכלים סטנדרטיים לניתוח מערכים עם משתנה תלוי רציף, תוך התייחסות לכך שהמשתנה התלוי תחום בחלקו החיובי של הציר הממשי \R^+.

## מושגי יסוד בניתוח הישרדות
* אירוע או מאורע (Event): התרחשות מוגדרת היטב של התופעה בה אנו מתעניינים. 
* תקופת המעקב (Follow-up period): התקופה בה נערך מעקב אחרי הפרט או אוכלוסייה של פרטים כדי לבדוק האם מתרחש אירוע. לדוגמה, בניסוי קליני להערכת יעילות טיפול מסוים לטרשת נפוצה, מתעניינים במשך הזמן העובר מתחילת הטיפול עד להתרחשות התקף. בדרך כלל מגבילים את משך הניסוי הקליני בזמן, למשל מחליטים לעקוב אחרי החולים במשך שנתיים בלבד לאחר תחילת הטיפול, כלומר תקופת המעקב היא שנתיים.
* זמן שרידות אמיתי (True survival time), לעיתים נקרא גם זמן לכישלון (Failure time) או זמן המאורע (Event time): משך הזמן העובר מנקודת ייחוס בזמן ועד להתרחשות המאורע עבור תצפית מסוימת. לרוב יסומן זמן זה באמצעות T_i.
* צינזור (Censoring): אם תקופת המעקב הסתיימה עבור פרט מסוים לפני התרחשות האירוע, נאמר כי התצפית היא מצונזרת (Censored observation). במקרה של תצפית מצונזרת, אנו יודעים בוודאות כי האירוע לא התרחש בזמן תקופת התצפית, אך לא יודעים אם הוא יתרחש בעתיד ואם כן מתי. הזמן הנמדד עבור תצפית זו הוא הזמן מתחילת המעקב עד זמן הצנזור. לרוב זמן זה עבור התצפית ה-i מיוצג על ידי C_i, והוא נקרא "זמן צנזור אמיתי" (True censoring Time). אם התצפית היא מצונזרת נהוג להוסיף סימן + בצידו הימני של הזמן-למאורע על מנת לסמן זאת. לדוגמה:7^+ מייצג תצפית מצונזרת, שעבורה ידוע שהזמן-למאורע גדול מ-7.
* תצפיות שבסיכון (At risk observations): כלומר כל התצפיות עד זמן y שטרם עברו צינזור, או ש-y קטן מהזמן-לאירוע שלהן. בכתיב מתמטי כל התצפיות i', שעבורן i':y_{i}\geq y . מושג הצנזור קרוב למושג קיטום (truncation), המתאם מדגם שבאופן סיסטמטי לא מכיל תצפיות מסוימות. לדוגמה מדגם שנאסף מבתי חולים ציבוריים לא מכיל כלל מידע על תצפיות מבתי חולים פרטיים, ולכן בהתייחס אל תצפיות מבתי חולים פרטיים הוא קטום. 
* באיסוף התצפיות החוקר מבחין בזמן השרידות האמיתי T_i או בזמן הצנזורC_i, המוקדם מבניהם. דוגמה: חולה המשתתף בניסוי קליני שמשכו שנתיים יכול להחליט לפרוש מהניסוי בכל עת. אם פרש מהניסוי אחרי חצי שנה, למשל, וזאת בטרם התרחשות האירוע, אז הזמן של אותו חולה הוא חצי שנה. מתמטית נהוג להגדיר את הזמן הנמדד בפועל כמשתנה מקרי עבור התצפית ה-i כ: Y_i= \min (C_i,T_i). נהוג גם להגדיר פונקציית סטטוס \delta_i \in \{0,1\}, מחזירה \delta_i =1 אם תצפית מסוימת i עברה אירוע בזמן תקופת המעקב, כלומר אם T_i \leq C_i. כמו כן, הפונקציה מחזירה \delta_i =0, אם תצפית מסוימת i היא מצונזרת כלומר T_i > C_i. במדגם אופייני לניתוח הישרדות החוקר מחזיק ב-n צמדים של (Y,\delta), המציינים את האוסף (y_1,\delta_1), ..., (y_n,\delta_n).
* פונקציית ההישרדות S(t): פונקציה שמחזירה עבור כל זמן t את ההסתברות שזמן השרידות האמיתי T גדול מ-t. פורמלית, S(t)=P(T>t), פונקציית ההישרדות היא פונקציה מונוטונית יורדת של t, המחזירה ערכים שתחומים בין 1 ל-0.
* פונקציית הסיכון (hazard function), המסומנת בדרך כלל ב-h(t), מבטאת את הסיכון הנקודתי להתרחשות אירוע בזמן t, והיא מתקבלת כאשר מחשבים את ההסתברות המותנית להתרחשות אירוע במרווח זמן מסוים \Delta t כאשר נתון כי האירוע לא התרחש לפני כן, מחלקים בגודל מרווח הזמן, ומשאיפים גודל זה ל-0: \lambda(t) =\lim_{\Delta t \to 0} \frac {P(t t)} {\Delta t } = \frac {f(t)} {S(t)},
כאשר f(t) היא פונקציית הצפיפות של T.

יש לשים לב כי הסיכון אינו הסתברות, וערכו יכול להיות גדול מ-1.

## גישות נפוצות לתיאור ניתוח השרדות
### לוח חיים
לוח חיים (או לוח תמותה, תלוי בהקשר) הוא טבלה תיאורית של נתוני ההישרדות. לכל נקודת זמן או מרווח זמן נתונים:
* מספר הפרטים באוכלוסייה שנמצאו בסיכון באותה נקודת זמן. גודל זה הוא מספר הפרטים שלא נצפה כי הם חוו אירוע עד נקודת זמן זו, בין אם הם חוו אירוע בנקודת זמן מאוחרת יותר ובין אם התצפית נקטמה.
* מספר הפרטים שחוו את האירוע בנקודת זמן זו.
משני נתונים אלה אפשר לחשב אמדן לערך פונקציית ההישרדות באותה נקודת זמן. על ידי שימוש בהתפלגות הבינומית ניתן לחשב רווח סמך לאמדן זה.

### עקומת קפלן-מאייר
עקומת קפלן-מאייר היא קירוב אפרמטרי לפונקציית ההישרדות, המתחשב בדגימות עם צנזור באמצעות הכללתם במכנה הביטוי לכל זמן  t, כאשר C הוא פרק הזמן עד לצנזור. האמדן נתון על ידי הנוסחה:
 \widehat S(t) = \prod\limits_{i:\ t_i\le t} \left(1 - \frac{d_i}{n_i}\right)
כאשר
* t_i היא נקודת זמן בה התרחש לפחות אירוע אחד.
* d_i הוא מספר האירועים שהתרחשו בנקודת הזמן t_i.
* n_i הוא מספר הפרטים באוכלוסייה ששרדו (כלומר לא חוו אירוע) עד זמן t_i.

קפלן ומאייר הוכיחו כי אמד זה הוא אמד נראות מקסימלית עבור פונקציית ההישרדות.

## שיטות להסקה בניתוח הישרדות

### מבחן לוג הדרגות (log rank test)
זהו מבחן סטטיסטי א-פרמטרי לבדיקת ההשערה כי אין הבדל בין פונקציות ההשרדות של שתי אוכלוסיות בלתי תלויות. מבחן זה נקרא לעיתים בשם מבחן מאנטל-קוקס. התפלגותו האסימפטוטית של סטטיסטי המבחן היא התפלגות חי בריבוע כאשר מספר דרגות החופש שווה למספר האוכלוסיות פחות 1.

### רגרסיית הישרדות פרמטרית
זהו מודל רגרסיה פרמטרי המקשר בין משך הזמן (או פונקציה שלו) ובין משתנים מסבירים. פורמלית:
g(T) = \beta_0 +\beta_1 X_1 +...+\beta_k X_k + \sigma\epsilon
כאשר
* g(T) היא פונקציה מונוטונית של T
* X_1,...,X_k הם משתנים מסבירים.
* \beta_1,...,\beta_k הם מקדמי הרגרסיה.
* ל-\epsilon יש התפלגות המקבלת ערכים אי שליליים.
* \sigma הוא מקדם מישקול (scaling) התלוי בהתפלגות של \epsilon.

אמידת הפרמטרים נעשית על ידי שיטת הנראות המקסימלית.

### מודל הסיכונים הפרופורציונליים (מודל קוקס)
מודל הסיכונים הפרופורציונליים הוא מודל רגרסיה סמי-פרמטרי המקשר בין פונקציית הסיכון ובין משתנים מסבירים. בעוד שעקומת קפלן מאייר ומבחן לוג הדרגות מתאימים למצב בו יש משתנה קטגורי יחיד המשפיע על פונקציית הסיכון (כגון קבוצת טיפול וקבוצת ביקורת, או מעשנים ולא מעשנים), מודל קוקס מתאים גם למצבים בהם יש יותר ממשתנה אחד המשפיע על ערך פונקציית הסיכון, ומשתנים אלה אינם מוגבלים להיות משתנים קטגוריים, אלא יכולים להיות משתנים בכל סולם מדידה. 

המודל הוצג על ידי סיר דויד קוקס בשנת 1972, והמאמר בו הוא הוצג נמנה עם 100 המאמרים המצוטטים ביותר בספרות המדעית.

ההגדרה הפורמלית של המודל היא
h(t | X_1, ..., X_k) = h_0(t) \cdot \exp (\beta_0+\beta_1X_1+...+\beta_k X_k)
כאשר
* X_1,...,X_k הם משתנים מסבירים.
* h(t | X_1, ..., X_k) הוא הסיכון בנקודת הזמן t בהינתן X_1,...,X_k.
* h_0(t) הוא הסיכון הבסיסי (הבלתי מותנה) בנקודת הזמן t.
* \beta_1,...,\beta_k הם מקדמי הרגרסיה.
המודל אינו מניח הנחה התפלגותית לגבי הזמן T, אך מניח כי יחס הסיכונים הוא קבוע, ולא תלוי בזמן. הנחה זו מופרת כאשר משווים שתי קבוצות של נבדקים, כאשר עד זמן t פונקציית הסיכון של קבוצה אחת גבוהה יותר מפונקציית הסיכון של קבוצה שנייה, אך לאחר זמן t היא נמוכה יותר ביחס לקבוצה השנייה.  

בגרסתו הפשוטה המודל גם מניח כי המשתנים המסבירים אינם תלויים בזמן. עם זאת, ניתן בקלות להרחיב את המודל כך שיטפל גם במשתנים מסבירים מסוג זה.  

כמו כן, המודל אינו אומד ישירות את פונקציית הסיכון אלא את לוג יחסי הסיכונים (hazard ratio):\log \biggr( \frac {h(t | X_1, ..., X_k)} {h_0(t)} \biggr). כלומר, במודל קוקס \hat\beta_i הוא אמדן ללוג יחס הסיכונים (HR) של המשתנה המסביר X_i בהינתן שאר המשתנים המסבירים אל מול הסיכון בביסליין h_0. לפיכך, שימוש במודל אינו מתאים אם לחוקר יש עניין רב בהבנת פונקציית הסיכון בבייסליין, h_0. למשל במקרים בהם החוקר מעוניין לחזות את הסיכון לתמותה לאחר זמן t עבור נבדק מסוים. 

אמידת הפרמטרים \beta_1,...,\beta_k נעשית בשיטת הנראות המקסימלית החלקית (Partial Maximum Likelihood), אשר אינה דורשת אמידה של h_0(t).   

#### נראות מקסימלית חלקית במודל קוקס
נניח שהתצפית ה-i לא מצונזרת. כלומר עבורה פונקציית הסטטוס \delta_i = 1. לפיכך y_i= \min (C_i,T_i)=T_i, כלומר הזמן שנמדד במדגם עבור תצפית זו הוא זמן השרידות האמיתי. 

במקרה זה פונקציית הסיכון לתצפית ה-i היא: h(y_i|x_i)=h_0(y_i)\exp  \Bigl( \textstyle \sum_{j=1}^ p x_{ij}
\beta_j \displaystyle\Bigr), והסיכון הכולל בזמן y_i לכל התצפיות שבסיכון (at risk): \sum_{i':y_{i'}\geq y_i} h_0(y_i) \exp \Bigl( \textstyle \sum_{j=1}^ p x_{i'j}
\beta_j \displaystyle\Bigr) . לכן, ההסתברות שהתצפית ה-i היא זו שתיפול בזמן y_i, היא היחס:  

\frac{ h_0(y_i)\exp \Bigl(\textstyle \sum_{j=1}^ p x_{ij}
\beta_j \displaystyle\Bigr)}{\sum_{i':y_{i'}\geq y_i}  h_0(y_i) \exp \Bigl( \textstyle \sum_{j=1}^ p x_{i'j}
\beta_j \displaystyle\Bigr) }. ניתן לראות שהביטוי h_0(y_i) ניתן לצמצום מהמונה והמכנה.

הנראות המקסימלית החלקית מוגדרת כמכפלה מעבר לכל התצפיות הלא מצונזרות של ביטוי זה:

PL (\beta)= \prod_{i:\delta_i=1}\frac{\exp \Bigl( \textstyle \sum_{j=1}^ p x_{ij}
\beta_j \displaystyle\Bigr)}{\sum_{i':y_{i'}\geq y_i} \exp \Bigl( \textstyle \sum_{j=1}^ p x_{i'j}
\beta_j \displaystyle\Bigr) }

והבטאות של מודל קוקס הן אלו אשר ממקסמות אותו.

### גישות נוספות
גישות אפשריות נוספות לניתוח נתוני השרדות כוללות את מודל הסיכונים המצטברים, מודלים בייסיאניים ומודלים של למידת מכונה. כמו כן ישנן הרחבות למודלים שהוזכרו כאן למקרים בהם מופרות חלק מההנחות.

ישנם מצבים בהם אירוע יכול להתרחש בצורות שונות (למשל אדם מעשן נמצא בסיכון גבוה גם לסרטן וגם למחלת לב. הוא עלול למות כתוצאה מסרטן או ממחלת לב, אך סיבת המוות לא יכולה להיות שתי המחלות יחד). המודל המתאים לניתוח הנתונים במצבים כאלה הוא מודל הסיכונים המתחרים.

## דוגמה
בדוגמה זו נשתמש בקובץ הנתונים lung המצורף לחבילתsurvival של תוכנת R. הניתוחים המוצגים בוצעו על ידי פונקציות מחבילה זו.

יש להדגיש כי הדוגמה ממחישה את האופן בו מתבצעים הניתוחים השונים ומסבירה כיצד לפרט את התוצאות, אך אין זו דוגמה לתהליך ניתוח מסודר של נתוני הישרדות.

### הנתונים
בקובץ lung יש נתונים אודות 228 חולי סרטן הריאה. הנתונים כוללים את מין וגיל הנבדקים, וכן נתונים נוספים אודות מצבם התפקודי והתזונתי של החולים. כן נתון משך הזמן בימים שעבר ממועד איסוף הנתונים (תחילת התצפית) ועד התרחשות אירוע מוות או סיום התצפית, ומשתנה המציין האם התצפית היא מלאה או קטומה. לצורך הדוגמה נשתמש בנתוני המין והגיל בלבד. להלן חלק מהנתונים:

Sex=1 מציין גבר ו-sex=2 מציין אישה. Status=1 מציין תצפית קטומה ו-status=2 מציין התרחשות אירוע מוות.

מתוך 228 חולים, 165 נפטרו במהלך תקופת התצפית.

### חציון זמן ההישרדות
מכיוון שאנו איננו יודעים את משך הזמן עד המוות של החולים שלא נפטרו אי אפשר לאמוד את הזמן הממוצע עד המוות, אבל אפשר לאמוד את משך הזמן החציוני בעזרת אמדן קפלן-מאייר לפונקציית ההישרדות. בעזרת הפונקציה suvrfit של R עולה כי משך הזמן החציוני עד למוות באוכלוסייה שווה ל-310 ימים. חישוב נפרד לגברים ונשים מעלה כי משך הזמן החציוני עד למוות עבור גברים הוא 270 יום, ועבור נשים הוא 426 ימים.

### לוח חיים (או לוח תמותה)
מוצגות פקודות R ו-6 השורות הראשונות של לוח התמותה (נתוני התמותה חושבו לכל החולים יחד):

מהשורה הראשונה של הלוח ניתן ללמוד כי:
•	ביום החמישי של תקופת המעקב היו 228 חולים בסיכון
•	חולה אחד מת ביום החמישי
•	ערך פונקציית ההשרדות עבור t=5 שווה ל-0.996. הערך חושב על ידי חלוקת 227 (מספר השורדים מעבר ליום ה-5) ב-228 (מספר החולים הכולל)

### עקומת קפלן-מאייר
ניתן להפיק את העקומה ב-R בעזרת הפקודות הבאות:

העקומה המתקבלת היא:

מרכז

ניתן לראות כי ערך פונקציית ההשרדות של הנשים גבוה מזה של הגברים ברוב הזמן. לדוגמה, ההסתברות כי אישה תשרוד יותר מ-400 ימים קרובה ל-0.6, בעוד עבור הגברים הסתברות זו היא בערך 0.3.

### מבחן לוג הדרגות (log rank test)
ביצוע המבחן ב-R:

ערך \chi^2 שהתקבל הוא 10.3 וערך ה-p הוא 0.00131. בהנחה כי רמת המובהקות שנקבעה מראש היא \alpha=0.05, ניתן לומר כי יש הבדל מובהק בין פונקציות ההישרדות של הגברים ושל הנשים.

### רגרסיית השרדות פרמטרית
נבצע רגרסיית הסתברות פרמטרית כאשר המשתנה המוסבר הוא הלוג של משך הזמן והמשתנים המסבירים הם מין וגיל. אנו מניחים כי לטעות \epsilon יש התפלגות וויבול עם פרמטר scale השווה ל-1.

טיב ההתאמה של המודל לנתונים נבדק על ידי סטטיסטי \chi^2 של מבחן יחס הנראות (על פי משפט וילקס). ערכו הוא 7.15 עם שתי דרגות חופש, וערך ה-p הוא 0.028. בהנחה כי רמת המובהקות שנקבעה מראש היא \alpha=0.05, ניתן לומר כי יש הבדל מובהק בין המודל הכולל את המשתנים המסבירים והמודל הבסיסי שאינו כולל משתנים מסבירים. עם זאת, אנו רואים כי המקדם של משתנה המין שונה מ-0 באופן מובהק סטטיסטית, אך המקדם של משתנה הגיל אינו שונה מאפס באופן מובהק.

### מודל הסיכונים הפרופורציונליים (מודל קוקס)
נבצע הרצה של המודל עם משתנה מסביר אחד - משתנה המין:

בתחתית הפלט יש מספר מבחנים לבדיקת טיב ההתאמה של המודל לנתונים בהשוואה למודל ללא משתנים מסבירים: מבחן התאימות (Concordance), מבחן יחס הנראות, מבחן ואלד ומבחן לוג הדרגות. בכל המבחנים התקבלו תוצאות מובהקות (ברמת מובהקות של \alpha=0.05). כן נתון ערך המתאם R^2 והוא שווה ל-0.046. ערך זה נמוך למדי, ורומז על כך שייתכן כי ישנם עוד משתנים מסבירים היכולים להסביר את השונות בין זמני ההישרדות.
ערכו של מקדם משנה המין הוא 0.531-, וערך זה מובהק סטטיסטית ברמת מובהקות של \alpha=0.05. מאחר שגברים מקודדים בקובץ הנתונים על ידי המספר 1 ונשים על ידי המספר 2, נובע כי כשאר ערך המין משתנה מ-1 ל-2, כלומר כאשר הסיכון מחושב עבור אישה ולא עבור גבר, הסיכון הכללי קטן. כאשר מפעילים את הפונקציה המעריכית על מקדם זה, אנו מקבלים כי יחס הסיכונים הוא 0.588. פירוש הדבר הוא כי הסיכון למוות של נשים נמוך בג-40% מזה של הגברים בכל נקודת זמן. רווח סמך ברמת סמך של 95% ליחס הסיכונים הוא (0.4237-0.0816) ואינו מכיל בתוכו את 1, וזה מתיישב עם הערך המובהק של מקדם הגיל ברגרסיה.

ניתן להפיק גרף של פונקציות ההישרדות שנאמדו על ידי המודל בעזרת הפקודות הבאות ב-R:

העקומה המתקבלת היא:
מרכז
שוב ניתן לראות כי ערכי פונקציית ההשרדות של הנשים גבוהים מאלה של הגברים. החל מנקודה מסוימת שתי העקומות מקבילות זו לזו, וזאת כתוצאה של הנחת הסיכונים הפרופורציונליים.

## לקריאה נוספת

## קישורים חיצוניים

* 
* 

## הערות שוליים

קטגוריה:ערכים שבהם תבנית בריטניקה אינה מתאימה
קטגוריה:סטטיסטיקה
קטגוריה:מודלים סטטיסטיים
קטגוריה:הישרדות