From 32960e48b692cb26537e3e81e79ee5b822740051 Mon Sep 17 00:00:00 2001 From: Schick Date: Tue, 30 Apr 2024 17:01:16 +0200 Subject: [PATCH 01/20] - added three new functions (createPriorMix(), mvpostmix() and getPosteriorOutput()) to enable use of full covariance matrix for posteriors - added if-condition to getPosterior() to differenciate between se_matrix and se_vector input - added new attribute to getPosteriorI() to also include full covariance matrices in output --- R/posterior.R | 182 ++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 175 insertions(+), 7 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index fa2a218..b960c5d 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -47,11 +47,29 @@ getPosterior <- function( if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) { - posterior_list <- getPosteriorI( - prior_list = prior_list, - mu_hat = mu_hat, - se_hat = se_hat, - calc_ess = calc_ess) + if (dim(se_hat)[2] > 1 && dim(se_hat)[2] == length(prior_list)) { + + prior_mix <- createPriorMix(prior_list) + + posterior <- mvpostmix( + priormix = prior_mix, + mu_hat = mu_hat, + se_hat = se_hat) + + posterior_list <- getPosteriorOutput( + posterior_list = posterior, + prior_list = prior_list, + calc_ess = calc_ess) + + } else if (dim(se_hat)[2] == 1) { + + posterior_list <- getPosteriorI( + prior_list = prior_list, + mu_hat = mu_hat, + se_hat = se_hat, + calc_ess = calc_ess) + + } } else if (is.null(mu_hat) && is.null(se_hat) && !is.null(data)) { @@ -83,7 +101,7 @@ getPosteriorI <- function( calc_ess = FALSE ) { - + checkmate::check_data_frame(data_i, null.ok = TRUE) checkmate::check_list(prior_list, names = "named", any.missing = FALSE) checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE) @@ -96,7 +114,7 @@ getPosteriorI <- function( checkmate::check_data_frame(data_i, null.ok = FALSE) checkmate::assert_names(names(data_i), must.include = "response") checkmate::assert_names(names(data_i), must.include = "dose") - + anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1) mu_hat <- summary(anova_res)$coefficients[, 1] se_hat <- summary(anova_res)$coefficients[, 2] @@ -136,6 +154,9 @@ getPosteriorI <- function( } + attr(post_list, "full covariance matrices") <- replicate(length(prior_list)-1, diag(c(se_hat)), simplify = FALSE) + names(attr(post_list, "full covariance matrices")) <- c("comp1", "comp2", "comp3", "robust") + return (post_list) } @@ -176,4 +197,151 @@ getPostCombsI <- function ( return (post_combs) +} + +################################################################################ +# @all: content of createPriorMix() is from DoseFinding (see: bMCTtest.R lines 73-88) +# mvpostmix() has been migrated as is from DoseFinding package +# getPosteriorOutput() turns posterior object from mvpostmix() into RBesT type object +# new variables have been defined in setup.R for new tests +# new tests for all new functions have been added in test-posterior.R +# -> please check out code & outputs and make changes if necessary + +createPriorMix <- function(prior) { + + checkmate::check_list(prior, names = "named", any.missing = FALSE, null.ok = FALSE) + + k <- length(prior) + + n_comps <- unlist(lapply(prior, ncol)) + args <- lapply(1:k, function(x) 1:n_comps[x]) + comp_ind <- do.call("expand.grid", args) + + n_comps_prior <- nrow(comp_ind) + + prior_weight <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][1, comp_ind[y,x]])), nrow = n_comps_prior) + prior_weight <- apply(prior_weight, 1, prod) + prior_mean <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][2, comp_ind[y,x]])), nrow = n_comps_prior) + prior_sd <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][3, comp_ind[y,x]])), nrow = n_comps_prior) + + prior_weight <- as.list(prior_weight) + prior_mean <- asplit(prior_mean, 1) + prior_vc <- lapply(asplit(prior_sd^2, 1), diag) + prior_mix <- list(prior_weight, prior_mean, prior_vc) + + return(prior_mix) + +} + +mvpostmix <- function( + + priormix, + mu_hat, + se_hat + +) { + + checkmate::check_list(priormix, names = "named", any.missing = FALSE, null.ok = FALSE) + checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = FALSE) + checkmate::check_double(mu_hat, null.ok = FALSE, lower = -Inf, upper = Inf) + checkmate::check_matrix(se_hat, any.missing = FALSE, null.ok = FALSE) + checkmate::check_double(se_hat, null.ok = FALSE, lower = -Inf, upper = Inf) + + logSumExp <- function(lx){ + + lm <- max(lx) + lm + log(sum(exp(lx - lm))) + + } + + dataPrec <- solve(se_hat) + priorPrec <- lapply(priormix[[3]], solve) + postPrec <- lapply(priorPrec, function(x) x + dataPrec) + SigmaPred <- lapply(priormix[[3]], function(x) x + se_hat) + + lw <- numeric(length(priormix[[1]])) + postmix <- vector("list", 3) + names(postmix) <- c("weights", "mean", "covmat") + + for(i in 1:3) + postmix[[i]] <- vector("list", length(lw)) + + ## The posterior distribution is a mixture of multivariate normals with updated mixture weights. + ## Posterior weights are updated based on the prior predictive (marginal) probabilities of the data under each + ## component of the mixture. + ## In the case of a MVN likelihood with known covariance and MVN priors for the mean the + ## prior predictive distributions are MVN distribution with mean vectors equal to the prior components' mean vectors + ## and covariance matrices which are the sum of the prior components' covariance matrices and the "known" covariance + ## matrix of the data (for which S_hat is plugged in here) + for(i in 1:length(lw)){ + + lw[i] <- log(priormix[[1]][[i]]) +mvtnorm::dmvnorm(mu_hat, priormix[[2]][[i]], SigmaPred[[i]], log = TRUE) + postmix[[2]][[i]] <- solve(priorPrec[[i]] + dataPrec) %*% (priorPrec[[i]] %*% priormix[[2]][[i]] + dataPrec %*% mu_hat) + postmix[[3]][[i]] <- solve(priorPrec[[i]] + dataPrec) + + } + + postmix[[1]] <- as.list(exp(lw - logSumExp(lw))) + + names_comps <- c("comp1", "comp2", "comp3", "robust") + + for(i in 1:3) + names(postmix[[i]]) <- names_comps[1:length(lw)] + + postmix + +} + +getPosteriorOutput <- function( + + posterior_list, + prior_list, + calc_ess = FALSE + +) { + + checkmate::check_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + checkmate::check_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + + posterior_output <- vector("list", length(prior_list)) + names(posterior_output) <- names(prior_list) + + posterior_list_weights <- t(posterior_list$weights) + posterior_list_means <- t(sapply(posterior_list$mean, rbind)) + ctr_data <- rbind(rbind(posterior_list_weights, posterior_list_means[,1]), sapply(1:(length(posterior_list$covmat)), function(x) posterior_list$covmat[[x]][1,1])) + + pos_data_frame <- as.data.frame(ctr_data) + + posterior_ctr <- mixnorm(unlist(pos_data_frame$comp1), + unlist(pos_data_frame$comp2), + unlist(pos_data_frame$comp3), + unlist(pos_data_frame$robust), + sigma = sigma(prior_list$Ctr)) + + colnames(posterior_ctr)[4] <- "robust" + posterior_output[[1]] <- posterior_ctr + + posterior_output[2:length(prior_list)] <- lapply(2:length(posterior_output), function(x) { + + covmat_sum <- sum(sapply(1:4, function(y) posterior_list$covmat[[y]][x,x])) + mixnorm(comp1 = c(sum(unlist(posterior_list$weights)), sum(unlist(posterior_list_means[,x])), covmat_sum), sigma = sigma(prior_list$DG_1)) + + }) + + class(posterior_output) <- "postList" + + if (calc_ess) { + + attr(posterior_output, "ess") <- getESS(posterior_output) + + } else { + + attr(posterior_output, "ess") <- numeric(0) + + } + + attr(posterior_output, "full covariance matrices") <- posterior_list$covmat + + return(posterior_output) + } \ No newline at end of file From 6968e1f20710e31e7312f77b1be6058883cf42cb Mon Sep 17 00:00:00 2001 From: Schick Date: Tue, 30 Apr 2024 17:03:11 +0200 Subject: [PATCH 02/20] - added parameter definitions to setup.R for new tests - added new tests to test-posterior.R for new functions in posterior.R (full covariance matrix) --- tests/testthat/setup.R | 35 +++++++++++++++++++ tests/testthat/test-posterior.R | 60 +++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 2a7bda7..63e8753 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -154,3 +154,38 @@ crit_pval = getCritProb( # alpha_crit_val = alpha_crit_val, # simple = TRUE # ) + +# Create covmat test case +mu_hat <- c(10, 20, 30, 40, 50) +se_hat_vector <- matrix(c(2.45, 1.76, 2.19, 2.94, 1.66), nrow = 5, ncol = 1) + +se_hat_matrix <- matrix(c(2.45, 0.01, 0.01, 0.01, 0.01, + 0.01, 1.76, 0.01, 0.01, 0.01, + 0.01, 0.01, 2.19, 0.01, 0.01, + 0.01, 0.01, 0.01, 2.94, 0.01, + 0.01, 0.01, 0.01, 0.01, 1.66), nrow = 5, ncol = 5) + +data("metaData") + +dataset_covmat <- filter(as.data.frame(metaData), bname == "BRINTELLIX") + +histcontrol_covmat <- filter( + dataset_covmat, + dose == 0, + primtime == 8, + indication == "MAJOR DEPRESSIVE DISORDER", + protid != 5) + +hist_data_covmat <- data.frame( + trial = histcontrol_covmat$nctno, + est = histcontrol_covmat$rslt, + se = histcontrol_covmat$se, + sd = histcontrol_covmat$sd, + n = histcontrol_covmat$sampsize) + +dose_levels_covmat <- c(0, 2.5, 5, 7.5, 10) + +prior_list_covmat <- getPriorList( + hist_data = hist_data_covmat, + dose_levels = dose_levels_covmat, + robust_weight = 0.3) diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R index c26bf01..f2906bd 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -105,3 +105,63 @@ test_that("getPostCombsI returns an object with correct attributes", { expect_equal(result$weights, 4) expect_true(all(result$vars == c(4, 4))) }) + +test_that("getPriorMix works correctly", { + + # function call without parameters + expect_error(createPriorMix()) + + # test createPriorMix function + prior_mix <- createPriorMix(prior_list_covmat) + expect_type(prior_mix, "list") + expect_length(prior_mix, 3) + +}) + +test_that("mvpostmix works correctly", { + + prior_mix <- createPriorMix(prior_list_covmat) + + # test mvpostmix function + expect_error(mvpostmix(prior_mix, mu_hat, se_hat_vector)) + expect_no_error(mvpostmix(prior_mix, mu_hat, se_hat_matrix)) + posterior <- mvpostmix(prior_mix, mu_hat, se_hat_matrix) + expect_type(posterior, "list") + expect_length(posterior, 3) + +}) + +test_that("getPosteriorOutput works correctly", { + + # create posterior with matrix + prior_mix <- createPriorMix(prior_list_covmat) + + posterior <- mvpostmix( + priormix = prior_mix, + mu_hat = mu_hat, + se_hat = se_hat_matrix) + + # create posterior with vector + posterior_vector <- getPosteriorI( + prior_list = prior_list_covmat, + mu_hat = mu_hat, + se_hat = se_hat_vector, + calc_ess = FALSE + ) + + # test getPosteriorOutput function + posterior_list <- getPosteriorOutput(posterior, prior_list_covmat, calc_ess) + expect_type(posterior_list, "list") + expect_s3_class(posterior_list, "postList") + + lapply(1:(length(dose_levels_covmat)-1), function(x) { + expect_equal(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)], posterior$covmat$Comp1[row(posterior$covmat$Comp1) == col(posterior$covmat$Comp1) + x]) + expect_equal(length(which(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)] >= 0)), length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)])) + }) + + # compare posterior result object with matrix to object with vector + expect_length(posterior_list, length(posterior_vector)) + expect_type(posterior_list, typeof(posterior_vector)) + expect_s3_class(posterior_list, S3Class(posterior_vector)) + +}) \ No newline at end of file From 11f4e193c79e5a5c70a82b4f459768616bea3599 Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 2 May 2024 17:01:27 +0200 Subject: [PATCH 03/20] - getPosterior() is now calling mvpostmix() straight out of DoseFinding package --- R/posterior.R | 71 +++++---------------------------------------------- 1 file changed, 6 insertions(+), 65 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index b960c5d..3c06de0 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -51,10 +51,10 @@ getPosterior <- function( prior_mix <- createPriorMix(prior_list) - posterior <- mvpostmix( + posterior <- DoseFinding::mvpostmix( priormix = prior_mix, mu_hat = mu_hat, - se_hat = se_hat) + S_hat = se_hat) posterior_list <- getPosteriorOutput( posterior_list = posterior, @@ -232,65 +232,6 @@ createPriorMix <- function(prior) { return(prior_mix) } - -mvpostmix <- function( - - priormix, - mu_hat, - se_hat - -) { - - checkmate::check_list(priormix, names = "named", any.missing = FALSE, null.ok = FALSE) - checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = FALSE) - checkmate::check_double(mu_hat, null.ok = FALSE, lower = -Inf, upper = Inf) - checkmate::check_matrix(se_hat, any.missing = FALSE, null.ok = FALSE) - checkmate::check_double(se_hat, null.ok = FALSE, lower = -Inf, upper = Inf) - - logSumExp <- function(lx){ - - lm <- max(lx) - lm + log(sum(exp(lx - lm))) - - } - - dataPrec <- solve(se_hat) - priorPrec <- lapply(priormix[[3]], solve) - postPrec <- lapply(priorPrec, function(x) x + dataPrec) - SigmaPred <- lapply(priormix[[3]], function(x) x + se_hat) - - lw <- numeric(length(priormix[[1]])) - postmix <- vector("list", 3) - names(postmix) <- c("weights", "mean", "covmat") - - for(i in 1:3) - postmix[[i]] <- vector("list", length(lw)) - - ## The posterior distribution is a mixture of multivariate normals with updated mixture weights. - ## Posterior weights are updated based on the prior predictive (marginal) probabilities of the data under each - ## component of the mixture. - ## In the case of a MVN likelihood with known covariance and MVN priors for the mean the - ## prior predictive distributions are MVN distribution with mean vectors equal to the prior components' mean vectors - ## and covariance matrices which are the sum of the prior components' covariance matrices and the "known" covariance - ## matrix of the data (for which S_hat is plugged in here) - for(i in 1:length(lw)){ - - lw[i] <- log(priormix[[1]][[i]]) +mvtnorm::dmvnorm(mu_hat, priormix[[2]][[i]], SigmaPred[[i]], log = TRUE) - postmix[[2]][[i]] <- solve(priorPrec[[i]] + dataPrec) %*% (priorPrec[[i]] %*% priormix[[2]][[i]] + dataPrec %*% mu_hat) - postmix[[3]][[i]] <- solve(priorPrec[[i]] + dataPrec) - - } - - postmix[[1]] <- as.list(exp(lw - logSumExp(lw))) - - names_comps <- c("comp1", "comp2", "comp3", "robust") - - for(i in 1:3) - names(postmix[[i]]) <- names_comps[1:length(lw)] - - postmix - -} getPosteriorOutput <- function( @@ -312,10 +253,10 @@ getPosteriorOutput <- function( pos_data_frame <- as.data.frame(ctr_data) - posterior_ctr <- mixnorm(unlist(pos_data_frame$comp1), - unlist(pos_data_frame$comp2), - unlist(pos_data_frame$comp3), - unlist(pos_data_frame$robust), + posterior_ctr <- mixnorm(unlist(pos_data_frame$Comp1), + unlist(pos_data_frame$Comp2), + unlist(pos_data_frame$Comp3), + unlist(pos_data_frame$Comp4), sigma = sigma(prior_list$Ctr)) colnames(posterior_ctr)[4] <- "robust" From a3e68c03212b1d4ec63dc69476bbf75dbcd6a10b Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 2 May 2024 17:03:56 +0200 Subject: [PATCH 04/20] - added getPosterior() call for testing --- tests/testthat/setup.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 63e8753..a0937e0 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -189,3 +189,10 @@ prior_list_covmat <- getPriorList( hist_data = hist_data_covmat, dose_levels = dose_levels_covmat, robust_weight = 0.3) + +posterior_covmat <- getPosterior( + prior_list = prior_list_covmat, + mu_hat = mu_hat, + se_hat = se_hat_matrix, + calc_ess = FALSE +) \ No newline at end of file From 5a59a4bc6ad643049ff608db129826392da3f438 Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 16 May 2024 11:10:35 +0200 Subject: [PATCH 05/20] - renamed getPosteriorOutput() function to postmix2RBesT() and updated input parameters - updated postmix2RBesT() to implement more efficient coding methods - extracted ESS calculation steps into its own function calcESS() in order to avoid redundant code --- R/posterior.R | 134 ++++++++++++++++++++++++-------------------- man/getPosterior.Rd | 6 +- 2 files changed, 77 insertions(+), 63 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index 3c06de0..8ef2092 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -8,7 +8,7 @@ #' @param data dataframe containing the information of dose and response. Default NULL #' Also a simulateData object can be provided. #' @param mu_hat vector of estimated mean values (per dose group). -#' @param se_hat vector of estimated standard deviations (per dose group). +#' @param S_hat vector or matrix of estimated standard deviations (per dose group). #' @param calc_ess boolean variable, indicating whether effective sample size should be calculated. Default FALSE #' @return posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group #' @examples @@ -33,7 +33,7 @@ getPosterior <- function( prior_list, data = NULL, mu_hat = NULL, - se_hat = NULL, + S_hat = NULL, calc_ess = FALSE ) { @@ -42,36 +42,36 @@ getPosterior <- function( checkmate::check_list(prior_list, names = "named", any.missing = FALSE) checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf) - checkmate::check_vector(se_hat, any.missing = FALSE, null.ok = TRUE) - checkmate::check_double(se_hat, null.ok = TRUE, lower = 0, upper = Inf) + checkmate::check_vector(S_hat, any.missing = FALSE, null.ok = TRUE) + checkmate::check_double(S_hat, null.ok = TRUE, lower = 0, upper = Inf) - if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) { + if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) { - if (dim(se_hat)[2] > 1 && dim(se_hat)[2] == length(prior_list)) { - + if (dim(S_hat)[2] > 1 && dim(S_hat)[2] == length(prior_list)) { prior_mix <- createPriorMix(prior_list) posterior <- DoseFinding::mvpostmix( - priormix = prior_mix, - mu_hat = mu_hat, - S_hat = se_hat) + priormix = prior_mix, + mu_hat = mu_hat, + S_hat = S_hat) - posterior_list <- getPosteriorOutput( + posterior_list <- postmix2RBesT( posterior_list = posterior, prior_list = prior_list, + prior_mix = prior_mix, calc_ess = calc_ess) - } else if (dim(se_hat)[2] == 1) { + } else if (dim(S_hat)[2] == 1) { posterior_list <- getPosteriorI( prior_list = prior_list, mu_hat = mu_hat, - se_hat = se_hat, + se_hat = S_hat, calc_ess = calc_ess) } - } else if (is.null(mu_hat) && is.null(se_hat) && !is.null(data)) { + } else if (is.null(mu_hat) && is.null(S_hat) && !is.null(data)) { posterior_list <- lapply(split(data, data$simulation), getPosteriorI, prior_list = prior_list, calc_ess = calc_ess) @@ -144,18 +144,10 @@ getPosteriorI <- function( names(post_list) <- names(prior_list) class(post_list) <- "postList" - if (calc_ess) { - - attr(post_list, "ess") <- getESS(post_list) - - } else { - - attr(post_list, "ess") <- numeric(0) - - } + attr(post_list, "ess") <- calcEss(calc_ess, post_list) attr(post_list, "full covariance matrices") <- replicate(length(prior_list)-1, diag(c(se_hat)), simplify = FALSE) - names(attr(post_list, "full covariance matrices")) <- c("comp1", "comp2", "comp3", "robust") + #names(attr(post_list, "full covariance matrices")) <- c("comp1", "comp2", "comp3", "robust") return (post_list) @@ -199,17 +191,9 @@ getPostCombsI <- function ( } -################################################################################ -# @all: content of createPriorMix() is from DoseFinding (see: bMCTtest.R lines 73-88) -# mvpostmix() has been migrated as is from DoseFinding package -# getPosteriorOutput() turns posterior object from mvpostmix() into RBesT type object -# new variables have been defined in setup.R for new tests -# new tests for all new functions have been added in test-posterior.R -# -> please check out code & outputs and make changes if necessary - createPriorMix <- function(prior) { - checkmate::check_list(prior, names = "named", any.missing = FALSE, null.ok = FALSE) + checkmate::assert_list(prior, names = "named", any.missing = FALSE, null.ok = FALSE) k <- length(prior) @@ -233,56 +217,86 @@ createPriorMix <- function(prior) { } -getPosteriorOutput <- function( +postmix2RBesT <- function( posterior_list, prior_list, + prior_mix, calc_ess = FALSE ) { - checkmate::check_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) - checkmate::check_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + checkmate::assert_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + summed_values <- vector("list", length(prior_list)) posterior_output <- vector("list", length(prior_list)) - names(posterior_output) <- names(prior_list) - - posterior_list_weights <- t(posterior_list$weights) - posterior_list_means <- t(sapply(posterior_list$mean, rbind)) - ctr_data <- rbind(rbind(posterior_list_weights, posterior_list_means[,1]), sapply(1:(length(posterior_list$covmat)), function(x) posterior_list$covmat[[x]][1,1])) - - pos_data_frame <- as.data.frame(ctr_data) - posterior_ctr <- mixnorm(unlist(pos_data_frame$Comp1), - unlist(pos_data_frame$Comp2), - unlist(pos_data_frame$Comp3), - unlist(pos_data_frame$Comp4), - sigma = sigma(prior_list$Ctr)) - - colnames(posterior_ctr)[4] <- "robust" - posterior_output[[1]] <- posterior_ctr - - posterior_output[2:length(prior_list)] <- lapply(2:length(posterior_output), function(x) { + posterior_output <- lapply(seq_along(prior_list), function(i) { + + patterns <- lapply(prior_list[[i]]["m",], function(x) { + + matches <- sapply(seq_along(posterior_list$weights), function(j) { + + prior_mix[[2]][[j]][i] == x + + }) + + result <- which(matches) + + return(result) + + }) - covmat_sum <- sum(sapply(1:4, function(y) posterior_list$covmat[[y]][x,x])) - mixnorm(comp1 = c(sum(unlist(posterior_list$weights)), sum(unlist(posterior_list_means[,x])), covmat_sum), sigma = sigma(prior_list$DG_1)) + diagonals <- lapply(posterior_list$covmat, diag) + diagonals_bound <- sapply(diagonals, cbind) + + summed_values <- lapply(patterns, function(x) { + + c(sum(as.numeric(t(t(posterior_list$weights[x])))), + mean(as.numeric(t(sapply(posterior_list$mean, rbind))[x,i])), + mean(diagonals_bound[i,x]) + ) + + }) + + additional_args <- list(sigma = sigma(prior_list[[i]])) + + args <- c(summed_values, additional_args) + + result <- do.call(RBesT::mixnorm, args) + + return(result) }) + names(posterior_output) <- c("Ctr", paste0("DG_", seq_along(posterior_output[-1]))) + class(posterior_output) <- "postList" + attr(posterior_output, "ess") <- calcEss(calc_ess, posterior_output) + + attr(posterior_output, "full covariance matrices") <- posterior_list$covmat + + return(posterior_output) + +} + +calcEss <- function(calc_ess, posterior_output) { + + checkmate::assert_logical(calc_ess, null.ok = FALSE, len = 1) + checkmate::assert_list(posterior_output, names = "named", any.missing = FALSE, null.ok = FALSE) + if (calc_ess) { - attr(posterior_output, "ess") <- getESS(posterior_output) + post_ESS <- getESS(posterior_output) } else { - attr(posterior_output, "ess") <- numeric(0) + post_ESS <- numeric(0) } - attr(posterior_output, "full covariance matrices") <- posterior_list$covmat - - return(posterior_output) + return(post_ESS) } \ No newline at end of file diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 37e3bbc..8b82c0d 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -8,7 +8,7 @@ getPosterior( prior_list, data = NULL, mu_hat = NULL, - se_hat = NULL, + S_hat = NULL, calc_ess = FALSE ) } @@ -20,9 +20,9 @@ Also a simulateData object can be provided.} \item{mu_hat}{vector of estimated mean values (per dose group).} -\item{se_hat}{vector of estimated standard deviations (per dose group).} - \item{calc_ess}{boolean variable, indicating whether effective sample size should be calculated. Default FALSE} + +\item{se_hat}{vector of estimated standard deviations (per dose group).} } \value{ posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group From 393ac24e84968bddab19077786df9b589632e0d7 Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 16 May 2024 11:13:40 +0200 Subject: [PATCH 06/20] - added full definition of a prior list as an alternative to getPriorList() in setup.R - added more and better fitting unit tests for postmix2RBesT() to cover updates to the function --- tests/testthat/setup.R | 31 +++++++++++++++++++++++++------ tests/testthat/test-posterior.R | 27 ++++++++++++++++++++++----- 2 files changed, 47 insertions(+), 11 deletions(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index a0937e0..c4447db 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -156,13 +156,32 @@ crit_pval = getCritProb( # ) # Create covmat test case +mixnorm_test <- mixnorm(comp1=c(0.2,2,3), comp2=c(0.2,5,6), comp3=c(0.2,8,9), comp4=c(0.2,11,12), robust=c(0.2,14,15), sigma=9.734) + +mixnorm_DG1 <- mixnorm(comp1=c(0.5,2,3), comp2=c(0.5,3,6), sigma=9.651) + +mixnorm_DG2 <- mixnorm(comp1=c(1,8,2), sigma=9.932) + +mixnorm_DG3 <- mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) + +mixnorm_DG4 <- mixnorm(comp1=c(1/4,10,3), comp2=c(1/4,3,6), comp3=c(1/4,4,7), comp4=c(1/4,9,8), sigma=9.236) + +prior_list <- vector("list", 5) +prior_list[[1]] <- mixnorm_test +prior_list[[2]] <- mixnorm_DG1 +prior_list[[3]] <- mixnorm_DG2 +prior_list[[4]] <- mixnorm_DG3 +prior_list[[5]] <- mixnorm_DG4 + +names(prior_list) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") + mu_hat <- c(10, 20, 30, 40, 50) -se_hat_vector <- matrix(c(2.45, 1.76, 2.19, 2.94, 1.66), nrow = 5, ncol = 1) +se_hat_vector <- matrix(c(3.11, 1.76, 0.38, 0.93, 1.66), nrow = 5, ncol = 1) -se_hat_matrix <- matrix(c(2.45, 0.01, 0.01, 0.01, 0.01, +se_hat_matrix <- matrix(c(3.11, 0.01, 0.01, 0.01, 0.01, 0.01, 1.76, 0.01, 0.01, 0.01, - 0.01, 0.01, 2.19, 0.01, 0.01, - 0.01, 0.01, 0.01, 2.94, 0.01, + 0.01, 0.01, 0.38, 0.01, 0.01, + 0.01, 0.01, 0.01, 0.93, 0.01, 0.01, 0.01, 0.01, 0.01, 1.66), nrow = 5, ncol = 5) data("metaData") @@ -193,6 +212,6 @@ prior_list_covmat <- getPriorList( posterior_covmat <- getPosterior( prior_list = prior_list_covmat, mu_hat = mu_hat, - se_hat = se_hat_matrix, + S_hat = se_hat_matrix, calc_ess = FALSE -) \ No newline at end of file +) diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R index f2906bd..1fa0a49 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -136,10 +136,10 @@ test_that("getPosteriorOutput works correctly", { # create posterior with matrix prior_mix <- createPriorMix(prior_list_covmat) - posterior <- mvpostmix( + posterior <- DoseFinding::mvpostmix( priormix = prior_mix, mu_hat = mu_hat, - se_hat = se_hat_matrix) + S_hat = se_hat_matrix) # create posterior with vector posterior_vector <- getPosteriorI( @@ -150,15 +150,32 @@ test_that("getPosteriorOutput works correctly", { ) # test getPosteriorOutput function - posterior_list <- getPosteriorOutput(posterior, prior_list_covmat, calc_ess) + posterior_list <- postmix2RBesT(posterior, prior_list_covmat, calc_ess = FALSE) expect_type(posterior_list, "list") expect_s3_class(posterior_list, "postList") lapply(1:(length(dose_levels_covmat)-1), function(x) { - expect_equal(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)], posterior$covmat$Comp1[row(posterior$covmat$Comp1) == col(posterior$covmat$Comp1) + x]) - expect_equal(length(which(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)] >= 0)), length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)])) + + expect_equal(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)], + posterior$covmat$Comp1[row(posterior$covmat$Comp1) == col(posterior$covmat$Comp1) + x]) + + expect_equal(length(which(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)] >= 0)), + length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)])) + }) + lapply(seq_along(prior_list), function(i) { + + sapply(seq_along(posterior$weights), function(j) { + + expect_in(prior_list[[i]]["m",], prior_mix[[2]][[j]][i]) + + }) + + }) + + expect_equal(length(posterior$weights), prod(lengths(prior_list_covmat)/nrow(prior_list_covmat$Ctr))) + # compare posterior result object with matrix to object with vector expect_length(posterior_list, length(posterior_vector)) expect_type(posterior_list, typeof(posterior_vector)) From bbef1929ef98f30955be7721da3612eafd8ed0f8 Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 23 May 2024 16:54:49 +0200 Subject: [PATCH 07/20] - added new function createMapping() to avoid redundant code - priorList2priorMix(): migrated mapping code to separate function createMapping() - getPosteriorOutput(): completely revamped getPosteriorOutput() to now use already available mapping sheets from priorList2priorMix() for the creation of the posterior list -> also renamed to postMix2posteriorList() --- R/posterior.R | 116 +++++++++++++++++++++++++------------------------- 1 file changed, 57 insertions(+), 59 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index 8ef2092..c4b4f23 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -48,17 +48,17 @@ getPosterior <- function( if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) { if (dim(S_hat)[2] > 1 && dim(S_hat)[2] == length(prior_list)) { - prior_mix <- createPriorMix(prior_list) + + prior_mix <- priorList2priorMix(prior_list) posterior <- DoseFinding::mvpostmix( priormix = prior_mix, mu_hat = mu_hat, S_hat = S_hat) - posterior_list <- postmix2RBesT( + posterior_list <- postMix2posteriorList( posterior_list = posterior, prior_list = prior_list, - prior_mix = prior_mix, calc_ess = calc_ess) } else if (dim(S_hat)[2] == 1) { @@ -191,37 +191,39 @@ getPostCombsI <- function ( } -createPriorMix <- function(prior) { - - checkmate::assert_list(prior, names = "named", any.missing = FALSE, null.ok = FALSE) +priorList2priorMix <- function (prior_list) { - k <- length(prior) + checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) - n_comps <- unlist(lapply(prior, ncol)) - args <- lapply(1:k, function(x) 1:n_comps[x]) + # create mapping + args <- createMapping(prior_list) comp_ind <- do.call("expand.grid", args) - n_comps_prior <- nrow(comp_ind) - prior_weight <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][1, comp_ind[y,x]])), nrow = n_comps_prior) + # map information -> mapping function? + prior_weight <- matrix( + sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, + function (y) prior_list[[x]][1, comp_ind[y, x]])), nrow = n_comps_prior) + + prior_mean <- matrix(sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, function (y) prior_list[[x]][2, comp_ind[y, x]])), nrow = n_comps_prior) + prior_sd <- matrix(sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, function (y) prior_list[[x]][3, comp_ind[y, x]])), nrow = n_comps_prior) + prior_weight <- apply(prior_weight, 1, prod) - prior_mean <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][2, comp_ind[y,x]])), nrow = n_comps_prior) - prior_sd <- matrix(sapply(1:k, function(x) sapply(1:n_comps_prior, function(y) prior[[x]][3, comp_ind[y,x]])), nrow = n_comps_prior) + # create prior_mix object prior_weight <- as.list(prior_weight) - prior_mean <- asplit(prior_mean, 1) - prior_vc <- lapply(asplit(prior_sd^2, 1), diag) - prior_mix <- list(prior_weight, prior_mean, prior_vc) + prior_mean <- asplit(prior_mean, 1) + prior_vc <- lapply(asplit(prior_sd^2, 1), diag) + prior_mix <- list(prior_weight, prior_mean, prior_vc) - return(prior_mix) + return (prior_mix) } -postmix2RBesT <- function( +postMix2posteriorList <- function ( posterior_list, prior_list, - prior_mix, calc_ess = FALSE ) { @@ -229,48 +231,34 @@ postmix2RBesT <- function( checkmate::assert_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) - summed_values <- vector("list", length(prior_list)) - posterior_output <- vector("list", length(prior_list)) + # create mapping + args <- createMapping(prior_list) + comp_ind <- do.call("expand.grid", args) - posterior_output <- lapply(seq_along(prior_list), function(i) { - - patterns <- lapply(prior_list[[i]]["m",], function(x) { - - matches <- sapply(seq_along(posterior_list$weights), function(j) { - - prior_mix[[2]][[j]][i] == x - - }) - - result <- which(matches) - - return(result) - - }) - - diagonals <- lapply(posterior_list$covmat, diag) - diagonals_bound <- sapply(diagonals, cbind) - - summed_values <- lapply(patterns, function(x) { - - c(sum(as.numeric(t(t(posterior_list$weights[x])))), - mean(as.numeric(t(sapply(posterior_list$mean, rbind))[x,i])), - mean(diagonals_bound[i,x]) - ) - - }) - - additional_args <- list(sigma = sigma(prior_list[[i]])) - - args <- c(summed_values, additional_args) - - result <- do.call(RBesT::mixnorm, args) - - return(result) - - }) + # map posterior information + posterior_weight <- lapply(1:length(prior_list), + function (x) sapply(1:length(args[[x]]), + function (y) sum(unlist(posterior$weights[which(comp_ind[ ,x]==args[[x]][y])])))) + + posterior_mean <- lapply(1:length(prior_list), + function (x) sapply(1:length(args[[x]]), + function (y) (Reduce("+", posterior$mean[which(comp_ind[ ,x]==args[[x]][y])]) / length(posterior$mean[which(comp_ind[,x]==args[[x]][y])]))[x, ])) - names(posterior_output) <- c("Ctr", paste0("DG_", seq_along(posterior_output[-1]))) + posterior_sd <- lapply(1:length(prior_list), + function (x) sapply(1:length(args[[x]]), + function (y) (Reduce("+", lapply(posterior$covmat, diag)[which(comp_ind[ ,x]==args[[x]][y])]) / length(lapply(posterior$covmat, diag)[which(comp_ind[ ,x]==args[[x]][y])]))[x])) + + combined_vectors <- mapply(function(x, y, z) Map(c, x, y, z), posterior_weight, posterior_mean, posterior_sd, SIMPLIFY = FALSE) + + # create posterior list + posterior_output <- lapply(1:length(combined_vectors), + function(x) do.call(RBesT::mixnorm, c(combined_vectors[[x]], sigma = sigma(prior_list[[x]])))) + + # TODO: talk about + #colnames(posterior_output$Ctr)[max(length(colnames(posterior_output$Ctr)))] <- "robust" + + names(posterior_output) <- c("Ctr", paste0("DG_", + seq_along(posterior_output[-1]))) class(posterior_output) <- "postList" @@ -299,4 +287,14 @@ calcEss <- function(calc_ess, posterior_output) { return(post_ESS) +} + +createMapping <- function(prior_list) { + + n_dg <- length(prior_list) + n_comps <- unlist(lapply(prior_list, ncol)) + args <- lapply(1:n_dg, function (x) 1:n_comps[x]) + + return(args) + } \ No newline at end of file From 3495c2a6c7905dda00d722d075d1218dce2fcc06 Mon Sep 17 00:00:00 2001 From: Schick Date: Thu, 23 May 2024 16:59:58 +0200 Subject: [PATCH 08/20] - updated tests to new function names - added two new tests and updated other tests to account for new syntax in postMix2posteriorList() --- tests/testthat/test-posterior.R | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R index 1fa0a49..ddca125 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -111,8 +111,8 @@ test_that("getPriorMix works correctly", { # function call without parameters expect_error(createPriorMix()) - # test createPriorMix function - prior_mix <- createPriorMix(prior_list_covmat) + # test priorList2priorMix function + prior_mix <- priorList2priorMix(prior_list_covmat) expect_type(prior_mix, "list") expect_length(prior_mix, 3) @@ -120,12 +120,12 @@ test_that("getPriorMix works correctly", { test_that("mvpostmix works correctly", { - prior_mix <- createPriorMix(prior_list_covmat) + prior_mix <- priorList2priorMix(prior_list_covmat) # test mvpostmix function - expect_error(mvpostmix(prior_mix, mu_hat, se_hat_vector)) - expect_no_error(mvpostmix(prior_mix, mu_hat, se_hat_matrix)) - posterior <- mvpostmix(prior_mix, mu_hat, se_hat_matrix) + expect_error(DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_vector)) + expect_no_error(DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_matrix)) + posterior <- DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_matrix) expect_type(posterior, "list") expect_length(posterior, 3) @@ -134,7 +134,7 @@ test_that("mvpostmix works correctly", { test_that("getPosteriorOutput works correctly", { # create posterior with matrix - prior_mix <- createPriorMix(prior_list_covmat) + prior_mix <- priorList2priorMix(prior_list_covmat) posterior <- DoseFinding::mvpostmix( priormix = prior_mix, @@ -149,8 +149,8 @@ test_that("getPosteriorOutput works correctly", { calc_ess = FALSE ) - # test getPosteriorOutput function - posterior_list <- postmix2RBesT(posterior, prior_list_covmat, calc_ess = FALSE) + # test postMix2posteriorList function + posterior_list <- postMix2posteriorList(posterior, prior_list_covmat, calc_ess = FALSE) expect_type(posterior_list, "list") expect_s3_class(posterior_list, "postList") @@ -166,16 +166,17 @@ test_that("getPosteriorOutput works correctly", { lapply(seq_along(prior_list), function(i) { - sapply(seq_along(posterior$weights), function(j) { + expect_equal(length(posterior_weight[[i]]), length(posterior_mean[[i]])) + expect_equal(length(posterior_mean[[i]]), length(posterior_sd[[i]])) + + sapply(seq_along(posterior_weight[[i]]), function(j) { - expect_in(prior_list[[i]]["m",], prior_mix[[2]][[j]][i]) + expect_length(combined_vectors[[i]][[j]], 3) }) }) - expect_equal(length(posterior$weights), prod(lengths(prior_list_covmat)/nrow(prior_list_covmat$Ctr))) - # compare posterior result object with matrix to object with vector expect_length(posterior_list, length(posterior_vector)) expect_type(posterior_list, typeof(posterior_vector)) From 354ca694f7e88f32d122c0aa26306b2149b3f8d4 Mon Sep 17 00:00:00 2001 From: Schick Date: Tue, 4 Jun 2024 09:41:17 +0200 Subject: [PATCH 09/20] - turned redundant code snippets into their own small functions getIndx() and createMapping() within the function postMix2posteriorList() - tidied and restructured code --- R/posterior.R | 76 +++++++++++++++++++++++------------------- tests/testthat/setup.R | 18 +++++----- 2 files changed, 51 insertions(+), 43 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index c4b4f23..7f8e309 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -231,42 +231,59 @@ postMix2posteriorList <- function ( checkmate::assert_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) + getIndx <- function (x, y) which(comp_mat_ind[, x] == comp_indx[[x]][y]) + + createMapping <- function (prior_list) { + + n_comps <- unlist(lapply(prior_list, ncol)) + comp_indx <- lapply(seq_along(prior_list), function (x) seq_len(n_comps[x])) + + return (comp_indx) + + } + # create mapping - args <- createMapping(prior_list) - comp_ind <- do.call("expand.grid", args) + comp_indx <- createMapping(prior_list) + comp_mat_ind <- do.call("expand.grid", comp_indx) # map posterior information - posterior_weight <- lapply(1:length(prior_list), - function (x) sapply(1:length(args[[x]]), - function (y) sum(unlist(posterior$weights[which(comp_ind[ ,x]==args[[x]][y])])))) + posterior_weight <- lapply(seq_along(prior_list), function (x) + sapply(seq_along(comp_indx[[x]]), function (y) + sum(unlist(posterior_list$weights[getIndx(x, y)])))) - posterior_mean <- lapply(1:length(prior_list), - function (x) sapply(1:length(args[[x]]), - function (y) (Reduce("+", posterior$mean[which(comp_ind[ ,x]==args[[x]][y])]) / length(posterior$mean[which(comp_ind[,x]==args[[x]][y])]))[x, ])) + posterior_mean <- lapply(seq_along(prior_list), function (x) + sapply(seq_along(comp_indx[[x]]), function (y) + mean(do.call(cbind, posterior_list$mean[getIndx(x, y)])[x, ]))) - posterior_sd <- lapply(1:length(prior_list), - function (x) sapply(1:length(args[[x]]), - function (y) (Reduce("+", lapply(posterior$covmat, diag)[which(comp_ind[ ,x]==args[[x]][y])]) / length(lapply(posterior$covmat, diag)[which(comp_ind[ ,x]==args[[x]][y])]))[x])) + posterior_sd <- lapply(seq_along(prior_list), function (x) + sapply(seq_along(comp_indx[[x]]), function (y) + mean(do.call(rbind, lapply(posterior_list$covmat, diag)[getIndx(x, y)])[, x]))) - combined_vectors <- mapply(function(x, y, z) Map(c, x, y, z), posterior_weight, posterior_mean, posterior_sd, SIMPLIFY = FALSE) + combined_vectors <- mapply(function (x, y, z) + Map(c, x, y, z), posterior_weight, posterior_mean, posterior_sd, + SIMPLIFY = FALSE) # create posterior list - posterior_output <- lapply(1:length(combined_vectors), - function(x) do.call(RBesT::mixnorm, c(combined_vectors[[x]], sigma = sigma(prior_list[[x]])))) + posterior_list_RBesT <- lapply(seq_along(combined_vectors), function (x) + do.call(RBesT::mixnorm, + c(combined_vectors[[x]], sigma = sigma(prior_list[[x]])))) - # TODO: talk about - #colnames(posterior_output$Ctr)[max(length(colnames(posterior_output$Ctr)))] <- "robust" - - names(posterior_output) <- c("Ctr", paste0("DG_", - seq_along(posterior_output[-1]))) - - class(posterior_output) <- "postList" - - attr(posterior_output, "ess") <- calcEss(calc_ess, posterior_output) + ## fix component names + names(posterior_list_RBesT) <- names(prior_list) + comp_names <- lapply(prior_list, colnames) + for (i in seq_along(posterior_list_RBesT)) { + + colnames(posterior_list_RBesT[[i]]) <- comp_names[[i]] + + } + rm(i) - attr(posterior_output, "full covariance matrices") <- posterior_list$covmat + ## set attributes + class(posterior_list_RBesT) <- "postList" + attr(posterior_list_RBesT, "ess") <- calcEss(calc_ess, posterior_list_RBesT) + attr(posterior_list_RBesT, "covariance matrices") <- posterior_list$covmat - return(posterior_output) + return (posterior_list_RBesT) } @@ -289,12 +306,3 @@ calcEss <- function(calc_ess, posterior_output) { } -createMapping <- function(prior_list) { - - n_dg <- length(prior_list) - n_comps <- unlist(lapply(prior_list, ncol)) - args <- lapply(1:n_dg, function (x) 1:n_comps[x]) - - return(args) - -} \ No newline at end of file diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index c4447db..d1c43fc 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -158,13 +158,13 @@ crit_pval = getCritProb( # Create covmat test case mixnorm_test <- mixnorm(comp1=c(0.2,2,3), comp2=c(0.2,5,6), comp3=c(0.2,8,9), comp4=c(0.2,11,12), robust=c(0.2,14,15), sigma=9.734) -mixnorm_DG1 <- mixnorm(comp1=c(0.5,2,3), comp2=c(0.5,3,6), sigma=9.651) +mixnorm_DG1 <- mixnorm(comp1=c(0.5,2,3), comp2=c(0.5,3,6), sigma=9.651) -mixnorm_DG2 <- mixnorm(comp1=c(1,8,2), sigma=9.932) +mixnorm_DG2 <- mixnorm(comp1=c(1,8,2), sigma=9.932) -mixnorm_DG3 <- mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) +mixnorm_DG3 <- mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) -mixnorm_DG4 <- mixnorm(comp1=c(1/4,10,3), comp2=c(1/4,3,6), comp3=c(1/4,4,7), comp4=c(1/4,9,8), sigma=9.236) +mixnorm_DG4 <- mixnorm(comp1=c(1/4,10,3), comp2=c(1/4,3,6), comp3=c(1/4,4,7), comp4=c(1/4,9,8), sigma=9.236) prior_list <- vector("list", 5) prior_list[[1]] <- mixnorm_test @@ -178,11 +178,11 @@ names(prior_list) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") mu_hat <- c(10, 20, 30, 40, 50) se_hat_vector <- matrix(c(3.11, 1.76, 0.38, 0.93, 1.66), nrow = 5, ncol = 1) -se_hat_matrix <- matrix(c(3.11, 0.01, 0.01, 0.01, 0.01, - 0.01, 1.76, 0.01, 0.01, 0.01, - 0.01, 0.01, 0.38, 0.01, 0.01, - 0.01, 0.01, 0.01, 0.93, 0.01, - 0.01, 0.01, 0.01, 0.01, 1.66), nrow = 5, ncol = 5) +se_hat_matrix <- matrix(c(3.11, 0.00, 0.00, 0.00, 0.00, + 0.00, 1.76, 0.00, 0.00, 0.00, + 0.00, 0.00, 0.38, 0.00, 0.00, + 0.00, 0.00, 0.00, 0.93, 0.00, + 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) data("metaData") From 5af1a756345a7d3628b6cce85a6ba51275c934f4 Mon Sep 17 00:00:00 2001 From: Schick Date: Fri, 7 Jun 2024 11:58:15 +0200 Subject: [PATCH 10/20] - removed setup with clinDR data, since the prior list can now be defined manually --- tests/testthat/setup.R | 31 +++---------------------------- 1 file changed, 3 insertions(+), 28 deletions(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index d1c43fc..dc4c76d 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -164,7 +164,7 @@ mixnorm_DG2 <- mixnorm(comp1=c(1,8,2), sigma=9.932) mixnorm_DG3 <- mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) -mixnorm_DG4 <- mixnorm(comp1=c(1/4,10,3), comp2=c(1/4,3,6), comp3=c(1/4,4,7), comp4=c(1/4,9,8), sigma=9.236) +mixnorm_DG4 <- mixnorm(comp1=c(1/3,10,3), comp2=c(1/3,3,6), comp3=c(1/3,4,7), sigma=9.236) prior_list <- vector("list", 5) prior_list[[1]] <- mixnorm_test @@ -184,33 +184,8 @@ se_hat_matrix <- matrix(c(3.11, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.93, 0.00, 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) -data("metaData") - -dataset_covmat <- filter(as.data.frame(metaData), bname == "BRINTELLIX") - -histcontrol_covmat <- filter( - dataset_covmat, - dose == 0, - primtime == 8, - indication == "MAJOR DEPRESSIVE DISORDER", - protid != 5) - -hist_data_covmat <- data.frame( - trial = histcontrol_covmat$nctno, - est = histcontrol_covmat$rslt, - se = histcontrol_covmat$se, - sd = histcontrol_covmat$sd, - n = histcontrol_covmat$sampsize) - -dose_levels_covmat <- c(0, 2.5, 5, 7.5, 10) - -prior_list_covmat <- getPriorList( - hist_data = hist_data_covmat, - dose_levels = dose_levels_covmat, - robust_weight = 0.3) - -posterior_covmat <- getPosterior( - prior_list = prior_list_covmat, +posterior <- getPosterior( + prior_list = prior_list, mu_hat = mu_hat, S_hat = se_hat_matrix, calc_ess = FALSE From d6096ffcd6f799af04107bde5b71a4aba1248867 Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 7 Jun 2024 13:07:22 +0100 Subject: [PATCH 11/20] Additional comments around posterior --- tests/testthat/setup.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index dc4c76d..04e038d 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -177,7 +177,8 @@ names(prior_list) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") mu_hat <- c(10, 20, 30, 40, 50) se_hat_vector <- matrix(c(3.11, 1.76, 0.38, 0.93, 1.66), nrow = 5, ncol = 1) - +#se_hat_vector <- matrix(c(sqrt(3.11), sqrt(1.76), sqrt(0.38), sqrt(0.93), sqrt(1.66)), nrow = 5, ncol = 1) +#Please note that a match between the two approaches is only there in case we are using the sqrt. TBD whether this is what we want (and need a good documentation) or whether the sqrt (or ^2) is calculated internally se_hat_matrix <- matrix(c(3.11, 0.00, 0.00, 0.00, 0.00, 0.00, 1.76, 0.00, 0.00, 0.00, 0.00, 0.00, 0.38, 0.00, 0.00, From a20776381670fdac6592d2b2f114e70288a64f33 Mon Sep 17 00:00:00 2001 From: Bossert Date: Fri, 7 Jun 2024 13:19:30 +0100 Subject: [PATCH 12/20] Comment around mistake in posteriorI function. --- R/posterior.R | 2 +- tests/testthat/setup.R | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/R/posterior.R b/R/posterior.R index 7f8e309..7661b9c 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -145,7 +145,7 @@ getPosteriorI <- function( class(post_list) <- "postList" attr(post_list, "ess") <- calcEss(calc_ess, post_list) - + #SB: The following line is not correct and needs to be changed by using the se information from the post_list object attr(post_list, "full covariance matrices") <- replicate(length(prior_list)-1, diag(c(se_hat)), simplify = FALSE) #names(attr(post_list, "full covariance matrices")) <- c("comp1", "comp2", "comp3", "robust") diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 04e038d..eb85f1f 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -185,6 +185,12 @@ se_hat_matrix <- matrix(c(3.11, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.93, 0.00, 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) +#se_hat_matrix <- matrix(c(3.11, 0.10, 0.20, 0.00, 0.00, + # 0.10, 1.76, 0.00, 0.00, 0.00, + # 0.20, 0.00, 0.38, 0.00, 0.00, + # 0.00, 0.00, 0.00, 0.93, 0.00, + # 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) +#We also need to include tests for cases with of diagonal elements posterior <- getPosterior( prior_list = prior_list, mu_hat = mu_hat, From 833f5de44ef357f9cdfc57f8d05327e0bcb5db18 Mon Sep 17 00:00:00 2001 From: Schick Date: Wed, 26 Jun 2024 11:33:28 +0200 Subject: [PATCH 13/20] - reworked if conditions to check if input is matrix or vector in getPosterior() - updated how the covariance matrices are collected within getPosteriorI() and postMix2posteriorList() - createMapping() is now its own function to allow use in both getPosteriorI() and postMix2posteriorList() for mapping of covariance matrices --- R/posterior.R | 68 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 51 insertions(+), 17 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index 7661b9c..efe2226 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -47,7 +47,27 @@ getPosterior <- function( if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) { - if (dim(S_hat)[2] > 1 && dim(S_hat)[2] == length(prior_list)) { + if (is.matrix(S_hat)) { + + is_matrix_S_hat <- TRUE + + } else if (is.vector(S_hat)) { + + is_matrix_S_hat <- FALSE + + se_hat <- S_hat + rm(S_hat) + + } else { + + stop("S_hat has to be either a vector or matrix") + + } + + if (is_matrix_S_hat) { + + stopifnot("dim of S_hat must equal length of prior list" = + dim(S_hat)[2] == length(prior_list)) prior_mix <- priorList2priorMix(prior_list) @@ -61,12 +81,12 @@ getPosterior <- function( prior_list = prior_list, calc_ess = calc_ess) - } else if (dim(S_hat)[2] == 1) { + } else if (!is_matrix_S_hat) { posterior_list <- getPosteriorI( prior_list = prior_list, mu_hat = mu_hat, - se_hat = S_hat, + se_hat = se_hat, calc_ess = calc_ess) } @@ -144,10 +164,18 @@ getPosteriorI <- function( names(post_list) <- names(prior_list) class(post_list) <- "postList" + comp_indx <- createMapping(prior_list) + comp_mat_ind <- do.call("expand.grid", comp_indx) + attr(post_list, "ess") <- calcEss(calc_ess, post_list) - #SB: The following line is not correct and needs to be changed by using the se information from the post_list object - attr(post_list, "full covariance matrices") <- replicate(length(prior_list)-1, diag(c(se_hat)), simplify = FALSE) - #names(attr(post_list, "full covariance matrices")) <- c("comp1", "comp2", "comp3", "robust") + + diagonals <- lapply(seq_along(comp_mat_ind[, 1]), function(x) { + + lapply(seq_along(comp_mat_ind[x, ]), function(y) return(post_list[[y]]["s", comp_mat_ind[x,y]])) + + }) + + attr(post_list, "full covariance matrices") <- lapply(seq_along(diagonals), function(x) diag(diagonals[[x]])) return (post_list) @@ -233,15 +261,6 @@ postMix2posteriorList <- function ( getIndx <- function (x, y) which(comp_mat_ind[, x] == comp_indx[[x]][y]) - createMapping <- function (prior_list) { - - n_comps <- unlist(lapply(prior_list, ncol)) - comp_indx <- lapply(seq_along(prior_list), function (x) seq_len(n_comps[x])) - - return (comp_indx) - - } - # create mapping comp_indx <- createMapping(prior_list) comp_mat_ind <- do.call("expand.grid", comp_indx) @@ -260,7 +279,7 @@ postMix2posteriorList <- function ( mean(do.call(rbind, lapply(posterior_list$covmat, diag)[getIndx(x, y)])[, x]))) combined_vectors <- mapply(function (x, y, z) - Map(c, x, y, z), posterior_weight, posterior_mean, posterior_sd, + Map(c, x, y, z), posterior_weight, posterior_mean, lapply(posterior_sd, sqrt), SIMPLIFY = FALSE) # create posterior list @@ -281,7 +300,14 @@ postMix2posteriorList <- function ( ## set attributes class(posterior_list_RBesT) <- "postList" attr(posterior_list_RBesT, "ess") <- calcEss(calc_ess, posterior_list_RBesT) - attr(posterior_list_RBesT, "covariance matrices") <- posterior_list$covmat + + diagonals <- lapply(seq_along(comp_mat_ind[, 1]), function(x) { + + lapply(seq_along(comp_mat_ind[x, ]), function(y) return(posterior_list_RBesT[[y]]["s", comp_mat_ind[x,y]])) + + }) + + attr(posterior_list_RBesT, "covariance matrices") <- lapply(seq_along(diagonals), function(x) diag(diagonals[[x]])) return (posterior_list_RBesT) @@ -306,3 +332,11 @@ calcEss <- function(calc_ess, posterior_output) { } +createMapping <- function (prior_list) { + + n_comps <- unlist(lapply(prior_list, ncol)) + comp_indx <- lapply(seq_along(prior_list), function (x) seq_len(n_comps[x])) + + return (comp_indx) + +} From 3521e08d3c6c146cc5eb72065c5149d9034942bc Mon Sep 17 00:00:00 2001 From: Schick Date: Wed, 26 Jun 2024 11:38:24 +0200 Subject: [PATCH 14/20] - updated setup.R to include a defined covariance matrix with all values not zero - added new tests to test-posterior.R to cover scenarios with a matrix without zeros - added new tests to test-posterior.R for equal results for a scenario with se_hat input as vector and S_hat input as matrix - created new function test_covmatrix_symmetry() in test-posterior.R to avoid redundant code --- tests/testthat/setup.R | 37 +++++---- tests/testthat/test-posterior.R | 132 +++++++++++++++++++++----------- 2 files changed, 108 insertions(+), 61 deletions(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index eb85f1f..9b1d1b5 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -176,24 +176,31 @@ prior_list[[5]] <- mixnorm_DG4 names(prior_list) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") mu_hat <- c(10, 20, 30, 40, 50) -se_hat_vector <- matrix(c(3.11, 1.76, 0.38, 0.93, 1.66), nrow = 5, ncol = 1) -#se_hat_vector <- matrix(c(sqrt(3.11), sqrt(1.76), sqrt(0.38), sqrt(0.93), sqrt(1.66)), nrow = 5, ncol = 1) -#Please note that a match between the two approaches is only there in case we are using the sqrt. TBD whether this is what we want (and need a good documentation) or whether the sqrt (or ^2) is calculated internally -se_hat_matrix <- matrix(c(3.11, 0.00, 0.00, 0.00, 0.00, - 0.00, 1.76, 0.00, 0.00, 0.00, - 0.00, 0.00, 0.38, 0.00, 0.00, - 0.00, 0.00, 0.00, 0.93, 0.00, - 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) - -#se_hat_matrix <- matrix(c(3.11, 0.10, 0.20, 0.00, 0.00, - # 0.10, 1.76, 0.00, 0.00, 0.00, - # 0.20, 0.00, 0.38, 0.00, 0.00, - # 0.00, 0.00, 0.00, 0.93, 0.00, - # 0.00, 0.00, 0.00, 0.00, 1.66), nrow = 5, ncol = 5) -#We also need to include tests for cases with of diagonal elements +se_hat_vector <- c(1.0, 3.0, 5.0, 9.0, 6.0) +se_hat_vector_sqrt <- c(sqrt(1), sqrt(3), sqrt(5), sqrt(9), sqrt(6)) + +se_hat_matrix <- matrix(c(1.00, 0.00, 0.00, 0.00, 0.00, + 0.00, 3.00, 0.00, 0.00, 0.00, + 0.00, 0.00, 5.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 9.00, 0.00, + 0.00, 0.00, 0.00, 0.00, 6.00), nrow = 5, ncol = 5) + +se_hat_matrix2 <- matrix(c(1.00, 0.10, 0.20, 0.30, 0.40, + 0.10, 3.00, 0.10, 0.20, 0.30, + 0.20, 0.10, 5.00, 0.10, 0.20, + 0.30, 0.20, 0.10, 9.00, 0.10, + 0.40, 0.30, 0.20, 0.10, 6.00), nrow = 5, ncol = 5) + posterior <- getPosterior( prior_list = prior_list, mu_hat = mu_hat, S_hat = se_hat_matrix, calc_ess = FALSE ) + +posterior_noZero <- getPosterior( + prior_list = prior_list, + mu_hat = mu_hat, + S_hat = se_hat_matrix2, + calc_ess = FALSE +) diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R index ddca125..3ece016 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -1,3 +1,23 @@ +# Function to test if covariance matrix is symmerical +test_covmatrix_symmetry <- function ( + + posterior, + mu_hat + + ) { + + lapply(1:(length(mu_hat)-1), function(x) { + + expect_equal(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)], + posterior$covmat$Comp1[row(posterior$covmat$Comp1) == col(posterior$covmat$Comp1) + x]) + + expect_equal(length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)]), + length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)])) + + }) + +} + test_that("getPosterior works correctly", { dummy_data <- getModelData(data, names(mods)[1]) @@ -106,35 +126,22 @@ test_that("getPostCombsI returns an object with correct attributes", { expect_true(all(result$vars == c(4, 4))) }) -test_that("getPriorMix works correctly", { +test_that("priorList2priorMix works correctly", { # function call without parameters - expect_error(createPriorMix()) + expect_error(priorList2priorMix()) # test priorList2priorMix function - prior_mix <- priorList2priorMix(prior_list_covmat) + prior_mix <- priorList2priorMix(prior_list) expect_type(prior_mix, "list") expect_length(prior_mix, 3) }) -test_that("mvpostmix works correctly", { - - prior_mix <- priorList2priorMix(prior_list_covmat) - - # test mvpostmix function - expect_error(DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_vector)) - expect_no_error(DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_matrix)) - posterior <- DoseFinding::mvpostmix(prior_mix, mu_hat, se_hat_matrix) - expect_type(posterior, "list") - expect_length(posterior, 3) - -}) - -test_that("getPosteriorOutput works correctly", { +test_that("postMix2posteriorList works correctly", { # create posterior with matrix - prior_mix <- priorList2priorMix(prior_list_covmat) + prior_mix <- priorList2priorMix(prior_list) posterior <- DoseFinding::mvpostmix( priormix = prior_mix, @@ -142,44 +149,77 @@ test_that("getPosteriorOutput works correctly", { S_hat = se_hat_matrix) # create posterior with vector - posterior_vector <- getPosteriorI( - prior_list = prior_list_covmat, + posterior_vector <- getPosterior( + prior_list = prior_list, mu_hat = mu_hat, - se_hat = se_hat_vector, + S_hat = se_hat_vector, calc_ess = FALSE ) - # test postMix2posteriorList function - posterior_list <- postMix2posteriorList(posterior, prior_list_covmat, calc_ess = FALSE) + ### test postMix2posteriorList function - se_hat_matrix only zeros + posterior_list <- postMix2posteriorList(posterior, prior_list, calc_ess = FALSE) expect_type(posterior_list, "list") expect_s3_class(posterior_list, "postList") - lapply(1:(length(dose_levels_covmat)-1), function(x) { - - expect_equal(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)], - posterior$covmat$Comp1[row(posterior$covmat$Comp1) == col(posterior$covmat$Comp1) + x]) - - expect_equal(length(which(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)] >= 0)), - length(posterior$covmat$Comp1[row(posterior$covmat$Comp1) + x == col(posterior$covmat$Comp1)])) - - }) - - lapply(seq_along(prior_list), function(i) { - - expect_equal(length(posterior_weight[[i]]), length(posterior_mean[[i]])) - expect_equal(length(posterior_mean[[i]]), length(posterior_sd[[i]])) - - sapply(seq_along(posterior_weight[[i]]), function(j) { - - expect_length(combined_vectors[[i]][[j]], 3) - - }) - - }) + test_covmatrix_symmetry(posterior, mu_hat) # compare posterior result object with matrix to object with vector expect_length(posterior_list, length(posterior_vector)) + expect_length(attr(posterior_list, "covariance matrices"), length(attr(posterior_vector, "full covariance matrices"))) expect_type(posterior_list, typeof(posterior_vector)) expect_s3_class(posterior_list, S3Class(posterior_vector)) -}) \ No newline at end of file + + ### test postMix2posteriorList function - se_hat_matrix not only zeros + mvpostmix_noZero <- DoseFinding::mvpostmix( + priormix = prior_mix, + mu_hat = mu_hat, + S_hat = se_hat_matrix2) + + posterior_noZero <- getPosterior( + prior_list = prior_list, + mu_hat = mu_hat, + S_hat = se_hat_matrix2, + calc_ess = FALSE + ) + + expect_type(posterior_noZero, "list") + expect_s3_class(posterior_noZero, "postList") + + test_covmatrix_symmetry(mvpostmix_noZero, mu_hat) + + # compare posterior result object with matrix to object with vector + expect_length(posterior_noZero, length(posterior_vector)) + expect_length(attr(posterior_noZero, "covariance matrices"), length(attr(posterior_vector, "full covariance matrices"))) + expect_type(posterior_noZero, typeof(posterior_vector)) + expect_s3_class(posterior_noZero, S3Class(posterior_vector)) + + + ### test similarity of results for posterior from a matrix compared to posterior from a vector with square roots + se_hat <- c(1, 2, 3, 4, 5) + S_hat <- diag(se_hat) + + posterior_matrix_S <- getPosterior( + prior_list = prior_list, + mu_hat = mu_hat, + S_hat = S_hat, + calc_ess = FALSE + ) + + posterior_vector_se <- getPosterior( + prior_list = prior_list, + mu_hat = mu_hat, + S_hat = sqrt(se_hat), + calc_ess = FALSE + ) + + expect_equal(posterior_matrix_S$Ctr, posterior_vector_se$Ctr) + expect_equal(posterior_matrix_S$DG_1, posterior_vector_se$DG_1) + expect_equal(posterior_matrix_S$DG_2, posterior_vector_se$DG_2) + expect_equal(posterior_matrix_S$DG_3, posterior_vector_se$DG_3) + expect_equal(posterior_matrix_S$DG_4, posterior_vector_se$DG_4) + + + + +}) From 6480b733070a68ece120510bfe4b52d4a3e11310 Mon Sep 17 00:00:00 2001 From: Schick Date: Mon, 8 Jul 2024 09:15:02 +0200 Subject: [PATCH 15/20] - updated analysis_normal vignette to fit new input variables of getPosterior() --- vignettes/analysis_normal.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd index af6159f..aea87fb 100644 --- a/vignettes/analysis_normal.Rmd +++ b/vignettes/analysis_normal.Rmd @@ -162,7 +162,7 @@ As outlined in Fleischer et al. [Bayesian MCPMod. Pharmaceutical Statistics. 202 posterior <- getPosterior( prior = prior_list, mu_hat = new_trial$rslt, - se_hat = new_trial$se, + S_hat = new_trial$se, calc_ess = TRUE) summary(posterior) From 5cf8f919f906e2d59d3aa17ebfcaf124ee2f4072 Mon Sep 17 00:00:00 2001 From: Schick Date: Tue, 16 Jul 2024 11:11:28 +0200 Subject: [PATCH 16/20] - updated pkgdown.yml - changed input variable of getPosterior() in analyis_normal_noninformative vignette from se_hat to S_hat --- _pkgdown.yml | 1 + vignettes/analysis_normal_noninformative.qmd | 413 +++++++++++++++++++ 2 files changed, 414 insertions(+) create mode 100644 vignettes/analysis_normal_noninformative.qmd diff --git a/_pkgdown.yml b/_pkgdown.yml index 1b0d774..934cedf 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -12,4 +12,5 @@ articles: navbar: ~ contents: - analysis_normal + - analysis_normal_noninformative - Simulation_Example diff --git a/vignettes/analysis_normal_noninformative.qmd b/vignettes/analysis_normal_noninformative.qmd new file mode 100644 index 0000000..89c4488 --- /dev/null +++ b/vignettes/analysis_normal_noninformative.qmd @@ -0,0 +1,413 @@ +--- +title: "Analysis Example of Bayesian MCPMod for Continuous Data with non-informative prior" +format: + html: + fig-height: 3.5 + self-contained: true + toc: true + number-sections: true + bibliography: references.bib +vignette: > + %\VignetteIndexEntry{Analysis Example of Bayesian MCPMod for Continuous Data with non-informative prior} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{quarto::html} +--- + +```{r setup} +library(BayesianMCPMod) +library(RBesT) +library(clinDR) +library(dplyr) +library(tibble) +library(reactable) + +set.seed(7015) +``` + +```{r} +#| code-summary: setup +#| code-fold: true +#| message: false +#| warning: false + +#' Display Parameters Table +#' +#' This function generates a markdown table displaying the names and values of parameters +#' from a named list. +#' +#' @param named_list A named list where each name represents a parameter name and the list +#' element represents the parameter value. Date values in the list are automatically +#' converted to character strings for display purposes. +#' +#' @return Prints a markdown table with two columns: "Parameter Name" and "Parameter Values". +#' The function does not return a value but displays the table directly to the output. +#' +#' @importFrom knitr kable +#' @examples +#' params <- list("Start Date" = as.Date("2020-01-01"), +#' "End Date" = as.Date("2020-12-31"), +#' "Threshold" = 10) +#' display_params_table(params) +#' +#' @export +display_params_table <- function(named_list) { + display_table <- data.frame() + value_names <- data.frame() + for (i in 1:length(named_list)) { + # dates will display as numeric by default, so convert to char first + if (class(named_list[[i]]) == "Date") { + named_list[[i]] = as.character(named_list[[i]]) + } + if (!is.null(names(named_list[[i]]))) { + value_names <- rbind(value_names, paste(names(named_list[[i]]), collapse = ', ')) + } + values <- data.frame(I(list(named_list[[i]]))) + display_table <- rbind(display_table, values) + } + + round_numeric <- function(x, digits = 3) { + if (is.numeric(x)) { + return(round(x, digits)) + } else { + return(x) + } + } + + display_table[1] <- lapply(display_table[1], function(sublist) { + lapply(sublist, round_numeric) + }) + + class(display_table[[1]]) <- "list" + + if (nrow(value_names) == 0) { + knitr::kable( + cbind(names(named_list), display_table), + col.names = c("Name", "Value") + ) + } else { + knitr::kable( + cbind(names(named_list), value_names, display_table), + col.names = c("Name", "Value Labels", "Value") + ) + } +} +``` + +# Introduction + +This vignette demonstrates the application of the {BayesianMCPMod} package for +analyzing a phase 2 dose-finding trial using the Bayesian MCPMod approach. + +# Prior Specification + +Ideally, priors are grounded in historical data. This approach allows for the synthesis of +prior knowledge with current data, enhancing the accuracy of trial evaluations. + +The focus of this vignette is more generic, however. We specify weakly-informative +priors across all dose groups to allow the trial data to have a stronger influence on the analysis. + +```{r} +dose_levels <- c(0, 2.5, 5, 10) + +prior_list <- lapply(dose_levels, function(dose_group) { + RBesT::mixnorm(weak = c(w = 1, m = 0, s = 200), sigma = 10) +}) + +names(prior_list) <- c("Ctr", paste0("DG_", dose_levels[-1])) +``` + +# Dose-Response Model Shapes + +Candidate models are specified using the {DoseFinding} package. Models can be +parameterized using guesstimates or by directly providing distribution parameters. +Note that the linear candidate model does not require parameterization. + +[**Note:** The LinLog model is rarely used and not currently supported by `{BayesianMCPMod}`.]{.aside} + +In the code below, the models are "guesstimated" using the `DoseFinding::guesst` function. +The `d` option usually takes a single value (a dose level), and the corresponding `p` +for the maximum effect achieved at `d`. + + +```{r} +# Guesstimate estimation +exp_guesst <- DoseFinding::guesst( + model = "exponential", + d = 5, p = 0.2, Maxd = max(dose_levels) +) +emax_guesst <- DoseFinding::guesst( + model = "emax", + d = 2.5, p = 0.9 +) +sigEmax_guesst <- DoseFinding::guesst( + model = "sigEmax", + d = c(2.5, 5), p = c(0.5, 0.95) +) +logistic_guesst <- DoseFinding::guesst( + model = "logistic", + d = c(5, 10), p = c(0.1, 0.85) +) +``` + +In some cases, you need to provide more information. For instance, `sigEmax` +requires a pair of `d` and `p` values, and `exponential` requires the specification of +the maximum dose for the trial (`Maxd`). + +[See the help files for model specifications by typing `?DoseFinding::guesst` in your console]{.aside} + + + +Of course, you can also specify the models directly on the parameter scale (without using `DoseFinding::guesst`). + +For example, you can get a betaMod model by specifying `delta1` and `delta2` +parameters (`scale` is assumed to be `1.2` of the maximum dose), or a quadratic model with the `delta2` parameter. + +```{r} +betaMod_params <- c(delta1 = 1, delta2 = 1) +quadratic_params <- c(delta2 = -0.1) +``` + +Now, we can go ahead and create a `Mods` object, which will be used in the remainder +of the vignette. + +```{r} +mods <- DoseFinding::Mods( + linear = NULL, + # guesstimate scale + exponential = exp_guesst, + emax = emax_guesst, + sigEmax = sigEmax_guesst, + logistic = logistic_guesst, + # parameter scale + betaMod = betaMod_params, + quadratic = quadratic_params, + # Options for all models + doses = dose_levels, + maxEff = -1, + placEff = -12.8 +) + +plot(mods) +``` + +The `mods` object we just created above contains the full model parameters, which can be helpful for +understanding how the guesstimates are translated onto the parameter scale. + +```{r} +display_params_table(mods) +``` + +And we can see the assumed treatment effects for the specified dose groups below: + +```{r} +knitr::kable(DoseFinding::getResp(mods, doses = dose_levels)) +``` + +## Trial Data + +We will use the trial with ct.gov number NCT00735709 as our phase 2 trial data, +available in the `{clinDR}` package [@nct00735709_2024a]. + +```{r} +data("metaData") + +trial_data <- dplyr::filter( + dplyr::filter(tibble::tibble(metaData), bname == "BRINTELLIX"), + primtime == 8, + indication == "MAJOR DEPRESSIVE DISORDER", + protid == 5 +) + +n_patients <- c(128, 124, 129, 122) +``` + +# Posterior Calculation + +In the first step of Bayesian MCPMod, the posterior is calculated by combining +the prior information with the estimated results of the trial [@fleischer_2022]. + +```{r} +posterior <- getPosterior( + prior_list = prior_list, + mu_hat = trial_data$rslt, + S_hat = trial_data$se, + calc_ess = TRUE +) + +knitr::kable(summary(posterior)) +``` + +# Bayesian MCPMod Test Step + +The testing step of Bayesian MCPMod is executed using a critical value on the +probability scale and a pseudo-optimal contrast matrix. + +The critical value is calculated using (re-estimated) contrasts for frequentist +MCPMod to ensure error control when using weakly-informative priors. + +A pseudo-optimal contrast matrix is generated based on the variability of the +posterior distribution (see [@fleischer_2022] for more details). + +```{r} +crit_pval <- getCritProb( + mods = mods, + dose_levels = dose_levels, + se_new_trial = trial_data$se, + alpha_crit_val = 0.05 +) + +contr_mat <- getContr( + mods = mods, + dose_levels = dose_levels, + sd_posterior = summary(posterior)[, 2] +) +``` + +Please note that there are different ways to derive the contrasts. +The following code shows the implementation of some of these ways but it is not +executed and the contrast specification above is used. + +```{r} +#| eval: false + +# i) the frequentist contrast +contr_mat_prior <- getContr( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + prior_list = prior_list) +# ii) re-estimated frequentist contrasts +contr_mat_prior <- getContr( + mods = mods, + dose_levels = dose_levels, + se_new_trial = trial_data$se) +# iii) Bayesian approach using number of patients for new trial and prior distribution +contr_mat_prior <- getContr( + mods = mods, + dose_levels = dose_levels, + dose_weights = n_patients, + prior_list = prior_list) +``` + +The Bayesian MCP testing step is then executed: + +```{r} +BMCP_result <- performBayesianMCP( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval) +``` + +Summary information: + +```{r} +BMCP_result +``` + +The testing step is significant, indicating a non-flat dose-response shape. +All models are significant, with the `emax` model indicating the greatest deviation +from the null hypothesis. + +# Model Fitting and Visualization + +In the model fitting step the posterior distribution is used as basis. + +Both simplified and full fitting are performed. + +For the simplified fit, the multivariate normal distribution of the control group +is approximated and reduced by a one-dimensional normal distribution. + +The actual fit (on this approximated posterior distribution) is then performed +using generalized least squares criterion. In contrast, for the full fit, the +non-linear optimization problem is addressed via the Nelder-Mead algorithm +[@neldermead_2024a] implemented by the `{nloptr}` package. + +The output of the fit includes information about the predicted effects for the +included dose levels, the generalized AIC, and the corresponding weights. + +For the considered case, the simplified and the full fit are very similar, so we +present the full fit. + +```{r} +# If simple = TRUE, uses approx posterior +# Here we use complete posterior distribution +fit <- getModelFits( + models = mods, + dose_levels = dose_levels, + posterior = posterior, + simple = FALSE) +``` + +Estimates for dose levels not included in the trial: + +```{r} +display_params_table(stats::predict(fit, doses = c(0, 2.5, 4, 5, 7, 10))) +``` + +Plots of fitted dose-response models and an AIC-based average model: + +```{r} +plot(fit) +``` + +To assess the uncertainty, one can additionally visualize credible bands +(orange shaded areas, default levels are 50% and 95%). + +These credible bands are calculated with a bootstrap method as follows: + +- Samples from the posterior distribution are drawn and for every sample the +simplified fitting step and a prediction is performed. + +- These predictions are then used to identify and visualize the specified quantiles. + +```{r} +plot(fit, cr_bands = TRUE) +``` + +The bootstrap based quantiles can also be directly calculated via the +`getBootstrapQuantiles()` function. + +For this example, only 6 quantiles are bootstrapped for each model fit. + +```{r} +bootstrap_quantiles <- getBootstrapQuantiles( + model_fits = fit, + quantiles = c(0.025, 0.5, 0.975), + doses = c(0, 2.5, 4, 5, 7, 10), + n_samples = 6 +) +``` + +```{r} +#| code-fold: true +reactable::reactable( + data = bootstrap_quantiles, + groupBy = "models", + columns = list( + doses = colDef(aggregate = "count", format = list(aggregated = colFormat(suffix = " doses"))), + "2.5%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 4))), + "50%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 4))), + "97.5%" = colDef(aggregate = "mean", format = list(aggregated = colFormat(prefix = "mean = ", digits = 2), cell = colFormat(digits = 4))) + ) +) +``` + +**Technical note:** The median quantile of the bootstrap based procedure is not +necessary similar to the main model fit, as they are derived via different procedures. + +The main fit (black line) minimizes residuals for the posterior distribution, +while the bootstrap median is the median fit of random sampling. + +# Additional note + +Testing and modeling can also be combined via `performBayesianMCPMod()`, +but this is not run here. + +```{r} +#| eval: false +performBayesianMCPMod( + posterior_list = posterior, + contr = contr_mat, + crit_prob_adj = crit_pval, + simple = FALSE) +``` From 4e4f001f24587e7d3f307da6594aedb05bb91482 Mon Sep 17 00:00:00 2001 From: Schick Date: Tue, 16 Jul 2024 11:25:35 +0200 Subject: [PATCH 17/20] - added package prefix --- tests/testthat/setup.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 9b1d1b5..bc28087 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -156,15 +156,15 @@ crit_pval = getCritProb( # ) # Create covmat test case -mixnorm_test <- mixnorm(comp1=c(0.2,2,3), comp2=c(0.2,5,6), comp3=c(0.2,8,9), comp4=c(0.2,11,12), robust=c(0.2,14,15), sigma=9.734) +mixnorm_test <- RBesT::mixnorm(comp1=c(0.2,2,3), comp2=c(0.2,5,6), comp3=c(0.2,8,9), comp4=c(0.2,11,12), robust=c(0.2,14,15), sigma=9.734) -mixnorm_DG1 <- mixnorm(comp1=c(0.5,2,3), comp2=c(0.5,3,6), sigma=9.651) +mixnorm_DG1 <- RBesT::mixnorm(comp1=c(0.5,2,3), comp2=c(0.5,3,6), sigma=9.651) -mixnorm_DG2 <- mixnorm(comp1=c(1,8,2), sigma=9.932) +mixnorm_DG2 <- RBesT::mixnorm(comp1=c(1,8,2), sigma=9.932) -mixnorm_DG3 <- mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) +mixnorm_DG3 <- RBesT::mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0.5,1), sigma=9.134) -mixnorm_DG4 <- mixnorm(comp1=c(1/3,10,3), comp2=c(1/3,3,6), comp3=c(1/3,4,7), sigma=9.236) +mixnorm_DG4 <- RBesT::mixnorm(comp1=c(1/3,10,3), comp2=c(1/3,3,6), comp3=c(1/3,4,7), sigma=9.236) prior_list <- vector("list", 5) prior_list[[1]] <- mixnorm_test From 0fe144f013b09b3bf9f50424a81be79281d9c5c2 Mon Sep 17 00:00:00 2001 From: Schick Date: Mon, 22 Jul 2024 14:26:34 +0200 Subject: [PATCH 18/20] - updated tests to only allow input of prior_lists as objects of RBesT package - renamed prior_list in setup.R --- tests/testthat/setup.R | 18 +++++++------- tests/testthat/test-posterior.R | 42 ++++++++++++++++++++------------- 2 files changed, 35 insertions(+), 25 deletions(-) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index bc28087..3b7ae13 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -166,14 +166,14 @@ mixnorm_DG3 <- RBesT::mixnorm(comp1=c(1/3,6,8), comp2=c(1/3,7,9), comp3=c(1/3,0 mixnorm_DG4 <- RBesT::mixnorm(comp1=c(1/3,10,3), comp2=c(1/3,3,6), comp3=c(1/3,4,7), sigma=9.236) -prior_list <- vector("list", 5) -prior_list[[1]] <- mixnorm_test -prior_list[[2]] <- mixnorm_DG1 -prior_list[[3]] <- mixnorm_DG2 -prior_list[[4]] <- mixnorm_DG3 -prior_list[[5]] <- mixnorm_DG4 +prior_list_matrix <- vector("list", 5) +prior_list_matrix[[1]] <- mixnorm_test +prior_list_matrix[[2]] <- mixnorm_DG1 +prior_list_matrix[[3]] <- mixnorm_DG2 +prior_list_matrix[[4]] <- mixnorm_DG3 +prior_list_matrix[[5]] <- mixnorm_DG4 -names(prior_list) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") +names(prior_list_matrix) <- c("Ctr","DG_1","DG_2","DG_3","DG_4") mu_hat <- c(10, 20, 30, 40, 50) se_hat_vector <- c(1.0, 3.0, 5.0, 9.0, 6.0) @@ -192,14 +192,14 @@ se_hat_matrix2 <- matrix(c(1.00, 0.10, 0.20, 0.30, 0.40, 0.40, 0.30, 0.20, 0.10, 6.00), nrow = 5, ncol = 5) posterior <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = se_hat_matrix, calc_ess = FALSE ) posterior_noZero <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = se_hat_matrix2, calc_ess = FALSE diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R index 3ece016..7145097 100644 --- a/tests/testthat/test-posterior.R +++ b/tests/testthat/test-posterior.R @@ -20,11 +20,21 @@ test_covmatrix_symmetry <- function ( test_that("getPosterior works correctly", { dummy_data <- getModelData(data, names(mods)[1]) + prior_list_noRBesT <- list(1,2,3,4) + + expect_error(getPosterior( + prior_list = prior_list_noRBesT, + mu_hat = mu_hat, + S_hat = se_hat_matrix, + calc_ess = FALSE + )) # Test getPosterior function posterior_list <- getPosterior( - data = getModelData(data, names(mods)[1]), - prior_list = prior_list + prior_list = prior_list_matrix, + mu_hat = mu_hat, + S_hat = se_hat_vector, + calc_ess = FALSE ) expect_type(posterior_list, "list") expect_s3_class(posterior_list, "postList") @@ -78,21 +88,21 @@ test_that("getPriorList input parameters do work as intented", { test_that("getPosteriorI works correctly", { # Prepare test data and parameters data_i <- data.frame( - dose = c(0, 1, 2, 3), - response = c(10, 20, 30, 40) + dose = c(0, 1, 2, 3, 4), + response = c(10, 20, 30, 40, 50) ) - prior_list <- list(1, 2, 3, 4) - mu_hat <- c(10, 20, 30, 40) - se_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1) + #prior_list <- list(1, 2, 3, 4) + #mu_hat <- c(10, 20, 30, 40) + #se_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1) # Test getPosteriorI function - post_list <- getPosteriorI(data_i, prior_list, mu_hat, se_hat) + post_list <- getPosteriorI(data_i = data_i, prior_list = prior_list_matrix, mu_hat = mu_hat, se_hat = se_hat_vector) expect_type(post_list, "list") expect_s3_class(post_list, "postList") # Test mu_hat and sd_hat both null branch - post_list <- getPosteriorI(data_i, prior_list, NULL, NULL) + post_list <- getPosteriorI(data_i, prior_list_matrix, NULL, NULL) expect_type(post_list, "list") expect_s3_class(post_list, "postList") }) @@ -132,7 +142,7 @@ test_that("priorList2priorMix works correctly", { expect_error(priorList2priorMix()) # test priorList2priorMix function - prior_mix <- priorList2priorMix(prior_list) + prior_mix <- priorList2priorMix(prior_list_matrix) expect_type(prior_mix, "list") expect_length(prior_mix, 3) @@ -141,7 +151,7 @@ test_that("priorList2priorMix works correctly", { test_that("postMix2posteriorList works correctly", { # create posterior with matrix - prior_mix <- priorList2priorMix(prior_list) + prior_mix <- priorList2priorMix(prior_list_matrix) posterior <- DoseFinding::mvpostmix( priormix = prior_mix, @@ -150,14 +160,14 @@ test_that("postMix2posteriorList works correctly", { # create posterior with vector posterior_vector <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = se_hat_vector, calc_ess = FALSE ) ### test postMix2posteriorList function - se_hat_matrix only zeros - posterior_list <- postMix2posteriorList(posterior, prior_list, calc_ess = FALSE) + posterior_list <- postMix2posteriorList(posterior, prior_list_matrix, calc_ess = FALSE) expect_type(posterior_list, "list") expect_s3_class(posterior_list, "postList") @@ -177,7 +187,7 @@ test_that("postMix2posteriorList works correctly", { S_hat = se_hat_matrix2) posterior_noZero <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = se_hat_matrix2, calc_ess = FALSE @@ -200,14 +210,14 @@ test_that("postMix2posteriorList works correctly", { S_hat <- diag(se_hat) posterior_matrix_S <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = S_hat, calc_ess = FALSE ) posterior_vector_se <- getPosterior( - prior_list = prior_list, + prior_list = prior_list_matrix, mu_hat = mu_hat, S_hat = sqrt(se_hat), calc_ess = FALSE From 820b22bc20c1ff5b0d76aa35a51f72d1c88e72e5 Mon Sep 17 00:00:00 2001 From: Schick Date: Mon, 22 Jul 2024 14:27:26 +0200 Subject: [PATCH 19/20] - added new if-clause to getPosterior() to check if prior_list is object of RBesT package --- R/posterior.R | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/R/posterior.R b/R/posterior.R index efe2226..5aa6d2d 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -45,6 +45,11 @@ getPosterior <- function( checkmate::check_vector(S_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(S_hat, null.ok = TRUE, lower = 0, upper = Inf) + is_matrix_S_hat <- FALSE + + stopifnot("prior_list must be an object of RBesT package" = + all(sapply(prior_list, function(x) is(x, "normMix") | is(x, "betaMix") | is(x, "mix")))) + if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) { if (is.matrix(S_hat)) { @@ -142,9 +147,10 @@ getPosteriorI <- function( } else if (!is.null(mu_hat) && !is.null(se_hat)) { stopifnot("m_hat length must match number of dose levels" = - length(prior_list) == length(mu_hat), - "se_hat length must match number of dose levels" = - length(prior_list) == length(se_hat)) + length(prior_list) == length(mu_hat)) + # , + # "se_hat length must match number of dose levels" = + # length(prior_list) == length(se_hat)) } else { From aa5881de0a81d2847e29d7711251714f8c18db47 Mon Sep 17 00:00:00 2001 From: Andersen Date: Tue, 23 Jul 2024 13:16:33 +0200 Subject: [PATCH 20/20] adaptations regarding the new S_hat functionality as well as addind Jonas as ctb and adjustign docu --- DESCRIPTION | 3 +- R/BMCPMod.R | 344 +++++++++++++++++------------------ R/posterior.R | 248 ++++++++++++------------- man/assessDesign.Rd | 2 +- man/getContr.Rd | 2 +- man/getCritProb.Rd | 2 +- man/getPosterior.Rd | 12 +- man/performBayesianMCP.Rd | 8 +- man/performBayesianMCPMod.Rd | 6 +- 9 files changed, 314 insertions(+), 313 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index c9ddc6b..cfbd7d0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,7 +7,8 @@ Authors@R: c( person("Stephan", "Wojciekowski", , "stephan.wojciekowski@boehringer-ingelheim.com", role = c("aut", "cre")), person("Lars", "Andersen", , "lars.andersen@boehringer-ingelheim.com", role = "aut"), person("Sebastian", "Bossert", , "sebastian.bossert@boehringer-ingelheim.com", role = "aut"), - person("Steven", "Brooks", , "steven.brooks@boehringer-ingelheim.com", role = "ctb") + person("Steven", "Brooks", , "steven.brooks@boehringer-ingelheim.com", role = "ctb"), + person("Jonas", "Schick", , "jonas.schick@boehringer-ingelheim.com", role = "ctb") ) Description: Bayesian MCPMod (Fleischer et al. (2022) ) is an innovative method that improves the diff --git a/R/BMCPMod.R b/R/BMCPMod.R index 41f00e1..67db71c 100644 --- a/R/BMCPMod.R +++ b/R/BMCPMod.R @@ -4,20 +4,20 @@ #' #' @param n_patients Vector specifying the planned number of patients per dose group #' @param mods An object of class "Mods" as specified in the DoseFinding package. -#' @param prior_list A prior_list object specifying the utilized prior for the different dose groups -#' @param sd A positive value, specification of assumed sd +#' @param prior_list A prior_list object specifying the utilized prior for the different dose groups +#' @param sd A positive value, specification of assumed sd #' @param n_sim Number of simulations to be performed #' @param alpha_crit_val (Un-adjusted) Critical value to be used for the MCP testing step. Passed to the getCritProb() function for the calculation of adjusted critical values (on the probability scale). Default is 0.05. #' @param simple Boolean variable defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE. #' @param reestimate Boolean variable defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE #' @param contr An object of class 'optContr' as created by the getContr() function. Allows specification of a fixed contrasts matrix. Default NULL #' @param dr_means A vector, allows specification of individual (not model based) assumed effects per dose group. Default NULL -#' +#' #' @return Returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size -#' +#' #' @examples #' if (interactive()) { # takes typically > 5 seconds -#' +#' #' mods <- DoseFinding::Mods(linear = NULL, #' linlog = NULL, #' emax = c(0.5, 1.2), @@ -27,41 +27,41 @@ #' sd <- 12 #' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), #' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), -#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , +#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , #' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , #' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) #' n_patients <- c(40, 60, 60, 60, 60) -#' +#' #' success_probabilities <- assessDesign( #' n_patients = n_patients, #' mods = mods, #' prior_list = prior_list, #' sd = sd, #' n_sim = 1e2) # speed up exammple run time -#' +#' #' success_probabilities -#' +#' #' } -#' +#' #' @export assessDesign <- function ( - + n_patients, mods, prior_list, - + sd, - + n_sim = 1e3, alpha_crit_val = 0.05, simple = TRUE, reestimate = FALSE, - + contr = NULL, dr_means = NULL - + ) { - + checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE) checkmate::check_class(mods, classes = "Mods") checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE) @@ -70,66 +70,66 @@ assessDesign <- function ( checkmate::check_double(alpha_crit_val, lower = 0, upper = 1) checkmate::check_logical(simple) # TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric # this is not applicable at the moment - + dose_levels <- attr(mods, "doses") - + data <- simulateData( n_patients = n_patients, dose_levels = dose_levels, sd = sd, mods = mods, - n_sim = n_sim, + n_sim = n_sim, dr_means = dr_means) - + model_names <- colnames(data)[-c(1:3)] - + crit_prob_adj <- getCritProb( mods = mods, dose_levels = dose_levels, dose_weights = n_patients, alpha_crit_val = alpha_crit_val) - + if (!reestimate & is.null(contr)) { - + contr <- getContr( mods = mods, dose_levels = dose_levels, dose_weights = n_patients, prior_list = prior_list) - + } - + eval_design <- lapply(model_names, function (model_name) { - + posterior_list <- getPosterior( data = getModelData(data, model_name), prior_list = prior_list) - + if (reestimate & is.null(contr)) { - + post_sds <- sapply(posterior_list, function (post) summary(post)[, 2]) - + contr <- apply(post_sds, 2, function (post_sd) getContr( mods = mods, dose_levels = dose_levels, sd_posterior = post_sd)) - - } - + + } + b_mcp_mod <- performBayesianMCPMod( posterior_list = posterior_list, contr = contr, crit_prob_adj = crit_prob_adj, simple = simple) - + }) - + avg_success_rate <- mean(sapply(eval_design, function (bmcpmod) { attr(bmcpmod$BayesianMCP, "successRate") })) - + names(eval_design) <- model_names - + attr(eval_design, "avgSuccessRate") <- avg_success_rate attr(eval_design, "placEff") <- ifelse(test = is.null(dr_means), yes = attr(mods, "placEff"), @@ -139,26 +139,26 @@ assessDesign <- function ( no = diff(range(dr_means))) attr(eval_design, "sampleSize") <- n_patients attr(eval_design, "priorESS") <- round(getESS(prior_list), 1) - + return (eval_design) - + } #' @title getContr -#' +#' #' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr() of the DoseFinding package. #' Hereby 4 different options can be distinguished that are automatically executed based on the input that is provided -#' 1) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated. +#' 1) Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated. #' In detail, in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for #' each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on #' pseudo-optimal contrasts via regular MCPMod are calculated from the -#' regular MCPMod for these specific weights +#' regular MCPMod for these specific weights #' 2) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the #' regular MCPMod for these specific weights #' 3) Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated #' 4) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the -#' regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances -#' +#' regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances +#' #' @param mods An object of class 'Mods' as created by the function 'DoseFinding::Mods()' #' @param dose_levels Vector containing the different dosage levels. #' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option 1) is planned these two inputs should be provided on the same scale (e.g. patient numbers). Default NULL @@ -166,7 +166,7 @@ assessDesign <- function ( #' @param sd_posterior A vector of positive values with information about the variability of the posterior distribution, only required for Option 3. Default NULL #' @param se_new_trial A vector of positive values with information about the observed variability, only required for Option 4. Default NULL -#' +#' #' @examples #' dose_levels <- c(0, 0.5, 2, 4, 8) #' mods <- DoseFinding::Mods( @@ -177,98 +177,98 @@ assessDesign <- function ( #' doses = dose_levels, #' maxEff = 6) #' sd_posterior <- c(2.8, 3, 2.5, 3.5, 4) -#' +#' #' contr_mat <- getContr( #' mods = mods, #' dose_levels = dose_levels, -#' sd_posterior = sd_posterior) -#' +#' sd_posterior = sd_posterior) +#' #' @return An object of class 'optContr' as provided by the function 'DoseFinding::optContr()'. -#' +#' #' @export getContr <- function ( - + mods, dose_levels, dose_weights = NULL, prior_list = NULL, sd_posterior = NULL, se_new_trial = NULL - + ) { - + checkmate::check_class(mods, classes = "Mods") checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(attr(mods, "doses"))) checkmate::check_double(dose_weights, any.missing = FALSE, len = length(attr(mods, "doses"))) checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE) - + # frequentist & re-estimation - if (!is.null(se_new_trial) & + if (!is.null(se_new_trial) & is.null(dose_weights) & is.null(prior_list) & is.null(sd_posterior)) { - + w <- NULL S <- diag((se_new_trial)^2) - + # frequentist & no re-estimation - } else if (!is.null(dose_weights) & + } else if (!is.null(dose_weights) & is.null(se_new_trial) & is.null(prior_list) & is.null(sd_posterior)) { - + w <- dose_weights S <- NULL - + # Bayesian & re-estimation - } else if (!is.null(sd_posterior) & + } else if (!is.null(sd_posterior) & is.null(se_new_trial) & is.null(prior_list) & is.null(dose_weights)) { - + w <- NULL S <- diag((sd_posterior)^2) - + # Bayesian & no re-estimation - } else if (!is.null(dose_weights) & !is.null(prior_list) & + } else if (!is.null(dose_weights) & !is.null(prior_list) & is.null(se_new_trial) & is.null(sd_posterior)) { - + w <- dose_weights + suppressMessages(round(unlist(lapply(prior_list, RBesT::ess)))) S <- NULL - + } else { - + stop (paste("Provided combiations of 'se_new_trial',", "'dose_weights', 'prior_list', 'sd_posterior' not allowed.", "See ?getContr for allowed combinations.")) - + } - + if (is.null(w)) { - + contr <- DoseFinding::optContr( models = mods, doses = dose_levels, S = S) - + } else { - + contr <- DoseFinding::optContr( models = mods, doses = dose_levels, w = w) - + } - + return (contr) - + } #' @title getCritProb -#' +#' #' @description This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that -#' when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). +#' when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). #' Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished #' 1) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the #' regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package. #' 2) Frequentist approach + re-estimation: If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the #' regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal() function of the DoseFinding package. -#' +#' #' @param mods An object of class "Mods" as specified in the DoseFinding package. #' @param dose_levels Vector containing the different dosage levels. #' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL @@ -286,44 +286,44 @@ getContr <- function ( #' mods = mods, #' dose_weights = c(50,50,50,50,50), #reflecting the planned sample size #' dose_levels = dose_levels, -#' alpha_crit_val = 0.05) +#' alpha_crit_val = 0.05) #' @return Multiplicity adjusted critical value on the probability scale. -#' +#' #' @export getCritProb <- function ( - + mods, dose_levels, dose_weights = NULL, se_new_trial = NULL, alpha_crit_val = 0.025 - + ) { - + checkmate::check_class(mods, classes = "Mods") checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(dose_weights)) checkmate::check_double(dose_weights, any.missing = FALSE, len = length(dose_levels)) checkmate::check_double(alpha_crit_val, lower = 0, upper = 1) - + contr <- getContr(mods = mods, dose_levels = dose_levels , dose_weights = dose_weights, se_new_trial = se_new_trial) - + crit_prob <- stats::pnorm(DoseFinding::critVal( corMat = contr$corMat, alpha = alpha_crit_val, df = 0, alternative = "one.sided")) - + return (crit_prob) - + } #' @title performBayesianMCPMod -#' +#' #' @description Performs Bayesian MCP Test step and modeling in a combined fashion. See performBayesianMCP() function for MCP Test step and getModelFits() for the modeling step -#' +#' #' @param posterior_list An object of class 'postList' as created by getPosterior() containing information about the (mixture) posterior distribution per dose group #' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step. #' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale). @@ -347,124 +347,124 @@ getCritProb <- function ( #' alpha_crit_val = 0.05) #' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), #' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), -#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , +#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , #' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , #' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) #' mu <- c(0, 1, 1.5, 2, 2.5) -#' se <- c(5, 4, 6, 7, 8) +#' S_hat <- c(5, 4, 6, 7, 8) #' posterior_list <- getPosterior( #' prior_list = prior_list, #' mu_hat = mu, -#' se_hat = se) +#' S_hat = S_hat) #' performBayesianMCPMod(posterior_list = posterior_list, #' contr = contr_mat, #' crit_prob_adj = critVal, #' simple = FALSE) -#' +#' #' @return Bayesian MCP test result as well as modeling result. -#' +#' #' @export performBayesianMCPMod <- function ( - + posterior_list, contr, crit_prob_adj, simple = FALSE - + ) { - + checkmate::check_class(posterior_list, "postList") checkmate::check_class(contr, "optContr") checkmate::check_class(crit_prob_adj, "numeric") checkmate::check_logical(simple) - + if (inherits(posterior_list, "postList")) { - + posterior_list <- list(posterior_list) - + } - + if (inherits(contr, "optContr")) { - + model_shapes <- colnames(contr$contMat) dose_levels <- as.numeric(rownames(contr$contMat)) - + } else if (length(contr) == length(posterior_list)) { - + model_shapes <- colnames(contr[[1]]$contMat) dose_levels <- as.numeric(rownames(contr[[1]]$contMat)) - + } else { - + stop ("Argument 'contr' must be of type 'optContr'") - + } b_mcp <- performBayesianMCP( posterior_list = posterior_list, contr = contr, crit_prob_adj = crit_prob_adj) - + fits_list <- lapply(seq_along(posterior_list), function (i) { - + if (b_mcp[i, 1]) { sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj") - + model_fits <- getModelFits( models = model_shapes, dose_levels = dose_levels, posterior = posterior_list[[i]], simple = simple) - + model_fits <- addSignificance(model_fits, sign_models) - + } else { - + NULL - + } - + }) - + bmcpmod <- list(BayesianMCP = b_mcp, Mod = fits_list) class(bmcpmod) <- "BayesianMCPMod" - + return (bmcpmod) - + } addSignificance <- function ( - + model_fits, sign_models - + ) { - + names(sign_models) <- NULL - + model_fits_out <- lapply(seq_along(model_fits), function (i) { - + c(model_fits[[i]], significant = sign_models[i]) - + }) - + attributes(model_fits_out) <- attributes(model_fits) - + return (model_fits_out) - + } #' @title performBayesianMCP -#' +#' #' @description Performs Bayesian MCP Test step, as described in Fleischer et al. (2022). #' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups). #' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale). #' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. 2022. Bayesian MCPMod. Pharmaceutical Statistics. 21(3): 654-670. doi:10.1002/pst.2193 -#' @param posterior_list An object derived with getPosterior with information about the (mixture) posterior distribution per dose group +#' @param posterior_list An object derived with getPosterior with information about the (mixture) posterior distribution per dose group #' @param contr An object of class 'optContr' as created by the getContr() function. It contains the contrast matrix to be used for the testing step. #' @param crit_prob_adj A getCritProb object, specifying the critical value to be used for the testing (on the probability scale) -#' +#' #' @examples #' mods <- DoseFinding::Mods(linear = NULL, #' linlog = NULL, @@ -484,52 +484,52 @@ addSignificance <- function ( #' alpha_crit_val = 0.05) #' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), #' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), -#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , +#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , #' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , #' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) #' mu <- c(0, 1, 1.5, 2, 2.5) -#' se <- c(5, 4, 6, 7, 8) +#' S_hat <- c(5, 4, 6, 7, 8) #' posterior_list <- getPosterior( #' prior_list = prior_list, #' mu_hat = mu, -#' se_hat = se) -#' +#' S_hat = S_hat) +#' #' performBayesianMCP(posterior_list = posterior_list, #' contr = contr_mat, #' crit_prob_adj = critVal) -#' -#' @return Bayesian MCP test result, with information about p-values for the individual dose-response shapes and overall significance -#' +#' +#' @return Bayesian MCP test result, with information about p-values for the individual dose-response shapes and overall significance +#' #' @export performBayesianMCP <- function( - + posterior_list, contr, crit_prob_adj - + ) { - + checkmate::check_class(posterior_list, "postList") checkmate::check_class(contr, "optContr") checkmate::check_class(crit_prob_adj, "numeric") checkmate::check_numeric(crit_prob_adj, lower = 0, upper = Inf) - + if (inherits(posterior_list, "postList")) { - + posterior_list <- list(posterior_list) - + } - + if (inherits(contr, "optContr")) { - + b_mcp <- t(sapply(posterior_list, BayesMCPi, contr, crit_prob_adj)) - + } else { - + b_mcp <- t(mapply(BayesMCPi, posterior_list, contr, crit_prob_adj)) - + } - + class(b_mcp) <- "BayesianMCP" attr(b_mcp, "critProbAdj") <- crit_prob_adj attr(b_mcp, "successRate") <- mean(b_mcp[, 1]) @@ -537,70 +537,70 @@ performBayesianMCP <- function( test = is.na(attr(posterior_list[[1]], "ess")), yes = numeric(0), no = rowMeans(sapply(posterior_list, function (posteriors) { - + attr(posteriors, "ess") - + }))) - + return (b_mcp) - + } getModelSuccesses <- function (b_mcp) { - + stopifnot(inherits(b_mcp, "BayesianMCP")) - + model_indices <- grepl("post_probs.", colnames(b_mcp)) model_names <- colnames(b_mcp)[model_indices] |> sub(pattern = "post_probs.", replacement = "", x = _) - + model_successes <- colMeans(b_mcp[, model_indices] > b_mcp[, "crit_prob_adj"]) - + names(model_successes) <- model_names - + return (model_successes) - + } BayesMCPi <- function ( - + posterior_i, contr, crit_prob_adj - + ) { - + getPostProb <- function ( - + contr_j, # j: dose level post_combs_i # i: simulation outcome - + ) { - + ## Test statistic = sum over all components of ## posterior weight * normal probability distribution of ## critical values for doses * estimated mean / sqrt(product of critical values for doses) - + ## Calculation for each component of the posterior contr_theta <- apply(post_combs_i$means, 1, `%*%`, contr_j) contr_var <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2) contr_weights <- post_combs_i$weights - + ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j) post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var))) - + return (post_probs) - + } - + post_combs_i <- getPostCombsI(posterior_i) post_probs <- apply(contr$contMat, 2, getPostProb, post_combs_i) - + res <- c(sign = ifelse(max(post_probs) > crit_prob_adj, 1, 0), crit_prob_adj = crit_prob_adj, max_post_prob = max(post_probs), post_probs = post_probs) - + return (res) - -} \ No newline at end of file + +} diff --git a/R/posterior.R b/R/posterior.R index 5aa6d2d..1e905c4 100644 --- a/R/posterior.R +++ b/R/posterior.R @@ -1,10 +1,10 @@ #' @title getPosterior -#' +#' #' @description Either the patient level data or both mu_hat as well as sd_hat must to be provided. #' If patient level data is provided mu_hat and se_hat are calculated within the function using a linear model. #' This function calculates the posterior for every dose group independently via the RBesT function postmix(). -#' -#' @param prior_list a prior list with information about the prior to be used for every dose group +#' +#' @param prior_list a prior list with information about the prior to be used for every dose group #' @param data dataframe containing the information of dose and response. Default NULL #' Also a simulateData object can be provided. #' @param mu_hat vector of estimated mean values (per dose group). @@ -14,280 +14,280 @@ #' @examples #' prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), #' DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), -#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , +#' DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , #' DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , #' DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) #' mu <- c(0, 1, 1.5, 2, 2.5) -#' se <- c(5, 4, 6, 7, 8) -#' +#' S_hat <- c(5, 4, 6, 7, 8) +#' #' posterior_list <- getPosterior( #' prior_list = prior_list, #' mu_hat = mu, -#' se_hat = se) -#' +#' S_hat = S_hat) +#' #' summary(posterior_list) -#' +#' #' @export getPosterior <- function( - + prior_list, data = NULL, mu_hat = NULL, S_hat = NULL, calc_ess = FALSE - + ) { - + checkmate::check_data_frame(data, null.ok = TRUE) checkmate::check_list(prior_list, names = "named", any.missing = FALSE) checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf) checkmate::check_vector(S_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(S_hat, null.ok = TRUE, lower = 0, upper = Inf) - + is_matrix_S_hat <- FALSE - + stopifnot("prior_list must be an object of RBesT package" = all(sapply(prior_list, function(x) is(x, "normMix") | is(x, "betaMix") | is(x, "mix")))) - + if (!is.null(mu_hat) && !is.null(S_hat) && is.null(data)) { - + if (is.matrix(S_hat)) { - + is_matrix_S_hat <- TRUE - + } else if (is.vector(S_hat)) { - + is_matrix_S_hat <- FALSE - + se_hat <- S_hat rm(S_hat) - + } else { - + stop("S_hat has to be either a vector or matrix") - + } - + if (is_matrix_S_hat) { - + stopifnot("dim of S_hat must equal length of prior list" = dim(S_hat)[2] == length(prior_list)) - + prior_mix <- priorList2priorMix(prior_list) - + posterior <- DoseFinding::mvpostmix( priormix = prior_mix, mu_hat = mu_hat, S_hat = S_hat) - + posterior_list <- postMix2posteriorList( - posterior_list = posterior, - prior_list = prior_list, + posterior_list = posterior, + prior_list = prior_list, calc_ess = calc_ess) - + } else if (!is_matrix_S_hat) { - + posterior_list <- getPosteriorI( prior_list = prior_list, mu_hat = mu_hat, se_hat = se_hat, calc_ess = calc_ess) - + } - + } else if (is.null(mu_hat) && is.null(S_hat) && !is.null(data)) { - + posterior_list <- lapply(split(data, data$simulation), getPosteriorI, prior_list = prior_list, calc_ess = calc_ess) - + } else { - + stop ("Either 'data' or 'mu_hat' and 'se_hat' must not be NULL.") - + } - + if (length(posterior_list) == 1) { - + posterior_list <- posterior_list[[1]] - + } - + return (posterior_list) - + } getPosteriorI <- function( - + data_i = NULL, prior_list, mu_hat = NULL, se_hat = NULL, calc_ess = FALSE - + ) { - + checkmate::check_data_frame(data_i, null.ok = TRUE) checkmate::check_list(prior_list, names = "named", any.missing = FALSE) checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf) checkmate::check_vector(se_hat, any.missing = FALSE, null.ok = TRUE) checkmate::check_double(se_hat, null.ok = TRUE, lower = 0, upper = Inf) - + if (is.null(mu_hat) && is.null(se_hat)) { - + checkmate::check_data_frame(data_i, null.ok = FALSE) checkmate::assert_names(names(data_i), must.include = "response") checkmate::assert_names(names(data_i), must.include = "dose") - + anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1) mu_hat <- summary(anova_res)$coefficients[, 1] se_hat <- summary(anova_res)$coefficients[, 2] - + } else if (!is.null(mu_hat) && !is.null(se_hat)) { - - stopifnot("m_hat length must match number of dose levels" = + + stopifnot("m_hat length must match number of dose levels" = length(prior_list) == length(mu_hat)) # , - # "se_hat length must match number of dose levels" = + # "se_hat length must match number of dose levels" = # length(prior_list) == length(se_hat)) - + } else { - + stop ("Both mu_hat and se_hat must be provided.") - + } - + post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat, SIMPLIFY = FALSE) - + if (is.null(names(prior_list))) { - + names(prior_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1]))) - + } - + names(post_list) <- names(prior_list) class(post_list) <- "postList" - + comp_indx <- createMapping(prior_list) comp_mat_ind <- do.call("expand.grid", comp_indx) - + attr(post_list, "ess") <- calcEss(calc_ess, post_list) - + diagonals <- lapply(seq_along(comp_mat_ind[, 1]), function(x) { - + lapply(seq_along(comp_mat_ind[x, ]), function(y) return(post_list[[y]]["s", comp_mat_ind[x,y]])) - + }) - + attr(post_list, "full covariance matrices") <- lapply(seq_along(diagonals), function(x) diag(diagonals[[x]])) - + return (post_list) - + } #' @title getESS -#' +#' #' @description This function calculates the effective sample size for every dose group via the RBesT function ess(). -#' +#' #' @param post_list A posterior list object, for which the effective sample size for each dose group should be calculated #' #' @return A vector of the effective sample sizes for each dose group -#' +#' #' @export getESS <- function ( - + post_list - + ) { - + # make s3 method for postList object suppressMessages(round(sapply(post_list, RBesT::ess), 1)) - + } getPostCombsI <- function ( - + posterior_i - + ) { - + post_params <- list( weights = lapply(posterior_i, function (x) x[1, ]), means = lapply(posterior_i, function (x) x[2, ]), vars = lapply(posterior_i, function (x) x[3, ]^2)) - + post_combs <- lapply(post_params, expand.grid) post_combs$weights <- apply(post_combs$weights, 1, prod) - + return (post_combs) - + } priorList2priorMix <- function (prior_list) { - + checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) - + # create mapping args <- createMapping(prior_list) - comp_ind <- do.call("expand.grid", args) + comp_ind <- do.call("expand.grid", args) n_comps_prior <- nrow(comp_ind) - + # map information -> mapping function? prior_weight <- matrix( - sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, + sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, function (y) prior_list[[x]][1, comp_ind[y, x]])), nrow = n_comps_prior) - + prior_mean <- matrix(sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, function (y) prior_list[[x]][2, comp_ind[y, x]])), nrow = n_comps_prior) prior_sd <- matrix(sapply(1:length(prior_list), function (x) sapply(1:n_comps_prior, function (y) prior_list[[x]][3, comp_ind[y, x]])), nrow = n_comps_prior) - + prior_weight <- apply(prior_weight, 1, prod) - + # create prior_mix object prior_weight <- as.list(prior_weight) prior_mean <- asplit(prior_mean, 1) prior_vc <- lapply(asplit(prior_sd^2, 1), diag) prior_mix <- list(prior_weight, prior_mean, prior_vc) - + return (prior_mix) - + } - + postMix2posteriorList <- function ( - + posterior_list, prior_list, calc_ess = FALSE - + ) { - + checkmate::assert_list(posterior_list, names = "named", any.missing = FALSE, null.ok = FALSE) checkmate::assert_list(prior_list, names = "named", any.missing = FALSE, null.ok = FALSE) - + getIndx <- function (x, y) which(comp_mat_ind[, x] == comp_indx[[x]][y]) - + # create mapping comp_indx <- createMapping(prior_list) comp_mat_ind <- do.call("expand.grid", comp_indx) - + # map posterior information posterior_weight <- lapply(seq_along(prior_list), function (x) sapply(seq_along(comp_indx[[x]]), function (y) - sum(unlist(posterior_list$weights[getIndx(x, y)])))) - + sum(unlist(posterior_list$weights[getIndx(x, y)])))) + posterior_mean <- lapply(seq_along(prior_list), function (x) sapply(seq_along(comp_indx[[x]]), function (y) mean(do.call(cbind, posterior_list$mean[getIndx(x, y)])[x, ]))) - + posterior_sd <- lapply(seq_along(prior_list), function (x) sapply(seq_along(comp_indx[[x]]), function (y) mean(do.call(rbind, lapply(posterior_list$covmat, diag)[getIndx(x, y)])[, x]))) - + combined_vectors <- mapply(function (x, y, z) Map(c, x, y, z), posterior_weight, posterior_mean, lapply(posterior_sd, sqrt), SIMPLIFY = FALSE) - + # create posterior list posterior_list_RBesT <- lapply(seq_along(combined_vectors), function (x) do.call(RBesT::mixnorm, @@ -297,52 +297,52 @@ postMix2posteriorList <- function ( names(posterior_list_RBesT) <- names(prior_list) comp_names <- lapply(prior_list, colnames) for (i in seq_along(posterior_list_RBesT)) { - + colnames(posterior_list_RBesT[[i]]) <- comp_names[[i]] - + } rm(i) - + ## set attributes class(posterior_list_RBesT) <- "postList" attr(posterior_list_RBesT, "ess") <- calcEss(calc_ess, posterior_list_RBesT) - + diagonals <- lapply(seq_along(comp_mat_ind[, 1]), function(x) { - + lapply(seq_along(comp_mat_ind[x, ]), function(y) return(posterior_list_RBesT[[y]]["s", comp_mat_ind[x,y]])) - + }) - + attr(posterior_list_RBesT, "covariance matrices") <- lapply(seq_along(diagonals), function(x) diag(diagonals[[x]])) - + return (posterior_list_RBesT) - + } calcEss <- function(calc_ess, posterior_output) { - + checkmate::assert_logical(calc_ess, null.ok = FALSE, len = 1) checkmate::assert_list(posterior_output, names = "named", any.missing = FALSE, null.ok = FALSE) - + if (calc_ess) { - + post_ESS <- getESS(posterior_output) - + } else { - + post_ESS <- numeric(0) - + } - + return(post_ESS) - + } createMapping <- function (prior_list) { - + n_comps <- unlist(lapply(prior_list, ncol)) comp_indx <- lapply(seq_along(prior_list), function (x) seq_len(n_comps[x])) - + return (comp_indx) - + } diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd index e018b7e..7b67845 100644 --- a/man/assessDesign.Rd +++ b/man/assessDesign.Rd @@ -57,7 +57,7 @@ mods <- DoseFinding::Mods(linear = NULL, sd <- 12 prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2), DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), - DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , + DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) n_patients <- c(40, 60, 60, 60, 60) diff --git a/man/getContr.Rd b/man/getContr.Rd index 9232ac3..e18ec0e 100644 --- a/man/getContr.Rd +++ b/man/getContr.Rd @@ -59,6 +59,6 @@ sd_posterior <- c(2.8, 3, 2.5, 3.5, 4) contr_mat <- getContr( mods = mods, dose_levels = dose_levels, - sd_posterior = sd_posterior) + sd_posterior = sd_posterior) } diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd index 187debd..68978b5 100644 --- a/man/getCritProb.Rd +++ b/man/getCritProb.Rd @@ -48,5 +48,5 @@ critVal <- getCritProb( mods = mods, dose_weights = c(50,50,50,50,50), #reflecting the planned sample size dose_levels = dose_levels, - alpha_crit_val = 0.05) + alpha_crit_val = 0.05) } diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd index 8b82c0d..b904a28 100644 --- a/man/getPosterior.Rd +++ b/man/getPosterior.Rd @@ -20,9 +20,9 @@ Also a simulateData object can be provided.} \item{mu_hat}{vector of estimated mean values (per dose group).} -\item{calc_ess}{boolean variable, indicating whether effective sample size should be calculated. Default FALSE} +\item{S_hat}{vector or matrix of estimated standard deviations (per dose group).} -\item{se_hat}{vector of estimated standard deviations (per dose group).} +\item{calc_ess}{boolean variable, indicating whether effective sample size should be calculated. Default FALSE} } \value{ posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group @@ -35,17 +35,17 @@ This function calculates the posterior for every dose group independently via th \examples{ prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), - DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , + DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) mu <- c(0, 1, 1.5, 2, 2.5) -se <- c(5, 4, 6, 7, 8) +S_hat <- c(5, 4, 6, 7, 8) posterior_list <- getPosterior( prior_list = prior_list, mu_hat = mu, - se_hat = se) - + S_hat = S_hat) + summary(posterior_list) } diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd index 21105fd..8957fcb 100644 --- a/man/performBayesianMCP.Rd +++ b/man/performBayesianMCP.Rd @@ -40,16 +40,16 @@ critVal <- getCritProb( alpha_crit_val = 0.05) prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), - DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , + DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) mu <- c(0, 1, 1.5, 2, 2.5) -se <- c(5, 4, 6, 7, 8) +S_hat <- c(5, 4, 6, 7, 8) posterior_list <- getPosterior( prior_list = prior_list, mu_hat = mu, - se_hat = se) - + S_hat = S_hat) + performBayesianMCP(posterior_list = posterior_list, contr = contr_mat, crit_prob_adj = critVal) diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd index d2fb698..ef8709c 100644 --- a/man/performBayesianMCPMod.Rd +++ b/man/performBayesianMCPMod.Rd @@ -40,15 +40,15 @@ critVal <- getCritProb( alpha_crit_val = 0.05) prior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2), DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2), - DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , + DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) , DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) , DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13), sigma = 2)) mu <- c(0, 1, 1.5, 2, 2.5) -se <- c(5, 4, 6, 7, 8) +S_hat <- c(5, 4, 6, 7, 8) posterior_list <- getPosterior( prior_list = prior_list, mu_hat = mu, - se_hat = se) + S_hat = S_hat) performBayesianMCPMod(posterior_list = posterior_list, contr = contr_mat, crit_prob_adj = critVal,