#' Compute Custom Network Statistics
#'
#' This function allows for computing custom network statistics for
#' weighted adjacency matrices (partial correlations). The statistics are computed for
#' each of the sampled matrices, resulting in a distribution.
#'
#' @param object An object of class \code{estimate}.
#'
#' @param FUN A custom function for computing the statistic. The first argument must be
#' a partial correlation matrix.
#'
#' @param iter Number of iterations (posterior samples; defaults to the number in the object).
#'
#' @param select Logical. Should the graph be selected ? The default is currently \code{FALSE}.
#'
#' @param cred Numeric. Credible interval between 0 and 1 (default is 0.95) that is used for selecting the graph.
#'
#' @param progress Logical. Should a progress bar be included (defaults to \code{TRUE}) ?
#'
#' @param ... Arguments passed to the function.
#'
#' @return An object defined by \code{FUN}.
#'
#' @details
#'
#' The user has complete control of this function. Hence, care must be taken as to what \code{FUN}
#' returns and in what format. The function should return a single number (one for the entire GGM)
#' or a vector (one for each node). This ensures that the print and \code{\link{plot.roll_your_own}}
#' will work.
#'
#'
#' When \code{select = TRUE}, the graph is selected and then the network statistics are computed based on
#' the weigthed adjacency matrix. This is accomplished internally by multiplying each of the sampled
#' partial correlation matrices by the adjacency matrix.
#'
#' @examples
#' \donttest{
#' ####################################
#' ###### example 1: assortment #######
#' ####################################
#' # assortment
#' library(assortnet)
#'
#' Y <- BGGM::bfi[,1:10]
#' membership <- c(rep("a", 5), rep("c", 5))
#'
#' # fit model
#' fit <- estimate(Y = Y, iter = 250,
#' progress = FALSE)
#'
#' # membership
#' membership <- c(rep("a", 5), rep("c", 5))
#'
#' # define function
#' f <- function(x,...){
#' assortment.discrete(x, ...)$r
#'}
#'
#'
#' net_stat <- roll_your_own(object = fit,
#' FUN = f,
#' types = membership,
#' weighted = TRUE,
#' SE = FALSE, M = 1,
#' progress = FALSE)
#'
#' # print
#' net_stat
#'
#'
#' ############################################
#' ###### example 2: expected influence #######
#' ############################################
#' # expected influence from this package
#' library(networktools)
#'
#' # data
#' Y <- depression
#'
#' # fit model
#' fit <- estimate(Y = Y, iter = 250)
#'
#' # define function
#' f <- function(x,...){
#' expectedInf(x,...)$step1
#' }
#'
#' # compute
#' net_stat <- roll_your_own(object = fit,
#' FUN = f,
#' progress = FALSE)
#'
#' #######################################
#' ### example 3: mixed data & bridge ####
#' #######################################
#' # bridge from this package
#' library(networktools)
#'
#' # data
#' Y <- ptsd[,1:7]
#'
#' fit <- estimate(Y,
#' type = "mixed",
#' iter = 250)
#'
#' # clusters
#' communities <- substring(colnames(Y), 1, 1)
#'
#' # function is slow
#' f <- function(x, ...){
#' bridge(x, ...)$`Bridge Strength`
#' }
#'
#' net_stat <- roll_your_own(fit,
#' FUN = f,
#' select = TRUE,
#' communities = communities,
#' progress = FALSE)
#'
#' }
#'
#' @export
roll_your_own <- function(object,
FUN,
iter = NULL,
select = FALSE,
cred = 0.95,
progress = TRUE,
...){
if(! all( c("estimate", "default") %in% class(object)) ){
stop("class must be 'estimate'")
}
if(!isFALSE(select)){
sel <- select(object, cred = cred)
adj <- sel$adj
} else {
p <- ncol(object$pcor_mat)
adj <- matrix(1, p, p)
}
if(is.null(iter)){
iter <- object$iter
}
pcors <- object$post_samp$pcors[, , 51:(iter + 50)]
if(isTRUE(progress)){
pb <- utils::txtProgressBar(min = 0, max = iter, style = 3)
}
results <- sapply(1:iter, function(x) {
pcors_s <- pcors[, , x] * adj
est <- FUN(pcors_s, ...)
if(isTRUE(progress)){
utils::setTxtProgressBar(pb, x)
}
est
})
returned_object <- list(results = results, iter = iter)
class(returned_object) <- c("BGGM",
"roll_your_own")
return(returned_object)
}
print_roll_your_own <- function(x, cred = 0.95, ...) {
cat("BGGM: Bayesian Gaussian Graphical Models \n")
cat("--- \n")
cat("Network Stats: Roll Your Own\n")
cat("Posterior Samples:", x$iter, "\n")
cat("--- \n")
cat("Estimates: \n\n")
lb <- (1-cred) / 2
ub <- 1 - lb
dims <- dim(x$results)
if(is.null(dims)){
mu <- mean(x$results)
scale <- sd(x$results)
res <- data.frame(Post.mean = round(mean(x$results), 3),
Post.sd = round(sd(x$results), 3),
Cred.lb = round(quantile(x$results, probs = lb), 3),
Cred.ub = round(quantile(x$results, probs = ub), 3) )
} else {
mu <- apply( x$results, 1, mean)
p <- length(mu)
scale <- apply( x$results, 1, sd)
ci_lb <- apply( x$results, 1, quantile, lb)
ci_ub <- apply( x$results, 1, quantile, ub)
res<- data.frame(Node = 1:p,
Post.mean = round(mu, 3),
Post.sd = round(scale, 3),
Cred.lb = round(ci_lb, 3),
Cred.ub = round(ci_ub, 3))
}
print(res, row.names = FALSE)
cat("--- \n")
}
#' Plot \code{roll_your_own} Objects
#'
#' @name plot.roll_your_own
#'
#' @param x An object of class \code{roll_your_own}
#'
#' @param fill Character string specifying the color for the ridges.
#'
#' @param alpha Numeric. Transparancey of the ridges
#'
#' @param ... Currently ignored
#'
#' @return An object of class \code{ggplot}
#'
#' @importFrom ggridges geom_density_ridges
#'
#' @examples
#' \donttest{
#' ####################################
#' ###### example 1: assortment #######
#' ####################################
#' # assortment
#' library(assortnet)
#'
#' Y <- BGGM::bfi[,1:10]
#' membership <- c(rep("a", 5), rep("c", 5))
#'
#' # fit model
#' fit <- estimate(Y = Y, iter = 250,
#' progress = FALSE)
#'
#' # membership
#' membership <- c(rep("a", 5), rep("c", 5))
#'
#' # define function
#' f <- function(x,...){
#' assortment.discrete(x, ...)$r
#'}
#'
#' net_stat <- roll_your_own(object = fit,
#' FUN = f,
#' types = membership,
#' weighted = TRUE,
#' SE = FALSE, M = 1,
#' progress = FALSE)
#'
#' # plot
#' plot(net_stat)
#'
#' }
#' @export
plot.roll_your_own <- function(x, fill = "#CC79A7", alpha = 0.5, ...){
dims <- dim(x$results)
if(is.null(dims)){
dat <- data.frame(x= x$results, y = 1)
plt <- ggplot(dat, aes(x = x,
y = as.factor(y))) +
geom_density_ridges(fill = fill,
alpha = alpha)
} else {
dat <- reshape::melt(t(x$results))
mus <- rowMeans(x$results)
dat$order <- factor(dat$X2, levels = unique(dat$X2)[order(mus)],
labels = unique(dat$X2)[order(mus)] )
zeros <- which(with(dat, tapply(value, X2, sum)) == 0)
if(length(zeros) > 0){
zeros_dat <- data.frame(X1 = 1, X2= names(zeros), value = 0, order = names(zeros))
dat_new <- subset(dat, X2 != names(zeros)[1])
for(i in 2:length(zeros)){
dat_new <- subset(dat_new, X2 != names(zeros)[i])
}
dat <- rbind.data.frame(dat_new, zeros_dat)
}
plt <- ggplot(dat, aes(x = value,
group = as.factor(order),
y = as.factor(order))) +
geom_density_ridges(fill = fill,
alpha = alpha,
rel_min_height = 0.001) +
ylab("")
}
plt
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.