Skip to content

Commit

Permalink
update plotEstimator() when including B_0
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhizuio
zhizuio committed Aug 5, 2024
1 parent 78e2b38 commit 6909650
Show file tree
Hide file tree
Showing 10 changed files with 718 additions and 228 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ Suggests: R.rsp, knitr, rmarkdown, BDgraph, data.table, plyr, scrime
LazyData: true
NeedsCompilation: yes
SystemRequirements: C++17
Packaged: 2024-08-05 09:11:35 UTC; zhiz
Packaged: 2024-08-05 13:10:39 UTC; zhiz
Author: Marco Banterle [aut],
Zhi Zhao [aut, cre],
Leonardo Bottolo [ctb],
Expand Down
1 change: 1 addition & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

* Fixed a bug about HRR model introduced in v2.2-0 by bringing back the control of `a_omega` and `b_omega` in function `BayesSUR()`
* Fixed a bug in function `SUR_Chain::stepWGibbs()` for extracting a submatrix
* Simplified the calculation in function `SUR_Chain::logPGamma`

# BayesSUR 2.2-1

Expand Down
238 changes: 140 additions & 98 deletions R/plotEstimator.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,15 +74,15 @@
#' }
#'
#' @export
plotEstimator <- function(x, estimator = NULL,
colorScale.gamma = grey((100:0) / 100),
colorScale.beta = c("blue", "white", "red"),
plotEstimator <- function(x, estimator = NULL,
colorScale.gamma = grey((100:0) / 100),
colorScale.beta = c("blue", "white", "red"),
legend.cex.axis = 1, name.responses = NA,
name.predictors = NA, xlab = "", ylab = "",
fig.tex = FALSE, output = "ParamEstimator",
header = "", header.cex = 2, tick = FALSE,
mgp = c(2.5, 1, 0), cex.main = 1.5,
title.beta = NA, title.gamma = NA, title.Gy = NA,
name.predictors = NA, xlab = "", ylab = "",
fig.tex = FALSE, output = "ParamEstimator",
header = "", header.cex = 2, tick = FALSE,
mgp = c(2.5, 1, 0), cex.main = 1.5,
title.beta = NA, title.gamma = NA, title.Gy = NA,
beta.type = "marginal", Pmax = 0, ...) {
if (!inherits(x, "BayesSUR")) {
stop("Use only with a \"BayesSUR\" object")
Expand All @@ -91,28 +91,33 @@ plotEstimator <- function(x, estimator = NULL,
stop("Please specify correct argument estimator!")
}

beta_hat <- getEstimator(x, estimator = "beta",
beta.type = beta.type, Pmax = Pmax)
beta_hat <- getEstimator(x,
estimator = "beta",
beta.type = beta.type, Pmax = Pmax
)
x$output[-1] <- paste0(x$output$outFilePath, x$output[-1])
gamma_hat <- as.matrix(read.table(x$output$gamma))
response_idx <- seq_len(ncol(gamma_hat))
predictor_idx <- seq_len(nrow(gamma_hat))
nonpen <- nrow(beta_hat) - nrow(gamma_hat)
predictor_idx0 <- c(seq_len(nonpen), nonpen + predictor_idx)
if (nonpen > 0) {
rownames(beta_hat) <- c(colnames(read.table(x$output$X0, header = TRUE)),
colnames(read.table(x$output$X, header = TRUE)))
rownames(beta_hat) <- c(
colnames(read.table(x$output$X0, header = TRUE)),
colnames(read.table(x$output$X, header = TRUE))
)
} else {
rownames(beta_hat) <- colnames(read.table(x$output$X, header = TRUE))
}
colnames(beta_hat) <- colnames(read.table(x$output$Y, header = TRUE))

covariancePrior <- x$input$covariancePrior
if ((covariancePrior == "HIW") && ("Gy" %in% estimator)) {
Gy_hat <- as.matrix(read.table(x$output$Gy))
colnames(Gy_hat) <- rownames(Gy_hat) <-
colnames(Gy_hat) <- rownames(Gy_hat) <-
colnames(read.table(x$output$Y, header = TRUE))
}

## BUG TO BE FIXED!!!
# specify the labels of axes
if (is.na(name.responses)[1]) name.responses <- response_idx
Expand All @@ -122,59 +127,67 @@ plotEstimator <- function(x, estimator = NULL,
stop("The length of the given response names are not consistent with the data!")
}
}

if (is.na(name.predictors)) name.predictors <- predictor_idx
if (name.predictors[1] == "auto") name.predictors <- rownames(beta_hat)
if (is.character(name.predictors)) {
if (length(name.predictors) != nrow(beta_hat)) {
stop("The length of the given predictor names are not consistent with the data!")
}
}

opar <- par(no.readonly = TRUE)
on.exit(par(opar))
if (!fig.tex) {
par(mar = c(6, 6, 5.1, 4.1))
par(mfrow = c(1, sum(estimator %in% c("beta", "gamma", "Gy"))))

if ("beta" %in% estimator) {
# floor(100*constant)+100-1 colours that your want in the legend bar which has the white middle color
if (stats::sd(beta_hat) == 0) {
colorbar <- colorRampPalette(c(colorScale.beta[2], colorScale.beta[2]))(1000)
} else {
colorbar <- c(colorRampPalette(c(colorScale.beta[1], colorScale.beta[2]))(
floor(1000 / (-(max(beta_hat) - min(beta_hat)) / min(beta_hat) - 1))),
colorRampPalette(c(colorScale.beta[2], colorScale.beta[3]))(1000)[-1])
colorbar <- c(
colorRampPalette(c(colorScale.beta[1], colorScale.beta[2]))(
floor(1000 / (-(max(beta_hat) - min(beta_hat)) / min(beta_hat) - 1))),
colorRampPalette(c(colorScale.beta[2], colorScale.beta[3]))(1000)[-1]
)
}
if (is.na(title.beta)) title.beta <- expression(hat(bold(B)))

image(
z = beta_hat, x = predictor_idx, y = response_idx,
z = beta_hat, x = predictor_idx0, y = response_idx,
col = colorbar, mgp = mgp,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = xlab, ylab = ylab, main = title.beta,
cex.main = cex.main, cex.lab = 1.5, ...)
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = xlab, ylab = ylab, main = title.beta,
cex.main = cex.main, cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorbar,
zlim = c(min(beta_hat), max(beta_hat)),
legend.cex.axis = legend.cex.axis)
vertical.image.legend(
color = colorbar,
zlim = c(min(beta_hat), max(beta_hat)),
legend.cex.axis = legend.cex.axis
)
if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = response_idx, labels = name.responses, tick = tick)
axis(1, at = predictor_idx, labels = name.predictors, tick = tick)
axis(1, at = predictor_idx0, labels = name.predictors, tick = tick)
}
}
if ("gamma" %in% estimator) {
if (is.na(title.gamma)) title.gamma <- expression(hat(bold(Gamma)))
image(
z = gamma_hat, x = predictor_idx, y = response_idx,
z = gamma_hat, x = predictor_idx, y = response_idx,
col = colorScale.gamma, mgp = mgp,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = xlab, ylab = ylab, main = title.gamma,
cex.main = cex.main, cex.lab = 1.5, ...)
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = xlab, ylab = ylab, main = title.gamma,
cex.main = cex.main, cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorScale.gamma, zlim = c(0, 1),
legend.cex.axis = legend.cex.axis)
vertical.image.legend(
color = colorScale.gamma, zlim = c(0, 1),
legend.cex.axis = legend.cex.axis
)
if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = response_idx, labels = name.responses, tick = tick)
Expand All @@ -184,20 +197,23 @@ plotEstimator <- function(x, estimator = NULL,
axis(1, at = predictor_idx, labels = name.predictors, tick = tick)
}
}

