SAS – симуляционный расчет мощности исследований биоэквивалентности

Цель данной работы – разработать в SAS University edition код для расчета мощности исследований биоэквивалености при известном CV и числе субъектов исследования.

За основу алгоритма была взята статья Определение объема выборки для исследований биоэквивалентности с использованием компьютерного моделирования. Для алгоритма использовался принцип полной генерации данных (т.е. данных всех субъектов), анализ выполнялся с помощью PROC GLM, а мощность определялась как доля успешных симуляций относительно общего количества.

Прототип был написан в SAS Studio и выполнялся в фоновом режиме. В случае успеха можно было бы использовать код для более масштабных симуляций… Но… Судя по логам симуляция на 1000 итераций выполнялась 1 минуту и 14 секунд… Что делает подобный подход полностью непригодным.

Работа выполнялась на процессоре Intel Core i7 7600U – да, не лучший вариант. И сама виртульная машина SAS работает в двухядерном режиме… Но даже в этом случае 1:14, это фиаско, Карл…

Не исключена вероятность, что код написан не оптимально и для симуляций надо использовать какие-либо другие подходы… (Update: а они есть…) Тем не менее код ниже…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

ODS select none; *Выключаем вывод
filename junk dummy; *Выключаем логи
proc printto log=junk; run;

proc datasets lib=work nolist; *Удаляем старые таблицы
delete list;
run;

%macro prog; *Макрос prog циклично генерирует дадасет и подвергает его анализу
%DO c = 0 %TO 1000; *1000 итераций

DATA a (keep = x p sq sb f); *Таблица a - генерируемые данные каждой итерации (x - параметр, p - период, sq - последовательность, sb - субъект, f - препарат)
DO i = 0 to 1; * i - последовательность
DO i1 = 1 to 24; * i1 - субъекты
s = 1 + RAND ('NORMAL')*0.1; * s - генерация среднего для субъекта
DO i2 = 0 to 1; * i2 - период
x = s + RAND ('NORMAL')*0.3; * в каждом периоде добавляем субхъекту внутрииндивидуальную вариацию
p = i2;
sq = i;
sb = i1+i*24;
* Присваем значение "препарат" субъекту в зависимости от периода и последовательности
if p = 0 and sq = 0 then f = 0; else if p = 1 and sq = 0 then f = 1; else if p = 0 and sq = 1 then f = 1; else f = 0;
OUTPUT a; * записываем
END;
END;
END;
RUN;

PROC GLM data=a; * используем GLM для расчетов
CLASS p sq f sb;
MODEL x = sq sb(sq) p f;
LSMEANS f / PDIFF CL ALPHA=0.1; * получем доверительный интервал для LSMEANS f
ods output LSMeanDiffCL=LSMeanDiffCL; * выводим в таблицу LSMeanDiffCL
RUN;

proc datasets lib=work nolist; * удаляем данные субъектов
delete a;
run;

DATA list; * добавляем из таблицы LSMeanDiffCL данные в таблицу List
set list work.lsmeandiffcl;
RUN;

%END;
%mend prog;
%prog;

DATA list (keep=a); *отмечаем в таблице успешные симуляции если доверительный интервал соответствует конвенционному
set list;
if LowerCL > -0.22314355131 and UpperCL < 0.22314355131 then a=1; else a=0; output;
RUN;

ODS select all; *Включаем вывод
PROC FREQ data=list; *Анализируем таблицу частот.
tables a;
RUN;

*Готово

Update:

Проблема в том, что SAS не эффективно работает с макросами. Каждый цикл макроса запускает PROC и DATA процедуры, запуск которых сам по себе требует больших ресурсов. Код будет работать быстрее в том случае, если минимизировать количество вхождений в PROC и DATA процедуры. Т.е. все итеративные процедуры надо переместить внутрь PROC и DATA. В итоге мы генерируем не 1000 таблиц с данными, мы генерируем одну, но с дополнительным полем для идентификации каждого набора. А анализ выполняем с помощью выражения BY по признаку идентификатора группы данных. В итоге… 1000 итераций заняли 30 секунд… Результат лучше, но скажем так… не позволяет выполнять серьезные симуляции. К счастью 100000 симуляций заняло не 30000 секунд, а всего 229 (3:49). Что позволяет надеяться но что, что при больших объемах процесс оптимизирован. Можно ли считать такую скорость нормальной для выполнения больших симуляций? Думаю нет. К этому надо добавить, что большие симуляции потребуют создания огромных таблиц, которые по своему размеру могут превышать размеры RAM, что существенно замедлит процесс.

В итоге. По быстродействию SAS не может составить хоть какую конкуренцию с кодом написанным на C/C++ или Java, но разработка симуляций на C/C++/Java требует совершенно других ресурсов.

Оптимизированный код:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

ods graphics off;
ods exclude all;
ods results off;
options nonotes;

filename junk dummy;
proc printto log=junk; run;

proc datasets lib=work nolist;
delete list;
run;

DATA a (keep = id x p sq sb f);
v=0;

DO scnt = 1 to 10000;

DO i = 0 to 1;
DO i1 = 1 to 24;
s = 1 + RAND ('NORMAL')*0.1;
DO i2 = 0 to 1;
x = s + RAND ('NORMAL')*0.3;
p = i2;
sq = i;
sb = i1+i*24;
if p = 0 and sq = 0 then
f = 0;
else if p = 1 and sq = 0 then
f = 1;
else if p = 0 and sq = 1 then
f = 1;
else f = 0;
id = scnt;
OUTPUT a;
END;
END;
END;
END;
RUN;

PROC GLM data=a;
BY id;
CLASS p sq f sb;
MODEL x = sq sb(sq) p f;
LSMEANS f / PDIFF CL ALPHA=0.1;
ods output LSMeanDiffCL=LSMeanDiffCL;
RUN;

DATA LSMeanDiffCL (keep=a);
set LSMeanDiffCL;
if LowerCL > -0.22314355131 and UpperCL < 0.22314355131 then a=1; else a=0; output;
RUN;

ods graphics on;
ods exclude none;
ods results on;
options notes;

PROC FREQ data=LSMeanDiffCL;
tables a;
RUN;

Поправка: sb = i1+i*24;