Feature Request: Calculate Nework Stats on Posteriors

[mjbsp]

Thanks for the hardwork on the really nice package.

I think it would be useful to provide a sensible and flexible way to calculate network stats on posterior estimates, so that it's possible to get estimates of variation/uncertainty of these measures. Although it might be possible to include various measures in the package itself, it might be more flexible in the long run to provide a way for people to use their own functions and functions from other packages.

Here is an example calculating assortivity; however, in an ideal world it'd be possible to swap out assortivity with anything that took a weight matrix.

library(BGGM)
library(assortnet)
library(tidyverse)

Y <- BGGM::bfi[,1:10]

# fit model
fit <- estimate(Y = Y, 
                analytic = FALSE, 
                iter = 1000)

# calculate assortivity

## select pcors
pcors <- dplyr::select(fit$posterior_samples, contains("pcors"))

## list of columns belowinging in each group
## e.g., first 5 are "a", last 5 are "c"

membership <- c(rep("a", 5), rep("c", 5))

## To store results
assortr <- list()

## calculate assortivity for each posterior
for(i in 1:nrow(pcors)){
  
  # turn posteriors into matrix
  m <- matrix(unlist(pcors[i,]), nrow = 10, ncol = 10)
  diag(m) <- 0
  
  # calculator assortivity
  assortr[[i]] <- assortment.discrete(m, types = membership, weighted = TRUE, SE = FALSE, M = 1)$r
  
}

## mean assortivity across posterior
assortivity.estimate <- mean(unlist(assortr))

## CI
assortivity.ci <- quantile(unlist(assortr), probs = c(.025, .975))

[donaldRwilliams]

This wont be a problem to implement. Do you have any suggestion for a general name.

Note that the function I add will essentially be like FUN from the apply family. Hence, the user will be entirely responsible for passing the appropriate arguments.

so perhaps

roll_your_own <- function(object, FUN, ...)

where the ... the takes the arguments that you passed to assortment.discrete.

Just a thought

[donaldRwilliams]

See here. Note it will work for any network stat, but gotta be careful with what the function returns.
For example, you could make a custom function to compute whatever you want and provide it to
FUN.

You will also want the developmental version which is still a work in progress but most will not change in any meaningful way. So essentially version 2.0.0

if (!requireNamespace("remotes")) {
  install.packages("remotes")
}
remotes::install_github("donaldRwilliams/BGGM")

Y <- BGGM::bfi[,1:10]

membership <- c(rep("a", 5), rep("c", 5))

# fit model
fit <- estimate(Y = Y,
                analytic = FALSE,
                iter = 1000)

# list of columns belowinging in each group
# e.g., first 5 are "a", last 5 are "c"

membership <- c(rep("a", 5), rep("c", 5))

f <- function(x,...){
  
  assortment.discrete(x, ...)$r
  
  }

net_stat <- roll_your_own(object = fit,
                          FUN = f,
                          types = membership,
                          weighted = TRUE,
                          SE = FALSE, M = 1)

hist(net_stat)

[donaldRwilliams]

To follow up, here is another example with expected influence.

if (!requireNamespace("remotes")) {
  install.packages("remotes")
}
remotes::install_github("donaldRwilliams/BGGM")

# need this pacakges
library(networktools)
library(ggplot2)
library(reshape)
library(BGGM)

# data from net tools
Y <- depression

# fit model
fit <- estimate(Y = Y, iter = 5000)

# define function
f <- function(x,...){
  expectedInf(x,...)$step1
}

# compute
net_stat <- roll_your_own(object = fit,
                          FUN = f)

# colmeans
colMeans(t(net_stat))

# full distribution
hist(net_stat[1,])

# reshape
results <-  reshape::melt(t(net_stat)) 

# plot
ggplot(results, aes(x =as.factor(X2), 
                    y = value, 
                    group = as.factor(X2))) + 
  geom_boxplot() +
  scale_x_discrete(labels = colnames(Y)) +
  theme(axis.text.x = element_text(angle = 90, 
                                   vjust = 0.5, 
                                   hjust=1))

[donaldRwilliams]

Added print function, say, for the net_stat object..

# from previous examples

# print
net_stat

# output

BGGM: Bayesian Gaussian Graphical Models 
--- 
Network Stats: Roll Your Own
Posterior Samples: 5000 
--- 
Estimates: 

 Node Post.mean Post.sd Cred.lb Cred.ub
    1     0.995   0.043   0.910   1.079
    2     0.991   0.045   0.901   1.077
    3     0.362   0.032   0.299   0.425
    4     0.727   0.042   0.644   0.809
    5     0.820   0.046   0.730   0.913
    6     0.834   0.043   0.752   0.919
    7     0.968   0.045   0.878   1.054
    8     0.765   0.043   0.682   0.850
    9     1.165   0.046   1.076   1.258
--- 

[donaldRwilliams]

The class roll_your_own is now fully supported, including summary output (above) and plot. If this is what you had in mind, then we can close this issue.

plot(net_stat)

[mjbsp]

This is great! Never expected that I'd go to bed and wake up with access to the feature!

[donaldRwilliams]

@mjbsp No problem. It was a really good idea :-)

[donaldRwilliams]

Bridge strength

#######################################
### example 3: mixed data & bridge ####
#######################################

# bridge from this package
library(networkTools)

# data
Y <- ptsd

fit <- estimate(Y,
                type = "mixed",
                iter = 1000)

# clusters
communities <- substring(colnames(Y), 1, 1)

f <- function(x, ...){
 networktools::bridge(x, ...)$`Bridge Strength`
}

net_stat <- roll_your_own(fit,
                          FUN = f,
                          select = TRUE,
                          communities = communities)

# print
net_stat

#plot
plot(net_stat)