Nothing
#' Calculate Expected Value of Sample Information (EVSI)
#'
#' @inheritParams calc_evppi
#' @param n additional sample size
#' @param n0 initial sample size
#' @param n_by_param if \code{TRUE}, each parameter in the metamodel can have a unique
#' initial and additional sample size. \code{n} and \code{n0} must be numerical
#' vectors of equal length to \code{params}, with each value corresponding to the
#' initial and additional sample sizes for each parameter in the metamodel. By default,
#' \code{n_by_param = FALSE}, and each value of \code{n} and \code{n0} is shared by
#' each parameter in the model. When \code{n_by_param = FALSE}, \code{n0} must be a single
#' numeric value, and \code{n} must be a numerical vector of additional sample sizes for
#' which EVSI is calculated from the metamodel.
#' @param progress \code{TRUE} or \code{FALSE} for whether or not function progress
#' should be displayed in console.
#'
#' @return A list containing 1) a data.frame with WTP thresholds, new prospective sample sizes
#' (if \code{n_by_param == FALSE}), and corresponding EVSIs
#' for the selected parameters and 2) a list of metamodels used to estimate EVSI for each
#' strategy at each willingness to pay threshold.
#'
#' @export
calc_evsi <- function(psa,
wtp,
params = NULL,
outcome = c("nhb", "nmb"),
k = -1,
n = 100,
n0 = 10,
n_by_param = FALSE,
pop = 1,
progress = TRUE) {
# define parameter values and make sure they correspond to a valid option
outcome <- match.arg(outcome)
# adjust outcome type
outcome <- paste0(outcome, "_loss_voi")
# number of wtp thresholds
n_wtps <- length(wtp)
# number of new sample sizes
n_n <- length(n)
if (n_by_param == TRUE) {
# vector to store evsi
evsi <- rep(0, n_wtps)
} else {
# matrix to store evsi
evsi <- matrix(rep(0, n_wtps * n_n), ncol = n_n)
}
mms_ls <- vector(mode = "list", length = n_wtps)
# calculate evppi at each wtp and new sample size
for (l in seq_len(n_wtps)) {
# run the metamodels
mms <- metamodel(analysis = "multiway",
psa = psa,
params = params,
outcome = outcome,
wtp = wtp[l],
type = "gam",
poly.order = 2,
k = k)
mms_ls[[l]] <- mms
if (progress == TRUE) {
if (l / (n_wtps / 10) == round(l / (n_wtps / 10), 0)) { # display progress every 10%
cat("\r", paste(l / n_wtps * 100, "% done", sep = " "))
}
}
if (n_by_param == TRUE) {
# predict from the regression models
predicted_loss_list <- lapply(mms$mods, function(m) predict_ga(m, n, n0))
# bind the columns to get a dataframe
predicted_loss_df <- bind_cols(predicted_loss_list)
# calculate the evsi as the average of the row maxima
row_maxes <- apply(predicted_loss_df, 1, max)
evsi[l] <- mean(row_maxes) * pop
} else {
for (i in seq_len(n_n)) {
# predict from the regression models
predicted_loss_list <- lapply(mms$mods, function(m) predict_ga(m, n[i], n0))
# bind the columns to get a dataframe
predicted_loss_df <- bind_cols(predicted_loss_list)
# calculate the evsi as the average of the row maxima
row_maxes <- apply(predicted_loss_df, 1, max)
evsi[l, i] <- mean(row_maxes) * pop
}
}
}
if (n_by_param == TRUE) {
df_evsi <- data.frame("WTP" = wtp, "EVSI" = evsi)
} else {
# data.frame to store EVPPI for each WTP threshold
df_evsi <- data.frame("WTP" = rep(wtp, n_n),
"n" = rep(n, each = n_wtps),
"EVSI" = c(evsi))
}
evsi_ls <- list(df_evsi = df_evsi, metamodel_ls = mms_ls)
class(evsi_ls) <- c("evsi")
return(evsi_ls)
}
#' Function to compute the preposterior for each of the
#' basis functions of the GAM model.
#'
#' @keywords internal
#'
#' @param object gam object
#' @param n scalar or vector of new sample size to compute evsi on
#' @param n0 scalar or vector of effective prior sample size
#' @return the conditional loss associated with a gam, effective prior sample size,
#' and new sample size
#' @importFrom stats coef
predict_ga <- function(object, n, n0) {
# Name of parameters
param_names <- colnames(object$model)
# Create dataframe with parameter values
param_vals <- data.frame(object$model[, -1])
# Name columns of dataframe
colnames(param_vals) <- param_names[-1]
# Number of parameters
n_params <- ncol(param_vals)
# Sanity checks
if (!(length(n) == 1 | length(n) == n_params)) {
stop("Variable 'n' should be either a scalar or a vector
the same size as the number of parameters")
}
if (!(length(n0) == 1 | length(n0) == n_params)) {
stop("Variable 'n0' should be either a scalar or a vector
the same size as the number of parameters")
}
# Make n & n0 consistent with the number of parameters
if (length(n) == 1) {
n <- rep(n, n_params)
}
if (length(n0) == 1) {
n0 <- rep(n0, n_params)
}
# Compute variance reduction factor
vrf <- sqrt(n / (n + n0))
# Number of smoothers
n_smooth <- length(object$smooth)
# Number of total basis functions
n_col_x <- length(object$coefficients)
# Number of observations
n_row_x <- nrow(object$model)
# Initialize matrix for preposterior of total basis functions
x <- matrix(NA, n_row_x, n_col_x)
x[, 1] <- 1
for (k in 1:n_smooth) {
klab <- substr(object$smooth[[k]]$label, 1, 1)
if (klab == "s") {
x_frag <- predict_smooth_ga(object$smooth[[k]], param_vals, vrf[k])
} else {
x_frag <- predict_matrix_tensor_smooth_ga(object$smooth[[k]], param_vals, vrf)
}
x[, object$smooth[[k]]$first.para:object$smooth[[k]]$last.para] <- x_frag
}
# Coefficients of GAM model
beta <- coef(object)
# Compute conditional Loss
l_tilde <- x %*% beta
return(l_tilde)
}
#' Function to compute the preposterior for each of the
#' basis functions of a smooth for one parameter
#'
#' @keywords internal
#' @return Pre-posterior for each basis function
#' @importFrom mgcv PredictMat
predict_smooth_ga <- function(object, param_vals, vrf = 1) {
# Produce basis functions for one parameter
x <- PredictMat(object, param_vals)
# Number of observations
n_obs <- nrow(x)
# Apply variance reduction to compute the preposterior
# for each of the basis functions
# Vector of ones
ones <- matrix(1, n_obs, 1)
# Compute phi on each of the basis function
x <- vrf * x + (1 - vrf) * (ones %*% colMeans(x))
return(x)
}
#' Predict matrix tensor smooth (GA)
#'
#' @description
#' Function to compute the preposterior for each of the
#' basis functions for one or more parameters and calculates
#' the tensor product if more than one parameter is selected
#' (Heavily based on function Predict.matrix.tensor.smooth from
#' mgcv package)
#' @return tensor product based on list of pre-posterior distributions
#' @keywords internal
#' @importFrom mgcv tensor.prod.model.matrix Predict.matrix PredictMat
predict_matrix_tensor_smooth_ga <- function(object,
param_vals,
vrf = rep(1, ncol(param_vals))) {
m <- length(object$margin)
x <- list()
for (i in seq_len(m)) {
term <- object$margin[[i]]$term
dat <- list()
for (j in seq_len(length(term))) {
dat[[term[j]]] <- param_vals[[term[j]]]
}
x[[i]] <- if (!is.null(object$mc[i])) {
PredictMat(object$margin[[i]], dat, n = length(dat[[1]]))
} else {
Predict.matrix(object$margin[[i]], dat)
}
n_obs <- nrow(x[[i]])
}
mxp <- length(object$XP)
if (mxp > 0) {
for (i in seq_len(mxp)) if (!is.null(object$XP[[i]])) {
x[[i]] <- x[[i]] %*% object$XP[[i]]
}
}
# Apply variance reduction to compute the preposterior
# for each of the basis functions
# Vector of ones
ones <- matrix(1, n_obs, 1)
# Initialize and fill list with preposterior of basis functions
# for each parameter
x_ga <- list()
for (i in seq_len(m)) {
x_ga[[i]] <- vrf[i] * x[[i]] + (1 - vrf[i]) * (ones %*% colMeans(x[[i]]))
}
# Compute tensor product
t_ga <- tensor.prod.model.matrix(x_ga)
return(t_ga)
}
#' Plot of Expected Value of Sample Information (EVSI)
#'
#' @description
#' Plots the \code{evsi} object created by \code{\link{calc_evsi}}.
#' EVSI is either plotted as a function of additional sample size
#' for each willingness-to-pay threshold provided, or
#' as a function of each willingness-to-pay
#' threshold, depending upon the usage of \code{\link{calc_evsi}}
#' used to create the \code{evsi} object.
#'
#' @param x object of class \code{evsi}, produced by function
#' \code{\link{calc_evsi}}
#' @param currency string with currency used in the cost-effectiveness analysis (CEA).
#' Default: $, but it could be any currency symbol or word (e.g., £, €, peso)
#' @param effect_units units of effectiveness. Default: QALY
#' @inheritParams add_common_aes
#' @keywords expected value of sample information
#' @return A \code{ggplot2} plot with the EVSI
#' @seealso \code{\link{calc_evsi}}
#' @import ggplot2
#' @importFrom scales comma
#' @export
plot.evsi <- function(x,
txtsize = 12,
currency = "$",
effect_units = "QALY",
n_y_ticks = 8,
n_x_ticks = 20,
xbreaks = NULL,
ybreaks = NULL,
xlim = c(0, NA),
ylim = NULL,
col = c("full", "bw"),
...) {
# Select EVSI data.frame, discarding metamodel list
x <- x[[1]]
x$WTP_thou <- x$WTP / 1000
col <- match.arg(col)
if (length(unique(x$WTP)) == 1) {
col <- "bw"
}
scale_text <- paste("Willingness to Pay\n(Thousand ", currency, "/", effect_units, ")", sep = "")
if (length(unique(x$WTP)) == 1 & "n" %in% names(x)) {
g <- ggplot(data = x,
aes_(x = as.name("n"), y = as.name("EVSI"))) +
xlab("Additional Sample Size")
} else if (!("n" %in% names(x)) | ("n" %in% names(x) & length(unique(x$n)) == 1)) {
g <- ggplot(data = x,
aes_(x = as.name("WTP_thou"), y = as.name("EVSI"))) +
xlab(scale_text)
} else if (length(unique(x$WTP > 1)) & length(unique(x$n)) > 1) {
x$WTP_thou <- as.factor(x$WTP_thou)
g <- ggplot(data = x,
aes_(x = as.name("n"), y = as.name("EVSI"), color = as.name("WTP_thou"))) +
xlab("Additional Sample Size")
} else {
stop("Insufficient data to plot")
}
g <- g +
geom_line() +
ylab(paste("EVSI (", currency, ")", sep = ""))
add_common_aes(g, txtsize, continuous = c("x", "y"),
n_x_ticks = n_x_ticks, n_y_ticks = n_y_ticks,
xbreaks = xbreaks, ybreaks = ybreaks,
xlim = xlim, ylim = ylim, scale_name = scale_text,
col = col, ...)
}
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.