Анализ медицинских данных в R и SAS Опции

[hot_assay]

Здравствуйте!
Буду благодарен за рекомендации по использованию R и SAS в анализе данных клинических исследований (хотя бы минимальные, общего характера, с учётом того, что пользователю, привыкшему к Graphical user interface, нужно въехать в специфику рассматриваемо ПО )) . Интересуют преимущества данных систем по сравнению с другими (например, Statsoft STATISTICA). Какие модули необходимо иметь для решения медицинских задач? Заменяют ли R и SAS друг друга в плане фукционала? Отдельное спасибо за ссылки с практическими примерами по рассматриваемому вопросу.

[TheThing]

Не с целью продолжения холивара на тему языков программирования, а просто обыденный случай.

Решил глянуть, как реализован пакет MDR в R - все хорошо, практически ничем по фичам не уступает аналогу GUI на JAVA http://sourceforge.net/projects/mdr/. Но ВРЕМЯ выполнения!! Массив из 250 наблюдений 25 предикторов - время анализа в R: 221 секунда, MDR Java: 1.2 секунд; массив 250 наблюдений 50 предикторов - R: 1950 секунд!!!, Java: 3,9 секунды. Ну и так далее...50 предикторов в генетических исследованиях - это очень мало, следовательно MDR в R пока сложно пользоваться (мне хотелось еще прикрутить алгоритм в shiny, там вообще браузер зависнет навека (((. Подобные пакеты,как mbMDR и другие работают также..

Решил взять другой пример, вычислить массивный полином (100) 500 000 раз. Код в R (время выполнения 4 минуты)

trpol2 <- function(n,x) {
mu <<- 10.0
pu <<- 0.0
pol <<- 1:100
tp1 <<- 2.0
tm1 <<- 1/2.0
for (i in 1:n) {
for (j in 1:100) {
mu <<- (mu + tp1) * tm1
pol[j] <<- mu
}
s <<- 0.0;
for (j in 1:100) {
s <<- x*s + pol[j];
}
pu <- s+pu;
}
print(pu)
}

trpol2(500000,0.2)

Код на Java, время выполнения в браузере! http://www.browxy.com/ - 2 секунды.

public class HelloWorld {
static public void main(String argv[]) {
float mu = (float)10.0;
float x,s;
float pu = (float)0.0;
int su, i, j, n;
float pol[] = new float[100];

n = 500000;
x = (float)0.2;
for(i=0; i<n; i++) {
for (j=0; j<100; j++) {
mu = (mu + (float)2.0) / (float)2.0;
pol[j] = mu;
}
s = (float)0.0;
for (j=0; j<100; j++) {
s = x*s + pol[j];
}
pu += s;
}
System.out.println(pu);
}
}

Самое смешное, код на Python, время выполнения 6 секунд:

n = 500000
x = 0.2

def t(x):
mu = 10.0
pu = 0.0
pol =[0] * 100
r = range(0,100)

for i in range(0,n):
for j in r:
pol[j] = mu = (mu + 2.0) / 2.0
su = 0.0
for j in r:
su = x * su + pol[j]
pu = pu + su
print pu

t(x)

Немного расстроился..Ладно с этими полиномами, но мне реально хотелось работать с многофакторным уменьшением размерности в R, видимо, не судьба..

Сообщение отредактировал TheThing - 15.09.2014 - 19:00

[p2004r]

Не с целью продолжения холивара на тему языков программирования, а просто обыденный случай.

да да, конечно

1) если так непосильно писать векторизированный код, а писать портянки с рекурсивными форами "массивного полинома" чувствуете себя в силе, то в Вашем распоряжении Rcpp и практически без изменений "код на ява" (в силу весьма ограниченного объема использования возможностей языка ) можно вставить в инлайн вставку кода c++


2) никакого примера о тормозах MDR в письме нет -- и комментировать собственно нечего.

[TheThing]

1) если так непосильно писать векторизированный код, а писать портянки с рекурсивными форами "массивного полинома" чувствуете себя в силе, то в Вашем распоряжении Rcpp и практически без изменений "код на ява" (в силу весьма ограниченного объема использования возможностей языка ) можно вставить в инлайн вставку кода c++

Значит, чтобы решить задачу с полиномом (а таких примеров множество), человек должен знать как R так и с++? В вышеприведенном примере я использовал базовые возможности языка Java, Python и т.д., не пользуясь никакими сторонними библиотеками. Хотелось бы увидеть Вашу PRO-версию векторизированного родного R кода, которая бы могла составить конкуренцию хотя бы Python по скорости выполнения..

2) никакого примера о тормозах MDR в письме нет -- и комментировать собственно нечего.

Прокомментируйте, пожалуйста и запаситесь терпением:

library(MDR)
data(mdr2)

fit.cv<-mdr.cv(data = mdr2, K = 3, cv = 10, ratio = NULL, equal = "HR", genotype = c(0, 1, 2))

Эта база имеет всего лишь 100 предикторов, у меня в реальной работе их тысячи. На Java и в SAS у меня это занимает секунды, а сколько у Вас в R?

Сообщение отредактировал TheThing - 23.10.2014 - 09:35

[p2004r]

Прокомментируйте, пожалуйста и запаситесь терпением:

library(MDR)
data(mdr2)

fit.cv<-mdr.cv(data = mdr2, K = 3, cv = 10, ratio = NULL, equal = "HR", genotype = c(0, 1, 2))

Эта база имеет всего лишь 100 предикторов, у меня в реальной работе их тысячи. На Java и в SAS у меня это занимает секунды, а сколько у Вас в R?

