R/utilities.R
In MaStatLab/MultiFit: Multiscale Fisher's Independence Test for Multivariate Dependence
Defines functions
multiSummary
multiTree
Documented in
multiSummary
multiTree
#' @title Summary of significant tests
#' @description Provide a post-hoc summary of significant tests. See vignettes for further examples.
#'
#' @param xy A list, whose first element corresponds to the matrix x as below, and
#' its second element corresponds to the matrix y as below.
#' if \code{xy} is not specified, \code{x} and \code{y} need to be assigned.
#' @param x A matrix, number of columns = dimension of random vector,
#' number of rows = number of observations.
#' @param y A matrix, number of columns = dimension of random vector,
#' number of rows = number of observations.
#' @param fit An object generated by \code{multiFit}.
#' @param alpha Numeric, only tests with adjusted \code{p}-values less than \code{alpha}
#' are presented in the output.
#' @param only.rk Positive integer vector. Show only tests that are ranked according to
#' \code{only.rk} and have adjusted \code{p}-value below \code{alpha}. If left as \code{NULL},
#' all tests with adjusted \code{p}-values less than \code{alpha}
#' are presented in the output.
#' @param use.pval String, choose between \code{"H"} (for Holm), \code{"Hcorrected"}
#' (for Holm on corrected \code{p}-values) or \code{"MH"} for modified Holm.
#' If left \code{NULL}, the order of preference is \code{"MH"}, \code{"Hcorrected"} and
#' then \code{"H"}, according to which is present in the object \code{fit}.
#' @param plot.tests Logical, plot the marginal scatter plots that are associated with
#' the presented significant tests.
#' @param pch Point style for plots. If left as \code{NULL}, a default combination of
#' crosses and bullets is applied.
#' @param rd Numeric, number of figures to round to when presenting ranges of variables.
#' @param plot.margin Logical, plot the marginal scatter plot of the margins that
#' are associated with each significant test, without highlighting which points
#' are conditioned on and are in the discretized 2x2 contingency table.
#' @return List whose elements are \code{significant.tests}, a data frame that summarizes
#' the main features of the tests and their overall ranking by \code{p}-value and
#' \code{original.scale.cuboids}, a list whose number of elements is equal to the number of
#' significant tests (the same number of rows of the data frame \code{significant.tests}). Each
#' element corresponds to a test and is a list whose elements are the marginal ranges of
#' the associated cuboid.
#' @examples
#' set.seed(1)
#' n = 300
#' Dx = Dy = 2
#' x = matrix(0, nrow=n, ncol=Dx)
#' y = matrix(0, nrow=n, ncol=Dy)
#' x[,1] = rnorm(n)
#' x[,2] = runif(n)
#' y[,1] = rnorm(n)
#' y[,2] = sin(5*pi*x[,2]) + 1/5*rnorm(n)
#' fit = multiFit(x=x, y=y, verbose=TRUE)
#' w = multiSummary(x=x, y=y, fit=fit, alpha=0.0001)
#' @export
multiSummary = function(xy, x=NULL, y=NULL, fit, alpha=0.05, only.rk=NULL,
use.pval=NULL, plot.tests=TRUE, pch=NULL, rd=2, plot.margin=FALSE) {
if (!missing(xy)) {
x=xy[[1]]
y=xy[[2]]
}
if (is.vector(x)) x=matrix(x,ncol=1)
if (is.vector(y)) y=matrix(y,ncol=1)
Dx = ncol(x)
Dy = ncol(y)
pvs=fit$pvs
a = attributes(pvs)$parameters
if(Dx!=a$Dx | Dy!=a$Dy) { stop("Mismatch: dimensions of variables differ from p-values data frame dimensions.") }
if(nrow(x)!=a$n | nrow(y)!=a$n) {
stop("Mismatch: number of observations in variables differ from p-values data frame number of observations.")
} else {
n = a$n
}
kx.header = paste0("kx.",1:Dx)
ky.header = paste0("ky.",1:Dy)
k.header = c(kx.header, ky.header)
lx.header = paste0("lx.",1:Dx)
ly.header = paste0("ly.",1:Dy)
l.header = c(lx.header, ly.header)
kl.header = c(k.header, l.header)
all.cuboids = do.call(rbind, lapply(fit$all, `[[`, "cuboids"))
all.cuboids = all.cuboids[!duplicated(all.cuboids[,kl.header]),]
all.tables = do.call(rbind, lapply(fit$all, `[[`, "tables"))
if (is.null(use.pval)) {
if (is.null(a$p.adjust.methods) & !a$correct) { use.pval = "pv" }
if (is.null(a$p.adjust.methods) & a$correct) { use.pval = "pv.correct" }
if (!is.null(a$p.adjust.methods)) { use.pval = a$p.adjust.methods[length(a$p.adjust.methods)] }
}
if (a$cutoff<alpha & (use.pval=="MH")) {
warning("P-values greater than ",a$cutoff," were not adjusted, but the requested significance level is ",alpha,".")
}
if (use.pval=="pv") {
warning("Showing tests for which p-value is less than ",alpha,". However, individual p-values were not adjusted for multiple testing.")
}
if (use.pval=="pv.correct") {
warning("Showing tests for which corrected p-value is less than ",alpha,". However, individual p-values were not adjusted for multiple testing.")
}
if (a$rank.transform) {
x.tran=apply(x,2,function(x) ecdf(x)(x)-1/n)
x.sort = apply(x,2,sort)
e_cdf.x = apply(x, 2, function(x) { x.sort=sort(x);
1:length(x.sort) / length(x.sort)})
y.tran=apply(y,2,function(y) ecdf(y)(y)-1/n)
e_cdf.y = apply(y, 2, function(y) { y.sort=sort(y);
1:length(y.sort) / length(y.sort)})
y.sort = apply(y,2,sort)
invECDF.x = function(x.tr, col)
{ sapply(col,function(c, xt, xs, ex) {ifelse(xt[c]==1,max(xs[,c]), xs[which(as.numeric(ex[,c]) >= as.numeric(xt[c]+1/n))[1],c])},xt=x.tr, xs=x.sort, ex=e_cdf.x) }
invECDF.y = function(y.tr, col)
{ sapply(col,function(c, yt, ys, ey) {ifelse(yt[c]==1,max(ys[,c]),ys[which(as.numeric(ey[,c]) >= as.numeric(yt[c]+1/n))[1],c])},yt=y.tr, ys=y.sort, ey=e_cdf.y) }
}
if (!a$rank.transform) {
scale01 = function(z){(1-10^(-7))*(z-min(z))/(max(z)-min(z))}
x.tran = apply(x,2,scale01)
y.tran = apply(y,2,scale01)
min.max.x = rbind(apply(x, 2, min),apply(x, 2, max))
invScale01.x=function(x.tr, col) {
sapply(col, function(c, xt, mmx) {ifelse(as.numeric(xt[c])==1,as.numeric(mmx[2,c]),
(as.numeric(mmx[2,c])-as.numeric(mmx[1,c]))*as.numeric(xt[c])/(1-10^(-7))+as.numeric(mmx[1,c]))},xt=x.tr, mmx=min.max.x) }
min.max.y = rbind(apply(y, 2, min),apply(y, 2, max))
invScale01.y=function(y.tr, col) {
sapply(col, function(c, yt, mmy) {ifelse(as.numeric(yt[c])==1,as.numeric(mmy[2,c]),
(as.numeric(mmy[2,c])-as.numeric(mmy[1,c]))*as.numeric(yt[c])/(1-10^(-7))+as.numeric(mmy[1,c]))},yt=y.tr, mmy=min.max.y) }
}
significant.pvs.mask = (pvs[,use.pval]<alpha & !is.na(pvs$pv))
if (sum(significant.pvs.mask)==0) { warning(paste0("No adjusted p-values are below ",alpha),
call. = FALSE); return(NULL) }
significant.pvs = pvs[significant.pvs.mask,use.pval]
significant.tests = all.tables[significant.pvs.mask,,drop=FALSE]
# if(!is.matrix(significant.tests)) significant.tests = matrix(significant.tests, nrow=1)
significant.cuboids = NULL
cuboids = NULL
order.p=order(significant.pvs)
if (is.null(only.rk)) {
only.rk = 1:nrow(significant.tests)
} else {
if (max(only.rk)>nrow(significant.tests)) only.rk = only.rk[only.rk<=nrow(significant.tests)]
}
if(nrow(significant.tests)==0) {
cat("\nThere are no tests with a p-value of less than ",alpha," (ordered from most- to less-significant).\n",sep="")
} else {
cat("\nThe following tests had a p-value of less than ",alpha,":\n",sep="")
for (t in only.rk) {
s=significant.tests[order.p[t],]
w=all.cuboids[all.cuboids$cuboid.idx==s["cuboid.idx"],]
significant.cuboids = rbind(significant.cuboids, w)
max.x.tran = ((2*w[lx.header]-0)*2^(-w[kx.header]-1))
if (a$rank.transform) max.x = invECDF.x(max.x.tran, 1:Dx)
if (!a$rank.transform) max.x = invScale01.x(max.x.tran, 1:Dx)
min.x.tran = ((2*w[lx.header]-2)*2^(-w[kx.header]-1))
if (a$rank.transform) min.x = invECDF.x(min.x.tran, 1:Dx)
if (!a$rank.transform) min.x = invScale01.x(min.x.tran, 1:Dx)
non.triv.idx.x = which(max.x.tran-min.x.tran<1)
max.y.tran = ((2*w[ly.header]-0)*2^(-w[ky.header]-1))
if (a$rank.transform) max.y = invECDF.y(max.y.tran, 1:Dy)
if (!a$rank.transform) max.y = invScale01.y(max.y.tran, 1:Dy)
min.y.tran = ((2*w[ly.header]-2)*2^(-w[ky.header]-1))
if (a$rank.transform) min.y = invECDF.y(min.y.tran, 1:Dy)
if (!a$rank.transform) min.y = invScale01.y(min.y.tran, 1:Dy)
non.triv.idx.y = which(max.y.tran-min.y.tran<1)
cuboid.x = lapply(1:Dx,
function(i) {
b = t(c(min.x[i], max.x[i]));
rownames(b)=i; colnames(b)=c("min","max"); b})
names(cuboid.x) = paste0("x",1:Dx)
cuboid.y = lapply(1:Dy,
function(j) {
b = t(c(min.y[j], max.y[j]));
rownames(b)=j; colnames(b)=c("min","max"); b})
names(cuboid.y) = paste0("y",1:Dy)
cuboid.add = list(c(cuboid.x, cuboid.y))
names(cuboid.add) = t
cuboids = c(cuboids, cuboid.add)
cat("Ranked #",t,", Test ",s["test.idx"],": x",as.integer(s["i"])," and y",as.integer(s["j"]),sep="")
if (length(non.triv.idx.x)>0 | length(non.triv.idx.y)>0) cat(" | ")
cat(sapply(non.triv.idx.x,function(idx) paste0(round(cuboid.x[[idx]][,"min"],rd),"<=","x",idx,"<",round(cuboid.x[[idx]][,"max"],rd))),sep=", ")
if (length(non.triv.idx.x)>0 & length(non.triv.idx.y)>0) cat(", ")
cat(sapply(non.triv.idx.y,function(idx) paste0(round(cuboid.y[[idx]][,"min"],rd),"<=","y",idx,"<",round(cuboid.y[[idx]][,"max"],rd))),sep=", ")
cat(" (p-value=",as.numeric(significant.pvs[order.p[t]]),")\n",sep="")
if (plot.tests) {
i=as.integer(s["i"]); j=as.integer(s["j"])
kx=as.integer(w[kx.header]); lx=as.integer(w[lx.header])
ky=as.integer(w[ky.header]); ly=as.integer(w[ly.header])
xl = (lx-1)*2^(-kx)
x.low = matrix(rep(xl, n), nrow=n, byrow=TRUE)
xh = lx*2^(-kx)
x.high = matrix(rep(xh, n), nrow=n, byrow=TRUE)
x.mask = x.tran>=x.low & x.tran<x.high
yl = (ly-1)*2^(-ky)
y.low = matrix(rep(yl, n), nrow=n, byrow=TRUE)
yh = ly*2^(-ky)
y.high = matrix(rep(yh, n),nrow=n, byrow=TRUE)
y.mask = y.tran>=y.low & y.tran<y.high
mask = apply(x.mask, 1, all) & apply(y.mask, 1, all)
x.loc = x[,i][mask]
y.loc = y[,j][mask]
xl.cond = (lx[-i]-1)*2^(-kx[-i])
x.low.cond = matrix(rep(xl.cond, n), nrow=n, byrow=TRUE)
xh.cond = lx[-i]*2^(-kx[-i])
x.high.cond = matrix(rep(xh.cond, n), nrow=n, byrow=TRUE)
x.mask.cond = x.tran[,-i]>=x.low.cond & x.tran[,-i]<x.high.cond
yl.cond = (ly[-j]-1)*2^(-ky[-j])
y.low.cond = matrix(rep(yl.cond, n), nrow=n, byrow=TRUE)
yh.cond = ly[-j]*2^(-ky[-j])
y.high.cond = matrix(rep(yh.cond, n),nrow=n, byrow=TRUE)
y.mask.cond = y.tran[,-j]>=y.low.cond & y.tran[,-j]<y.high.cond
mask.cond = apply(x.mask.cond, 1, all) & apply(y.mask.cond, 1, all)
if (is.null(pch)) { pch0="x"; pch1=20 } else { pch0=pch1=pch }
if (plot.margin) {
plot(x[,i],y[,j], col="grey", pch=pch0, xlab=paste0("x",i),
ylab=paste0("y",j), main="Margins")
}
plot(x[,i],y[,j], col="grey", pch=pch0, xlab=paste0("x",i),
ylab=paste0("y",j), main=paste0("Ranked #",t,": Test ",s["test.idx"]))
points(x[,i][mask.cond], y[,j][mask.cond], col="orange", pch=pch1)
points(x.loc, y.loc, col="red", pch=pch1)
x.i.min = cuboids[[which(t==only.rk)]][[paste0("x",i)]][,"min"]
x.i.max = cuboids[[which(t==only.rk)]][[paste0("x",i)]][,"max"]
y.j.min = cuboids[[which(t==only.rk)]][[paste0("y",j)]][,"min"]
y.j.max = cuboids[[which(t==only.rk)]][[paste0("y",j)]][,"max"]
abline(v=x.i.min, col="blue",lty=2)
abline(v=x.i.max, col="blue",lty=2)
abline(h=y.j.min, col="blue",lty=2)
abline(h=y.j.max, col="blue",lty=2)
segments(x0=x.i.min,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.max,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=(x.i.max+x.i.min)/2,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.min,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=(y.j.max+y.j.min)/2,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.max,col="blue")
}
}
}
s0 = matrix(only.rk,ncol=1)
s1 = matrix(significant.tests[order.p[only.rk],], nrow=length(only.rk))
colnames(s1)=colnames(significant.tests)
s2 = significant.pvs[order.p[only.rk]]
s3 = significant.cuboids[,kl.header]
S = cbind(s0,s1,s2,s3)
colnames(S)[1] = "overall.rank"
colnames(S)[10] = use.pval
rownames(S) = NULL
return(invisible(list(significant.tests=S,
original.scale.cuboids = cuboids)))
}
# ==============================================================================
#' @title Plot tree structure of tests on 2x2 contingency tables
#' @description Plot a post-hoc tree of all tests or all significant tests on 2x2 discretized
#' contingency tables. See vignettes for examples.
#'
#' @param xy A list (optional), whose first element corresponds to the matrix x as below, and
#' its second element corresponds to the matrix y as below.
#' if \code{xy} is not specified, \code{x} and \code{y} need to be assigned. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param x A matrix (optional), number of columns = dimension of random vector,
#' number of rows = number of observations. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param y A matrix (optional), number of columns = dimension of random vector,
#' number of rows = number of observations. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param fit An object generated by \code{multiFit}.
#' @param show.all Logical. If \code{TRUE}, all tests are shown. If \code{FALSE}
#' only tests who were ranked in each resolution amongst the top \code{M} ranking tests
#' are shown. See \code{?multiFit} for an explanation about the parameter \code{M} and see
#' documentation for further information.
#' @param max.node.size Numeric. Maximal node size. All nodes are scaled between \code{min.node.size} and
#' \code{max.node.size}, where larger nodes are associated smaller \code{p}-values of the corresponding tests
#' on 2x2 contingency tables.
#' @param min.node.size Numeric. Minimal node size. All nodes are scaled between \code{min.node.size} and
#' \code{max.node.size}, where larger nodes are associated smaller \code{p}-values of the corresponding tests
#' on 2x2 contingency tables.
#' @param use.pval String, choose between \code{"H"} (for Holm), \code{"Hcorrected"}
#' (for Holm on corrected \code{p}-values) or \code{"MH"} for modified Holm.
#' If left \code{NULL}, the order of preference is \code{"MH"}, \code{"Hcorrected"} and
#' then \code{"H"}, according to which is present in the object \code{fit}.
#' @param images.path String, path to save \code{png} images of nodes to. If not
#' specified, images are saved to \code{tempdir()}.
#' @param node.name String, prefix for file names for nodes \code{png}s.
#' @param filename String, file name for tree output. If left \code{NULL}, file name
#' is prefixed by \code{multiTree} and ends with system time. See documentation of
#' \code{qgraph::qgraph} for further information.
#' @param filetype String, default is \code{pdf}, See documentation of
#' \code{qgraph::qgraph} for further information.
#' @return The main output of multiTree is a \code{pdf} file with the directed acyclic graph
#' showing tests as nodes.
#' @return In addition, the function returns a list. Its elements are:
#' \code{qgraph.object}, the graphical object generated by the \code{qgraph} function. See
#' the \code{qgraph} package documentation for further details.
#' \code{qgraph.call}, the call for the tree generating function. Arguments for
#' the call: \code{adj}, the adjacency matrix, \code{nodes.size}, a numeric vector with the
#' scaled sizes of the nodes, \code{images}, the file names of the nodes images (may be
#' \code{NULL}), \code{filename} as passed to \code{multiTree} and passed over to \code{qgraph},
#' and \code{filetype} as passed to \code{multiTree} and passed over to \code{qgraph}.
#' @return Other elements of the returned list are \code{pvs.attributes}, the attributes summarizing the
#' data and the tests performed as stored in \code{fit}, and \code{n.nodes}, the number of nodes.
#' @export
multiTree = function(xy, x=NULL, y=NULL, fit, show.all=FALSE,
max.node.size=5, min.node.size=2.5, use.pval=NULL,
images.path=NULL, node.name="node", filename=NULL,
filetype="pdf") {
if (!requireNamespace("qgraph", quietly = TRUE)) {
stop("Package \"qgraph\" is required.",
call. = FALSE)
}
if (!missing(xy)) {
x=xy[[1]]
y=xy[[2]]
}
pvs=fit$pvs
a = attributes(pvs)$parameters
if (!is.null(x) & !is.null(y)) {
if (is.vector(x)) x=matrix(x,ncol=1)
if (is.vector(y)) y=matrix(y,ncol=1)
Dx = ncol(x)
Dy = ncol(y)
if(Dx!=a$Dx | Dy!=a$Dy) { stop("Mismatch: dimensions of variables differ from p-values data frame dimensions.") }
if(nrow(x)!=a$n | nrow(y)!=a$n) {
stop("Mismatch: number of observations in variables differ from p-values data frame number of observations.")
} else {
n = a$n
}
} else {
Dx=a$Dx
Dy=a$Dy
n = a$n
}
kx.header = paste0("kx.",1:Dx)
ky.header = paste0("ky.",1:Dy)
k.header = c(kx.header, ky.header)
lx.header = paste0("lx.",1:Dx)
ly.header = paste0("ly.",1:Dy)
l.header = c(lx.header, ly.header)
kl.header = c(k.header, l.header)
all.cuboids = do.call(rbind, lapply(fit$all, `[[`, "cuboids"))
all.tables = do.call(rbind, lapply(fit$all, `[[`, "tables"))
if (is.null(use.pval)) {
if (is.null(a$p.adjust.methods) & !a$correct) { use.pval = "pv" }
if (is.null(a$p.adjust.methods) & a$correct) { use.pval = "pv.correct" }
if (!is.null(a$p.adjust.methods)) { use.pval = a$p.adjust.methods[length(a$p.adjust.methods)] }
}
tree.struct = NULL
tree.struct = split(all.cuboids$child.of.test, all.cuboids$cuboid.idx)
# suppressMessages(library(png))
# suppressMessages(library(qgraph))
if (!is.null(x) & !is.null(y)) {
if (a$rank.transform) {
x.tran=apply(x,2,function(x) ecdf(x)(x)-1/n)
x.sort = apply(x,2,sort)
e_cdf.x = apply(x, 2, function(x) { x.sort=sort(x);
1:length(x.sort) / length(x.sort)})
y.tran=apply(y,2,function(y) ecdf(y)(y)-1/n)
e_cdf.y = apply(y, 2, function(y) { y.sort=sort(y);
1:length(y.sort) / length(y.sort)})
y.sort = apply(y,2,sort)
invECDF.x = function(x.tr, col)
{ sapply(col,function(c, xt, xs, ex) {ifelse(xt[c]==1,max(xs[,c]), xs[which(as.numeric(ex[,c]) >= as.numeric(xt[c]+1/n))[1],c])},xt=x.tr, xs=x.sort, ex=e_cdf.x) }
invECDF.y = function(y.tr, col)
{ sapply(col,function(c, yt, ys, ey) {ifelse(yt[c]==1,max(ys[,c]),ys[which(as.numeric(ey[,c]) >= as.numeric(yt[c]+1/n))[1],c])},yt=y.tr, ys=y.sort, ey=e_cdf.y) }
}
if (!a$rank.transform) {
scale01 = function(z){(1-10^(-7))*(z-min(z))/(max(z)-min(z))}
x.tran = apply(x,2,scale01)
y.tran = apply(y,2,scale01)
min.max.x = rbind(apply(x, 2, min),apply(x, 2, max))
invScale01.x=function(x.tr, col) {
sapply(col, function(c, xt, mmx) {ifelse(as.numeric(xt[c])==1,as.numeric(mmx[2,c]),
(as.numeric(mmx[2,c])-as.numeric(mmx[1,c]))*as.numeric(xt[c])/(1-10^(-7))+as.numeric(mmx[1,c]))},xt=x.tr, mmx=min.max.x) }
min.max.y = rbind(apply(y, 2, min),apply(y, 2, max))
invScale01.y=function(y.tr, col) {
sapply(col, function(c, yt, mmy) {ifelse(as.numeric(yt[c])==1,as.numeric(mmy[2,c]),
(as.numeric(mmy[2,c])-as.numeric(mmy[1,c]))*as.numeric(yt[c])/(1-10^(-7))+as.numeric(mmy[1,c]))},yt=y.tr, mmy=min.max.y) }
}
}
# Add network notation:
if (!show.all) {
ranked.tests = do.call(c, lapply(fit$all, `[[`, "rank.tests"))
} else {
ranked.tests = rep(TRUE,sum(!is.na(pvs$pv)))
}
nodes = all.tables[!is.na(pvs$pv), "test.idx"][ranked.tests]
n.nodes = length(nodes)
# Initialize adjacency matrix
adj = matrix(FALSE, nrow=n.nodes, ncol=n.nodes)
colnames(adj) = rownames(adj) = nodes
# Fill it!
for (s in 1:n.nodes) {
w = as.character(all.tables[all.tables[,"test.idx"]==nodes[s],"cuboid.idx"])
if (w>1) {
adj[as.character(tree.struct[[w]]),as.character(nodes[s])] = TRUE
}
}
# If data is given, plot the nodes, each png in low quality:
images=NULL
if (!is.null(x) & !is.null(y)) {
if (!requireNamespace("png", quietly = TRUE)) {
stop("Package \"png\" is required.",
call. = FALSE)
}
if (is.null(images.path)) images.path = tempdir()
W = all.cuboids[!duplicated(all.cuboids$cuboid.idx),]
rownames(W) = W$cuboid.idx
rownames(all.tables) = all.tables[,"test.idx"]
for (m in 1:n.nodes) {
w = as.character(all.tables[as.character(nodes[m]),"cuboid.idx"])
kx = unlist(W[w,kx.header])
ky = unlist(W[w,ky.header])
lx = unlist(W[w,lx.header])
ly = unlist(W[w,ly.header])
i = all.tables[nodes[m],"i"]
j = all.tables[nodes[m],"j"]
xl.cond = (lx[-i]-1)*2^(-kx[-i])
x.low.cond = matrix(rep(xl.cond, n), nrow=n, byrow=TRUE)
xh.cond = lx[-i]*2^(-kx[-i])
x.high.cond = matrix(rep(xh.cond, n), nrow=n, byrow=TRUE)
x.mask.cond = x.tran[,-i]>=x.low.cond & x.tran[,-i]<x.high.cond
yl.cond = (ly[-j]-1)*2^(-ky[-j])
y.low.cond = matrix(rep(yl.cond, n), nrow=n, byrow=TRUE)
yh.cond = ly[-j]*2^(-ky[-j])
y.high.cond = matrix(rep(yh.cond, n),nrow=n, byrow=TRUE)
y.mask.cond = y.tran[,-j]>=y.low.cond & y.tran[,-j]<y.high.cond
mask.cond = apply(x.mask.cond, 1, all) & apply(y.mask.cond, 1, all)
xl = (lx[i]-1)*2^(-kx[i])
xh = lx[i]*2^(-kx[i])
x.mask = x.tran[,i]>=xl & x.tran[,i]<xh
yl = (ly[j]-1)*2^(-ky[j])
yh = ly[j]*2^(-ky[j])
y.mask = y.tran[,j]>=yl & y.tran[,j]<yh
mask = as.logical(x.mask*y.mask*mask.cond)
x.loc = x[,i][mask]
y.loc = y[,j][mask]
png(file.path(images.path, paste0(node.name, m,".png")), width=80, height=60)
par(mar=(c(1,0.8,0,0)))
plot(x[,i],y[,j], col="grey", pch=".", xlab=NA, ylab=NA, xaxt='n', yaxt='n')
mtext(side=1, text=paste("x",i,sep=""), line=0.1, cex=0.7)
mtext(side=2, text=paste("y",j,sep=""), line=0.1, cex=0.7)
points(x[,i][mask.cond], y[,j][mask.cond], col="orange", pch=".")
points(x.loc, y.loc, col="red", pch=".")
max.x.tran = ((2*lx-0)*2^(-kx-1))
if (a$rank.transform) max.x = invECDF.x(max.x.tran, 1:Dx)
if (!a$rank.transform) max.x = invScale01.x(max.x.tran, 1:Dx)
min.x.tran = ((2*lx-2)*2^(-kx-1))
if (a$rank.transform) min.x = invECDF.x(min.x.tran, 1:Dx)
if (!a$rank.transform) min.x = invScale01.x(min.x.tran, 1:Dx)
max.y.tran = ((2*ly-0)*2^(-ky-1))
if (a$rank.transform) max.y = invECDF.y(max.y.tran, 1:Dy)
if (!a$rank.transform) max.y = invScale01.y(max.y.tran, 1:Dy)
min.y.tran = ((2*ly-2)*2^(-ky-1))
if (a$rank.transform) min.y = invECDF.y(min.y.tran, 1:Dy)
if (!a$rank.transform) min.y = invScale01.y(min.y.tran, 1:Dy)
x.i.min = min.x[i]
x.i.max = max.x[i]
y.j.min = min.y[j]
y.j.max = max.y[j]
abline(v=x.i.min, col="blue",lty=2)
abline(v=x.i.max, col="blue",lty=2)
abline(h=y.j.min, col="blue",lty=2)
abline(h=y.j.max, col="blue",lty=2)
segments(x0=x.i.min,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.max,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=(x.i.max+x.i.min)/2,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.min,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=(y.j.max+y.j.min)/2,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.max,col="blue")
dev.off()
}
images = paste0(images.path,"/node",1:n.nodes,".png")
}
nodes.pvs = pvs[!is.na(pvs$pv), use.pval][ranked.tests]
nodes.size = 1.5*(log(-log(nodes.pvs)+1)+.6)
if (max(nodes.size)>max.node.size) nodes.size = max.node.size*nodes.size/max(nodes.size)
if (max(nodes.size)!=min(nodes.size)) {
nodes.size = min.node.size+(max.node.size-min.node.size)*(nodes.size-min(nodes.size))/(max(nodes.size)-min(nodes.size))
}
if(is.null(filename)) filename=paste0("multiTree_",format(Sys.time(), "%Y-%m-%d_%H:%M:%S"))
if(is.null(filetype)) filetype = "pdf"
qg = qgraph::qgraph(input=adj, images=images, shape='square', vTrans=200, vsize=nodes.size,
border.width = 0, label.cex=.6, filename=filename, filetype=filetype)
return(invisible(list(qgraph.object=qg, qgraph.call="qgraph(input=adj, images=images, shape='square', vTrans=200, vsize=nodes.size, border.width=0, label.cex=.6, filename=filename, filetype=filetype)",
adj=adj, n.nodes = n.nodes, nodes.pvs=nodes.pvs, nodes.size=nodes.size,
images=images, pvs.attributes=a, filename=filename, filetype=filetype)))
}
MaStatLab/MultiFit documentation built on Jan. 11, 2020, 5:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions multiSummary multiTree
Documented in multiSummary multiTree
#' @title Summary of significant tests
#' @description Provide a post-hoc summary of significant tests. See vignettes for further examples.
#'
#' @param xy A list, whose first element corresponds to the matrix x as below, and
#' its second element corresponds to the matrix y as below.
#' if \code{xy} is not specified, \code{x} and \code{y} need to be assigned.
#' @param x A matrix, number of columns = dimension of random vector,
#' number of rows = number of observations.
#' @param y A matrix, number of columns = dimension of random vector,
#' number of rows = number of observations.
#' @param fit An object generated by \code{multiFit}.
#' @param alpha Numeric, only tests with adjusted \code{p}-values less than \code{alpha}
#' are presented in the output.
#' @param only.rk Positive integer vector. Show only tests that are ranked according to
#' \code{only.rk} and have adjusted \code{p}-value below \code{alpha}. If left as \code{NULL},
#' all tests with adjusted \code{p}-values less than \code{alpha}
#' are presented in the output.
#' @param use.pval String, choose between \code{"H"} (for Holm), \code{"Hcorrected"}
#' (for Holm on corrected \code{p}-values) or \code{"MH"} for modified Holm.
#' If left \code{NULL}, the order of preference is \code{"MH"}, \code{"Hcorrected"} and
#' then \code{"H"}, according to which is present in the object \code{fit}.
#' @param plot.tests Logical, plot the marginal scatter plots that are associated with
#' the presented significant tests.
#' @param pch Point style for plots. If left as \code{NULL}, a default combination of
#' crosses and bullets is applied.
#' @param rd Numeric, number of figures to round to when presenting ranges of variables.
#' @param plot.margin Logical, plot the marginal scatter plot of the margins that
#' are associated with each significant test, without highlighting which points
#' are conditioned on and are in the discretized 2x2 contingency table.
#' @return List whose elements are \code{significant.tests}, a data frame that summarizes
#' the main features of the tests and their overall ranking by \code{p}-value and
#' \code{original.scale.cuboids}, a list whose number of elements is equal to the number of
#' significant tests (the same number of rows of the data frame \code{significant.tests}). Each
#' element corresponds to a test and is a list whose elements are the marginal ranges of
#' the associated cuboid.
#' @examples
#' set.seed(1)
#' n = 300
#' Dx = Dy = 2
#' x = matrix(0, nrow=n, ncol=Dx)
#' y = matrix(0, nrow=n, ncol=Dy)
#' x[,1] = rnorm(n)
#' x[,2] = runif(n)
#' y[,1] = rnorm(n)
#' y[,2] = sin(5*pi*x[,2]) + 1/5*rnorm(n)
#' fit = multiFit(x=x, y=y, verbose=TRUE)
#' w = multiSummary(x=x, y=y, fit=fit, alpha=0.0001)
#' @export
multiSummary = function(xy, x=NULL, y=NULL, fit, alpha=0.05, only.rk=NULL,
use.pval=NULL, plot.tests=TRUE, pch=NULL, rd=2, plot.margin=FALSE) {
if (!missing(xy)) {
x=xy[[1]]
y=xy[[2]]
}
if (is.vector(x)) x=matrix(x,ncol=1)
if (is.vector(y)) y=matrix(y,ncol=1)
Dx = ncol(x)
Dy = ncol(y)
pvs=fit$pvs
a = attributes(pvs)$parameters
if(Dx!=a$Dx | Dy!=a$Dy) { stop("Mismatch: dimensions of variables differ from p-values data frame dimensions.") }
if(nrow(x)!=a$n | nrow(y)!=a$n) {
stop("Mismatch: number of observations in variables differ from p-values data frame number of observations.")
} else {
n = a$n
}
kx.header = paste0("kx.",1:Dx)
ky.header = paste0("ky.",1:Dy)
k.header = c(kx.header, ky.header)
lx.header = paste0("lx.",1:Dx)
ly.header = paste0("ly.",1:Dy)
l.header = c(lx.header, ly.header)
kl.header = c(k.header, l.header)
all.cuboids = do.call(rbind, lapply(fit$all, `[[`, "cuboids"))
all.cuboids = all.cuboids[!duplicated(all.cuboids[,kl.header]),]
all.tables = do.call(rbind, lapply(fit$all, `[[`, "tables"))
if (is.null(use.pval)) {
if (is.null(a$p.adjust.methods) & !a$correct) { use.pval = "pv" }
if (is.null(a$p.adjust.methods) & a$correct) { use.pval = "pv.correct" }
if (!is.null(a$p.adjust.methods)) { use.pval = a$p.adjust.methods[length(a$p.adjust.methods)] }
}
if (a$cutoff<alpha & (use.pval=="MH")) {
warning("P-values greater than ",a$cutoff," were not adjusted, but the requested significance level is ",alpha,".")
}
if (use.pval=="pv") {
warning("Showing tests for which p-value is less than ",alpha,". However, individual p-values were not adjusted for multiple testing.")
}
if (use.pval=="pv.correct") {
warning("Showing tests for which corrected p-value is less than ",alpha,". However, individual p-values were not adjusted for multiple testing.")
}
if (a$rank.transform) {
x.tran=apply(x,2,function(x) ecdf(x)(x)-1/n)
x.sort = apply(x,2,sort)
e_cdf.x = apply(x, 2, function(x) { x.sort=sort(x);
1:length(x.sort) / length(x.sort)})
y.tran=apply(y,2,function(y) ecdf(y)(y)-1/n)
e_cdf.y = apply(y, 2, function(y) { y.sort=sort(y);
1:length(y.sort) / length(y.sort)})
y.sort = apply(y,2,sort)
invECDF.x = function(x.tr, col)
{ sapply(col,function(c, xt, xs, ex) {ifelse(xt[c]==1,max(xs[,c]), xs[which(as.numeric(ex[,c]) >= as.numeric(xt[c]+1/n))[1],c])},xt=x.tr, xs=x.sort, ex=e_cdf.x) }
invECDF.y = function(y.tr, col)
{ sapply(col,function(c, yt, ys, ey) {ifelse(yt[c]==1,max(ys[,c]),ys[which(as.numeric(ey[,c]) >= as.numeric(yt[c]+1/n))[1],c])},yt=y.tr, ys=y.sort, ey=e_cdf.y) }
}
if (!a$rank.transform) {
scale01 = function(z){(1-10^(-7))*(z-min(z))/(max(z)-min(z))}
x.tran = apply(x,2,scale01)
y.tran = apply(y,2,scale01)
min.max.x = rbind(apply(x, 2, min),apply(x, 2, max))
invScale01.x=function(x.tr, col) {
sapply(col, function(c, xt, mmx) {ifelse(as.numeric(xt[c])==1,as.numeric(mmx[2,c]),
(as.numeric(mmx[2,c])-as.numeric(mmx[1,c]))*as.numeric(xt[c])/(1-10^(-7))+as.numeric(mmx[1,c]))},xt=x.tr, mmx=min.max.x) }
min.max.y = rbind(apply(y, 2, min),apply(y, 2, max))
invScale01.y=function(y.tr, col) {
sapply(col, function(c, yt, mmy) {ifelse(as.numeric(yt[c])==1,as.numeric(mmy[2,c]),
(as.numeric(mmy[2,c])-as.numeric(mmy[1,c]))*as.numeric(yt[c])/(1-10^(-7))+as.numeric(mmy[1,c]))},yt=y.tr, mmy=min.max.y) }
}
significant.pvs.mask = (pvs[,use.pval]<alpha & !is.na(pvs$pv))
if (sum(significant.pvs.mask)==0) { warning(paste0("No adjusted p-values are below ",alpha),
call. = FALSE); return(NULL) }
significant.pvs = pvs[significant.pvs.mask,use.pval]
significant.tests = all.tables[significant.pvs.mask,,drop=FALSE]
# if(!is.matrix(significant.tests)) significant.tests = matrix(significant.tests, nrow=1)
significant.cuboids = NULL
cuboids = NULL
order.p=order(significant.pvs)
if (is.null(only.rk)) {
only.rk = 1:nrow(significant.tests)
} else {
if (max(only.rk)>nrow(significant.tests)) only.rk = only.rk[only.rk<=nrow(significant.tests)]
}
if(nrow(significant.tests)==0) {
cat("\nThere are no tests with a p-value of less than ",alpha," (ordered from most- to less-significant).\n",sep="")
} else {
cat("\nThe following tests had a p-value of less than ",alpha,":\n",sep="")
for (t in only.rk) {
s=significant.tests[order.p[t],]
w=all.cuboids[all.cuboids$cuboid.idx==s["cuboid.idx"],]
significant.cuboids = rbind(significant.cuboids, w)
max.x.tran = ((2*w[lx.header]-0)*2^(-w[kx.header]-1))
if (a$rank.transform) max.x = invECDF.x(max.x.tran, 1:Dx)
if (!a$rank.transform) max.x = invScale01.x(max.x.tran, 1:Dx)
min.x.tran = ((2*w[lx.header]-2)*2^(-w[kx.header]-1))
if (a$rank.transform) min.x = invECDF.x(min.x.tran, 1:Dx)
if (!a$rank.transform) min.x = invScale01.x(min.x.tran, 1:Dx)
non.triv.idx.x = which(max.x.tran-min.x.tran<1)
max.y.tran = ((2*w[ly.header]-0)*2^(-w[ky.header]-1))
if (a$rank.transform) max.y = invECDF.y(max.y.tran, 1:Dy)
if (!a$rank.transform) max.y = invScale01.y(max.y.tran, 1:Dy)
min.y.tran = ((2*w[ly.header]-2)*2^(-w[ky.header]-1))
if (a$rank.transform) min.y = invECDF.y(min.y.tran, 1:Dy)
if (!a$rank.transform) min.y = invScale01.y(min.y.tran, 1:Dy)
non.triv.idx.y = which(max.y.tran-min.y.tran<1)
cuboid.x = lapply(1:Dx,
function(i) {
b = t(c(min.x[i], max.x[i]));
rownames(b)=i; colnames(b)=c("min","max"); b})
names(cuboid.x) = paste0("x",1:Dx)
cuboid.y = lapply(1:Dy,
function(j) {
b = t(c(min.y[j], max.y[j]));
rownames(b)=j; colnames(b)=c("min","max"); b})
names(cuboid.y) = paste0("y",1:Dy)
cuboid.add = list(c(cuboid.x, cuboid.y))
names(cuboid.add) = t
cuboids = c(cuboids, cuboid.add)
cat("Ranked #",t,", Test ",s["test.idx"],": x",as.integer(s["i"])," and y",as.integer(s["j"]),sep="")
if (length(non.triv.idx.x)>0 | length(non.triv.idx.y)>0) cat(" | ")
cat(sapply(non.triv.idx.x,function(idx) paste0(round(cuboid.x[[idx]][,"min"],rd),"<=","x",idx,"<",round(cuboid.x[[idx]][,"max"],rd))),sep=", ")
if (length(non.triv.idx.x)>0 & length(non.triv.idx.y)>0) cat(", ")
cat(sapply(non.triv.idx.y,function(idx) paste0(round(cuboid.y[[idx]][,"min"],rd),"<=","y",idx,"<",round(cuboid.y[[idx]][,"max"],rd))),sep=", ")
cat(" (p-value=",as.numeric(significant.pvs[order.p[t]]),")\n",sep="")
if (plot.tests) {
i=as.integer(s["i"]); j=as.integer(s["j"])
kx=as.integer(w[kx.header]); lx=as.integer(w[lx.header])
ky=as.integer(w[ky.header]); ly=as.integer(w[ly.header])
xl = (lx-1)*2^(-kx)
x.low = matrix(rep(xl, n), nrow=n, byrow=TRUE)
xh = lx*2^(-kx)
x.high = matrix(rep(xh, n), nrow=n, byrow=TRUE)
x.mask = x.tran>=x.low & x.tran<x.high
yl = (ly-1)*2^(-ky)
y.low = matrix(rep(yl, n), nrow=n, byrow=TRUE)
yh = ly*2^(-ky)
y.high = matrix(rep(yh, n),nrow=n, byrow=TRUE)
y.mask = y.tran>=y.low & y.tran<y.high
mask = apply(x.mask, 1, all) & apply(y.mask, 1, all)
x.loc = x[,i][mask]
y.loc = y[,j][mask]
xl.cond = (lx[-i]-1)*2^(-kx[-i])
x.low.cond = matrix(rep(xl.cond, n), nrow=n, byrow=TRUE)
xh.cond = lx[-i]*2^(-kx[-i])
x.high.cond = matrix(rep(xh.cond, n), nrow=n, byrow=TRUE)
x.mask.cond = x.tran[,-i]>=x.low.cond & x.tran[,-i]<x.high.cond
yl.cond = (ly[-j]-1)*2^(-ky[-j])
y.low.cond = matrix(rep(yl.cond, n), nrow=n, byrow=TRUE)
yh.cond = ly[-j]*2^(-ky[-j])
y.high.cond = matrix(rep(yh.cond, n),nrow=n, byrow=TRUE)
y.mask.cond = y.tran[,-j]>=y.low.cond & y.tran[,-j]<y.high.cond
mask.cond = apply(x.mask.cond, 1, all) & apply(y.mask.cond, 1, all)
if (is.null(pch)) { pch0="x"; pch1=20 } else { pch0=pch1=pch }
if (plot.margin) {
plot(x[,i],y[,j], col="grey", pch=pch0, xlab=paste0("x",i),
ylab=paste0("y",j), main="Margins")
}
plot(x[,i],y[,j], col="grey", pch=pch0, xlab=paste0("x",i),
ylab=paste0("y",j), main=paste0("Ranked #",t,": Test ",s["test.idx"]))
points(x[,i][mask.cond], y[,j][mask.cond], col="orange", pch=pch1)
points(x.loc, y.loc, col="red", pch=pch1)
x.i.min = cuboids[[which(t==only.rk)]][[paste0("x",i)]][,"min"]
x.i.max = cuboids[[which(t==only.rk)]][[paste0("x",i)]][,"max"]
y.j.min = cuboids[[which(t==only.rk)]][[paste0("y",j)]][,"min"]
y.j.max = cuboids[[which(t==only.rk)]][[paste0("y",j)]][,"max"]
abline(v=x.i.min, col="blue",lty=2)
abline(v=x.i.max, col="blue",lty=2)
abline(h=y.j.min, col="blue",lty=2)
abline(h=y.j.max, col="blue",lty=2)
segments(x0=x.i.min,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.max,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=(x.i.max+x.i.min)/2,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.min,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=(y.j.max+y.j.min)/2,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.max,col="blue")
}
}
}
s0 = matrix(only.rk,ncol=1)
s1 = matrix(significant.tests[order.p[only.rk],], nrow=length(only.rk))
colnames(s1)=colnames(significant.tests)
s2 = significant.pvs[order.p[only.rk]]
s3 = significant.cuboids[,kl.header]
S = cbind(s0,s1,s2,s3)
colnames(S)[1] = "overall.rank"
colnames(S)[10] = use.pval
rownames(S) = NULL
return(invisible(list(significant.tests=S,
original.scale.cuboids = cuboids)))
}
# ==============================================================================
#' @title Plot tree structure of tests on 2x2 contingency tables
#' @description Plot a post-hoc tree of all tests or all significant tests on 2x2 discretized
#' contingency tables. See vignettes for examples.
#'
#' @param xy A list (optional), whose first element corresponds to the matrix x as below, and
#' its second element corresponds to the matrix y as below.
#' if \code{xy} is not specified, \code{x} and \code{y} need to be assigned. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param x A matrix (optional), number of columns = dimension of random vector,
#' number of rows = number of observations. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param y A matrix (optional), number of columns = dimension of random vector,
#' number of rows = number of observations. If \code{xy},
#' \code{x} and \code{y} are missing or \code{NULL}, the tree nodes are blank. If
#' \code{xy} or \code{x} and \code{y} are provided, nodes are \code{png} images of the marginal
#' scatter plots that are associated with each test.
#' @param fit An object generated by \code{multiFit}.
#' @param show.all Logical. If \code{TRUE}, all tests are shown. If \code{FALSE}
#' only tests who were ranked in each resolution amongst the top \code{M} ranking tests
#' are shown. See \code{?multiFit} for an explanation about the parameter \code{M} and see
#' documentation for further information.
#' @param max.node.size Numeric. Maximal node size. All nodes are scaled between \code{min.node.size} and
#' \code{max.node.size}, where larger nodes are associated smaller \code{p}-values of the corresponding tests
#' on 2x2 contingency tables.
#' @param min.node.size Numeric. Minimal node size. All nodes are scaled between \code{min.node.size} and
#' \code{max.node.size}, where larger nodes are associated smaller \code{p}-values of the corresponding tests
#' on 2x2 contingency tables.
#' @param use.pval String, choose between \code{"H"} (for Holm), \code{"Hcorrected"}
#' (for Holm on corrected \code{p}-values) or \code{"MH"} for modified Holm.
#' If left \code{NULL}, the order of preference is \code{"MH"}, \code{"Hcorrected"} and
#' then \code{"H"}, according to which is present in the object \code{fit}.
#' @param images.path String, path to save \code{png} images of nodes to. If not
#' specified, images are saved to \code{tempdir()}.
#' @param node.name String, prefix for file names for nodes \code{png}s.
#' @param filename String, file name for tree output. If left \code{NULL}, file name
#' is prefixed by \code{multiTree} and ends with system time. See documentation of
#' \code{qgraph::qgraph} for further information.
#' @param filetype String, default is \code{pdf}, See documentation of
#' \code{qgraph::qgraph} for further information.
#' @return The main output of multiTree is a \code{pdf} file with the directed acyclic graph
#' showing tests as nodes.
#' @return In addition, the function returns a list. Its elements are:
#' \code{qgraph.object}, the graphical object generated by the \code{qgraph} function. See
#' the \code{qgraph} package documentation for further details.
#' \code{qgraph.call}, the call for the tree generating function. Arguments for
#' the call: \code{adj}, the adjacency matrix, \code{nodes.size}, a numeric vector with the
#' scaled sizes of the nodes, \code{images}, the file names of the nodes images (may be
#' \code{NULL}), \code{filename} as passed to \code{multiTree} and passed over to \code{qgraph},
#' and \code{filetype} as passed to \code{multiTree} and passed over to \code{qgraph}.
#' @return Other elements of the returned list are \code{pvs.attributes}, the attributes summarizing the
#' data and the tests performed as stored in \code{fit}, and \code{n.nodes}, the number of nodes.
#' @export
multiTree = function(xy, x=NULL, y=NULL, fit, show.all=FALSE,
max.node.size=5, min.node.size=2.5, use.pval=NULL,
images.path=NULL, node.name="node", filename=NULL,
filetype="pdf") {
if (!requireNamespace("qgraph", quietly = TRUE)) {
stop("Package \"qgraph\" is required.",
call. = FALSE)
}
if (!missing(xy)) {
x=xy[[1]]
y=xy[[2]]
}
pvs=fit$pvs
a = attributes(pvs)$parameters
if (!is.null(x) & !is.null(y)) {
if (is.vector(x)) x=matrix(x,ncol=1)
if (is.vector(y)) y=matrix(y,ncol=1)
Dx = ncol(x)
Dy = ncol(y)
if(Dx!=a$Dx | Dy!=a$Dy) { stop("Mismatch: dimensions of variables differ from p-values data frame dimensions.") }
if(nrow(x)!=a$n | nrow(y)!=a$n) {
stop("Mismatch: number of observations in variables differ from p-values data frame number of observations.")
} else {
n = a$n
}
} else {
Dx=a$Dx
Dy=a$Dy
n = a$n
}
kx.header = paste0("kx.",1:Dx)
ky.header = paste0("ky.",1:Dy)
k.header = c(kx.header, ky.header)
lx.header = paste0("lx.",1:Dx)
ly.header = paste0("ly.",1:Dy)
l.header = c(lx.header, ly.header)
kl.header = c(k.header, l.header)
all.cuboids = do.call(rbind, lapply(fit$all, `[[`, "cuboids"))
all.tables = do.call(rbind, lapply(fit$all, `[[`, "tables"))
if (is.null(use.pval)) {
if (is.null(a$p.adjust.methods) & !a$correct) { use.pval = "pv" }
if (is.null(a$p.adjust.methods) & a$correct) { use.pval = "pv.correct" }
if (!is.null(a$p.adjust.methods)) { use.pval = a$p.adjust.methods[length(a$p.adjust.methods)] }
}
tree.struct = NULL
tree.struct = split(all.cuboids$child.of.test, all.cuboids$cuboid.idx)
# suppressMessages(library(png))
# suppressMessages(library(qgraph))
if (!is.null(x) & !is.null(y)) {
if (a$rank.transform) {
x.tran=apply(x,2,function(x) ecdf(x)(x)-1/n)
x.sort = apply(x,2,sort)
e_cdf.x = apply(x, 2, function(x) { x.sort=sort(x);
1:length(x.sort) / length(x.sort)})
y.tran=apply(y,2,function(y) ecdf(y)(y)-1/n)
e_cdf.y = apply(y, 2, function(y) { y.sort=sort(y);
1:length(y.sort) / length(y.sort)})
y.sort = apply(y,2,sort)
invECDF.x = function(x.tr, col)
{ sapply(col,function(c, xt, xs, ex) {ifelse(xt[c]==1,max(xs[,c]), xs[which(as.numeric(ex[,c]) >= as.numeric(xt[c]+1/n))[1],c])},xt=x.tr, xs=x.sort, ex=e_cdf.x) }
invECDF.y = function(y.tr, col)
{ sapply(col,function(c, yt, ys, ey) {ifelse(yt[c]==1,max(ys[,c]),ys[which(as.numeric(ey[,c]) >= as.numeric(yt[c]+1/n))[1],c])},yt=y.tr, ys=y.sort, ey=e_cdf.y) }
}
if (!a$rank.transform) {
scale01 = function(z){(1-10^(-7))*(z-min(z))/(max(z)-min(z))}
x.tran = apply(x,2,scale01)
y.tran = apply(y,2,scale01)
min.max.x = rbind(apply(x, 2, min),apply(x, 2, max))
invScale01.x=function(x.tr, col) {
sapply(col, function(c, xt, mmx) {ifelse(as.numeric(xt[c])==1,as.numeric(mmx[2,c]),
(as.numeric(mmx[2,c])-as.numeric(mmx[1,c]))*as.numeric(xt[c])/(1-10^(-7))+as.numeric(mmx[1,c]))},xt=x.tr, mmx=min.max.x) }
min.max.y = rbind(apply(y, 2, min),apply(y, 2, max))
invScale01.y=function(y.tr, col) {
sapply(col, function(c, yt, mmy) {ifelse(as.numeric(yt[c])==1,as.numeric(mmy[2,c]),
(as.numeric(mmy[2,c])-as.numeric(mmy[1,c]))*as.numeric(yt[c])/(1-10^(-7))+as.numeric(mmy[1,c]))},yt=y.tr, mmy=min.max.y) }
}
}
# Add network notation:
if (!show.all) {
ranked.tests = do.call(c, lapply(fit$all, `[[`, "rank.tests"))
} else {
ranked.tests = rep(TRUE,sum(!is.na(pvs$pv)))
}
nodes = all.tables[!is.na(pvs$pv), "test.idx"][ranked.tests]
n.nodes = length(nodes)
# Initialize adjacency matrix
adj = matrix(FALSE, nrow=n.nodes, ncol=n.nodes)
colnames(adj) = rownames(adj) = nodes
# Fill it!
for (s in 1:n.nodes) {
w = as.character(all.tables[all.tables[,"test.idx"]==nodes[s],"cuboid.idx"])
if (w>1) {
adj[as.character(tree.struct[[w]]),as.character(nodes[s])] = TRUE
}
}
# If data is given, plot the nodes, each png in low quality:
images=NULL
if (!is.null(x) & !is.null(y)) {
if (!requireNamespace("png", quietly = TRUE)) {
stop("Package \"png\" is required.",
call. = FALSE)
}
if (is.null(images.path)) images.path = tempdir()
W = all.cuboids[!duplicated(all.cuboids$cuboid.idx),]
rownames(W) = W$cuboid.idx
rownames(all.tables) = all.tables[,"test.idx"]
for (m in 1:n.nodes) {
w = as.character(all.tables[as.character(nodes[m]),"cuboid.idx"])
kx = unlist(W[w,kx.header])
ky = unlist(W[w,ky.header])
lx = unlist(W[w,lx.header])
ly = unlist(W[w,ly.header])
i = all.tables[nodes[m],"i"]
j = all.tables[nodes[m],"j"]
xl.cond = (lx[-i]-1)*2^(-kx[-i])
x.low.cond = matrix(rep(xl.cond, n), nrow=n, byrow=TRUE)
xh.cond = lx[-i]*2^(-kx[-i])
x.high.cond = matrix(rep(xh.cond, n), nrow=n, byrow=TRUE)
x.mask.cond = x.tran[,-i]>=x.low.cond & x.tran[,-i]<x.high.cond
yl.cond = (ly[-j]-1)*2^(-ky[-j])
y.low.cond = matrix(rep(yl.cond, n), nrow=n, byrow=TRUE)
yh.cond = ly[-j]*2^(-ky[-j])
y.high.cond = matrix(rep(yh.cond, n),nrow=n, byrow=TRUE)
y.mask.cond = y.tran[,-j]>=y.low.cond & y.tran[,-j]<y.high.cond
mask.cond = apply(x.mask.cond, 1, all) & apply(y.mask.cond, 1, all)
xl = (lx[i]-1)*2^(-kx[i])
xh = lx[i]*2^(-kx[i])
x.mask = x.tran[,i]>=xl & x.tran[,i]<xh
yl = (ly[j]-1)*2^(-ky[j])
yh = ly[j]*2^(-ky[j])
y.mask = y.tran[,j]>=yl & y.tran[,j]<yh
mask = as.logical(x.mask*y.mask*mask.cond)
x.loc = x[,i][mask]
y.loc = y[,j][mask]
png(file.path(images.path, paste0(node.name, m,".png")), width=80, height=60)
par(mar=(c(1,0.8,0,0)))
plot(x[,i],y[,j], col="grey", pch=".", xlab=NA, ylab=NA, xaxt='n', yaxt='n')
mtext(side=1, text=paste("x",i,sep=""), line=0.1, cex=0.7)
mtext(side=2, text=paste("y",j,sep=""), line=0.1, cex=0.7)
points(x[,i][mask.cond], y[,j][mask.cond], col="orange", pch=".")
points(x.loc, y.loc, col="red", pch=".")
max.x.tran = ((2*lx-0)*2^(-kx-1))
if (a$rank.transform) max.x = invECDF.x(max.x.tran, 1:Dx)
if (!a$rank.transform) max.x = invScale01.x(max.x.tran, 1:Dx)
min.x.tran = ((2*lx-2)*2^(-kx-1))
if (a$rank.transform) min.x = invECDF.x(min.x.tran, 1:Dx)
if (!a$rank.transform) min.x = invScale01.x(min.x.tran, 1:Dx)
max.y.tran = ((2*ly-0)*2^(-ky-1))
if (a$rank.transform) max.y = invECDF.y(max.y.tran, 1:Dy)
if (!a$rank.transform) max.y = invScale01.y(max.y.tran, 1:Dy)
min.y.tran = ((2*ly-2)*2^(-ky-1))
if (a$rank.transform) min.y = invECDF.y(min.y.tran, 1:Dy)
if (!a$rank.transform) min.y = invScale01.y(min.y.tran, 1:Dy)
x.i.min = min.x[i]
x.i.max = max.x[i]
y.j.min = min.y[j]
y.j.max = max.y[j]
abline(v=x.i.min, col="blue",lty=2)
abline(v=x.i.max, col="blue",lty=2)
abline(h=y.j.min, col="blue",lty=2)
abline(h=y.j.max, col="blue",lty=2)
segments(x0=x.i.min,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.max,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=(x.i.max+x.i.min)/2,y0=y.j.min,y1=y.j.max,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.min,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=(y.j.max+y.j.min)/2,col="blue")
segments(x0=x.i.min,x1=x.i.max, y0=y.j.max,col="blue")
dev.off()
}
images = paste0(images.path,"/node",1:n.nodes,".png")
}
nodes.pvs = pvs[!is.na(pvs$pv), use.pval][ranked.tests]
nodes.size = 1.5*(log(-log(nodes.pvs)+1)+.6)
if (max(nodes.size)>max.node.size) nodes.size = max.node.size*nodes.size/max(nodes.size)
if (max(nodes.size)!=min(nodes.size)) {
nodes.size = min.node.size+(max.node.size-min.node.size)*(nodes.size-min(nodes.size))/(max(nodes.size)-min(nodes.size))
}
if(is.null(filename)) filename=paste0("multiTree_",format(Sys.time(), "%Y-%m-%d_%H:%M:%S"))
if(is.null(filetype)) filetype = "pdf"
qg = qgraph::qgraph(input=adj, images=images, shape='square', vTrans=200, vsize=nodes.size,
border.width = 0, label.cex=.6, filename=filename, filetype=filetype)
return(invisible(list(qgraph.object=qg, qgraph.call="qgraph(input=adj, images=images, shape='square', vTrans=200, vsize=nodes.size, border.width=0, label.cex=.6, filename=filename, filetype=filetype)",
adj=adj, n.nodes = n.nodes, nodes.pvs=nodes.pvs, nodes.size=nodes.size,
images=images, pvs.attributes=a, filename=filename, filetype=filetype)))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.