Блог

Не могли бы вы мне помочь?

Я пишу код R для автоматизации анализа нулевой модели нескольких сетей. Во-первых, код считывает несколько текстовых matrices файлов в R. Во-вторых, он вычисляет топологическую метрику для каждой сети. В-третьих, он рандомизирует каждую сеть N раз, используя нулевую модель. В-четвертых, он вычисляет одну и ту же топологическую метрику для всех рандомизированных версий исходных матриц.

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

Вот примеры data frames того, что необходимо проанализировать:

 networks <- paste("network", rep(1:3), sep = "")
randomizations <- seq(1:10)

observed.ex <- data.frame(network = networks,
                          observed = runif(3, min = 0, max = 1))

randomized.ex <- data.frame(network = sort(rep(networks, 10)),
                            randomization = rep(randomizations, 3),
                            randomized = rnorm(length(networks)*
                                                   length(randomizations),
                                               mean = 0.5, sd = 0.1))
  

На первом этапе окончательного анализа код оценивает P-значения, выполняя простые подсчеты. Как вы видите, мне нужно сделать копии вызова вычисления для каждой сети:

 randomized.network1 <- subset(randomized.ex, network == "network1")
sum(randomized.network1$randomized >= observed.ex$observed[1]) /
    length(randomized.network1$randomized)
sum(randomized.network1$randomized <= observed.ex$observed[1]) /
    length(randomized.network1$randomized)

randomized.network2 <- subset(randomized.ex, network == "network2")
sum(randomized.network2$randomized >= observed.ex$observed[2]) /
    length(randomized.network2$randomized)
sum(randomized.network2$randomized <= observed.ex$observed[2]) /
    length(randomized.network2$randomized)

randomized.network3 <- subset(randomized.ex, network == "network3")
sum(randomized.network3$randomized >= observed.ex$observed[3]) /
    length(randomized.network3$randomized)
sum(randomized.network3$randomized <= observed.ex$observed[3]) /
    length(randomized.network3$randomized)
  

На втором этапе окончательного анализа код создает графики плотности. Как вы видите, мне нужно сделать копии вызова вертикальной линии для каждой сети:

 ggplot(randomized.ex, aes(randomized))  
    geom_density()  
    facet_grid(network~.)  
    geom_vline(data=filter(randomized.ex, network == "network1"),
               aes(xintercept = observed.ex$observed[1]), colour = "red")   
    geom_vline(data=filter(randomized.ex, network == "network2"),
               aes(xintercept = observed.ex$observed[2]), colour = "red")   
    geom_vline(data=filter(randomized.ex, network == "network3"),
               aes(xintercept = observed.ex$observed[3]), colour = "red") 
  

Есть ли способ сделать этот окончательный анализ более общим, чтобы он всегда выполнял одни и те же вычисления и графики, независимо от того, сколько сетей считывается в начале?

Большое вам спасибо!

Похоже, это можно аккуратно обернуть в lapply цикл, который повторяется по каждому файлу. Как для вас работает приведенное ниже? Вы также можете передавать имена файлов, а не количество файлов (в настоящее время 1: 3), и первая строка «читается» в ваших текстовых матрицах.

 library(dplyr) #For %>%, group_by, and summarize
output <- lapply(1:3, function(network_num){
  network <- paste0("network", network_num)
  n_randomizations <- 10
  observed.ex <- runif(1)
  randomized.ex <- rnorm(n_randomizations, mean = 0.5, sd = 0.1)

  return(data.frame(network=network, observed=observed.ex, randomized=randomized.ex))
}) %>% do.call(what = rbind)

output %>%
  group_by(network) %>%
  summarize(p_value=mean(observed>=randomized))

ggplot(output)  
  geom_density(aes(randomized))  
  facet_grid(network~.)  
  geom_vline(aes(xintercept = observed), col="red")