95%-е доверительные интервалы на графике выживаемости

* Тема: Доверительные интервалы доли выживших на графике функции выживаемости.
* Ключевые слова: анализ выживаемости, графики, доверительные интервалы.
* Опубликован: 21.03.2008.
* Автор: Marta Garcia-Granero.
* Перевод: А. Балабанов.
* Размещение: http://www.spsstools.ru/Syntax/SurvivalAnalysis/Show95pcCIforFailurePointsOnSurvivalPlot.txt (.sps).

(Вопрос) Можно ли отобразить 95% доверительный интервал
 на графике функции выживаемости?

**************************
(Ответ) Для этого потребуется сначала подсчитать значения функции выживаемости,
 стандартные ошибки и границы доверительного интервала, а затем все это нарисовать
 на одном и том же графике.
* Автор решения: Marta Garcia-Granero, 06.18.2004.

* (1) Решение для одной группы наблюдений *.

* Пример данных *.
DATA LIST FREE /surtime(F8.0) status(F4.0).
BEGIN DATA
65 1 156 0 100 1 134 1 16 1 108 1 121 1 4 1
69 1 143 0 56 1 26 1 22 1 1 1 1 1 5 1 65 1
65 1 17 1 7 1 16 1 22 1 3 1 4 1 2 1 3 1 8 1
4 1 3 1 30 1 4 1 43 1 56 1
END DATA.
VALUE LABELS status 0 'Цензурированное' 1 'Событие'.

* Процедура KM с опцией /SAVE создаст 2 переменные: значения функции выживаемости (sur_1)
  и её стандартную ошибку (se_1). Для ориентира выведем также график функции*.
KM surtime /STATUS=status(1)
  /PRINT TABLE MEAN
  /PLOT SURVIVAL
  /SAVE SURVIVAL SE .

* Для цензурированных данных значения функции и её стандартная ошибка не
  вычисляются. Мы припишем им значения из ближайших предыдущих наблюдений *.
SORT CASES BY surtime(A).
IF (MISSING(sur_1)) sur_1 = LAG(sur_1) .
IF (MISSING(se_1)) se_1 = LAG(se_1) .

* Вычисляем нижние и верхние границы интервала, заменяя 
  границы, выходящие за пределы 0 или 1 на теоретический минимум/максимум *.
COMPUTE lower = MAX((sur_1-1.96*se_1),0) .
COMPUTE upper = MIN((sur_1+1.96*se_1),1) .
EXECUTE.

* Теперь рисуем доверительные уровни вместе с функцией выживаемости
  (сравните форму графика с формой графика из процедуры KV) *.
GRAPH
 /SCATTERPLOT(OVERLAY)=surtime WITH sur_1 lower upper.

* Далее следует "ручное" вмешательство (настройка графика).
* Потребуется: ввести интерполирующую линию, определив ей
  характер ступенчатой функции (Left Step), а также, по желанию,
  убрать маркеры.
  
* Кроме того, вы можете добавить заголовки к осям, поменять цвета и т.д.

* Можете также сохранить шаблон графика, чтобы в будущем меньше
  возиться с ручной настройкой. Использовать этот шаблон следует
  так (пример):
* GRAPH
   /SCATTERPLOT(OVERLAY)=surtime WITH sur_1 lower upper
   /TEMPLATE='C:\\temp\\surv.sgt'.


**************************
* (2) Решение для двух групп наблюдений. Несколько более сложное, т.к.
  процедура SCATTERPLOT(OVERLAY) не принимает группировку через ключевое поле BY *.

