Здравствуйте, коллеги! Прошу помощи в анализе данных.
Задача исследования - оценить связь между генотипом (15 SNP), "промежуточным фенотипом" (параметры биохимии, иммунологии и др.), исходом (ЗНО есть/нет).
Существующий сервис "SNPstats" (http://bioinfo.iconcologia.net/SNPstats_web) выдает отношения шансов, "adjustet by фактор1+фактор2+...", используется при этом "logistic regression models" (то есть, логит-регрессию?). Хотелось бы поточнее узнать, что значит "adjusted by".
Кроме того, в данной програме остается "за кадром", какой из факторов является ведущим.
Возможно, есть какие то альтернативные методы анализа, позволяющие оценить вклад конкретных факторов?
Посоветуйте, пожалуйста. Заранее благодарен.

И есть ещё один вопрос:
По разным SNP имеется разное количество генотипированных, как и разное количество известных значений по каждому из "промежуточных фенотипов" и исходов.
То есть, грубо говоря, выборки по каждому из SNP перекрываются только отчасти. Нужно ли в этом случае рассматривать проблему множественных сравнений?
Спасибо!

1) Что означает adjusted..Вы включаете в модель один фактор_1, рассчитываете отношения шансов для этого одного "снипа"..затем создаете новую модель, в которую включаете фактор_2, получаете отношения шансов для фактора_2. Теперь, если Вы включите фактор_1 и фактор_2 в одну модель, то получите скорректированные (adjusted) отношения шансов для фактора_1 и фактора_2, поскольку в модель включаются оба фактора, поэтому OR фактора_1 корректируются относительно присутствия фактора_2 и наоборот. То есть adjusted в данном случае - это просто расчет OR относительно количества предикторов, включенных в модель..

2) Как Вы знаете, снипы (SNPs) сами по себе обладают очень маленькими величинами эффектов, то есть сам по себе 1 полиморфизм не может определить предрасположенность к заболеванию (если Вы только не исследуете моногенные заболевания, в чем я сомневаюсь). Следовательно, риск развития наиболее распространенных полигенных заболеваний, таких как артериальная гипертензия, инфаркт миокарда, сахарный диабет и др., определяется совокупностью влияний снипов, а также факторами внешней среды. Поэтому, изучая вклад каждого снипа без учета взаимодействия (строя аддитивные модели) Вы никогда не докопаетесь до истины. Для этого необходимо изучать межгенные взаимодействия с целью определить, например, усиливает ли Снип1 влияние Снипа 2, какая связь между Снипом1 и Снипом2? Для решения подобных задач был разработан целый спектр подходов и алгоритмов, основные из них приведены на рисунке ниже.
Нажмите для просмотра прикрепленного файла

Из этого списка очень хорошие результаты показывают Multifactorial Dimensionality reduction и Random Forest с модификацией Catherine Strobl. Про первых подход можно ознакомиться на сайте MDR:
http://www.multifactordimensionalityreduction.org/
Статьи про MDR, например:
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1226028/
Вторую прикрепил здесь.
Нажмите для просмотра прикрепленного файла

Метод хороший и работает как для подхода ген-кандидат (как в Вашем случае), так и для полногеномных исследований, когда количество снипов будет от 100 тысяч (здесь логистическая регрессия загнется в результате "curse of dimensionality").
3) Какой из снипов делает наибольший вклад в риск развития заболевания позволит Вам определить, например, Random Forest с модификацией Catherine Strobl ("тетя" из Германии, которая всю жизнь занимается лишь случайными лесами, ей можно доверять ). Погуглите Strobl C. An Introduction to Recursive Partitioning: Rationale, Application and Characteristics of Classification and Regression Trees, Bagging and Random Forests. Должны быть статьи в свободном доступе, если не найдете, вышлю Вам. Также этот метод реализовал в скрипте R, если работаете в R (или пока не работаете) - тоже могу выслать.

Для статьи в хорошем журнале или для кандидатской (докторской) хорошего уровня, расчет лишь OR c р-values - явно недостаточно. Любая хорошая конференция, посвященная проблемам биоинформатики и определения SNP's, включает вышеприведенные методы.

Если данных достаточное количество, то стоило бы попробовать построить Bayesian networks http://www.bnlearn.com/

Благодарю, коллеги, за ответы. Буду разбираться...

Добрый день, коллеги!
Изучение вопроса привело меня к необходимости использовать random forests, а именно - модификацию, доступную в пакете party для R (cforest), где для построения леса используются "безошибочные" деревья, т.к. в этом случае адекватно оценивается важность переменных и нет преференций для количественных данных (как в "стандартном" random forests)[Strobl, 2007]. Кроме того, необходимо вычислять "кондиционную" значимость, которая адекватнее в случае коррелированных переменных [Strobl, 2009]. И тут возникает следующее препятствие: missing values, которые в моих данных для некоторых параметров составляют более 50%. Однако, и на этот счет у той же "тёти" Strobl нашлось решение [Hapfelmeier, 2013].
На данном этапе моего "мастерства" в R хватает только на импорт таблицы, вычисления значимости переменных и построения гистограммы значимости:

