1

Preciso comparar meu valor vindo de um bootstrap de matriz, com o valor único original, mas não sei exatamente como fazer isso com meus dados. Aqui eu tenho os dados de exemplo

###DADOSEXEMPLO##

  valores <- c(1, 5, 3, 8, 2, 9, 3, 2, 3 ,2, 5, 3, 6, 5, 1, 5, 3, 4, 2, 2, 5, 6, 5, 7, 1, 2, 
 8, 12, 5, 1, 5, 3, 6, 5, 9, 3,3,4,5,7,8,2,7,8,4,3,1,4,4 )
matrixnet<-matrix(valores , 7, 7)
matcor <- cor(matrixnet, method = "spearman")

  matrixnetwork = graph.adjacency(matcor, mode="undirected", 
                        weighted 
                        = TRUE, 
                        add.colnames=NULL, diag=FALSE)
 plot(matrixnetwork)

Aqui calculo o valor único do gráfico, neste caso para a medida entre

   centr_betw(matrixnetwork, directed = FALSE, normalized = TRUE)$res

Aqui aplico a permutação "alterar valor em vertes" em conjunto com a função específica que desejo aplicar, neste caso Betw

set.seed(123
bootcentr_Betw <- lapply(seq_len(Nperm), function(x){  
randomnet <- rewire(matrixnetwork, with=each_edge(0.5)) #rewire vertices with constant probability
  E(randomnet)$weight <- sample(E(matrixnetwork)$weight) #shuffle initial weights and assign 
 them randomly to edges
 return(centr_betw(randomnet)$res)
 })

Agora preciso comparar o valor original com o bootstrap do gráfico aleatório.

Abaixo tentei traçar o valor único original em relação ao bootstrap, replicando o código que obtive de outro lugar

# 4. Plote o valor observado em rela??o ? distribui??o de valores aleat?rios
 niveis<-row.names(bootstrength )
for(k in niveis)
{
  if(any(is.na(bootstrength [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(bootstrength [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(bootstrength [k,]),sd(bootstrength [k,])), 
file=paste(nome.arq,sep=""), 
                sep=" ",row.names=TRUE,col.names=FALSE)
  }
}

Confesso que não sei onde está o arquivo que ele gerou.

Depois disso tentei replicar as estatísticas que encontrei aqui no fórum, mas não funcionou.

 #5. Estime o valor P (signific?ncia)

 Estou tentando comparar um valor observado e sua permutação e gerar um valor p pra isso.

Eu fiz o teste t sugerido no código em uma resposta abaixo. Porém vi que algumas pessoas utilizam uma fórmula chama de teste de permutação. Que é Número de valores (estatísticas) que são maiores ou iguais ao valor observado e dividem pelo número de valores. No código, pval = soma(s) >= s0)/N; Alguém pode me ajudar a aplicar essa fórmula aos dados acima? Pois é diferente do teste que eu apliquei

1 Resposta 1

2

Primeiro vou refazer os dados, para ter um grafo com vértices de centralidades diferentes.

library(igraph)
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union

###DADOSEXEMPLO##

valores <- c(1, 5, 3, 8, 2, 9, 3, 2, 3 ,2, 5, 3, 6, 5, 1, 5, 3, 4, 2, 2, 5, 6, 5, 7, 1, 2, 
             8, 12, 5, 1, 5, 3, 6, 5, 9, 3,3,4,5,7,8,2,7,8,4,3,1,4,4 )
matrixnet <- matrix(valores, 7, 7)
matcor <- cor(matrixnet, method = "spearman")

Created on 2023-12-08

O grafo gerado por esta matrix de correlações tem todos os vértices com o mesmo grau 6. E a função centr_betw dá tudo zeros. Portanto vou mudar algumas correlações. As mais pequenas passam a ser iguais a zero.

cor_pequena <- 0.1
matcor[] <- ifelse(matcor < cor_pequena, 0, matcor)

Isto é suficiente para ter um grafo mais realista.

matrixnetwork = graph.adjacency(matcor, mode = "undirected", 
                                weighted = TRUE, 
                                add.colnames = NULL, diag = FALSE)
degree(matrixnetwork)
#> [1] 3 4 3 3 1 2 2

plot(matrixnetwork)

Created on 2023-12-08

Agora o bootstrap.

  • Primeiro calcular com centr_betw a estatística para o grafo original;
  • depois simular Mperm grafos e calcular a mesma estatística;
  • Traçar os gráficos das estatísticas bootstrap;
  • testes t para a média dão os p-valores pedidos.

orig <- centr_betw(matrixnetwork, directed = FALSE, normalized = TRUE)$res

Nperm <- 1e3
set.seed(123)
bootcentr_Betw <- replicate(Nperm, {
  randomnet <- rewire(matrixnetwork, with = each_edge(0.5))
  E(randomnet)$weight <- sample(E(matrixnetwork)$weight)
  centr_betw(randomnet)$res
})

op <- par(mfrow = c(2, 4))
for(k in seq_along(orig)) {
  main <- sprintf("Vértice: %d", k)
  bootcentr_Betw[k,] |> density() |> plot(main = main, sub = "Observed vs. randomized")
  abline(v = orig[k], col = "red")
}
par(op)


pval_eq <- numeric(length(orig))
# pval_gt <- numeric(length(orig))
# pval_le <- numeric(length(orig))
for(k in seq_along(orig)) {
  rand <- bootcentr_Betw[k, ]
  obs <- orig[k]
  pval_eq[k] <- t.test(rand, mu = obs)$p.value
  # pval_gt[k] <- t.test(rand, mu = obs, alternative = "greater")$p.value
  # pval_le[k] <- t.test(rand, mu = obs, alternative = "less")$p.value
}
setNames(
  pval_eq,
  # pval_gt,
  # pval_le,
  seq_along(orig)
)
#>             1             2             3             4             5 
#>  2.786037e-57 1.458862e-210  6.921656e-01 6.997004e-225  9.001682e-87 
#>             6             7 
#> 9.445648e-114  9.104672e-38

Created on 2023-12-08


Edição

O seguinte código calcula modularidades para um vetor de resoluções entre 0 e 2. Isso é feito pela função modularity_by_resolution.

modularity_by_resolution <- function(x, from, to, by) {
  resolution_vec <- seq(from, to, by = by)
  mod_vec <- sapply(resolution_vec, \(r) {
    x |>
      cluster_louvain(resolution = r) |>
      membership() |>
      modularity(x, membership = _, resolution = r)
  })
  list(resolution = resolution_vec, modularity = mod_vec)
}

res_0.01 <- modularity_by_resolution(matrixnetwork, 0, 2, by = 0.01)
main <- "Modularity by resolution\nResolution step: 0.01"
plot(res_0.01$resolution, res_0.01$modularity, type = "l", main = main)
abline(h = 0, lty = "dashed")


# com vários valores para o passo da resolução
res_vec <- c(0.01, 0.001)
# duas colunas, uma para cada valor da resolução
old_par <- par(mfrow = c(1L, 2L))
lapply(res_vec, \(r) {
  mod <- modularity_by_resolution(matrixnetwork, 0, 2, by = r)
  main <- paste0("Resolution step: ", r)
  with(mod, {
    plot(resolution, modularity, type = "l", main = main)
    abline(h = 0, lty = "dashed")
  })
})

#> [[1]]
#> NULL
#> 
#> [[2]]
#> NULL
# repor o parâmetro gráfico
par(old_par)

Created on 2023-12-11 with reprex v2.0.2


Edição 2

Eu tentei dessa forma para aplicar a mudança no valor de resolution, a cada calculo de número de grupos, na função cluster_louvain.

group_by_resolution <- function(x, from, to, by) {
  resolution_vec <- seq(from, to, by = by)
  mod_vec <- sapply(resolution_vec, \(r) {
    cluster_louvain(matrixnetwork, resolution = r)
  })
  list(resolution = resolution_vec, groups = mod_vec)
}

res_0.01 <- group_by_resolution(matrixnetwork, 0, 2, by = 0.01)
main <- "group by resolution\nResolution step: 0.01"
plot(res_0.01$resolution, res_0.01$groups, type = "l", main = main)
abline(h = 0, lty = "dashed")
res_vec <- c(0.01, 0.001)

Porém, na hora de gerar o gráfico, diz que meu objeto é duplo e não pode ser coercido. Pelo que entendi ele está trazendo em conjunto tanto o valor de quantos grupos tem, quanto quem pertence a cada grupo, mas eu queria apenas o valor de quantos grupos tem.

Também tentei de outra forma:

    modularity_by_resolution <- function(x, from, to, by) {
      resolution_vec <- seq(from, to, by = by)
      mod_vec <- sapply(resolution_vec, \(r) {
        x |>
          cluster_louvain(resolution = r) |>
          membership() |>
          groups(x, membership = _, resolution = r)
      })
      list(resolution = resolution_vec, groups = mod_vec)
    }
    res_0.01 <- modularity_by_resolution(matrixnetwork, 0, 2, by = 0.01)

mas o seguinte erro é exposto:  

Error in groups(x, membership = membership(cluster_louvain(x, resolution = r)), :
unused arguments (membership = membership(cluster_louvain(x, resolution = r)), resolution = r)


Edição 3

Solução do problema descrito na edição 2.

n_clusters_by_resolution <- function(x, from, to, by) {
  resolution_vec <- seq(from, to, by = by)
  mod_vec <- sapply(resolution_vec, \(r) {
    cluster_louvain(matrixnetwork, resolution = r)
  })
  list(resolution = resolution_vec, nclusters = lengths(mod_vec))
}

clust_0.01 <- n_clusters_by_resolution(matrixnetwork, 0, 2, by = 0.01)
main <- "group by resolution\nResolution step: 0.01"
plot(clust_0.01$resolution, clust_0.01$nclusters, type = "s", main = main)
1
  • A conversa se estendeu, então os comentários foram movidos para o chat.
    – Syzoth
    Commented 11/12/2023 às 19:56

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .