Как переставить сеть в igraph для R? - PullRequest
Новая
       20

Как переставить сеть в igraph для R?

1 голос
/ 04 апреля 2019

Я пытаюсь написать код для процедуры Монте-Карло в R. Моя цель - оценить значимость метрики, рассчитанной для взвешенной, однонаправленной, неориентированной сети, отформатированной для пакета igraph.

Пока я включил в код следующие шаги:

1. Создайте взвешенную, единственную, неориентированную сеть и рассчитайте наблюдаемую модульность Лувена 

nodes <- read.delim("nodes.txt")
links <- read.delim("links.txt")
anurosnet <- graph_from_data_frame(d=links, vertices=nodes, directed=F) 
anurosnet

modularity1 = cluster_louvain(anurosnet)
modularity1$modularity #observed value 

obs=modularity1$modularity
obs
real<-data.frame(obs)
real

2. Создать пустой вектор 

Nperm = 9 #I am starting with a low n, but intend to use at least 1000 permutations
randomized.modularity=matrix(nrow=length(obs),ncol=Nperm+1)
row.names(randomized.modularity)=names(obs)
randomized.modularity[,1]=obs 
randomized.modularity

3. Перестановка исходной сети с сохранением ее характеристик, расчет модульности Лувена для всех рандомизированных сетей и компиляция результатов в векторе 

i<-1
while(i<=Nperm){ 

  randomnet <- rewire(anurosnet, with=each_edge(0.5)) #rewire vertices with constant probability
  E(randomnet)$weight <- sample(E(anurosnet)$weight) #shuffle initial weights and assign them randomly to edges

  mod<-(cluster_louvain(randomnet))

  mod$modularity

  linha = mod$modularity

  randomized.modularity[,i+1]=linha
  print(i)
  i=i+1
}
randomized.modularity #Here the result is not as expected

4. График наблюдаемого значения против распределения рандомизированных значений 

niveis<-row.names(randomized.modularity)
for(k in niveis)
{
  if(any(is.na(randomized.modularity[k,]) == TRUE))
  {
    print(c(k, "metrica tem NA"))
  } else {
    nome.arq<- paste("modularity",k,".png", sep="")
    png(filename= nome.arq, res= 300, height= 15, width=21, unit="cm")
    plot(density(randomized.modularity[k,]), main="Observed vs. randomized",)
    abline(v=obs[k], col="red", lwd=2, xlab="")
    dev.off()
    print(k)
    nome.arq<- paste("Patefield_Null_mean_sd_",k,".txt", sep="")
    write.table(cbind(mean(randomized.modularity[k,]),sd(randomized.modularity[k,])), file=paste(nome.arq,sep=""), 
                sep=" ",row.names=TRUE,col.names=FALSE)
  }
}

5. Оценить P-значение (значимость) 

significance=matrix(nrow=nrow(randomized.modularity),ncol=3)
row.names(significance)=row.names(randomized.modularity)
colnames(significance)=c("p (rand <= obs)", "p (rand >= obs)", "p (rand=obs)")

signif.sup=function(x) sum(x>=x[1])/length(x)
signif.inf=function(x) sum(x<=x[1])/length(x)
signif.two=function(x) ifelse(min(x)*2>1,1,min(x)*2)

significance[,1]=apply(randomized.modularity,1,signif.inf)
significance[,2]=apply(randomized.modularity,1,signif.sup)
significance[,3]=apply(significance[,-3],1,signif.two)

significance

Что-то пошло не так в шаге 3. Я ожидал, что вектор будет заполнен 10 значениями, но по некоторым причинам он через некоторое время останавливается.

Слот "mod $ modularity" неожиданно получает 2 значения вместо 1.

Два файла TXT, упомянутые в начале кода, можно скачать здесь:

https://1drv.ms/t/s!AmcVKrxj94WClv8yQyqyl4IWk5mNvQ https://1drv.ms/t/s!AmcVKrxj94WClv8z_Pow5Tg2U7mjLw

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

1 Ответ

1 голос
/ 05 апреля 2019

Ваша ошибка связана с несоответствием размеров с вашей матрицей randomized.modularity и некоторыми из ваших рандомизированных результатов модульности. В вашем примере ваша матрица в итоге становится [1 x Nperm], однако иногда 2 оценки модульности возвращаются при перестановках. Чтобы это исправить, я просто храню результаты в виде списка. Остальную часть вашего анализа необходимо будет скорректировать, поскольку у вас несоответствие показателей модульности. 

library(igraph)

nodes <- read.delim("nodes.txt")
links <- read.delim("links.txt")
anurosnet <- graph_from_data_frame(d=links, vertices=nodes, directed=F) 
anurosnet

modularity1 = cluster_louvain(anurosnet)
modularity1$modularity #observed value 

obs <- modularity1$modularity
obs
real<-data.frame(obs)
real

Nperm = 100 #I am starting with a low n, but intend to use at least 1000 permutations
#randomized.modularity <- matrix(nrow=length(obs),ncol=Nperm+1)
#row.names(randomized.modularity) <- names(obs)
randomized.modularity <- list()
randomized.modularity[1] <- obs 
randomized.modularity


for(i in 1:Nperm){ 

  randomnet <- rewire(anurosnet, with=each_edge(0.5)) #rewire vertices with constant probability
  E(randomnet)$weight <- sample(E(anurosnet)$weight) #shuffle initial weights and assign them randomly to edges

  mod <- (cluster_louvain(randomnet))      
  mod$modularity      
  linha = mod$modularity     
  randomized.modularity <- c(randomized.modularity, list(linha))

}

randomized.modularity

Лучший способ написать цикл 

randomized.modularity <- lapply(seq_len(Nperm), function(x){  
                                            randomnet <- rewire(anurosnet, with=each_edge(0.5)) #rewire vertices with constant probability
                                            E(randomnet)$weight <- sample(E(anurosnet)$weight) #shuffle initial weights and assign them randomly to edges
                                            return(cluster_louvain(randomnet)$modularity)
                                          })
...