read.table("gen2.txt", h=T)->gen2
set.seed(100)
cfresults<-cforest(cr~.,data=gen2,controls=cforest_unbiased(mtry=10,ntree=100))
varimp(cfresults)->vi1
par(las=2)
barplot(sort(vi1),main="Variable importance",space=0.5,cex.names=0.6)

Прошу помощи у бывалых:
1) Каким должен быть код для вычисления кондиционной важности переменных при наличии пропусков данных? В приложении к [Hapfelmeier, 2013] есть код, но там другие задачи исследования и я не могу разобрать что к чему.
2) Как нарисовать индивидуальные деревья в cforest и правильно интерпретировать их?

Может быть, кто - то соизволит поковыряться с моими данными, на этот случай прилагаю таблицу.
Также прилагается Hapfelmeier, 2013

Попробуйте прочитать о важности переменных здесь http://gewissta.livejournal.com/ / Метод случайного леса.../
Код для вычисления обычной важности переменных там точно есть.
Я немного работал с платформой Salford Predictive Modeler
http://lib.giveawayoftheday.com/b/predictive_analytics/
Стандартные вычисления в среде этой платформы очень легко выполнять, просто используя интерфейс.
Но и с кодами также можно работать.

Коды, генерируемые программой, можно преобразовать в коды R или SAS и в другие форматы.

anserovtv
Большое Вам спасибо за первую ссылочку (до второй пока не добрался, рабочая сеть не позволяет)! Наглядно и доступно представлен алгоритм "классического" random forests.
Однако, автор верно заметил, что классический вариант RF "не переваривает" анализ смешанных количественных и категориальных данных, в этом случае необходим cforest.
Коллеги, если кто-то натыкался на пример использования cforest, будьте добры, поделитесь!

C кодом cforest ответил Вам в личку.
Не совсем понятно, зачем Вам рассматривать каждое дерево по отдельности и что-то интерпретировать, ведь каждое дерево имеет свою индивидуальную структуру и свой прогноз, а конечный результат сводится к голосованию за результат, который дает большинство деревьев. Можно вывести структуру каждого из 100 (как в Вашем случае) деревьев, но потом все это интерпретировать будет заданием очень непростым и вряд ли полезным.

После того, как cforest выдает Вам наиболее значимые предикторы на уровне СНИПов, Вам интересно ведь посмотреть какой генотип ассоциируется с риском развития заболевания или имеет наоборот протективный эффект. Для этого очень удобно использовать комбинацию методов: случайный лес, MDR, логистическая регрессия. Случайный лес - в качестве фильтра наиболее важных СНИПов, MDR - оценка взаимодействия СНИПов и оценка каждого из генотипов, а также их комбинаций, нахождения типа связи между СНИПами. Лог. регрессия - финальный этап, создание модели, которая включает лишь факторы, которые определились как важные из первых двух методов, валидация результатов.

P.S. В Вашей базе - что прогнозируется? nat или cr? Если cr - получается, что исследуете очень редкое заболевание?

Прогнозируется cr (некое заболевание).
Это урезанная база, здесь cr 24 из 165, в оригинальной базе 55 из 439. Но в оригинальной базе 45 предикторов, включая 13 снипов, и громадное количество пропусков в некоторых из них (более 50%). Поэтому пришлось несколько сузить масштабы анализа, убрав часть предикторов. Возможно, это неверно. Поправьте меня, если я не прав.

Мне сложно Вас поправить, поскольку я практически не знаю дизайна исследования, целей, задач и т.д. - если Вы решили урезать массив, значит на то были причины.

Посмотрел немного Вашу базу, прикольная (в том смысле, что собрано столько материала!!). Однако, например, первые три СНИПа IL1b_rs1143634, IL2_rs2069762, IL4_rs2070874, а точнее распределение их генотипов в группах больных и здоровых (cr) не позволяет их нормально анализировать, поскольку наблюдается или perfect separation (то есть классы больных и здоровых идеально расходятся и большинство алгоритмов сходит с ума поэтому) или алгоритм просто не сходится (did not converge), поскольку присутствуют не все комбинации генотипов в базе данных (может быть потому что Вы ее сократили) - это беда всех параметрических методов, в том числе и логит-регрессии(curse of dimensionality).

Что-то нужно думать с базой..

