Nothing
#' Plot elicited data, fitted marginals or model output
#'
#' @param Z list object that contains matrix \code{theta} of elicitations,
#' character \code{link} and character \code{target} as initialised by
#' \code{\link{designLink}} and updated by \code{\link{elicitPt}}
#' @param X design matrix (can be \code{NULL}, unless modelled output is
#' requested)
#' @param design.pt single integer that denotes design point of interest
#' @param elicited.fractiles logical, plot vertical lines for elicited
#' fractiles?
#' @param fitted.fractiles logical, plot vertical lines for fitted conditional
#' mean prior fractiles for this design point? Alternatively, a numeric vector of arbitrary fractiles to be
#' plotted from the fitted elicitation distribution. If \code{TRUE} then the
#' fractiles corresponding to the median, upper and lower level central CI
#' are plotted
#' @param fitted.curve, logical plot fitted conditional mean prior density for this design point?
#' @param CI.prob numeric scalar, locally specified probability assigned to the
#' elicited central credible interval of the current design point. Defaults to
#' \code{NULL} in which case the global value initially assigned by
#' \code{\link{designLink}} or as updated by \code{\link{elicitPt}} is used
#' @param estimated.probs numeric vector of values for which estimated
#' probabilities are to be estimated from the fitted elicitation
#' distribution for the target theta. Default is \code{NULL}.
#' The result is output to the console.
#' @param modelled.fractiles logical, plot vertical lines for modelled
#' fractiles from the conditional mean prior distribution fit to
#' all design points? This option requires a design matrix \code{X} of full column rank.
#' @param modelled.curve logical, plot modelled conditional mean prior density for
#' the entire model? This option requires a design matrix \code{X} of full column rank.
#' @param cumul.prob.bounds numeric vector of length two, giving plot bounds by
#' cumulative probability. This argument is ignored if there is not enough data
#' to fit a parametric distribution or if \code{theta.bounds} is not \code{NULL}
#' @param theta.bounds numeric vector giving support of response for plotting
#' purposes (can be \code{NULL}). This will overwrite \code{cumul.prob.bounds},
#' if applicable
#' @param ylim.max numeric maximum value of y-axis (can be \code{NULL})
#' @param xlog logical log x-axis
#' @param design.table logical include design dataframe, elicited fractiles and
#' modelled or fitted fractiles
#' @param n.pts numeric giving number of point to evalate density curve (if
#' plotted)
#' @return a plot to the current device. See \code{dev.cur()} to check.
#' @examples
#' # design matrix: two scenarios
#' X <- matrix(c(1, 1, 0, 1), nrow = 2)
#' rownames(X) <- c("scenario1", "scenario2")
#' colnames(X) <- c("covariate1", "covariate2")
#'
#' # logit link
#' # central credible intervals with probability = 1/2
#' Z <- designLink(design = X, link = "logit", CI.prob = 0.5)
#'
#' # 1st design point
#' # no elicited fractiles
#' indirect::plotDesignPoint(Z, design.pt = 1)
#' # elicited median
#' Z <- indirect::elicitPt(Z, design.pt = 1,
#' lower.CI.bound = NA,
#' median = 0.4,
#' upper.CI.bound = NA,
#' CI.prob = NULL)
#' indirect::plotDesignPoint(Z, design.pt = 1,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1))
#' # lower and upper quartiles and median
#' Z <- indirect::elicitPt(Z, design.pt = 1,
#' lower.CI.bound = 0.2,
#' median = 0.4,
#' upper.CI.bound = 0.6,
#' comment = "Completed.")
#' indirect::plotDesignPoint(Z, design.pt = 1,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1),
#' fitted.fractiles = TRUE, fitted.curve = TRUE)
#' indirect::plotDesignPoint(Z, design.pt = 1,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1),
#' fitted.fractiles = c(1/10, 1/4, 1/2, 3/4, 9/10),
#' fitted.curve = TRUE)
#'
#' # second design point
#' # central credible intervals with probability = 1/3
#' # elicit upper and lower tertiles
#' Z <- elicitPt(Z, design.pt = 2,
#' lower.CI.bound = 0.1,
#' upper.CI.bound = 0.3,
#' CI.prob = 1/3,
#' comment = "Switched to tertiles.")
#' indirect::plotDesignPoint(Z, design.pt = 2,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1))
#' indirect::plotDesignPoint(Z, design.pt = 2,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1),
#' fitted.fractiles = TRUE, fitted.curve = TRUE)
#' indirect::plotDesignPoint(Z, design.pt = 2,
#' elicited.fractiles = TRUE, theta.bounds = c(0, 1),
#' fitted.fractiles = c(1/10, 1/3, 1/2, 2/3, 9/10),
#' fitted.curve = TRUE)
plotDesignPoint <- function(Z, X = NULL, design.pt = NULL,
elicited.fractiles = TRUE, fitted.fractiles = FALSE, fitted.curve = FALSE,
CI.prob = NULL, estimated.probs = NULL,
modelled.fractiles = FALSE, modelled.curve = FALSE,
cumul.prob.bounds = c(0.05, 0.95), theta.bounds = NULL, ylim.max = NULL, xlog = FALSE,
design.table = TRUE, n.pts = 101) {
if (is.null(design.pt)) stop("specify design point")
theta <- Z$theta
link <- Z$link
design <- Z$design
fit.method <- Z$fit.method
# elicited data could be missing but design point must be specified
dat <- theta[design.pt, 1:(ncol(theta) - 1)]
# find fractiles
ci.prob <- theta[design.pt, "CI_prob"]
cumul.probs <- c((1 - ci.prob)/2, 1/2, 1 - (1 - ci.prob)/2)
# if specific fitted fractiles aren't requested then use fractiles
# for lower and upper central CI and median
if (is.logical(fitted.fractiles)) {
cumul.probs.fit <- cumul.probs
} else if (is.numeric(fitted.fractiles)) {
cumul.probs.fit <- fitted.fractiles
fitted.fractiles <- TRUE
} else {
stop("fitted.fractiles should be either a logical value or numeric")
}
if (fitted.fractiles || fitted.curve || modelled.fractiles || modelled.curve) {
x.axis.mf <- matrix(NA, nrow = 2, ncol = 2,
dimnames = list(c("min", "max"), c("fitted", "modelled")))
# transform elicitations to moments
mV.ls <- mV(Z, fit.method = fit.method)
fitted.mean <- mV.ls$m[design.pt]
fitted.sd <- sqrt(mV.ls$V[design.pt, design.pt])
if (!is.null(cumul.prob.bounds)) {
x.axis.mf[ , "fitted"] <- qnorm(cumul.prob.bounds, fitted.mean, fitted.sd)
} else {
x.axis.mf[ , "fitted"] <- qnorm(cumul.probs.fit, fitted.mean, fitted.sd)
}
if (modelled.fractiles || modelled.curve) {
if (is.null(X)) stop("design matrix X required")
# estimate mu and Sigma for given elicitations
muSigma.ls <- muSigma(Z, X, fit.method = fit.method, wls.method = "default")
# predictions for g(theta(design.pt))
modelled.mean <- X[design.pt, ]%*%muSigma.ls$mu
modelled.sd <- sqrt(X[design.pt, , drop = FALSE]%*%muSigma.ls$Sigma%*%t(X[design.pt, , drop = FALSE]))
if (!is.null(cumul.prob.bounds)) {
x.axis.mf[ , "modelled"] <- qnorm(cumul.prob.bounds, modelled.mean, modelled.sd)
} else {
x.axis.mf[ , "modelled"] <- qnorm(cumul.probs.fit, modelled.mean, modelled.sd)
}
}
x.axis.mf <- c(min(x.axis.mf["min", ], na.rm = TRUE), max(x.axis.mf["max", ], na.rm = TRUE))
# limits for plot's x axis
# theta.bounds takes priority over cumul.prob.bounds if not NULL
x.lim.cumul.prob <- FALSE
if (!is.null(theta.bounds)) {
x.axis.lim <- theta.bounds
} else if (!is.null(x.axis.mf)) {
x.axis.lim <- x.axis.mf
x.lim.cumul.prob <- TRUE
} else {
stop("not enough information provided for plotting fitted values")
}
# transform by link function
switch(link,
log = {
if (x.lim.cumul.prob) x.axis.lim <- exp(x.axis.lim)
if (fitted.fractiles) dat.fit <- exp(qnorm(cumul.probs.fit, fitted.mean, fitted.sd))
if (modelled.fractiles) {
dat.mod <- exp(qnorm(cumul.probs.fit, modelled.mean, modelled.sd))
}
if (fitted.curve || modelled.curve) {
if (xlog) {
x <- seq(from = log(x.axis.lim[1]), to = log(x.axis.lim[2]), length.out = n.pts)
x <- exp(x)
} else {
x <- seq(from = x.axis.lim[1], to = x.axis.lim[2], length.out = n.pts)
}
if (fitted.curve) dens.fit <- dlnorm(x, fitted.mean, fitted.sd)
if (modelled.curve) {
dens.mod <- dlnorm(x, modelled.mean, modelled.sd)
}
}
},
cloglog = {
if (x.lim.cumul.prob) x.axis.lim <- 1 - exp(-exp(x.axis.lim))
if (fitted.fractiles) dat.fit <- 1 - exp(-exp(qnorm(cumul.probs.fit, fitted.mean, fitted.sd)))
if (modelled.fractiles) {
dat.mod <- 1 - exp(-exp(qnorm(cumul.probs.fit, modelled.mean, modelled.sd)))
}
if (fitted.curve || modelled.curve) {
if (xlog) {
x <- seq(from = log(x.axis.lim[1]), to = log(x.axis.lim[2]), length.out = n.pts)
x <- exp(x)
} else {
x <- seq(from = x.axis.lim[1], to = x.axis.lim[2], length.out = n.pts)
}
if (fitted.curve) dens.fit <- dGompertzNorm(x, fitted.mean, fitted.sd)
if (modelled.curve) {
dens.mod <- dGompertzNorm(x, modelled.mean, modelled.sd)
}
}
},
logit = {
if (x.lim.cumul.prob) x.axis.lim <- exp(x.axis.lim)/(1 + exp(x.axis.lim))
if (fitted.fractiles) dat.fit <- exp(qnorm(cumul.probs.fit, fitted.mean, fitted.sd))/(1 + exp(qnorm(cumul.probs.fit, fitted.mean, fitted.sd)))
if (modelled.fractiles) {
dat.mod <- exp(qnorm(cumul.probs.fit, modelled.mean, modelled.sd))/(1 + exp(qnorm(cumul.probs.fit, modelled.mean, modelled.sd)))
}
if (fitted.curve || modelled.curve) {
if (xlog) {
x <- seq(from = log(x.axis.lim[1]), to = log(x.axis.lim[2]), length.out = n.pts)
x <- exp(x)
} else {
x <- seq(from = x.axis.lim[1], to = x.axis.lim[2], length.out = n.pts)
}
if (fitted.curve) dens.fit <- dLogitNorm(x, fitted.mean, fitted.sd)
if (modelled.curve) {
dens.mod <- dLogitNorm(x, modelled.mean, modelled.sd)
}
}
},
identity = {
if (fitted.fractiles) dat.fit <- qnorm(cumul.probs.fit, fitted.mean, fitted.sd)
if (modelled.fractiles) {
dat.mod <- qnorm(cumul.probs.fit, modelled.mean, modelled.sd)
}
if (fitted.curve || modelled.curve) {
if (xlog) {
x <- seq(from = log(x.axis.lim[1]), to = log(x.axis.lim[2]), length.out = n.pts)
x <- exp(x)
} else {
x <- seq(from = x.axis.lim[1], to = x.axis.lim[2], length.out = n.pts)
}
if (fitted.curve) dens.fit <- dnorm(x, fitted.mean, fitted.sd)
if (modelled.curve) {
dens.mod <- dnorm(x, modelled.mean, modelled.sd)
}
}
},
stop("link unspecified")
)
if (!fitted.fractiles) dat.fit <- rep(NA, 3)
if (is.null(ylim.max)) {
if (!modelled.curve && fitted.curve) max.y <- max(dens.fit, na.rm = TRUE)
if (modelled.curve && !fitted.curve) max.y <- max(dens.mod, na.rm = TRUE)
if (modelled.curve && fitted.curve) max.y <- max(dens.fit, dens.mod, na.rm = TRUE)
if (!modelled.curve && !fitted.curve) max.y <- 1 # arbitrary value
y.axis.lim <- c(0, max.y)
} else {
y.axis.lim <- c(0, ylim.max)
}
} else {
# elicited data only - no fitting
# info for design table, if requested: no additional fitted fractiles
cumul.probs.fit <- cumul.probs
dat.fit <- rep(NA, 3)
# plotting limits
if (!is.null(theta.bounds)) {
x.axis.lim <- theta.bounds
} else if (all(!is.na(theta[design.pt, ]))) {
x.axis.lim <- c(theta[design.pt, c("lower", "upper")])
x.axis.lim[1] <- x.axis.lim[1] - x.axis.lim[1]/10
x.axis.lim[2] <- x.axis.lim[2] + x.axis.lim[2]/10
} else {
x.axis.lim <- c(-10^6, 10^6)
}
if (is.null(ylim.max)) {
y.axis.lim <- c(0, 1)
} else {
y.axis.lim <- c(0, ylim.max)
}
}
# create plot
if (design.table) layout(matrix(c(1, 1, 2, 3), nrow = 2, ncol = 2, byrow = TRUE))
par(mar = c(4, 4, 2, 0) + 0.1)
if (xlog) {
xlog <- "x"
} else {
xlog <- ""
}
plot(x.axis.lim, y.axis.lim, type = "n", xlab = "", ylab = "Density", main = Z$target.caption, log = xlog, axes = FALSE)
axis(1, cex.axis = 2)
axis(2)
box()
if (elicited.fractiles) {
abline(v = dat, lty = 2)
for (i in 1:3) {
f <- MASS::fractions(cumul.probs[i])
f.num <- strsplit(as.character(f), split = "/")[[1]][1]
f.den <- strsplit(as.character(f), split = "/")[[1]][2]
text(dat[i], 0.98*max(y.axis.lim), substitute(f[z/zz], list(z = f.num, zz = f.den)), pos = 4)
}
}
if (fitted.curve) lines(x, dens.fit, col = "blue")
if (fitted.fractiles) {
abline(v = dat.fit, lty = 2, col = "blue")
for (i in 1:length(dat.fit)) {
f <- MASS::fractions(cumul.probs.fit[i])
f.num <- strsplit(as.character(f), split = "/")[[1]][1]
f.den <- strsplit(as.character(f), split = "/")[[1]][2]
text(dat.fit[i], 0, substitute(f[z/zz], list(z = f.num, zz = f.den)), pos = 4, col = "blue")
}
}
if (modelled.curve) lines(x, dens.mod, col = "orange")
if (modelled.fractiles) {
abline(v = dat.mod, lty = 2, col = "orange")
for (i in 1:length(dat.mod)) {
f <- MASS::fractions(cumul.probs.fit[i])
f.num <- strsplit(as.character(f), split = "/")[[1]][1]
f.den <- strsplit(as.character(f), split = "/")[[1]][2]
text(dat.mod[i], 0, substitute(f[z/zz], list(z = f.num, zz = f.den)), pos = 4, col = "orange")
}
}
# add design table if specified
if (design.table) {
designpoint <- matrix(signif(design[design.pt, ], 4), nrow = length(design[design.pt, ]), ncol = 1,
dimnames = list(Covariate = colnames(design), "Value"))
gplots::textplot(designpoint, valign = "top", mar = c(1, 1, 2, 1))
title(paste("Scenario:", row.names(design)[design.pt]))
box("figure")
cprobs <- c(cumul.probs, cumul.probs.fit)
cprobs.fracs <- as.character(MASS::fractions(cprobs))
all.cumul.probs <- as.character(MASS::fractions(sort(cprobs[!duplicated(cprobs.fracs)])))
elicit.ind <- match(as.character(MASS::fractions(cumul.probs)), all.cumul.probs)
fitted.ind <- match(as.character(MASS::fractions(cumul.probs.fit)), all.cumul.probs)
fractiles <- matrix(NA, nrow = length(all.cumul.probs), ncol = 2)
fractiles[elicit.ind, 1] <- signif(dat, 3)
if (modelled.fractiles) {
fractiles[fitted.ind, 2] <- signif(dat.mod, 3)
dimnames(fractiles) <- list(Fractiles = all.cumul.probs, c("Elicited", "Modelled"))
} else {
fractiles[fitted.ind, 2] <- signif(dat.fit, 3)
dimnames(fractiles) <- list(Fractiles = all.cumul.probs, c("Elicited", "Fitted"))
}
gplots::textplot(fractiles, valign = "top", mar = c(1, 1, 2, 1))
title("Fractiles")
box("figure")
}
if (!is.null(estimated.probs)) print(pdist(estimated.probs, Z, design.pt, fit.method = fit.method))
}
#' density for Gompertz transformed univariate Gaussian
#'
#' @param x, numeric real
#' @param mu, numeric real
#' @param sigma, numeric real positive
#' @return tranformed density on support (0, 1)
#' @examples
#' mu <- -1
#' sigma <- 1
#' z <- rnorm(10000, mu, sigma)
#' hist(1 - exp(-exp(z)), freq = FALSE)
#' curve(dGompertzNorm(x, mu = mu, sigma = sigma), col = 'red', add = TRUE, from = 0.01, to = 0.99)
#' integrate(function(x) dGompertzNorm(x, mu = mu, sigma = sigma), lower = 0, upper = 1) # equals 1
dGompertzNorm <- function(x, mu, sigma) {
dnorm(log(-log(1 - x)), mu, sigma)/((x - 1)*log(1 - x))
}
#' density for logit transformed univariate Gaussian
#'
#' @param x, numeric real
#' @param mu, numeric real
#' @param sigma, numeric real positive
#' @return tranformed density on support (0, 1)
#' @examples
#' mu <- -1
#' sigma <- 1
#' z <- rnorm(10000, mu, sigma)
#' hist(exp(z)/(1 + exp(z)), freq = FALSE)
#' curve(dLogitNorm(x, mu = mu, sigma = sigma), col = 'red', add = TRUE, from = 0.01, to = 0.99)
#' integrate(function(x) dLogitNorm(x, mu = mu, sigma = sigma), lower = 0, upper = 1) # equals 1
dLogitNorm <- function(x, mu, sigma) {
dnorm(log(x/(1 - x)), mu, sigma)/(x*(1 - x))
}
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.