## Partial EVPI functions ## Author: Mark Strong ## Dec 2012 ## Requires packages MASS and mgcv ## We assume that the PSA sample is held in two objects. ## The input parameter samples are held in a N x r matrix (or data frame), where N is the number of PSA runs ## and r is the number of parameters. ## The corresponding net benefits are held in a N x D matrix (or data frame), ## where D is the number of decision options. ## NOTE ## Bug fix 25th Jan 2019 ## There was a bug in the code for computing the upward bias of the GP derived EVPPI estimate ## This would have had an unpredictable effect on the upward bias (including making it negative) ## Thanks to Howard Thom and Wei Fang for reporting the bug ########################################################################################################## ## This is a generic function for estimating partial EVPI via GAM ## NB is the matrix (or data frame) of net benefits ## input.parameters is the matrix (or data frame) of input parameters ## regression.model is the specification of the regression model (must be in quotes) ## For parameters that are expected to interact use e.g. for two parameters "te(x1,x2)" ## For parameters that are not expected to interact use e.g. "s(x1)+s(x2)" ## If there are more than three parameters that are expected to interact set a maximum basis dimension of 4 ## This will avoid the model trying to estimate too many coefficients "te(x1,x2,x3,x4,k=4)" ## If there are more than six parameters that are expected to interact, the GP approach may be better. ## In some scenarios it may be possible to group parameters together that will interact into separate sets ## e.g. regression.model<-"te(x1,x2,x3)+te(x4,x5,x6)" ########################################################################################################### evpi.gam <- function(NB, input.parameters, regression.model) { library(mgcv) if (!is.null(dim(NB))) { NB <- NB - NB[, 1] } else { NB <- cbind(0, NB) } D <- ncol(NB) N <- nrow(NB) g.hat <- vector("list", D) g.hat[[1]] <- rep(0, N) for (d in 2:D) { print(paste("estimating g.hat for incremental NB for option", d, "versus 1")) f <- update(formula(NB[, d] ~ .), formula(paste(".~", regression.model))) model <- gam(f, data = data.frame(input.parameters)) g.hat[[d]] <- model$fitted } perfect.info <- mean(do.call(pmax, g.hat)) baseline <- max(unlist(lapply(g.hat, mean))) rm(g.hat) gc() partial.evpi <- perfect.info - baseline ## estimate EVPI return(partial.evpi) } ############################################################################## ## This is function includes estimation of SE and upward bias for EVPI via GAM ## S is the simulation size for the Monte Carlo computation of SE and bias ############################################################################## evpi.gam.SE.bias <- function(NB, input.parameters, regression.model, S = 1000) { library(mgcv) library(MASS) if (!is.null(dim(NB))) { NB <- NB - NB[, 1] } else { NB <- cbind(0, NB) } D <- ncol(NB) N <- nrow(NB) g.hat <- beta.hat <- Xstar <- V <- tilde.g <- vector("list", D) g.hat[[1]] <- rep(0, N) for (d in 2:D) { print(paste("estimating g.hat for incremental NB for option", d, "versus 1")) f <- update(formula(NB[, d] ~ .), formula(paste(".~", regression.model))) model <- gam(f, data = data.frame(input.parameters)) g.hat[[d]] <- model$fitted beta.hat[[d]] <- model$coef Xstar[[d]] <- predict(model, type = "lpmatrix") V[[d]] <- model$Vp } perfect.info <- mean(do.call(pmax, g.hat)) baseline <- max(unlist(lapply(g.hat, mean))) partial.evpi <- perfect.info - baseline ## estimate EVPI rm(g.hat) gc() print("computing standard error and upward bias via Monte Carlo") for (d in 2:D) { sampled.coef <- mvrnorm(S, beta.hat[[d]], V[[d]]) tilde.g[[d]] <- sampled.coef %*% t(Xstar[[d]]) } tilde.g[[1]] <- matrix(0, nrow = S, ncol = N) rm(V, beta.hat, Xstar, sampled.coef) gc() sampled.perfect.info <- rowMeans(do.call(pmax, tilde.g)) sampled.baseline <- do.call(pmax, lapply(tilde.g, rowMeans)) rm(tilde.g) gc() sampled.partial.evpi <- sampled.perfect.info - sampled.baseline SE <- sd(sampled.partial.evpi) upward.bias <- mean(sampled.partial.evpi) - partial.evpi return(list(partial.evpi = partial.evpi, SE = SE, upward.bias = upward.bias)) } ############################################################################# ## This function computes the log posterior density of the GP hyperparameters ############################################################################# post.density <- function(hyperparams, NB, input.parameters, inputs.of.interest) { dinvgamma <- function(x, alpha, beta) { (beta^alpha)/gamma(alpha) * x^(-alpha - 1) * exp(-beta/x) } input.matrix <- as.matrix(input.parameters[, inputs.of.interest, drop = FALSE]) colmin <- apply(input.matrix, 2, min) colmax <- apply(input.matrix, 2, max) colrange <- colmax - colmin input.matrix <- sweep(input.matrix, 2, colmin, "-") input.matrix <- sweep(input.matrix, 2, colrange, "/") N <- nrow(input.matrix) p <- ncol(input.matrix) H <- cbind(1, input.matrix) q <- ncol(H) a.sigma <- 0.001 b.sigma <- 0.001 ## hyperparameters for IG prior for sigma^2 a.nu <- 0.001 b.nu <- 1 ## hyperparameters for IG prior for nu delta <- exp(hyperparams)[1:p] nu <- exp(hyperparams)[p + 1] A <- exp(-(as.matrix(dist(t(t(input.matrix)/delta), upper = TRUE, diag = TRUE))^2)) Astar <- A + nu * diag(N) T <- chol(Astar) y <- backsolve(t(T), NB, upper.tri = FALSE) x <- backsolve(t(T), H, upper.tri = FALSE) tHAstarinvH <- t(x) %*% (x) betahat <- solve(tHAstarinvH) %*% t(x) %*% y residSS <- y %*% y - t(y) %*% x %*% betahat - t(betahat) %*% t(x) %*% y + t(betahat) %*% tHAstarinvH %*% betahat prior <- prod(dnorm(log(delta), 0, sqrt(1e+05))) * dinvgamma(nu, a.nu, b.nu) l <- -sum(log(diag(T))) - 1/2 * log(det(tHAstarinvH)) - (N - q + 2 * a.sigma)/2 * log(residSS/2 + b.sigma) + log(prior) names(l) <- NULL return(l) } ######################################################################### ## This function estimates the GP hyperparameters numerically using optim ## m is the number of PSA samples used to estimate the hyperparameters ######################################################################### estimate.hyperparameters <- function(NB, input.parameters, inputs.of.interest, m = 500) { if (!is.null(dim(NB))) { NB <- NB - NB[, 1] } else { NB <- cbind(0, NB) } NB <- as.matrix(NB) p <- length(inputs.of.interest) D <- ncol(NB) hyperparameters <- vector("list", D) hyperparameters[[1]] <- NA for (d in 2:D) { initial.values <- rep(0, p + 1) repeat { print(paste("calling optim function for net benefit", d)) log.hyperparameters <- optim(initial.values, fn = post.density, NB = NB[1:m, d], input.parameters = input.parameters[1:m, ], inputs.of.interest = inputs.of.interest, method = "Nelder-Mead", control = list(fnscale = -1, maxit = 10000, trace = 0))$par if (sum(abs(initial.values - log.hyperparameters)) < 0.01) { hyperparameters[[d]] <- exp(log.hyperparameters) break } initial.values <- log.hyperparameters } } return(hyperparameters) } ########################################################################## ## This function estimates the EVPI via GP ## The function takes the hyperparameters estimated above as an argument ## The function also reports standard error and upward bias ## S is the simulation size for the Monte Carlo computation of SE and bias ########################################################################## evpi.GP <- function(NB, input.parameters, inputs.of.interest, hyperparameters, compute.SE = TRUE, S = 500) { if (!is.null(dim(NB))) { NB <- NB - NB[, 1] } else { NB <- cbind(0, NB) } NB <- as.matrix(NB) p <- length(inputs.of.interest) D <- ncol(NB) input.matrix <- as.matrix(input.parameters[, inputs.of.interest, drop = FALSE]) colmin <- apply(input.matrix, 2, min) colmax <- apply(input.matrix, 2, max) colrange <- colmax - colmin input.matrix <- sweep(input.matrix, 2, colmin, "-") input.matrix <- sweep(input.matrix, 2, colrange, "/") N <- nrow(input.matrix) p <- ncol(input.matrix) H <- cbind(1, input.matrix) q <- ncol(H) V <- g.hat <- vector("list", D) g.hat[[1]] <- rep(0, N) for (d in 2:D) { print(paste("estimating g.hat for incremental NB for option", d, "versus 1")) delta.hat <- hyperparameters[[d]][1:p] nu.hat <- hyperparameters[[d]][p + 1] A <- exp(-(as.matrix(dist(t(t(input.matrix)/delta.hat), upper = TRUE, diag = TRUE))^2)) Astar <- A + nu.hat * diag(N) Astarinv <- chol2inv(chol(Astar)) rm(Astar) gc() AstarinvY <- crossprod(Astarinv, NB[, d]) tHAstarinv <- crossprod(H, Astarinv) tHAHinv <- solve(tHAstarinv %*% H) betahat <- tHAHinv %*% (tHAstarinv %*% NB[, d]) Hbetahat <- H %*% betahat resid <- NB[, d] - Hbetahat g.hat[[d]] <- Hbetahat + crossprod(A, crossprod(Astarinv, resid)) AAstarinvH <- A %*% t(tHAstarinv) sigmasqhat <- as.numeric(t(resid) %*% Astarinv %*% resid)/(N - q - 2) V[[d]] <- sigmasqhat * (nu.hat * diag(N) - nu.hat^2 * Astarinv + (H - AAstarinvH) %*% (tHAHinv %*% t(H - AAstarinvH))) rm(A, Astarinv, AstarinvY, tHAstarinv, tHAHinv, betahat, Hbetahat, resid, sigmasqhat) gc() } perfect.info <- mean(do.call(pmax, g.hat)) baseline <- max(unlist(lapply(g.hat, mean))) partial.evpi <- perfect.info - baseline if (compute.SE) { print("computing standard error and upward bias via Monte Carlo") tilde.g <- vector("list", D) tilde.g[[1]] <- matrix(0, nrow = S, ncol = min(N, 1000)) for (d in 2:D) { tilde.g[[d]] <- mvrnorm(S, g.hat[[d]][1:(min(N, 1000))], V[[d]][1:(min(N, 1000)), 1:(min(N, 1000))]) } sampled.perfect.info <- rowMeans(do.call(pmax, tilde.g)) sampled.baseline <- do.call(pmax, lapply(tilde.g, rowMeans)) rm(tilde.g) gc() sampled.partial.evpi <- sampled.perfect.info - sampled.baseline SE <- sd(sampled.partial.evpi) g.hat.short <- lapply(g.hat, function(x) x[1:(min(N, 1000))]) mean.partial.evpi <- mean(do.call(pmax, g.hat.short)) - max(unlist(lapply(g.hat.short, mean))) upward.bias <- mean(sampled.partial.evpi) - mean.partial.evpi rm(V, g.hat) gc() return(list(partial.evpi = partial.evpi, SE = SE, upward.bias = upward.bias)) } else { rm(V, g.hat) gc() return(list(partial.evpi = partial.evpi)) } }