Представьте следующую картину. У нас есть большая база, порядка полутора тысяч чел., в которой собраны наблюдения по разным параметрам, но с разной степенью "перекрывания" данных между ними. Нам известны:
1) генотипы по нескольким снипам для 400 чел.,
2) параметр Х1 для 200 чел, из которых со снипами "перекрывается" 150,
3) параметр Х2 для 300 чел, из которых со снипами перекрывается 250, а с Х1 - только 15
4) параметр Х3 для 100 чел, из которых все перекрываются со снипами, но с Х2 и Х1 только 10-20.
5) и таких параметров еще десяток, с разной степенью "перекрывания" с другими.
6) для каждого из пациентов известен статус по болезни.
Чтобы лучше понять, о чем идет речь, приложу файл с большой первоначальной базой (don_base111.txt).

Причина урезания базы заключается в попытке отобрать тех пациентов, для которых известно как можно больше параметров. Делается это из предположения, что для адекватного анализа RF пропущенные значения конечно могут быть, но должны быть сведены к какому то разумному пределу (не более 30 %).

То есть, до урезания базу можно представить так (база1.png), после - как-то так (база2.png).

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

Пока гляну сегодня по урезанной базе, что с остальными снипами..

Да, я пытаюсь объединить в одном анализе сразу несколько показателей из разных баз.

Я пробовал с помощью SNPStats (http://bioinfo.iconcologia.net/SNPstats) анализировать только снипы, по одиночке. Алгоритм заключается в анализе Хиквадратом сводной таблицы (если без кофакторов) или лог-регрессией (с кофакторами). Там получилась ассоциация со снипами в TNFa и IL4. Если с поправкой Бонферрони, то только TNFa.


А если имеет место взаимодействие снипов с параметрами? То есть, сам по себе снип и среди других снипов не информативен, а среди тех, у кого повышен какой то параметр, он играет роль?

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

Классификационная способность на выборке для обучения: 78%
Классификационная способность на выборке для тестирования: 63,8% (ну уже что-то )

Наилучшая модель - 3-х локусная: IL6_G(-174)C,IL8_rs4073,IL12_rs3212227 кросс-валидация 10 из 10

Правила для комбинаций генотипов:
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0.
IF IL6_G(-174)C = 0 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1.

0 - AA
1 - Aa
2 - aa
Classify as 0 - Control
Classify as 1 - Case

IL6_G(-174)C,IL8_rs4073,IL12_rs3212227

Cross-validation Statistics:

Training Balanced accuracy: 0.789
Training Accuracy: 0.6788
Training Sensitivity: 0.9444
Training Specificity: 0.6336
Training Odds Ratio: 29.3935 (4.5023,191.8959)
Training Χ²: 24.9117 (p < 0.0001)
Training Precision: 0.3049
Training Kappa: 0.3091
Training F-Measure: 0.461

Testing Balanced accuracy: 0.6348
Testing Accuracy: 0.6121
Testing Sensitivity: 0.6667
Testing Specificity: 0.6028
Testing Odds Ratio: 3.0357 (0.169,54.5173)
Testing Χ²: 0.6057 (p = 0.4364)
Testing Precision: 0.2222
Testing Kappa: 0.1473
Testing F-Measure: 0.3333
Cross-validation Consistency: 10/10

Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла
Нажмите для просмотра прикрепленного файла

1) Воспользоваться любым методом импутации данных

А для отбора значимых предикторов можно использовать library(Boruta), на этом форуме уже есть пример её использования. Метод ищет значимые, а не просто имеющие большое значение с точки зрения рандомфореста предикторы.

Доброго времени суток!
Возникло ещё несколько вопросов.


1) Скажите, если исключить из этого перечня случайный лес, что мы потеряем? Ведь МДР тоже может выделить "важные" снипы?

2) Как Вы считаете, не будет ли ошибкой (с точки зрения множественных сравнений), если пытаться с помощью MDRa оценить связи SNP с каждым из анализируемых параметров, чтобы в итоге сказать "с цветом глаз ассоциированы вот эти снипы, с цветом волос - вот эти, а с длиной ресниц - ничего" ?

Благодарю за ответ!
Объясню ситуацию:
1) Доступное количество наблюдений по снипам и статусу case/control = 430
2) Доступное количество наблюдений по параметрам х1, х2, х3 = 250, 200, 180
То есть, если использовать методы импутации, получится что половина (в лучшем случае) данных окажется импутированной.
Позвольте вопрос: правомерно ли брать в анализ переменную, в которой половина данных импутирована?

В литературе по MDR речь идет ещё о вычислении permutation accuracy [Ritchie, Motsinger, 2005], однако я что то не вижу в программе (MDR v. 3.0.2) ничего похожего... Не знаете, куда оно делось?

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

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

2) Если модель допускает некоторые инварианты импутации данных и если Вы делаете импутацию множество раз, то вместо одной выборки получается набор выборок в котором представлено распределение (по которому можно вычислить распределение) интересующей статистики в зависимости от а) имеющихся данных и б) неопределенности предположения о значении отстутствующих данных. Этакий "импутационный бутстреп"

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

Вот тут приатачен файл с примером анализа пакетом Boruta

http://forum.disser.ru/index.php?showtopic...amp;#entry16317