Nothing
#' Structure dependence test
#'
#' Test whether or not a network estimates can be considered structurally dependent, i.e.
#' a function of the network structure. By rewiring the graph and calculating
#' a particular statistic \eqn{t}, the test compares the observed mean of \eqn{t}
#' against the empirical distribution of it obtained from rewiring the network.
#'
#' @param graph A \code{\link{diffnet}} graph.
#' @param statistic A function that returns either a scalar or a vector.
#' @param R Integer scalar. Number of repetitions.
#' @param rewire.args List. Arguments to be passed to \code{\link{rewire_graph}}
#' @param x A \code{diffnet_struct_test} class object.
#' @param ... Further arguments passed to the method (see details).
#' @param main Character scalar. Title of the histogram.
#' @param xlab Character scalar. x-axis label.
#' @param breaks Passed to \code{\link{hist}}.
#' @param annotated Logical scalar. When TRUE marks the observed data average and the simulated data average.
#' @return A list of class \code{diffnet_struct_test} containing the following:
#' \item{graph}{The graph passed to \code{struct_test}.}
#' \item{p.value}{The resulting p-value of the test (see details).}
#' \item{t0}{The observed value of the statistic.}
#' \item{mean_t}{The average value of the statistic applied to the simulated networks.}
#' \item{R}{Number of simulations.}
#' \item{statistic}{The function \code{statistic} passed to \code{struct_test}.}
#' \item{boot}{A \code{boot} class object as return from the call to \code{boot}.}
#' \item{rewire.args}{The list \code{rewire.args} passed to \code{struct_test}.}
#'
#' @details
#'
#' \code{struct_test} computes the test by generating the null distribution using
#' Monte Carlo simulations (rewiring). \code{struct_test_asymp} computes the
#' test using an asymptotic approximation. While available, we do not recommend
#' using the asymptotic approximation since it has not shown good results when
#' compared to the MC approximation. Furthermore, the asymptotic version has only
#' been implemented for \code{graph} as static graph.
#'
#' The output from the \code{hist} method is the same as \code{\link{hist.default}}.
#'
#' \code{struct_test} is a wrapper for the function \code{\link[boot:boot]{boot}} from the
#' \pkg{boot} package. Instead of resampling data--vertices or edges--in each iteration the function
#' rewires the original graph using \code{\link{rewire_graph}} and applies
#' the function defined by the user in \code{statistic}.
#'
#' The default values to \code{rewire_graph} via \code{rewire.args} are:
#' \tabular{ll}{
#' \code{p} \tab Number or Integer with default \code{n_rewires(graph)}. \cr
#' \code{undirected} \tab Logical scalar with default \code{getOption("diffnet.undirected", FALSE)}. \cr
#' \code{copy.first} \tab Logical scalar with \code{TRUE}. \cr
#' \code{algorithm} \tab Character scalar with default \code{"swap"}.
#' }
#'
#' In \code{struct_test} \code{\dots} are passed to \code{boot}, otherwise are passed
#' to the corresponding method (\code{\link{hist}} for instance).
#'
#' From the \code{print} method, p-value for the null of the statistic been
#' equal between graph and its rewired versions is computed as follows
#'
#'
#' \deqn{%
#' p(\tau)=2\times\min\left(\mbox{Pr}(t\leq\tau), \mbox{Pr}(t\geq\tau)\right) %
#' }{ %
#' p(tau) = 2*min[Pr(t<=tau), Pr(t>=tau)] %
#' }
#'
#' Where \eqn{\mbox{Pr}\{\cdot\}}{Pr(.)} is approximated using the
#' Empirical Distribution Function retrieved from the simulations.
#'
#' For the case of the asymptotic approximation, under the null we have
#'
#' \deqn{%
#' \sqrt{n}\left(\hat\beta(Y,G)-\mu_\beta\right)\sim^d\mbox{N}\left(0,\sigma_\beta^2\right)
#' }{%
#' sqrt(n)*[mean(beta) - E[beta]]/Var[beta] ~ N(0,1)
#' }
#'
#'
#' The test is actually on development by Vega Yon and Valente. A copy of the
#' working paper can be distributed upon request to \email{g.vegayon@gmail.com}.
#'
#' The function \code{n_rewires} proposes a vector of number of rewirings that
#' are performed in each iteration.
#'
#' @family Functions for inference
#' @references
#' Vega Yon, George G. and Valente, Thomas W. (On development).
#'
#' Davidson, R., & MacKinnon, J. G. (2004). Econometric Theory and Methods. New York:
#' Oxford University Press.
#' @examples
#' # Creating a random graph
#' set.seed(881)
#' diffnet <- rdiffnet(100, 5, seed.graph="small-world")
#'
#' # Testing structure-dependency of threshold
#' res <- struct_test(
#' diffnet,
#' function(g) mean(threshold(g), na.rm=TRUE),
#' R=100
#' )
#'
#' res
#' hist(res)
#'
#' # Adding a legend
#' legend("topright", bty="n",
#' legend=c(
#' expression(t[0]:~Baseline),
#' expression(t:~Rewired~average)
#' )
#' )
#'
#' # Concatenating results
#' c(res, res)
#'
#' # Running in parallel fashion
#' res <- struct_test(
#' diffnet, function(g) mean(threshold(g), na.rm=TRUE),
#' R=100, ncpus=2, parallel="multicore"
#' )
#'
#' res
#'
#' hist(res)
#' @author George G. Vega Yon
#' @name struct_test
NULL
#' @export
#' @param p Either a Numeric scalar or vector of length \code{nslices(graph)-1}
#' with the number of rewires per links.
#' @rdname struct_test
n_rewires <- function(graph, p=c(20L, rep(.1, nslices(graph) - 1))) {
as.integer(round(unlist(nlinks(graph))*p))
}
struct_test_pval <- function(meanobs, meansim) {
n <- length(meanobs)
ans <- vector("numeric",n)
for (i in 1:n)
ans[i] <- 2*min(mean(meansim[,i] < meanobs[i]),
mean(meansim[,i] > meanobs[i]))
ans
}
#' @rdname struct_test
#' @export
struct_test <- function(
graph,
statistic,
R,
rewire.args=list(),
...
) {
# Checking defaults
if (!length(rewire.args$p)) rewire.args$p <- n_rewires(graph)
if (!length(rewire.args$undirected))
rewire.args$undirected <- getOption("diffnet.undirected", FALSE)
if (!length(rewire.args$copy.first)) rewire.args$copy.first <- TRUE
if (!length(rewire.args$algorithm)) rewire.args$algorithm <- "swap"
# Preparing the call to boot
rewire.args$graph <- graph
statisticpll <- function(d, i, fn, rewire.args, ...) {
fn(do.call(rewire_graph, rewire.args),...)
}
# Calling boot
boot_res <- boot::boot(1, statisticpll, R=R, fn=statistic, rewire.args=rewire.args,
...)
# The t0 must be applied with no rewiring!
boot_res$t0 <- statistic(graph)
# Calc pval
# To be conservative, in a two tail test we use the min of the two
# So, following davidson & mckinnon Confidence intrval section,
# p(tau) = 2 * min[F(tau), 1-F(tau)]
p.value <- struct_test_pval(boot_res$t0, boot_res$t)
# Creating the object
out <- list(
graph = graph,
p.value = p.value,
t0 = boot_res$t0,
mean_t = colMeans(boot_res$t, na.rm = TRUE),
var_t = apply(boot_res$t, 2, var, na.rm = TRUE),
R = R,
statistic = statistic,
boot = boot_res,
rewire.args = rewire.args
)
return(structure(out, class="diffnet_struct_test"))
}
#' @export
#' @param recursive Ignored
#' @rdname struct_test
c.diffnet_struct_test <- function(..., recursive=FALSE) {
# Checking arguments names
args <- list(...)
nm <- lapply(args, names)
if (!all(sapply(nm, function(x) identical(x, nm[[1]]))))
stop("arguments are not all the same type of \"diffnet_struct_test\" object")
# Checking graph dim
res <- args[[1]]
res$boot <- do.call(c, lapply(args, "[[", "boot"))
# res$p.value <- with(res, 2*min(mean(boot$t < boot$t0),
# mean(boot$t > boot$t0)))
res$p.value <- struct_test_pval(res$boot$t0, res$boot$t)
res$mean_t <- colMeans(res$boot$t, na.rm=TRUE)
res$R <- res$boot$R
res
}
#' @export
#' @rdname struct_test
print.diffnet_struct_test <- function(x, ...) {
# Neat column printing
netcol <- function(obs, expe, pval) {
txt <- paste0(sprintf(" %10.4f %10.4f %10.4f", obs, expe, pval), collapse="\n")
paste(
sprintf(" %10s %10s %10s", "observed", "expected", "p.val"),
txt, sep = "\n"
)
}
with(x, {
nsim <- ifelse(!is.na(R), R, 0)
cat("\nStructure dependence test\n",
"# Simulations : ", formatC(nsim, digits = 0, format = "f", big.mark = ","),"\n",
"# nodes : ", formatC(nnodes(x$graph), digits = 0, format = "f", big.mark = ","),"\n",
"# of time periods : ", formatC(nslices(x$graph), digits = 0, format = "f", big.mark = ","),"\n",
paste(rep("-",80), collapse=""),"\n",
" H0: E[beta(Y,G)|G] - E[beta(Y,G)] = 0 (no structure dependency)\n",
netcol(t0, mean_t, p.value),
sep="")
})
invisible(x)
}
#' @export
#' @param b0 Character scalar. When \code{annotated=TRUE}, label for the value of \code{b0}.
#' @param b Character scalar. When \code{annotated=TRUE}, label for the value of \code{b}.
#' @param ask Logical scalar. When \code{TRUE}, asks the user to type \code{<Enter>} to see each plot (as
#' many as statistics where computed).
#' @rdname struct_test
hist.diffnet_struct_test <- function(
x,
main="Empirical Distribution of Statistic",
xlab=expression(Values~of~t),
breaks=20,
annotated=TRUE,
b0=expression(atop(plain("") %up% plain("")), t[0]),
b =expression(atop(plain("") %up% plain("")), t[]),
ask = TRUE,
...) {
# Par parameters
if (ask & interactive()) {
oldpar <- par(ask=ask)
on.exit(par(oldpar))
}
out <- vector("list", ncol(x$boot$t))
for (i in 1:ncol(x$boot$t)) {
# Computing the histogram and its range
out[[i]] <- hist(x$boot$t[,i], breaks=breaks, plot=FALSE)
ran <- range(out[[i]]$mids, na.rm = TRUE)
if (annotated) {
# Annotated version (pretty)
mt <- mean(x$boot$t[,i], na.rm=TRUE)
ran <- range(c(ran, mt, x$boot$t0[i]), na.rm=TRUE)
hist(x$boot$t[,i], breaks=breaks, main=main, xlab=xlab, xlim = ran, ...)
} else {
# Not annotated version
hist(x$boot$t[,i], breaks=breaks, main=main, xlab=xlab, ...)
}
# Adding margin note
if (annotated) {
mtext(b0, side = 1, at=x$boot$t0[i])
mtext(b, side = 1, at=mt)
}
}
invisible(out)
}
# #' @rdname struct_test
# boot_thr <- struct_test(net, function(x) {
# t <- threshold(x)
# cbind(mean=mean(t, na.rm=TRUE), sd=sd(t, na.rm = TRUE))
# }
# , nsim)
#
# ttest<-(boot_thr$t[1,] - mean(threshold(net), na.rm=TRUE))/boot_thr$t[2,]
#' @export
#' @rdname struct_test
#' @param Y Numeric vector of length \eqn{n}.
#' @param statistic_name Character scalar. Name of the metric to compute. Currently
#' this can be either \code{"distance"},\code{">"},\code{"<"},\code{"=="}, \code{">="},
#' or \code{"<="}.
struct_test_asymp <- function(
graph, Y,
statistic_name="distance", p=2.0, ...) {
# Distance metric
D <- vertex_covariate_compare(graph, Y, statistic_name)
# Computing observed mean
m_obs <- mean(Matrix::rowSums(D*graph)/(Matrix::rowSums(graph) + 1e-15), na.rm=TRUE)
# Checking NA
isna <- which(is.na(Y))
if (length(isna)) {
warning("There are some NA values in Y: ",paste(isna, collapse=", "))
Y <- Y[-isna]
}
# Computing theoreticals
m_null <- struct_test_mean(Y, statistic_name)
v_null <- struct_test_var(Y, statistic_name)
p.value <- pnorm(sqrt(nnodes(graph))*(m_obs-m_null)/sqrt(v_null))
p.value <- ifelse(p.value > .5, 1-p.value, p.value)*2
# Creating the object
out <- list(
graph = graph,
p.value = p.value,
t0 = m_obs,
mean_t = m_null,
var_t = v_null,
R = NA,
statistic = statistic_name,
boot = NA,
rewire.args = NA
)
return(structure(out, class="diffnet_struct_test"))
}
# library(mvtnorm)
#
# rm(list=ls())
# set.seed(123)
# n <- 400
# G <- rgraph_ws(n,8,.2)
# G <- G/(Matrix::rowSums(G) + 1e-15)
# X <- cbind(rchisq(n, 4))
#
# # Case 1
# ape::Moran.I(as.vector(X), as.matrix(G))
# struct_test_approx(G, X, "distance")
#
# # Case 2
# X <- as.vector(
# solve(Matrix::Diagonal(n) - .4*G) %*% matrix(rmvnorm(1, rep(1,n)), ncol=1)
# )
# ape::Moran.I(as.vector(X), as.matrix(G))
# struct_test_approx(G, X, "distance")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.