Я пытаюсь создать фигуру, которая показывает связь между вирусными копиями и освещением генома (GCC). Вот как выглядят мои данные:
Сначала я только построил линейную регрессию, но мои руководители сказали мне, что это неправильно, и попробовал сигмоидальную кривую. Поэтому я сделал это с помощью geom_smooth:
library(scales)
ggplot(scatter_plot_new, aes(x = Copies_per_uL, y = Genome_cov, colour = Virus)) +
geom_point() +
scale_x_continuous(trans = log10_trans(), breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x))) +
geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, size = 1) +
theme_bw() +
theme(legend.position = 'top', legend.text = element_text(size = 10), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size=12), axis.title.y = element_text(margin = margin (r = 10)), axis.title.x = element_text(margin = margin(t = 10))) +
labs(x = "Virus copies/µL", y = "GCC (%)") +
scale_y_continuous(breaks=c(25,50,75,100))
Тем не менее, мои руководители говорят, что это тоже неправильно, потому что из-за кривых похоже, что GCC может превысить 100%, а это невозможно.
Мой вопрос: как лучше показать связь между вирусными копиями и GCC? Я хочу прояснить, что А) низкий уровень копий вируса = низкий GCC, и что B) после того, как определенное количество вируса копирует плато GCC.
Я исследовал множество различных методов - GAM, LOESS, логистику, кусочно - но я не знаю, как определить, какой метод лучше всего подходит для моих данных.
РЕДАКТИРОВАТЬ: это данные:
>print(scatter_plot_new)
Subsample Virus Genome_cov Copies_per_uL
1 S1.1_RRAV RRAV 100 92500
2 S1.2_RRAV RRAV 100 95900
3 S1.3_RRAV RRAV 100 92900
4 S2.1_RRAV RRAV 100 4049.54
5 S2.2_RRAV RRAV 96.9935 3809
6 S2.3_RRAV RRAV 94.5054 3695.06
7 S3.1_RRAV RRAV 3.7235 86.37
8 S3.2_RRAV RRAV 11.8186 84.2
9 S3.3_RRAV RRAV 11.0929 95.2
10 S4.1_RRAV RRAV 0 2.12
11 S4.2_RRAV RRAV 5.0799 2.71
12 S4.3_RRAV RRAV 0 2.39
13 S5.1_RRAV RRAV 4.9503 0.16
14 S5.2_RRAV RRAV 0 0.08
15 S5.3_RRAV RRAV 4.4147 0.08
16 S1.1_UMAV UMAV 5.7666 1.38
17 S1.2_UMAV UMAV 26.0379 1.72
18 S1.3_UMAV UMAV 7.4128 2.52
19 S2.1_UMAV UMAV 21.172 31.06
20 S2.2_UMAV UMAV 16.1663 29.87
21 S2.3_UMAV UMAV 9.121 32.82
22 S3.1_UMAV UMAV 92.903 627.24
23 S3.2_UMAV UMAV 83.0314 615.36
24 S3.3_UMAV UMAV 90.3458 632.67
25 S4.1_UMAV UMAV 98.6696 11180
26 S4.2_UMAV UMAV 98.8405 12720
27 S4.3_UMAV UMAV 98.7939 8680
28 S5.1_UMAV UMAV 98.6489 318200
29 S5.2_UMAV UMAV 99.1303 346100
30 S5.3_UMAV UMAV 98.8767 345100
источник
method.args=list(family=quasibinomial))
аргументыgeom_smooth()
в исходный код ggplot.se=FALSE
. Всегда приятно показать людям, насколько велика эта неопределенность ...Ответы:
Еще один способ сделать это - использовать байесовскую формулировку. Начать с нее может быть немного тяжело, но это значительно упрощает формулировку специфики вашей проблемы, а также позволяет лучше понять, где находится «неопределенность». является
Stan - это сэмплер Монте-Карло с относительно простым в использовании программным интерфейсом, библиотеки доступны для R и других, но я использую Python здесь
мы используем сигмовидную кишку, как и все: у нее есть биохимические мотивы, а с ней очень удобно работать с математической точки зрения. Хорошая параметризация для этой задачи:
где
alpha
определяет среднюю точку сигмовидной кривой (т.е. где она пересекает 50%) иbeta
определяет наклон, значения ближе к нулю более плоскиечтобы показать, как это выглядит, мы можем получить ваши данные и нанести их на график:
где
raw_data.txt
содержит данные, которые вы дали, и я преобразовал освещение в нечто более полезное. коэффициенты 5.5 и 3 выглядят хорошо и дают график очень похожий на другие ответы:чтобы «подогнать» эту функцию с помощью Stan, нам нужно определить нашу модель, используя ее собственный язык, который представляет собой смесь между R и C ++. простая модель будет что-то вроде:
который, надеюсь, читает хорошо. у нас есть
data
блок, который определяет данные, которые мы ожидаем при выборке модели,parameters
определяем вещи, которые выбираем, иmodel
определяем функцию правдоподобия. Вы говорите Стэну «скомпилировать» модель, что занимает некоторое время, а затем вы можете взять из нее некоторые данные. например:arviz
упрощает создание хороших диагностических графиков, в то время как печать подгонки дает вам хорошую сводку параметров в стиле R:большое стандартное отклонение
beta
говорит о том, что данные действительно не предоставляют много информации об этом параметре. также некоторые ответы, дающие более 10 значащих цифр в подборах моделей, несколько преувеличиваютПоскольку в некоторых ответах отмечалось, что для каждого вируса могут потребоваться свои параметры, я расширил модель, чтобы разрешить ее
alpha
иbeta
изменить на «Вирус». все становится немного странно, но два вируса почти наверняка имеют разныеalpha
значения (т. е. вам нужно больше копий / мкл RRAV для того же покрытия), и график, показывающий это:данные такие же, как и раньше, но я нарисовал кривую для 40 образцов задних.
UMAV
кажется относительно хорошо определенным, хотяRRAV
может следовать за тем же уклоном и нуждаться в большем количестве копий, или иметь более крутой уклон и похожее количество копий. большая часть задней массы нуждается в большем количестве копий, но эта неопределенность может объяснить некоторые различия в других ответах, которые находят разные вещиЯ в основном используется ответив это упражнение , чтобы улучшить свои знания Стана, и я поставил Jupyter тетрадь это здесь в случае , если кто заинтересован / хочет повторить это.
источник
(Отредактировано с учетом комментариев ниже. Спасибо @BenBolker & @WeiwenNg за полезный вклад.)
Подгонка дробной логистической регрессии к данным. Он хорошо подходит для процентных данных, которые ограничены от 0 до 100% и теоретически обоснованы во многих областях биологии.
Обратите внимание, что вам может потребоваться разделить все значения на 100, чтобы соответствовать ему, поскольку программы часто ожидают, что данные будут находиться в диапазоне от 0 до 1. И, как рекомендует Бен Болкер, для решения возможных проблем, вызванных строгими предположениями биномиального распределения относительно дисперсии, используйте квазибиномиальное распределение вместо.
Я сделал несколько предположений, основанных на вашем коде, например, что вас интересуют 2 вируса, и они могут демонстрировать разные паттерны (то есть может быть взаимодействие между типом вируса и количеством копий).
Во-первых, модель подойдет:
Если вы доверяете p-значениям, вывод не предполагает, что два вируса отличаются друг от друга. Это противоречит приведенным ниже результатам @ NickCox, хотя мы использовали разные методы. Я не был бы очень уверен в любом случае с 30 точками данных.
Во-вторых, сюжет:
Нетрудно составить код для визуализации вывода самостоятельно, но, похоже, есть пакет ggPredict, который сделает большую часть работы за вас (не могу ручаться за это, я сам не пробовал). Код будет выглядеть примерно так:Обновление: я больше не рекомендую код или функцию ggPredict в более общем смысле. Попробовав это, я обнаружил, что нанесенные точки не совсем отражают входные данные, а вместо этого изменены по какой-то странной причине (некоторые из нанесенных точек были выше 1 и ниже 0). Поэтому я рекомендую кодировать это самостоятельно, хотя это больше работы.
источник
family=quasibinomial()
чтобы избежать предупреждения (и основных проблем со слишком строгими предположениями об отклонениях). Воспользуйтесь советом @ mkt по другой проблеме.Это не отличный ответ от @mkt, но, в частности, графики не поместятся в комментарии. Сначала я подгоняю логистическую кривую в Stata (после регистрации предиктора) ко всем данным и получаю этот график
Уравнение
100
invlogit
(-4,192654 + 1,880951log10
(Copies
))Теперь я подгоняю кривые отдельно для каждого вируса в простейшем сценарии, когда вирус определяет переменную индикатора. Вот для записи сценарий Stata:
Это сильно влияет на крошечный набор данных, но P-значение для вируса выглядит благоприятным для согласования двух кривых совместно.
источник
Попробуйте сигмовидную функцию. Есть много формулировок этой формы, включая логистическую кривую. Гиперболическая касательная является еще одним популярным выбором.
Учитывая графики, я не могу исключить и простую пошаговую функцию. Боюсь, вы не сможете различить ступенчатую функцию и любое количество сигмовидных характеристик. У вас нет никаких наблюдений, когда ваш процент находится в диапазоне 50%, поэтому простая пошаговая формулировка может быть самым экономным выбором, который работает не хуже, чем более сложные модели.
источник
Вот 4PL (логистика с четырьмя параметрами), как ограниченная, так и неограниченная, с уравнением по К. А. Гольштейну, М. Гриффину, Дж. Хонгу, П. Д. Сэмпсону, «Статистический метод определения и сравнения пределов обнаружения биоанализов», Anal , Химреагент 87 (2015) 9795-9801. Уравнение 4PL показано на обоих рисунках, а значения параметров следующие: a = нижняя асимптота, b = коэффициент наклона, c = точка перегиба и d = верхняя асимптота.
На рисунке 1 a равно 0%, а d равно 100%:
Рисунок 2 не имеет ограничений на 4 параметра в уравнении 4PL:
Это было весело, я не претендую на то, что знаю что-то биологическое, и будет интересно посмотреть, как все это устроится!
источник
Я извлек данные из вашей диаграммы рассеяния, и мой поиск по уравнению нашел трехпараметрическое уравнение логистического типа в качестве хорошего кандидата: "y = a / (1.0 + b * exp (-1.0 * c * x))", где " х "- логарифмическая основа 10 на ваш участок. Подогнанные параметры были a = 9.0005947126706630E + 01, b = 1.2831794858584102E + 07 и c = 6.6483431489473155E + 00 для моих извлеченных данных, подгонка (log 10 x) исходных данных должна дать аналогичные результаты, если вы повторно подгоните исходные данные, используя мои значения в качестве начальных оценок параметров. Мои значения параметров дают R-квадрат = 0,983 и RMSE = 5,625 на извлеченных данных.
РЕДАКТИРОВАТЬ: Теперь, когда вопрос был отредактирован для включения фактических данных, вот график с использованием вышеуказанного уравнения с 3 параметрами и начальных оценок параметров.
источник
Так как мне пришлось открыть свой большой рот о Хевисайде, вот результаты. Я установил точку перехода на log10 (viruscopies) = 2.5. Затем я вычислил стандартные отклонения двух половин набора данных, то есть Хевисайд предполагает, что данные с обеих сторон имеют все производные = 0.
Стандартное отклонение от
правой стороны = 4,76 стандартное отклонение от правой стороны = 7,72
Так как оказывается, что в каждой партии 15 выборок, среднее значение стандартного отклонения составляет 6,24.
Если предположить, что «RMSE», приведенная в других ответах, является «среднеквадратичной ошибкой», то функция Хевисайда, по-видимому, выполняет, по крайней мере, даже лучше, чем большинство «Z-кривой» (заимствовано из номенклатуры фотографического ответа). Вот.
редактировать
Бесполезный график, но запрошенный в комментариях:
источник