if ("Gy" %in% estimator) {
if (is.na(title.Gy)) title.Gy <- "Estimated graph of responses"
image(
z = Gy_hat + diag(ncol(Gy_hat)),
x = response_idx, y = response_idx,
z = Gy_hat + diag(ncol(Gy_hat)),
x = response_idx, y = response_idx,
col = colorScale.gamma, mgp = mgp,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = ylab, ylab = ylab, main = title.Gy,
cex.main = cex.main, cex.lab = 1.5, ...)
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = ylab, ylab = ylab, main = title.Gy,
cex.main = cex.main, cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorScale.gamma,
zlim = c(min(Gy_hat), max(Gy_hat)),
legend.cex.axis = legend.cex.axis)
vertical.image.legend(
color = colorScale.gamma,
zlim = c(min(Gy_hat), max(Gy_hat)),
legend.cex.axis = legend.cex.axis
)
if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = response_idx, labels = name.responses, tick = tick)
Expand All @@ -206,101 +222,127 @@ plotEstimator <- function(x, estimator = NULL,
}
title(paste0("\n", header), cex.main = header.cex, outer = TRUE)
} else {
options(tikzMetricPackages = c("\\usepackage{amsmath}", "\\usepackage{bm}",
"\\usetikzlibrary{calc}"))
options(tikzMetricPackages = c(
"\\usepackage{amsmath}", "\\usepackage{bm}",
"\\usetikzlibrary{calc}"
))
tikz(paste0(output, ".tex"),
width = 3.6 * sum(estimator %in% c("beta", "gamma", "Gy")),
height = 4,
standAlone = TRUE,
packages = c("\\usepackage{tikz}",
"\\usepackage{amsmath}",
"\\usepackage{bm}",
"\\usepackage[active,tightpage,psfixbb]{preview}",
"\\PreviewEnvironment{pgfpicture}"))

width = 3.6 * sum(estimator %in% c("beta", "gamma", "Gy")),
height = 4,
standAlone = TRUE,
packages = c(
"\\usepackage{tikz}",
"\\usepackage{amsmath}",
"\\usepackage{bm}",
"\\usepackage[active,tightpage,psfixbb]{preview}",
"\\PreviewEnvironment{pgfpicture}"
)
)

