Я пытаюсь смоделировать данные с помощью начальной загрузки, чтобы создать доверительные интервалы для моих реальных данных с помощью графика воронки. Я основываюсь на стратегии принятого ответа на предыдущий вопрос. Вместо использования одного распределения вероятностей для моделирования моих данных я хочу изменить его, чтобы использовать разные распределения вероятностей в зависимости от части моделируемых данных.
Я очень признателен всем, кто может помочь ответить на вопрос или помочь мне сформулировать вопрос более четко.
Моя проблема заключается в написании соответствующего кода R для более сложной формы моделирования данных.
Текущий код:
n <- 1e4
set.seed(42)
sims <- sapply(1:80,
function(k)
rowSums(
replicate(k, sample((1:7)/10, n, TRUE, ps))) / k)
Этот код имитирует данные, где каждая точка данных имеет значение, которое является средним значением между 1:80 наблюдениями. Например, когда значения точек данных являются средним значением 10 наблюдений (k = 10), он случайным образом выбирает 10 значений (которые могут быть 0,1,0,2,0,3, 0,4, 0,5,0,6 или 0,7) на основе распределения вероятностей. ps, который дает вероятность каждого значения (на основе всего эмпирического распределения).
ps выглядит так:
ps <- prop.table(table((DF$mean_score)[DF$total_number_snps == 1]))
# 0.1 0.2 0.3 0.4 0.5 0.6 0.7
#0.582089552 0.194029851 0.124378109 0.059701493 0.029850746 0.004975124 0.004975124
например, вероятность того, что значение наблюдения 0.1, равна 0.582089552.
Теперь вместо использования одного частотного распределения для всех симуляций я хотел бы использовать разные частотные распределения, условно в зависимости от количества наблюдений, лежащих в основе каждой точки данных.
Я сделал таблицу cond_probs, в которой есть строка для каждой из моих реальных точек данных. Есть столбец с total количеством наблюдений и столбец, в котором указана частота каждого из значений для каждого наблюдения.
Пример таблицы cond_probs:
gene_name 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 total
A1 0.664 0.319 0.018 0.000 0.000 0.000 0.000 0.000 0.000 113.000
A2 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000
Таким образом, для точки данных A2 существует только 1 наблюдение, которое имеет значение 0.1. Следовательно, частота 0.1 наблюдений 1. Для A1 имеется 113 наблюдений, и большинство из них (0.664) имеют значение 0.1. Идея состоит в том, что cond_probs похож на ps, но cond_probs имеет распределение вероятностей для каждой точки данных, а не одно для всех данных.
Я хотел бы изменить приведенный выше код, чтобы в выборке использовалось cond_probs вместо ps для частотного распределения. И использовать количество наблюдений k в качестве критерия при выборе строки в cond_probs для выборки. Итак, это будет работать так:
Для точек данных с k количеством наблюдений:
посмотрите в cond_probs таблицу и случайным образом выберите строку, в которой total количество наблюдений по размеру совпадает с k: 0.9k-1.1k. Если таких строк нет, продолжайте.
После выбора точки данных используйте распределение вероятностей из этой строки в cond_probs так же, как ps используется в исходном коде, для случайной выборки k количества наблюдений и вывода среднего значения этих наблюдений.
Для каждой из n итераций replicate произвольно выбирайте с заменой новую точку данных из cond_probs из всех строк, где значение total похоже на текущее значение k (0.9k-1.1k).
Идея состоит в том, что для этого набора данных необходимо определить, какое распределение вероятностей использовать на основе количества наблюдений, лежащих в основе точки данных. Это связано с тем, что в этом наборе данных вероятность наблюдения зависит от количества наблюдений (гены с большим количеством SNP, как правило, имеют более низкий балл за наблюдение из-за генетической связи и фонового отбора).
ОБНОВЛЕНИЕ С ИСПОЛЬЗОВАНИЕМ ОТВЕТА НИЖЕ:
Я попытался использовать ответ ниже, и он работает для смоделированных данных cond_probs в примере, но не для моего реального файла cond_probs. Я импортировал и преобразовал свой файл cond_probs в матрицу с
cond_probs <- read.table("cond_probs.txt", header = TRUE, check.names = FALSE)
cond_probs <- as.matrix(cond_probs)
а в первом примере десять строк (из ~ 20000 строк) выглядят так:
>cond_probs
total 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
[1,] 109 0.404 0.174 0.064 0.183 0.165 0.009 0.000 0.000 0.000 0.000
[2,] 181 0.564 0.221 0.144 0.066 0.006 0.000 0.000 0.000 0.000 0.000
[3,] 289 0.388 0.166 0.118 0.114 0.090 0.093 0.028 0.003 0.000 0.000
[4,] 388 0.601 0.214 0.139 0.039 0.008 0.000 0.000 0.000 0.000 0.000
[5,] 133 0.541 0.331 0.113 0.000 0.008 0.008 0.000 0.000 0.000 0.000
[6,] 221 0.525 0.376 0.068 0.032 0.000 0.000 0.000 0.000 0.000 0.000
[7,] 147 0.517 0.190 0.150 0.054 0.034 0.048 0.007 0.000 0.000 0.000
[8,] 107 0.458 0.196 0.252 0.084 0.009 0.000 0.000 0.000 0.000 0.000
[9,] 13 0.846 0.154 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Если я бегу:
sampleSize <- 20
set.seed(42)
#replace 1:80 with 1: max number of SNPs in gene in dataset
sims_test <- sapply( 1:50, simulateData, sampleSize )
и посмотрите на средние значения выборки с x количеством наблюдений, я получаю только один результат, когда их должно быть 20.
Например:
> sims_test[[31]]
[1] 0.1
И sims_test заказывается не так, как sims:
>sims_test
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.1 0.1 0.1666667 0.200 0.14 0.2666667 0.2000000 0.2375 0.1888889
[2,] 0.1 0.1 0.1333333 0.200 0.14 0.2333333 0.1571429 0.2625 0.1222222
[3,] 0.1 0.1 0.3333333 0.225 0.14 0.1833333 0.2285714 0.2125 0.1555556
[4,] 0.1 0.1 0.2666667 0.250 0.10 0.1500000 0.2000000 0.2625 0.2777778
[5,] 0.1 0.1 0.3000000 0.200 0.16 0.2000000 0.2428571 0.1750 0.1000000
[6,] 0.1 0.1 0.3666667 0.250 0.16 0.1666667 0.2142857 0.2500 0.2000000
[7,] 0.1 0.1 0.4000000 0.300 0.12 0.2166667 0.1857143 0.2375 0.1666667
[8,] 0.1 0.1 0.4000000 0.250 0.10 0.2500000 0.2714286 0.2375 0.2888889
[9,] 0.1 0.1 0.1333333 0.300 0.14 0.1666667 0.1714286 0.2750 0.2888889
ОБНОВЛЕНИЕ 2
Используя cond_probs <- head(cond_probs,n), я определил, что код работает до n = 517, а затем для всех размеров больше этого он дает тот же результат, что и выше. Я не уверен, проблема ли это в самом файле или проблема с памятью. Я обнаружил, что если я удалю строку 518 и продублирую строки до этого несколько раз, чтобы сделать файл большего размера, это сработает, предполагая, что сама строка вызывает проблему. Строка 518 выглядит так:
9.000 0.889 0.000 0.000 0.000 0.111 0.000 0.000 0.000 0.000 0.000
Я нашел еще 4 оскорбительные строки:
9.000 0.444 0.333 0.111 0.111 0.000 0.000 0.000 0.000 0.000 0.000
9.000 0.444 0.333 0.111 0.111 0.000 0.000 0.000 0.000 0.000 0.000
9.000 0.111 0.222 0.222 0.111 0.111 0.222 0.000 0.000 0.000 0.000
9.000 0.667 0.111 0.000 0.000 0.000 0.222 0.000 0.000 0.000 0.000
Я не замечаю в них ничего необычного. Всего у них 9 сайтов. Если я удалю эти строки и запустил файл cond_probs, содержащий только строки ПЕРЕД этим, код заработает. Но должны быть и другие проблемные строки, так как вся команда cond_probs по-прежнему не работает.
Я попытался поместить эти проблемные строки обратно в меньший файл cond_probs, и этот файл затем заработал, поэтому я очень смущен, поскольку не кажется, что строки изначально проблематичны. С другой стороны, тот факт, что все они имеют 9 сайтов, предполагает некую причинную закономерность.
Я был бы счастлив поделиться всем файлом в частном порядке, если это поможет, поскольку я не знаю, что делать дальше для устранения неполадок.
Возникает еще одна проблема: я не уверен, работает ли код должным образом. Я сделал фиктивный файл cond_probs, в котором есть две точки данных с «общим» наблюдением «1»:
total 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000
1.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Поэтому я ожидал бы, что они оба будут отобраны для точек данных с наблюдением «1» и, следовательно, получат примерно 50% наблюдений со средним значением «0,2» и 50% со средним значением «0,6». Однако среднее значение всегда равно 0,2:
sims_test[[1]]
[1] 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
Даже если я сделаю выборку 10000 раз, все наблюдения будут 0,2, а не 0,6. Насколько я понимаю, код должен случайным образом выбирать новую строку из cond_probs с аналогичным размером для каждого наблюдения, но в этом случае, похоже, этого не происходит. Я неправильно понимаю код или все еще проблема с тем, что я ввел неверный код?
Весь файл cond_probs можно найти по следующему адресу:
ОБНОВЛЕНИЕ 3
Изменение sapply на lapply при запуске моделирования устранило эту проблему.
Еще одна причина, по которой я думаю, что оставить cond_probs как есть и выбрать распределение sampleSize количество раз может быть лучшим решением: вероятность выбора распределения должна быть связана с его частотой в cond_probs. Если мы объединим распределения, шансы выбора распределения с total 9 или 10 больше не будут зависеть от количества наблюдений с этими итогами. Пример: Если есть 90 дистрибутивов с total=10 и 10 с total=9, должна быть 90% возможность выбрать дистрибутив с total=10. Если мы объединим распределения, не станут ли шансы 50/50 для выбора распределения с «общим» = 9 или 10 (что не было бы идеальным)?