* Пример данных из учебника "Statistics at Square One" (Swinscow&Campell) *.
DATA LIST FREE/treatmen(F8.0) surtime(F8.1) status(F8.0).
BEGIN DATA
1 1 0 1 5 0 1 6 1 1 6 1 1 9 0 1 10 1 1 10 1 1 10 0 1 12 1
1 12 1 1 12 1 1 12 1 1 12 0 1 13 0 1 15 0 1 16 0 1 20 0
1 24 1 1 24 0 1 27 0 1 32 1 1 34 0 1 36 0 1 36 0 1 44 0
0 3 0 0 6 1 0 6 1 0 6 1 0 6 1 0 8 1 0 8 1 0 12 1 0 12 1
0 12 0 0 15 0 0 16 0 0 18 0 0 18 0 0 20 1 0 22 0 0 24 1
0 28 0 0 28 0 0 28 0 0 30 1 0 30 0 0 33 0 0 42 1
END DATA.

VAR LABELS treatmen 'Характер лечения'
 /surtime 'Время жизни (мес.)'
 /status 'Статус'.
VALUE LABELS treatmen 1'Gamma-Linoleic' 0'контрольная группа'
  /status 0 'Жив' 1 'Мёртв'.

* Получим оценки функции выживаемости и её стандартные ошибки для обеих групп *.
TEMPORARY.
SELECT IF (treatmen EQ 0).
KM surtime
  /STATUS=status(1)
  /PRINT TABLE MEAN
  /SAVE SURVIVAL SE .
TEMPORARY.
SELECT IF (treatmen EQ 1).
KM surtime
  /STATUS=status(1)
  /PRINT TABLE MEAN
  /SAVE SURVIVAL SE .

* Вычисление пределов 95%-го ДИ *.
DO REPEAT A=lower1 lower2 /B=upper1 upper2 /C=sur_1 sur_2 /D=se_1 se_2.
- COMPUTE A = C-1.96*D .
- COMPUTE B = C+1.96*D .
END REPEAT.
EXECUTE.

* Коррекция границ интервалов, выходящих за теоретические пределы функции (0-1) *.
DO IF (lower1 LT 0) AND (treatmen EQ 0).
- COMPUTE lower1=0.
END IF.
DO IF upper1 GT 1 AND (treatmen EQ 0).
- COMPUTE upper1=1.
END IF.
DO IF lower2 LT 0 AND (treatmen EQ 1).
- COMPUTE lower2=0.
END IF.
DO IF upper2 GT 1 AND (treatmen EQ 1).
- COMPUTE upper2=1.
END IF.
EXECUTE.

* Заполнение пропущенных значений для цензурированных данных *.
SORT CASES BY treatmen(A) surtime(A) .
DO IF ($casenum EQ 1) AND (MISSING(sur_1)).
- COMPUTE sur_1=1.
- COMPUTE lower1=1.
- COMPUTE upper1=1.
END IF.
DO IF MISSING(sur_1) AND (treatmen EQ 0).
- COMPUTE sur_1=LAG(sur_1).
- COMPUTE lower1=LAG(lower1).
- COMPUTE upper1=LAG(upper1).
END IF.
SORT CASES BY treatmen(D) surtime(A) .
DO IF ($casenum EQ 1) AND (MISSING(sur_2)).
- COMPUTE sur_2=1.
- COMPUTE lower2=1.
- COMPUTE upper2=1.
END IF.
DO IF MISSING(sur_2) AND (treatmen EQ 1).
- COMPUTE sur_2=LAG(sur_2).
- COMPUTE lower2=LAG(lower2).
- COMPUTE upper2=LAG(upper2).
END IF.
EXEC.


VARIABLE LABELS
 sur_1 'Функция1'
 sur_2 'Функция2'
 lower1 'НГ95%ДИ/ф-я1'
 upper1 'ВГ95%ДИ/ф-я1'
 lower2 'НГ95%ДИ/ф-я2'
 upper2 'ВГ95%ДИ/ф-я2'.
GRAPH
 /SCATTERPLOT(OVERLAY)=surtime WITH sur_1 lower1 upper1 sur_2 lower2 upper2
 /MISSING=VARIABLEWISE.

* Далее проводим аналогичные "ручные" настройки графика (после двойного щелчка мышкой на нём), 
  которые делали в Решении №1.