par(mfrow = c(1, sum(estimator %in% c("beta", "gamma", "Gy"))))
par(mar = c(6, 6, 4, 4) + 0.1)
if ("beta" %in% estimator) {
# floor(100*constant)+100-1 colours that your want in the legend bar which has the white middle color
colorbar <-
c(colorRampPalette(c(colorScale.beta[1], colorScale.beta[2]))(
floor(1000 / (-(max(beta_hat) - min(beta_hat)) / min(beta_hat) - 1))),
colorRampPalette(c(colorScale.beta[2], colorScale.beta[3]))(1000)[-1])
if (is.na(title.beta))
colorbar <-
c(
colorRampPalette(c(colorScale.beta[1], colorScale.beta[2]))(
floor(1000 / (-(max(beta_hat) - min(beta_hat)) / min(beta_hat) - 1))),
colorRampPalette(c(colorScale.beta[2], colorScale.beta[3]))(1000)[-1]
)
if (is.na(title.beta)) {
title.beta <- paste("Estimator", "$\\hat{\\bm{B}}$")

}

image(
z = beta_hat, x = predictor_idx, y = response_idx, col = colorbar,
z = beta_hat, x = predictor_idx0, y = response_idx, col = colorbar,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE), mgp = mgp,
xlab = xlab, ylab = ylab, main = title.beta,
cex.main = cex.main, cex.lab = 1.5, ...)
xlab = xlab, ylab = ylab, main = title.beta,
cex.main = cex.main, cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorbar,
zlim = c(min(beta_hat), max(beta_hat)),
legend.cex.axis = legend.cex.axis)
vertical.image.legend(
color = colorbar,
zlim = c(min(beta_hat), max(beta_hat)),
legend.cex.axis = legend.cex.axis
)
if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = response_idx, labels = name.responses, tick = tick)
# opar <- par(cex.axis=1)
axis(1, at = predictor_idx, labels = name.predictors, tick = tick)
axis(1, at = predictor_idx0, labels = name.predictors, tick = tick)
}
}
if ("gamma" %in% estimator) {
if (is.na(title.gamma))
if (is.na(title.gamma)) {
title.gamma <- paste("Estimator", "$\\hat{\\mathbf{\\Gamma}}$")
}
image(
z = gamma_hat, x = predictor_idx, y = seq_len(ncol(gamma_hat)),
col = colorScale.gamma,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
z = gamma_hat, x = predictor_idx, y = seq_len(ncol(gamma_hat)),
col = colorScale.gamma,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
xlab = xlab, ylab = ylab, main = title.gamma, mgp = mgp,
cex.main = cex.main, cex.lab = 1.5, ...)
cex.main = cex.main, cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorScale.gamma,
zlim = c(0, 1), legend.cex.axis = legend.cex.axis)
vertical.image.legend(
color = colorScale.gamma,
zlim = c(0, 1), legend.cex.axis = legend.cex.axis
)
if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = seq_len(ncol(gamma_hat)),
labels = name.responses, tick = tick)
axis(2,
at = seq_len(ncol(gamma_hat)),
labels = name.responses, tick = tick
)
if (nonpen > 0) {
name.predictors <- name.predictors[-c(1:nonpen)]
}
axis(1, at = predictor_idx, labels = name.predictors, tick = tick)
}
}