внутри mdr.cv() как и любой другой кроссвалидации находится <s>неонка</s> цикл. В этом легко убедится набрав просто mdr.cv.

Данный цикл Вы должны заменить (в нем нету ни глобального присваивания, ни прочих императивных извратов), как истинный питонист на foreach(). Оная конструкция прекрасно параллелится и на кластеры, и на некластеры. Делать это за Вас мне тоже неинтересно.

[TheThing]

внутри mdr.cv() как и любой другой кроссвалидации находится <s>неонка</s> цикл. В этом легко убедится набрав просто mdr.cv.

Данный цикл Вы должны заменить (в нем нету ни глобального присваивания, ни прочих императивных извратов), как истинный питонист на foreach(). Оная конструкция прекрасно параллелится и на кластеры, и на некластеры. Делать это за Вас мне тоже неинтересно.

С чего Вы взяли, что я истинный питонист?
Вы же сказали, что я не привел никакого кода в R по поводу метода MDR, поэтому обсуждать собственно нечего. Я привел пример своего ужасного неоптимизированного кода и в результате получил ответ чисто теоретического характера, а также что Вам это неинтересно. Почему просто не помочь сообществу в решении абсолютно конкретной задачи конкретным решением?

Сообщение отредактировал TheThing - 23.10.2014 - 21:23

[p2004r]

С чего Вы взяли, что я истинный питонист?
Вы же сказали, что я не привел никакого кода в R по поводу метода MDR, поэтому обсуждать собственно нечего. Я привел пример своего ужасного неоптимизированного кода и в результате получил ответ чисто теоретического характера, а также что Вам это неинтересно. Почему просто не помочь сообществу в решении абсолютно конкретной задачи конкретным решением?

хватит троллить, у вас нет к этому призвания. вы получили исчерпывающий совет, переписать единственный цикл в доступной вам функции. тон которым вы "просите" "помочь сообществу", а вернее себе, не оставляет никаких шансов что кто то кинется вам помогать.

переписать цикл кроссвалидации в неустраивающей вас процедуре тривиальная операция, и её код намного меньше занимает места, чем ваши упражнения в остроумии на этом форуме.

[TheThing]

тон которым вы "просите" "помочь сообществу", а вернее себе, не оставляет никаких шансов что кто то кинется вам помогать.

Поверьте, что количество моих постов и личных сообщений, где я помог людям, значительно превосходят те случаи, когда я просил о помощи
В этом вопросе мне помощь не нужна, я пользуюсь прекрасной реализацией MDR на Java, был интерес, как пользоваться этим алгоритмом в R, но благодаря Вам - интерес пропал, видимо это у нас взаимно

[p2004r]

Поверьте, что количество моих постов и личных сообщений, где я помог людям, значительно превосходят те случаи, когда я просил о помощи
В этом вопросе мне помощь не нужна, я пользуюсь прекрасной реализацией MDR на Java, был интерес, как пользоваться этим алгоритмом в R, но благодаря Вам - интерес пропал, видимо это у нас взаимно

понятно. видимо крутость этих "личных сообщений" недостаточно очевидна для окружающих, что бы вы не могли удержаться от поливания грязью и провокаций на этом форуме.

PS

вот код в mdr() который вызывает проблему

counts[, 1] <- apply(case, 1, compare, mat = part, k = k)
counts[, 2] <- apply(ctrl, 1, compare, mat = part, k = k)

поскольку

case <- cbind(rep(1, g^k), expand.grid(rep(geno, k)))
ctrl <- cbind(rep(0, g^k), expand.grid(rep(geno, k)))

решить это очень просто заменив apply() на эквивалент работающий параллельно --- например parApply(cl, X, MARGIN, FUN, ...)

Или возможно эффективными окажуться
parRapply(cl, x, fun, ...)
parCapply(cl, x, fun, ...)

В результате сколько машин вы поднимете в кластере, столько и будет эффективноть ваших расчетов. В пределе она сойдется к скорости рассылки case и ctrl между нодами кластера.

PPS заметьте, что это не "личное сообщение" и действительно полезно всем.

PPPS

я рыдал (С)

Код

> compare
function (mat, vec, k)
{
    b <- 1
    match <- 1:dim(mat)[1]
    while (b <= (k + 1)) {
        match <- match[mat[match, b] == as.numeric(vec[b])]
        b <- b + 1
    }
    return(length(match))
}
<environment: namespace:MDR>

Сообщение отредактировал p2004r - 24.10.2014 - 19:14

[p2004r]

Код

function (mat, vec, k)
{
    b <- 1
    match <- 1:dim(mat)[1]
    while (b <= (k + 1)) {
        match <- match[mat[match, b] == as.numeric(vec[b])]
        if (length(math)==0) break   ### !!!!!!!
        b <- b + 1
    }
    return(length(match))
}

но на удивление эта фигня довольно оптимально работает

после профилирования кода получается, что реальное ускорение возможно получить только на clusterApplyLB(или иной разумный сплит массива на обработку) вставленных вместо обычных apply(), иначе высокие накладные расходы на пересылку.

PS Оказывается всё уже написали http://cran.r-project.org/web/packages/for...tes/foreach.pdf на стр.11 реализация apply() поверх %dopar% с настраиваемым числом (увеличивать до загрузки всех ядер) передаваемых в неё колонок (поэтому исходную матрицу транспонировать надо будет)



Сообщение отредактировал p2004r - 25.10.2014 - 20:06