if ("Gy" %in% estimator) {
if (is.na(title.Gy))
if (is.na(title.Gy)) {
title.Gy <- paste("Estimator", "$\\hat{\\mathcal{G}}$")
}
image(
z = Gy_hat + diag(ncol(Gy_hat)),
x = seq_len(nrow(Gy_hat)),
y = seq_len(nrow(Gy_hat)),
col = colorScale.gamma,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
z = Gy_hat + diag(ncol(Gy_hat)),
x = seq_len(nrow(Gy_hat)),
y = seq_len(nrow(Gy_hat)),
col = colorScale.gamma,
axes = ifelse(is.na(name.responses)[1], TRUE, FALSE),
mgp = mgp,
xlab = ylab,
ylab = ylab,
main = title.Gy,
cex.main = cex.main,
cex.lab = 1.5, ...)
xlab = ylab,
ylab = ylab,
main = title.Gy,
cex.main = cex.main,
cex.lab = 1.5, ...
)
box()
vertical.image.legend(color = colorScale.gamma,
zlim = c(min(Gy_hat), max(Gy_hat)),
legend.cex.axis = legend.cex.axis)

vertical.image.legend(
color = colorScale.gamma,
zlim = c(min(Gy_hat), max(Gy_hat)),
legend.cex.axis = legend.cex.axis
)

if (!is.na(name.responses)[1]) {
par(las = 2, cex.axis = 1)
axis(2, at = seq_len(ncol(Gy_hat)),
labels = name.responses, tick = tick)
axis(1, at = seq_len(nrow(Gy_hat)),
labels = name.responses, tick = tick)
axis(2,
at = seq_len(ncol(Gy_hat)),
labels = name.responses, tick = tick
)
axis(1,
at = seq_len(nrow(Gy_hat)),
labels = name.responses, tick = tick
)
}
}
title(paste0("\n", header), cex.main = header.cex, outer = TRUE)
dev.off()
tools::texi2pdf(paste0(output, ".tex"))
}
}

# the function vertical.image.legend() is adapted from the R package "aqfig"
vertical.image.legend <- function(zlim, color, legend.cex.axis = 1) {
starting.par.settings <- par(no.readonly = TRUE)
Expand Down
Loading

0 comments on commit 6909650

Please sign in to comment.