diff --git a/.Rbuildignore b/.Rbuildignore index b38ebc96..cd9c03b3 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,7 +8,6 @@ ^_pkgdown\.yml$ ^docs$ ^pkgdown$ -^vignettes$ ^\.httr-oauth$ ^CONDUCT\.md$ ^man-roxygen$ diff --git a/.github/workflows/check-standard.yaml b/.github/workflows/check-standard.yaml index 5bbd54af..2ffa55cb 100644 --- a/.github/workflows/check-standard.yaml +++ b/.github/workflows/check-standard.yaml @@ -34,7 +34,7 @@ jobs: - uses: r-lib/actions/setup-r@v2 with: - extra-repositories: 'https://inla.r-inla-download.org/R/stable' + extra-repositories: 'https://inla.r-inla-download.org/R/testing' r-version: ${{ matrix.config.r }} http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true @@ -43,7 +43,7 @@ jobs: - name: Add inla repo run: | - cat('options(repos = c(INLA = "https://inla.r-inla-download.org/R/stable/", getOption("repos")))', file = "~/.Rprofile", append = TRUE) + cat('options(repos = c(INLA = "https://inla.r-inla-download.org/R/testing/", getOption("repos")))', file = "~/.Rprofile", append = TRUE) # cat("\noptions(tinytex.verbose = TRUE)\n", file = "~/.Rprofile", append = TRUE) shell: Rscript {0} @@ -53,7 +53,7 @@ jobs: shell: Rscript {0} # - name: Install JAGS -# run: | +# run: | # sudo apt-get update -y # sudo apt-get install -y jags @@ -67,15 +67,13 @@ jobs: - name: Install jags (macOS-latest) if: runner.os == 'macOS' - run : | - rm '/usr/local/bin/gfortran' - brew install jags + run : brew install jags - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: - any::rcmdcheck, local::. - needs: check + dependencies: '"all"' + extra-packages: | + rcmdcheck - uses: r-lib/actions/check-r-package@v2 with: diff --git a/BCEA.Rproj b/BCEA.Rproj index 9f964990..6843888d 100644 --- a/BCEA.Rproj +++ b/BCEA.Rproj @@ -15,4 +15,4 @@ LaTeX: pdfLaTeX BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source -PackageRoxygenize: rd,collate,namespace +PackageRoxygenize: rd,collate,namespace,vignette diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index 5eddd3d7..adf2dc35 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 2.4.2 -Date: 2022-08-23 11:44:30 UTC -SHA: f0469ce213fba001efa2707da7fa79289b337105 +Version: 2.4.4 +Date: 2023-06-05 12:52:55 UTC +SHA: 770d3d2d910eca75ebaaccf4c290c20e4668a95a diff --git a/DESCRIPTION b/DESCRIPTION index 95793658..cb5d13b0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,11 @@ Package: BCEA Type: Package Title: Bayesian Cost Effectiveness Analysis -Version: 2.4.2.9000 +Version: 2.4.5 Authors@R: c( person("Gianluca", "Baio", - email = "gianluca@stats.ucl.ac.uk", - role = c("aut", "cre"), + email = "g.baio@ucl.ac.uk", + role = c("aut", "cre", "cph"), comment = c(ORCID = "0000-0003-4314-2570")), person("Andrea", "Berardi", email = "a.berardi@ucl.ac.uk", @@ -23,36 +23,38 @@ Imports: cli (>= 3.3.0), dplyr, ggplot2, - gridExtra, graphics, + gridExtra, MASS, Matrix, + MCMCvis, purrr, Rdpack, reshape2, rlang, - scales + scales, + voi Depends: R (>= 3.5.0) Suggests: coda, grid, - INLA, knitr, markdown, - MCMCvis, mgcv, plotly, + plotrix, + RColorBrewer, rjags, rmarkdown, rstan, - RColorBrewer, splancs, testthat (>= 2.1.0), vdiffr, withr RdMacros: Rdpack -Additional_repositories: https://inla.r-inla-download.org/R/stable/ +VignetteBuilder: knitr +Additional_repositories: https://inla.r-inla-download.org/R/testing/ Description: Produces an economic evaluation of a sample of suitable variables of cost and effectiveness / utility for two or more interventions, e.g. from a Bayesian model in the form of MCMC simulations. @@ -64,6 +66,9 @@ URL: https://gianluca.statistica.it/software/bcea/, https://gianluca.statistica.it/, https://github.com/giabaio/BCEA/ NeedsCompilation: no -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3.9000 Encoding: UTF-8 BugReports: https://github.com/n8thangreen/BCEA/issues/ +Remotes: + chjackson/voi +LazyData: false diff --git a/NAMESPACE b/NAMESPACE index 21791667..bbdfef24 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -94,9 +94,9 @@ export(validate_bcea) import(dplyr) import(ggplot2) import(grid) -import(purrr) import(reshape2) importFrom(MASS,kde2d) +importFrom(MCMCvis,MCMCchains) importFrom(Matrix,Matrix) importFrom(Matrix,expm) importFrom(Rdpack,reprompt) @@ -109,7 +109,6 @@ importFrom(dplyr,slice) importFrom(grDevices,colors) importFrom(grDevices,colours) importFrom(grDevices,dev.new) -importFrom(grDevices,dev.off) importFrom(grDevices,devAskNewPage) importFrom(grDevices,grey.colors) importFrom(grDevices,pdf.options) @@ -131,7 +130,10 @@ importFrom(graphics,text) importFrom(grid,unit) importFrom(gridExtra,grid.arrange) importFrom(purrr,keep) +importFrom(purrr,list_cbind) +importFrom(purrr,map) importFrom(purrr,map_int) +importFrom(purrr,map_lgl) importFrom(purrr,pluck) importFrom(reshape2,melt) importFrom(rlang,.data) @@ -140,10 +142,7 @@ importFrom(scales,label_dollar) importFrom(stats,as.formula) importFrom(stats,cov) importFrom(stats,density) -importFrom(stats,dist) -importFrom(stats,dnorm) importFrom(stats,lm) -importFrom(stats,optim) importFrom(stats,qnorm) importFrom(stats,qqline) importFrom(stats,qqnorm) @@ -153,9 +152,8 @@ importFrom(stats,rnorm) importFrom(stats,runif) importFrom(stats,sd) importFrom(stats,setNames) -importFrom(stats,update) importFrom(utils,methods) importFrom(utils,modifyList) -importFrom(utils,select.list) importFrom(utils,str) importFrom(utils,tail) +importFrom(voi,evppi) diff --git a/NEWS.md b/NEWS.md index 907810a0..8ecf9da3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,47 @@ -# BCEA 2.4.2.9000 +# BCEA 2.4.5 + +_October 2023_ + +Moved internal EVPPI calculation out of `BCEA` and now uses `voi` package instead. +Refactoring but retaining same interface and functionality. + +* `evppi()` tested against all use cases in BCEA book (1c1457d2) +* Select parameters by position (as well as name) in new `evppi()` (f2e4d005) +* Use single parameter case only like `voi` package for methods `sal` and `so` (#140) +* New `evppi()` matching output of old `evppi()` (1e2c5e7) +* Latest development version of `voi` needed when we use `check = TRUE` in `voi::evppi()` in order to access fitting data (6e436b5, 94f5fc5) +* No longer require `INLA` package to be available inside of `BCEA` so can remove direct dependency. This helps with passing CRAN checks and GitHub Actions () + +# BCEA 2.4.4 +_June 2023_ + +* Patch to fix a CRAN checks error. Suggested package `MCMCvis` wasn't used conditionally in unit test. Moved to Required packages in `DESCRIPTION`. + +# BCEA 2.4.3 +_May 2023_ + +## Bug fixes + +* Consistent colours across plots for each intervention for grid of plots in `plot.bcea()` (cf1ee43) +* `make.report()` change variable name (f940f2e) +* Fixed issue with summary table where names of interventions in the wrong order (6a006e3) +* `summary.bcea()` now only prints results for chosen comparisons and not always all of them. `kstar` and `best` in `bcea()` object were not updated with subset of interventions (#125) + +## Refactoring + +* `withr::with_par()` used in plotting function `plot.bcea()` to only temporarily change graphics parameters. (725c536) +* Using `@md` and markdown syntax in function documentation +* Update `psa.struct()` to add the absolute value in the formula to compute the weights (1cea278) +* Use `dplyr` piping new syntax from `.data$*` to simply using speech marks `"*"` (2b280ad) + +## Miscellaneous + +* Template added for GitHub Issues (0ea59fa) # BCEA 2.4.2 -August 2022 +_August 2022_ ## Bug fixes @@ -24,15 +62,14 @@ August 2022 ## New features * Can now specify what order the interventions labels are in the legend for ce plane (and contour plots) for base R and ggplot2 i.e. reference first or second with optional `ref_first` argument (cc38f07) -* Can specify currency for axes in `ce-plane.plot` and `ceac.plot` `ggplot2` versions (6808aa6) -* Argument added to `ceplane.plot` of `icer_annot` to annotate each of the ICER points with the text label of the intervention name. Only for `ggplot2` at the moment. (a7b4beb) +* Can specify currency for axes in `ceplane.plot()` and `ceac.plot()` `ggplot2` versions (6808aa6) +* Argument added to `ceplane.plot()` of `icer_annot` to annotate each of the ICER points with the text label of the intervention name. Only for `ggplot2` at the moment. (a7b4beb) * Added `pos` argument to `contour2()` so that its consistent with `contour()` and `ceplane.plot()`. (50f8f8b) * Allow passing `ref` argument by name as well as index in `bcea()`. (9eab459) - # BCEA 2.4.1.2 -April 2022 +_April 2022_ ## Bug fixes @@ -64,7 +101,7 @@ April 2022 # BCEA 2.4.1.1 -Oct 2021 +_Oct 2021_ ## Major refactoring @@ -156,9 +193,9 @@ Oct 2015 * New function for EVPPI using SPDE-INLA * Modifications to the EVPPI functions * Documentation updated -* Allows `xlim` & `ylim` in the `ceplane.plot`, `contour` and `contour2` functions +* Allows `xlim` & `ylim` in the `ceplane.plot()`, `contour()` and `contour2()` functions * It is now possible to run `bcea` for a scalar wtp -* Old `evppi` function and method has been renamed `evppi0`, which means there's also a new `plot.evppi0` method +* Old `evppi()` function and method has been renamed `evppi0`, which means there's also a new `plot.evppi0` method # BCEA 2.1-0 13 Jan 2015 diff --git a/R/BCEA-package.R b/R/BCEA-package.R index 6be356d3..4b0fd1ac 100644 --- a/R/BCEA-package.R +++ b/R/BCEA-package.R @@ -1,39 +1,23 @@ +#------------------------------------------------------------------------------# +# +# ___ ___ ___ +# _____ / /\ / /\ / /\ +# / /::\ / /:/ / /:/_ / /::\ +# / /:/\:\ / /:/ / /:/ /\ / /:/\:\ +# / /:/~/::\ / /:/ ___ / /:/ /:/_ / /:/~/::\ +# /__/:/ /:/\:| /__/:/ / /\ /__/:/ /:/ /\ /__/:/ /:/\:\ +# \ \:\/:/~/:/ \ \:\ / /:/ \ \:\/:/ /:/ \ \:\/:/__\/ +# \ \::/ /:/ \ \:\ /:/ \ \::/ /:/ \ \::/ +# \ \:\/:/ \ \:\/:/ \ \:\/:/ \ \:\ +# \ \::/ \ \::/ \ \::/ \ \:\ +# \__\/ \__\/ \__\/ \__\/ +# +# +#------------------------------------------------------------------------------# -#' BCEA: A package for Bayesian Cost-Effectiveness Analysis -#' -#' A package to post-process the results of a Bayesian health economic model -#' and produce standardised output for the analysis of the results. -#' -#' \tabular{ll}{ Package: \tab BCEA\cr -#' Type: \tab Package\cr -#' Version: \tab 2.4.2\cr -#' Date: \tab 2022-09-03\cr -#' License: \tab GPL2 \cr -#' LazyLoad: \tab Yes\cr } -#' -#' BCEA produces a health economic evaluation given a random sample of -#' suitable variables of costs and clinical benefits for two or more -#' interventions, e.g. using results of a Bayesian model (possibly based on -#' MCMC) in the form of simulations from the posterior distributions. -#' Compares one of the interventions (the "reference") to the others -#' ("comparators"). Produces many summaries and plots to analyse the results. -#' -#' @aliases BCEA-package BCEA -#' -#' @author Gianluca Baio, Andrea Berardi, Anna Heath, Nathan Green -#' -#' @references -#' \insertRef{Baio2011}{BCEA} -#' -#' \insertRef{Baio2013}{BCEA} -#' -#' \insertRef{Baio2017}{BCEA} -#' -#' @keywords package -#' -#' @docType package -#' @name BCEA-package -#' -#' @import dplyr ggplot2 purrr reshape2 -#' @importFrom Rdpack reprompt +#' @keywords internal +"_PACKAGE" + +## usethis namespace: start +## usethis namespace: end NULL diff --git a/R/CEriskav_plot_graph.R b/R/CEriskav_plot_graph.R index b213c3ef..389f9e67 100644 --- a/R/CEriskav_plot_graph.R +++ b/R/CEriskav_plot_graph.R @@ -88,8 +88,8 @@ CEriskav_plot_ggplot <- function(he, pos_legend) { eib_dat <- melt(he$eibr[, default_comp, , drop = FALSE], value.name = "eibr") %>% - rename(k = .data$Var1, - r = .data$Var3) %>% + rename(k = "Var1", + r = "Var3") %>% mutate(r = as.factor(.data$r)) eibr_plot <- @@ -122,8 +122,8 @@ CEriskav_plot_ggplot <- function(he, pos_legend) { evi_dat <- melt(he$evir, value.name = "evir") %>% - rename(r = .data$Var2, - k = .data$Var1) %>% + rename(r = "Var2", + k = "Var1") %>% mutate(r = as.factor(.data$r)) evir_plot <- diff --git a/R/CreateInputs.R b/R/CreateInputs.R index c5230942..9f42022d 100644 --- a/R/CreateInputs.R +++ b/R/CreateInputs.R @@ -14,6 +14,12 @@ createInputs.default <- function(inputs, print_is_linear_comb = TRUE) { + # remove NA columns + if (sum(is.na(inputs)) > 0) { + inputs <- inputs[ , colSums(is.na(inputs)) == 0] + message("Dropped any columns containing NAs") + } + if (!is.logical(print_is_linear_comb)) stop("print_is_linear_comb must be logical.", call. = FALSE) diff --git a/R/bcea.default.R b/R/bcea.default.R index 2d6ddd9e..21cd956a 100644 --- a/R/bcea.default.R +++ b/R/bcea.default.R @@ -26,7 +26,7 @@ bcea.default <- function(eff, if (is.null(ref)) { ref <- 1 - message("No reference selected. Defaulting to first intervention.") + message("No reference selected. Defaulting to first intervention.") } if (!is.null(k) && length(k) == 1) @@ -69,7 +69,7 @@ bcea.default <- function(eff, df_ce <- df_ce %>% select(-ref) %>% - rename(ref = .data$ints) %>% + rename(ref = "ints") %>% merge(df_ce, by = c("ref", "sim"), suffixes = c("0", "1"), @@ -94,6 +94,7 @@ bcea.default <- function(eff, #' @rdname bcea #' @param ... Additional arguments +#' @importFrom MCMCvis MCMCchains #' @export bcea.rjags <- function(eff, ...) { diff --git a/R/best_interv_given_k.R b/R/best_interv_given_k.R index afb8ef2c..974378e8 100644 --- a/R/best_interv_given_k.R +++ b/R/best_interv_given_k.R @@ -13,7 +13,7 @@ best_interv_given_k <- function(eib, ref, comp) { - + if (length(comp) == 1) { best <- rep(ref, NROW(eib)) diff --git a/R/ce_table.R b/R/ce_table.R index 0cf516da..2400c133 100644 --- a/R/ce_table.R +++ b/R/ce_table.R @@ -54,8 +54,7 @@ tabulate_means <- function(he, comp_label = NULL, ...) { - if (is.null(comp_label)) - comp_label <- 1:he$n_comparisons + comp_label <- comp_label %||% seq_len(he$n_comparisons) data.frame( lambda.e = vapply(1:he$n_comparisons, diff --git a/R/ceac.plot.R b/R/ceac.plot.R index 28154f27..57074c7f 100644 --- a/R/ceac.plot.R +++ b/R/ceac.plot.R @@ -64,7 +64,6 @@ ceac.plot.bcea <- function(he, pos = c(1, 0), graph = c("base", "ggplot2", "plotly"), ...) { - graph <- match.arg(graph) he <- setComparisons(he, comparison) diff --git a/R/ceac_plot_graph.R b/R/ceac_plot_graph.R index bbdc9510..b89d0b5f 100644 --- a/R/ceac_plot_graph.R +++ b/R/ceac_plot_graph.R @@ -122,21 +122,21 @@ ceac_ggplot <- function(he, ceac_dat <- he[[ceac]] n_lines <- ncol(ceac_dat) + len_k <- length(he$k) data_psa <- tibble(k = rep(he$k, times = n_lines), ceac = c(ceac_dat), - comparison = as.factor(rep(1:n_lines, - each = length(he$k)))) + comparison = as.factor(rep(1:n_lines, each = len_k))) graph_params <- helper_ggplot_params(he, graph_params) legend_params <- make_legend_ggplot(he, pos_legend) theme_add <- purrr::keep(extra_params, is.theme) - ggplot(data_psa, aes(.data$k, .data$ceac)) + + ggplot(data_psa, aes(x = .data$k, y = .data$ceac)) + geom_line(aes(linetype = .data$comparison, - size = .data$comparison, + size = factor(.data$comparison), colour = factor(.data$comparison))) + theme_ceac() + theme_add + # theme @@ -176,8 +176,7 @@ ceac_plot_plotly <- function(he, sapply(comparisons_label, function(x) rep(x, length(he$k))) ))) - if (is.null(graph_params$line$types)) - graph_params$line$type <- rep_len(1:6, he$n_comparisons) + graph_params$line$type <- graph_params$line$type %||% rep_len(1:6, he$n_comparisons) # opacities if (!is.null(graph_params$area$color)) diff --git a/R/ceaf.plot.R b/R/ceaf.plot.R index bf05ce95..1f3efb58 100644 --- a/R/ceaf.plot.R +++ b/R/ceaf.plot.R @@ -101,8 +101,8 @@ ceaf.plot.pairwise <- function(mce, } else { df <- data.frame(k = mce$k, ceaf = mce$ceaf) - ggceaf <- - ggplot(df, aes(x = .data$k, y = .data$ceaf)) + + + ggplot(df, aes(x = .data$k, y = .data$ceaf)) + theme_bw() + geom_line() + coord_cartesian(ylim = c(-0.05, 1.05)) + @@ -119,8 +119,6 @@ ceaf.plot.pairwise <- function(mce, face = "bold", size = 14.3, hjust = 0.5)) - - return(ggceaf) } } diff --git a/R/compute_xxx.R b/R/compute_xxx.R index 246991dc..e8bf2fb2 100644 --- a/R/compute_xxx.R +++ b/R/compute_xxx.R @@ -264,7 +264,7 @@ compute_IB <- function(df_ce, k) { df_ce <- df_ce %>% filter(ints != .data$ref) %>% - rename(comps = .data$ints) + rename(comps = "ints") ib_df <- data.frame(k = rep(k, each = nrow(df_ce)), @@ -304,7 +304,7 @@ compute_ICER <- function(df_ce) { summarise( ICER = mean(.data$delta_c)/mean(.data$delta_e)) %>% ungroup() %>% - select(.data$ICER) %>% # required to match current format + select("ICER") %>% # required to match current format unlist() %>% setNames(comp_names) } @@ -331,7 +331,7 @@ comp_names_from_ <- function(df_ce) { filter(.data$ref != .data$ints) %>% distinct() %>% arrange(.data$ints) %>% - select(.data$interv_names) %>% + select("interv_names") %>% unlist() } diff --git a/R/diag.evppi.R b/R/diag.evppi.R index 53d268c2..0471a9fe 100644 --- a/R/diag.evppi.R +++ b/R/diag.evppi.R @@ -15,7 +15,7 @@ #' 3) None of the residual stands out from the basic random pattern of residuals. #' This suggests that there are no outliers. #' -#' The second possible diagnostic is the qqplot for the fitted value. This is a +#' The second possible diagnostic is the Q-Q plot for the fitted value. This is a #' graphical method for comparing the fitted values distributions with the #' assumed underlying normal distribution by plotting their quantiles against #' each other. First, the set of intervals for the quantiles is chosen. A point @@ -28,7 +28,7 @@ #' on a \code{bcea} model. #' @template args-he #' @param plot_type The type of diagnostics to be performed. It can be the 'residual -#' plot' or the 'qqplot plot'. +#' plot' (\code{residuals}) or the Q-Q (quantile-quantile) plot (\code{qqplot}). #' @param interv Specifies the interventions for which diagnostic tests should be #' performed (if there are many options being compared) #' @return Plot @@ -65,7 +65,9 @@ diag.evppi <- function(evppi, } -#' +#' Residual Plot +#' @keywords internal hplot +#' evppi_residual_plot <- function(evppi, he, interv) { @@ -101,6 +103,9 @@ evppi_residual_plot <- function(evppi, } +#' Q-Q Plot +#' @keywords internal hplot +#' #' @importFrom graphics par #' @importFrom stats qqnorm qqline #' diff --git a/R/eib_plot_graph.R b/R/eib_plot_graph.R index 222759dd..20ae07e8 100644 --- a/R/eib_plot_graph.R +++ b/R/eib_plot_graph.R @@ -88,7 +88,7 @@ eib_plot_ggplot <- function(he, data = data.frame("kstar" = he$kstar), colour = "grey50", linetype = 2, - size = 0.5) + + linewidth = 0.5) + scale_linetype_manual( "", labels = graph_params$labels, @@ -166,8 +166,7 @@ eib_plot_plotly <- function(he, n_comp <- length(comparison) - if (is.null(plot_aes$line$types)) - plot_aes$line$types <- rep(1:6, ceiling(he$n_comparisons/6))[1:he$n_comparisons] + plot_aes$line$types <- plot_aes$line$types %||% rep(1:6, ceiling(he$n_comparisons/6))[1:he$n_comparisons] comparisons.label <- paste0(he$interventions[he$ref], " vs ", he$interventions[he$comp]) diff --git a/R/evi.plot.mixedAn.R b/R/evi.plot.mixedAn.R index 85c40194..1041eeb9 100644 --- a/R/evi.plot.mixedAn.R +++ b/R/evi.plot.mixedAn.R @@ -82,9 +82,7 @@ evi.plot.mixedAn <- function(he, alt_legend <- pos base.graphics <- all(pmatch(graph, c("base", "ggplot2")) != 2) - if (is.null(y.limits)){ - y.limits <- range(he$evi, he$evi.star) - } + y.limits <- y.limits %||% range(he$evi, he$evi.star) if (base.graphics) { diff --git a/R/evppi.R b/R/evppi.R index be5719f7..0b1d59ff 100644 --- a/R/evppi.R +++ b/R/evppi.R @@ -39,21 +39,21 @@ #' SPDE-INLA should be plotted. Default set to `FALSE`. #' @param residuals A logical value indicating whether the fitted values for #' the SPDE-INLA method should be outputted. Default set to `TRUE`. -#' @param ... Additional arguments. The default methods to compute the EVPPI -#' are: +#' @param method Character string to select which method to use. The default methods are recommended. +#' However, it is possible (mainly for backward compatibility) to use different methods. +#' @param ... Additional arguments. Details of the methods to compute the EVPPI and their additional arguments are: #' - For single-parameter: -#' GAM regression. +#' - Generalized additive model (GAM) (default). +#' - The method of Strong & Oakley use `method` as string `so`. +#' The user *needs* to also specify the number of "blocks" (e.g. `n.blocks=20`). +#' Note that the multi-parameter version for this method has been deprecated. +#' - The method of Sadatsafavi \emph{et al.} where `method` takes as value a string of either `sad` or `sal`. +#' It is then possible to also specify the number of "separators" (e.g. `n.seps=3`). +#' If none is specified, the default value `n.seps=1` is used. +#' Note that the multi-parameter version for this method has been deprecated. #' - For multi-parameter: -#' INLA/SPDE. However, it is possible (mainly for backward compatibility) to -#' use different methods. For single-parameter, the user can specify the method -#' of Sadatsafavi \emph{et al.} or the method of Strong & Oakley. In order to do so, it -#' is necessary to include the extra parameter \code{method} which takes as -#' value a string \code{"sad"} in the former case and a string \code{"so"} in -#' the latter. In case "sal" is selected, then it is possible to also specify -#' the number of "separators" (e.g. \code{n.seps=3}). If none is specified, the -#' default value \code{n.seps=1} is used. If \code{"so"} is used as method for -#' the calculation of the EVPPI, then the user *needs* to also specify the -#' number of "blocks" (e.g. \code{n.blocks=20}). +#' - INLA/SPDE (default). +#' - Gaussian process regression with `method` of `gp`. #' #' @section GAM regression: #' For multi-parameter, the user can select 3 possible methods. If @@ -141,6 +141,7 @@ #' \insertRef{Heath2016}{BCEA} #' #' @export +#' @md #' #' @examples #' # See Baio G., Dawid A.P. (2011) for a detailed description of the @@ -172,7 +173,8 @@ #' plot(EVPPI.sad) #' #' # Compute the EVPPI using INLA/SPDE -#' x_inla <- evppi(he = m, 39:40, input = inp$mat) +#' if (require("INLA")) +#' x_inla <- evppi(he = m, 39:40, input = inp$mat) #' #' # using GAM regression #' x_gam <- evppi(he = m, 39:40, input = inp$mat, method = "GAM") @@ -181,14 +183,16 @@ #' x_gp <- evppi(he = m, 39:40, input = inp$mat, method = "GP") #' #' # plot results -#' plot(x_inla) +#' if (require("INLA")) plot(x_inla) #' points(x_inla$k, x_inla$evppi, type = "l", lwd = 2, lty = 2) #' points(x_gam$k, x_gam$evppi, type = "l", col = "red") #' points(x_gp$k, x_gp$evppi, type = "l", col = "blue") #' -#' plot(x_inla$k, x_inla$evppi, type = "l", lwd = 2, lty = 2) -#' points(x_gam$k, x_gam$evppi, type = "l", col = "red") -#' points(x_gp$k, x_gp$evppi, type = "l", col = "blue") +#' if (require("INLA")) { +#' plot(x_inla$k, x_inla$evppi, type = "l", lwd = 2, lty = 2) +#' points(x_gam$k, x_gam$evppi, type = "l", col = "red") +#' points(x_gp$k, x_gp$evppi, type = "l", col = "blue") +#' } #' #' data(Smoking) #' treats <- c("No intervention", "Self-help", diff --git a/R/evppi.default.R b/R/evppi.default.R index ac89cd6e..613d4697 100644 --- a/R/evppi.default.R +++ b/R/evppi.default.R @@ -1,626 +1,143 @@ #' @rdname evppi -#' #' @export #' +evppi.default <- function(he, ...) { + stop("No method available", call. = FALSE) +} + + +#' @rdname evppi +#' +#' @importFrom voi evppi +#' @importFrom purrr list_cbind map +#' @examples +#' data(Vaccine, package = "BCEA") +#' treats <- c("Status quo", "Vaccination") +#' bcea_vacc <- bcea(e.pts, c.pts, ref = 2, interventions = treats) +#' inp <- createInputs(vaccine_mat) +#' evppi(bcea_vacc, c("beta.1.", "beta.2."), inp$mat) +#' @export +#' evppi.bcea <- function(he, - param_idx, + param_idx = NULL, input, - N = NULL, + N = NULL, plot = FALSE, - residuals = TRUE, ...) { + residuals = TRUE, + method = NULL, ...) { - if (is.null(colnames(input))) { - colnames(input) <- paste0("theta", seq_len(dim(input)[2])) - } - if (is.numeric(param_idx[1]) || is.integer(param_idx[1])) { - params <- colnames(input)[param_idx] - } else { - params <- param_idx - for (i in seq_along(params)) { - param_idx[i] <- which(colnames(input) == params[i]) - } - class(param_idx) <- "numeric" - } - if (is.null(N)) { - N <- he$n_sim - } + comp_ids <- c(he$comp, he$ref) + outputs <- list(e = he$e[, comp_ids], + c = he$c[, comp_ids], + k = he$k) - robust <- NULL - extra_args <- list(...) + if (!is.null(method)) + method <- tolower(method) - suppress.messages <- - if (!exists("suppress.messages", where = extra_args)) { - FALSE - } else { - extra_args$suppress.messages - } + # allow alternative name for Sadatsafavi method + if (length(method) > 0 && method == "sad") method <- "sal" - extra_args$select <- - if (!exists("select", where = extra_args)) { - if (N >= he$n_sim) { - seq_len(he$n_sim) - } else { - sample(seq_len(he$n_sim), size = N, replace = FALSE) - } + # replace column numbers with names + pars <- + if (all(is.numeric(param_idx))) { + if (any(!param_idx %in% 1:ncol(input))) + stop("Column number(s) not in available parameter inputs") + + param_idx <- names(input)[param_idx] + } else { + param_idx } - inputs <- data.frame(input)[extra_args$select, ] + n_sims <- nrow(input) + n_outputs <- nrow(outputs$c) - # Sets default for method of calculation - # If number of params <=4, then use GAM, if not defaults to INLA/SPDE - - if (!exists("method", where = extra_args)) { + # this way passes on error checking to voi::evppi + if (n_sims == n_outputs) { - extra_args$method <- - if (length(param_idx) <= 4) { - list(rep("GAM", he$n_comparators - 1), - rep("GAM", he$n_comparators - 1)) + # subset number of PSA samples; default all + row_idxs <- + if (is.null(N) || (N >= n_sims)) { + 1:n_sims } else { - list(rep("INLA", he$n_comparators - 1), - rep("INLA", he$n_comparators - 1)) + sample(1:n_sims, size = N, replace = FALSE) } - print(paste("method:", extra_args$method)) - + input <- data.frame(input[row_idxs, ]) + outputs$e <- outputs$e[row_idxs, ] + outputs$c <- outputs$c[row_idxs, ] } - if (!inherits(extra_args$method, "list")) { - - if (extra_args$method != "sad" && extra_args$method != "so") { - if (length(extra_args$method) > 1) { - extra_args$method <- list(extra_args$method, - extra_args$method) - } else { - extra_args$method <- - list(rep(extra_args$method, he$n_comparators - 1), - rep(extra_args$method, he$n_comparators - 1)) - } - } - } + res <- voi::evppi(outputs, inputs = input, pars = pars, method = method, + check = TRUE, plot_inla_mesh = plot, ...) - if (inherits(extra_args$method, "list")) { - - len_methods <- - length(extra_args$method[[1]]) + - length(extra_args$method[[2]]) - - if (len_methods != 2*(he$n_comparators - 1)) { - stop(paste("The argument 'method' must be a list of length 2 with", - he$n_comparators - 1, "elements each."), call. = FALSE) - } - } + voi_methods <- unname(attr(res, "methods")) + voi_models <- attr(res, "models") - if (!exists("int.ord", where = extra_args)) { - extra_args$int.ord <- - list(rep(1, he$n_comparators - 1), - rep(1, he$n_comparators - 1)) - } + # method name returned from evppi + method_nm <- voi_methods[1] - if (!inherits(extra_args$int.ord, "list")) { - - extra_args$int.ord <- - list( - rep(extra_args$int.ord[1], he$n_comparators - 1), - rep(extra_args$int.ord[2], he$n_comparators - 1) - ) - } + form <- + if (method_nm == "gam") { + paste("te(", paste(param_idx, ",", sep = "", + collapse = ""), "bs='cr')") + } else {NULL} - if (!inherits(extra_args$method, "list")) { - if (extra_args$method == "sal" || extra_args$method == "sad") { - method <- "Sadatsafavi et al" - n.blocks <- NULL - - n_seps <- - if (!exists("n_seps", where = extra_args)) { - 1 - } else { - extra_args$n_seps} - - if (length(params) == 1) { - d <- he$n_comparators - n <- he$n_sim - w <- params - param <- inputs[, w] - o <- order(param) - param <- param[o] - nSegs <- matrix(1, d, d) - nSegs[1, 2] <- n_seps - nSegs[2, 1] <- n_seps - res <- segPoints <- numeric() - - for (k in seq_along(he$k)) { - nbs <- he$U[, k, ] - nbs <- nbs[o, ] - for (i in seq_len(d - 1)) { - for (j in (i + 1):d) { - cm <- cumsum(nbs[, i] - nbs[, j])/n - - if (nSegs[i, j] == 1) { - l <- which.min(cm) - u <- which.max(cm) - if (cm[u] - max(cm[1], cm[n]) > min(cm[1], - cm[n]) - cm[l]) { - segPoint <- u - } else { - segPoint <- l - } - if (segPoint > 1 && segPoint < n) { - segPoints <- c(segPoints, segPoint) - } - } - - if (nSegs[i, j] == 2) { - distMaxMin <- 0 - distMinMax <- 0 - minL <- Inf - maxL <- -Inf - for (sims in seq_len(n)) { - if (cm[sims] > maxL) { - maxLP <- sims - maxL <- cm[sims] - } else { - if (maxL - cm[sims] > distMaxMin) { - distMaxMin <- maxL - cm[sims] - segMaxMinL <- maxLP - segMaxMinR <- sims - } - } - if (cm[sims] < minL) { - minLP <- sims - minL <- cm[sims] - } else { - if (cm[sims] - minL > distMinMax) { - distMinMax <- cm[sims] - minL - segMinMaxL <- minLP - segMinMaxR <- sims - } - } - } - siMaxMin <- cm[segMaxMinL] + distMaxMin + (cm[n] - cm[segMaxMinR]) - siMinMax <- -cm[segMaxMinL] + distMinMax - (cm[n] - cm[segMinMaxR]) - - if (siMaxMin > siMinMax) { - segPoint <- c(segMaxMinL, segMaxMinR) - } - else { - segPoint <- c(segMinMaxL, segMinMaxR) - } - if (segPoint[1] > 1 && segPoint[1] < n) { - segPoints <- c(segPoints, segPoint[1]) - } - if (segPoint[2] > 1 && segPoint[2] < n) { - segPoints <- c(segPoints, segPoint[2]) - } - } - } - } - if (length(segPoints) > 0) { - segPoints2 <- unique(c(0, segPoints[order(segPoints)], n)) - res[k] <- 0 - for (j in seq_len(length(segPoints2) - 1)) { - res[k] <- - res[k] + max(colSums( - matrix(nbs[(1 + segPoints2[j]):segPoints2[j + 1], ], - ncol = d)))/n - } - res[k] <- res[k] - max(colMeans(nbs)) - } - else { - res[k] <- 0 - } - } - } - if (length(params) > 1) { - res <- list() - for (lp in seq_along(params)) { - d <- he$n_comparators - n <- he$n_sim - w <- params[lp] - param <- inputs[, w] - o <- order(param) - param <- param[o] - nSegs <- matrix(1, d, d) - nSegs[1, 2] <- n_seps - nSegs[2, 1] <- n_seps - temp <- segPoints <- numeric() - - for (k in seq_along(he$k)) { - nbs <- he$U[, k, ] - nbs <- nbs[o, ] - - for (i in seq_len(d - 1)) { - for (j in (i + 1):d) { - cm <- cumsum(nbs[, i] - nbs[, j])/n - if (nSegs[i, j] == 1) { - l <- which.min(cm) - u <- which.max(cm) - if (cm[u] - max(cm[1], cm[n]) > min(cm[1], cm[n]) - cm[l]) { - segPoint <- u - } else { - segPoint <- l - } - if (segPoint > 1 && segPoint < n) { - segPoints <- c(segPoints, segPoint) - } - } - if (nSegs[i, j] == 2) { - distMaxMin <- 0 - distMinMax <- 0 - minL <- Inf - maxL <- -Inf - for (sims in seq_len(n)) { - if (cm[sims] > maxL) { - maxLP <- sims - maxL <- cm[sims] - } else { - if (maxL - cm[sims] > distMaxMin) { - distMaxMin <- maxL - cm[sims] - segMaxMinL <- maxLP - segMaxMinR <- sims - } - } - if (cm[sims] < minL) { - minLP <- sims - minL <- cm[sims] - } else { - if (cm[sims] - minL > distMinMax) { - distMinMax <- cm[sims] - minL - segMinMaxL <- minLP - segMinMaxR <- sims - } - } - } - siMaxMin <- cm[segMaxMinL] + distMaxMin + (cm[n] - cm[segMaxMinR]) - siMinMax <- -cm[segMaxMinL] + distMinMax - (cm[n] - cm[segMinMaxR]) - if (siMaxMin > siMinMax) { - segPoint <- c(segMaxMinL, segMaxMinR) - } else { - segPoint <- c(segMinMaxL, segMinMaxR) - } - if (segPoint[1] > 1 && segPoint[1] < n) { - segPoints <- c(segPoints, segPoint[1]) - } - if (segPoint[2] > 1 && segPoint[2] < n) { - segPoints <- c(segPoints, segPoint[2]) - } - } - } - } - if (length(segPoints) > 0) { - segPoints2 <- unique(c(0, segPoints[order(segPoints)], n)) - temp[k] <- 0 - for (j in seq_len(length(segPoints2) - 1)) { - temp[k] <- - temp[k] + max(colSums( - matrix(nbs[(1 + segPoints2[j]):segPoints2[j + 1], ], - ncol = d)))/n - } - temp[k] <- temp[k] - max(colMeans(nbs)) - } else { - temp[k] <- 0 - } - } - res[[lp]] <- temp - } - names(res) <- params - } - - res <- list(evppi = res, - index = params, - parameters = params, - k = he$k, - evi = he$evi, - method = method) - } - if (extra_args$method == "so") { - method <- "Strong & Oakley (univariate)" - n_seps <- NULL - if (!exists("n.blocks", where = extra_args)) { - stop("Please specify the param_idx 'n.blocks' to use the Strong and Oakley univariate method", - call. = FALSE) - } else { - n.blocks <- extra_args$n.blocks - } - S <- he$n_sim - J <- S/extra_args$n.blocks - check <- S%%extra_args$n.blocks - if (check > 0) { - stop("number of simulations/number of blocks must be an integer. - Please select a different value for n.blocks \n", call. = FALSE) - } - D <- he$n_comparators - if (length(param_idx) == 1) { - sort.order <- order(inputs[, params]) - sort.U <- array(NA, dim(he$U)) - evpi <- res <- numeric() - for (i in seq_along(he$k)) { - evpi[i] <- he$evi[i] - sort.U[, i, ] <- he$U[sort.order, i, ] - U.array <- array(sort.U[, i, ], - dim = c(J, extra_args$n.blocks, D)) - mean.k <- apply(U.array, c(2, 3), mean) - partial.info <- mean(apply(mean.k, 1, max)) - res[i] <- partial.info - max(apply(he$U[, i, ], 2, mean)) - } - } - if (length(param_idx) > 1) { - res <- list() - for (j in seq_along(param_idx)) { - sort.order <- order(inputs[, params[j]]) - sort.U <- array(NA, dim(he$U)) - evpi <- evppi.temp <- numeric() - - for (i in seq_along(he$k)) { - evpi[i] <- he$evi[i] - sort.U[, i, ] <- he$U[sort.order, i, ] - U.array <- array(sort.U[, i, ], - dim = c(J, n.blocks, D)) - mean.k <- apply(U.array, c(2, 3), mean) - partial.info <- mean(apply(mean.k, 1, max)) - evppi.temp[i] <- partial.info - max(apply(he$U[, i, ], 2, mean)) - } - res[[j]] <- evppi.temp - } - names(res) <- params - } - - res <- - list(evppi = res, - index = params, - parameters = params, - k = he$k, - evi = he$evi, - method = method) - } - } + # fitted values + get_fitted_values <- residuals && (method_nm %in% c("inla", "gp", "gam")) - if (inherits(extra_args$method, "list")) { - time <- list() - time[[1]] <- list() - time[[2]] <- list() - - fit.full <- vector("list") - fit.full[[1]] <- matrix( - data = 0, - nrow = length(extra_args$select), - ncol = he$n_comparators) - - fit.full[[2]] <- - matrix(data = 0, - nrow = length(extra_args$select), - ncol = he$n_comparators) - - for (k in 1:2) { - for (l in seq_len(he$n_comparisons)) { - - x <- - prep.x(he = he, - seq_rows = extra_args$select, - k = k, - l = l) - - method <- toupper(extra_args$method[[k]][l]) - - if (method == "GAM" || method == "G") { - method <- "GAM" - mesh <- robust <- NULL - if (!requireNamespace("mgcv", quietly = TRUE)) { - stop("You need to install the package 'mgcv'. Please run in your R terminal:\n - install.packages('mgcv')", call. = FALSE) - } - if (requireNamespace("mgcv", quietly = TRUE)) { - if (!suppress.messages) { - cat("\n") - cat("Calculating fitted values for the GAM regression \n") - } - - inp <- names(inputs)[param_idx] - - form <- - if (exists("formula", where = extra_args)) { - extra_args$formula - } else { - paste("te(", paste(inp, ",", sep = "", - collapse = ""), "bs='cr')") - } - - fit <- fit.gam(parameter = param_idx, - inputs = inputs, - x = x, - form = form) - } - } - if (method == "GP") { - mesh <- NULL - if (!suppress.messages) { - cat("\n") - cat("Calculating fitted values for the GP regression \n") - # If the number of simulations to be used to estimate the - # hyper-params is set then use that, else use N/2 - } - - n_sim <- - if (!exists("n_sim", where = extra_args)) { - N/2 - } else { - extra_args$n_sim - } - - fit <- fit.gp(parameter = param_idx, - inputs = inputs, - x = x, - n.sim = n_sim) - } - if (method == "INLA") { - if (!requireNamespace("INLA", quietly = TRUE)) { - stop("You need to install the packages 'INLA'. Please run in your R terminal:\n - install.packages('INLA',repos=c(getOption('repos'),INLA='https://inla.r-inla-download.org/R/stable/'), dep=TRUE)\n - If you are under MS Windows, you may need to install the most recent version\n - of 'Rtools' - see https://cran.r-project.org/bin/windows/Rtools/", call. = FALSE) - } - - if (requireNamespace("INLA", quietly = TRUE)) { - if (!is.element("INLA", (.packages()))) { - attachNamespace("INLA") - } - if (length(param_idx) < 2) { - stop("The INLA method can only be used with 2 or more params", call. = FALSE) - } - if (!suppress.messages) { - cat("\n") - cat("Finding projections \n") - } - projections <- make.proj(parameter = param_idx, - inputs = inputs, - x = x, - k = k, - l = l) - data <- projections$data - if (!suppress.messages) { - cat("Determining Mesh \n") - } - - cutoff <- - if (!exists("cutoff", where = extra_args)) { - 0.3 - } else { - extra_args$cutoff - } - - convex.inner <- - if (!exists("convex.inner", where = extra_args)) { - 0.4 - } else { - extra_args$convex.inner - } - - convex.outer <- - if (!exists("convex.outer", where = extra_args)) { - -0.7 - } else { - extra_args$convex.outer - } - - max.edge <- - if (!exists("max.edge", where = extra_args)) { - 0.7 - } else { - extra_args$max.edge - } - - mesh <- - make.mesh( - data = data, - convex.inner = convex.inner, - convex.outer = convex.outer, - cutoff = cutoff, - max.edge = max.edge - ) - plot.mesh(mesh = mesh$mesh, - data = data, - plot = plot) - if (!suppress.messages) { - cat("Calculating fitted values for the GP regression using INLA/SPDE \n") - } - if (exists("h.value", where = extra_args)) { - h.value <- extra_args$h.value - } - else { - h.value <- 5e-05 - } - if (exists("robust", where = extra_args)) { - if (extra_args$robust) { - family <- "T" - robust <- TRUE - } - else { - family <- "gaussian" - robust <- FALSE - } - } else { - family <- "gaussian" - robust <- FALSE - } - if (exists("int.ord", where = extra_args)) { - int.ord <- extra_args$int.ord[[k]][l] - } - else { - int.ord <- 1 - } - fit <- fit.inla( - parameter = param_idx, - inputs = inputs, - x = x, - mesh = mesh$mesh, - data.scale = data, - int.ord = int.ord, - convex.inner = convex.inner, - convex.outer = convex.outer, - cutoff = cutoff, - max.edge = max.edge, - h.value = h.value, - family = family - ) - } - } - fit.full[[k]][,l] <- fit$fitted - - # calculate time taken - if (method == "INLA") { - time. <- c(projections$time, mesh$time, fit$time) - time. <- sum(time.) - time[[k]][l] <- (time.) - } else { - time. <- fit$time - time[[k]][l] <- time. - } - } - } - if (!suppress.messages) cat("Calculating EVPPI \n") + # named differently in different methods + fitted_lup <- c(inla = "fitted", gp = "fitted", gam = "fitted.values") + + if (get_fitted_values) { + fitted_txt <- fitted_lup[method_nm] - comp <- compute.evppi(he = he, fit.full = fit.full) + fitted_c_ls <- + purrr::map(voi_models[[1]]$c, + ~unname(as.data.frame(.x[[fitted_txt]]))) |> + rev() - name <- prepare.output(parameters = params, inputs = inputs) + fitted_c <- + suppressMessages(purrr::list_cbind(fitted_c_ls)) |> + as.matrix() |> + cbind(0) - time[[3]] <- comp$time - names(time) <- c("Fitting for Effects", - "Fitting for Costs", - "Calculating EVPPI") - names(extra_args$method) <- c("Methods for Effects", "Methods for Costs") + fitted_e_ls <- + purrr::map(voi_models[[1]]$e, + ~unname(as.data.frame(.x[[fitted_txt]]))) |> + rev() - if (residuals) { - res <- list( - evppi = comp$EVPPI, - index = params, - k = he$k, - evi = he$evi, - parameters = name, - time = time, - method = extra_args$method, - fitted.costs = fit.full[[2]], - fitted.effects = fit.full[[1]], - select = extra_args$select) - } else { - res <- list( - evppi = comp$EVPPI, - index = params, - k = he$k, - evi = he$evi, - parameters = name, - time = time, - method = extra_args$method) - } + fitted_e <- + suppressMessages(purrr::list_cbind(fitted_e_ls)) |> + as.matrix() |> + cbind(0) + } else { + fitted_c <- NULL + fitted_e <- NULL } - structure(res, class = "evppi") -} - - -#' @rdname evppi -#' @export -#' -evppi.default <- function(he, ...) { - stop("No method available", call. = FALSE) + residuals_out <- + list(fitted.costs = fitted_c, + fitted.effects = fitted_e, + select = row_idxs) + + out <- c( + list( + evppi = res$evppi, + index = pars, + k = res$k, + evi = he$evi, + parameters = paste(pars, collapse = " and "), + time = list("Fitting for Effects" = NULL, + "Fitting for Costs" = NULL, + "Calculating EVPPI" = NULL), + method = list("Methods for Effects" = voi_methods, + "Methods for Costs" = voi_methods)), + residuals_out, + list(formula = form, + pars = res$pars, + res = res)) + + structure(out, class = c("evppi", class(out))) } diff --git a/R/evppi2stage.R b/R/evppi2stage.R deleted file mode 100644 index 9b853ac0..00000000 --- a/R/evppi2stage.R +++ /dev/null @@ -1,105 +0,0 @@ -# ##TODO:... -# -# #' Two-stage EVPPI -# #' -# #' Performs the (rather computationally expensive) computations for the EVPPI. Effectively, considers -# #' the vector of parameters theta=(phi,psi), where phi is the parameter of interest and psi is -# #' the rest of the parameters. First needs to run the Bayesian model assuming psi as random variables -# #' and taking phi as a fixed value. Then performs the economic analysis calling BCEA and then -# #' computes some relevant syntheses. -# #' -# #' @param phi a vector with simulations from the posterior distribution of the the parameter with respect -# #' to which the EVPPI is computed (phi), obtained by the full model -# #' @param phi.name a string containing the name of the parameter phi, as represented in the txt file -# #' of the full model for psi. In this case phi is considered as fixed. This is has to be -# #' the same as one of the elements of the list data! -# #' @param n.out number of simulations used to compute the average with respect to phi (default = 1000) -# #' @param m a BCEA object in which the full economic analysis is stored as a result of the call to bcea -# #' @param data a list with the data to be passed to JAGS to run the model for the vector of parameters psi -# #' (conditionally on each value of the fixed parameter phi) -# #' @param model.file a txt file containing the full model for psi (but considering the parameter phi as fixed. -# #' This is in general a modification of the original model file) -# #' @param params a collection of parameters to be monitored in the conditional model for psi (given phi) -# #' @param inits specifies the initial values for the MCMC procedure (default = NULL) -# #' @param n.iter number of iterations to be run -# #' @param n.burnin number of iterations to be discarded as burn-in -# #' @param n.thin number of iterations to be considered for thinning (default NULL) -# #' @param n.mc number of simulations retained to compute the MC estimations (default=1000) -# #' @param working.dir a string containing the working directory (default is NULL) -# #' @param ec.vars a string with the name of the R file that contains the commands required to -# #' compute the relevant variables of costs and benefits. These are usually -# #' obtained as functions of the parameters estimated by the full Bayesian model -# #' @param bcea.call a string containing the call to BCEA to perform the C/E evaluation. By default -# #' assumes no options and the existence of two vectors e,c in which suitable -# #' values for the measures of effectiveness and costs are stored -# #' -# #' @return List -# #' -# evppi2stage <- function(phi, -# phi.name, -# n.out = 1000, -# m, -# data, -# model.file, -# params, -# inits = NULL, -# n.iter, -# n.burnin, -# n.thin = NULL, -# n.mc = 1000, -# working.dir = NULL, -# ec.vars, -# bcea.call = "bcea(eff, cost)") { -# -# # Replicates the C/E analysis for a fixed value of the relevant parameter at each iteration -# Ustar.phi <- -# matrix(NA, n.out, length(m$k)) # matrix of conditional net benefits -# -# for (i in 1:n.out) { -# cmd <- paste(phi.name, " <- phi[i]", sep = "") -# eval(parse(text = cmd)) -# if (is.null(n.thin)) { -# n.thin <- floor((n.iter - n.burnin) / (n.mc / 2)) -# } -# -# model.evppi <- -# jags(data, -# inits, -# params, -# model.file = model.file, -# n.chains = 2, -# n.iter, -# n.burnin, -# n.thin, -# DIC = TRUE, -# working.directory = working.dir, -# progress.bar = "text") -# -# attach.bugs(model.evppi$BUGSoutput) -# -# # performs the economic analysis given the simulations for psi conditionally on phi -# cmd <- paste("source('", ec.vars, "')", sep = "") -# eval(parse(text = cmd)) -# -# # runs BCEA to perform the economic analysis conditionally on phi -# m.evppi <- eval(parse(text = bcea.call)) -# -# # now computes the maximum expected utility for the current iteration for each value of k -# Ustar.phi[i, ] <- apply(apply(m.evppi$U, c(2, 3), mean), 1, max) -# rm(m.evppi) -# detach.bugs() -# print(i) # prints a counter to show the progression in the analysis -# } -# -# # Computes the average value of the maximum expected utility for each value of k -# Vstar.phi <- apply(Ustar.phi, 2, mean) -# Umax <- apply(apply(m$U, c(2, 3), mean), 1, max) -# -# # Finally computes the EVPPI for each value of k -# EVPPI <- Vstar.phi - Umax -# -# list(evppi = EVPPI, -# Vstar.phi = Vstar.phi, -# Ustar.phi = Ustar.phi) -# } -# diff --git a/R/evppi_helpers.R b/R/evppi_helpers.R deleted file mode 100644 index 5ecd7f08..00000000 --- a/R/evppi_helpers.R +++ /dev/null @@ -1,574 +0,0 @@ -# evppi() helper functions ------------------------------------------------ - - -#' Prepare Delta arrays -#' -#' @template args-he -#' @param seq_rows Rows of (e,c) to keep -#' @param k e or c? 1 or 2. -#' @param l Columns of (e,c) to keep -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -prep.x <- function(he, - seq_rows, - k, - l){ - if (k == 1) { - x <- as.matrix(he$delta_e)[seq_rows, l] - } - if (k == 2) { - x <- as.matrix(he$delta_c)[seq_rows, l] - } - return(x) -} - - -#' Gaussian Additive Model Fitting -#' -#' @param parameter Parameter -#' @param inputs Inputs -#' @param x Response variable -#' @param form Formula -#' -#' @return List -#' -#' @importFrom stats update -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -fit.gam <- function(parameter, - inputs, - x, - form) { - tic <- proc.time() - model <- mgcv::gam(update(formula(x ~ .), - formula(paste(".~", form))), - data = data.frame(inputs)) - hat <- model$fitted - - ##TODO: should this be used instead of hat? - # N <- nrow(inputs) - # p <- length(parameter) - # fitted <- matrix(hat, nrow = N, ncol = p) - - formula <- form - toc <- proc.time() - tic - time <- toc[3] - names(time) <- "Time to fit GAM regression (seconds)" - - list(fitted = hat, - formula = formula, - fit = model, - time = time) -} - -#' Gaussian Process Fitting -#' -#' @param hyperparams Hyperparameters -#' @param parameter Parameters -#' @param x Response variable -#' @param input.matrix Input data matrix -#' -#' @importFrom stats dist dnorm -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -post.density <- function(hyperparams, - parameter, - x, - input.matrix) { - - dinvgamma <- function(x, alpha, beta) { - (beta^alpha)/gamma(alpha) * - x^(-alpha - 1) * - exp(-beta/x) - } - - N <- length(x) - p <- length(parameter) - H <- cbind(1, input.matrix) - q <- ncol(H) - a.sigma <- 0.001 - b.sigma <- 0.001 - a.nu <- 0.001 - b.nu <- 1 - delta <- exp(hyperparams)[1:p] - nu <- exp(hyperparams)[p + 1] - A <- exp(-(as.matrix(dist(t(t(input.matrix)/delta), - upper = TRUE, diag = TRUE))^2)) - Astar <- A + nu * diag(N) - chol_Astar <- chol(Astar) - y <- backsolve(t(chol_Astar),(x), upper.tri = FALSE) - x. <- backsolve(t(chol_Astar), H, upper.tri = FALSE) - tHAstarinvH <- t(x.) %*% (x.) - betahat <- solve(tHAstarinvH) %*% t(x.) %*% y - residSS <- y %*% y - t(y) %*% x. %*% betahat - t(betahat) %*% - t(x.) %*% y + t(betahat) %*% tHAstarinvH %*% betahat - prior <- prod(dnorm(log(delta), 0, sqrt(1e+05))) * - dinvgamma(nu, a.nu, b.nu) - l <- -sum(log(diag(chol_Astar))) - 1/2 * log(det(tHAstarinvH)) - - (N - q + 2 * a.sigma)/2 * log(residSS/2 + b.sigma) + log(prior) - names(l) <- NULL - - return(l) -} - - -#' Estimate hyperparameters -#' -#' @param x x -#' @param input.matrix Input matrix -#' @param parameter Parameters -#' @param n.sim Number of simulations -#' -#' @importFrom stats optim -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -estimate.hyperparams <- function(x, - input.matrix, - parameter, - n.sim) { - p <- length(parameter) - initial.values <- rep(0, p + 1) - repeat { - log.hyperparams <- - optim(initial.values, - fn = post.density, - parameter = parameter, - x = x[1:n.sim], - input.matrix = input.matrix[1:n.sim, ], - method = "Nelder-Mead", - control = list(fnscale = -1, - maxit = 10000, trace = 0))$par - - if (sum(abs(initial.values - log.hyperparams)) < 0.01) { - hyperparams <- exp(log.hyperparams) - break - } - initial.values <- log.hyperparams - } - - hyperparams -} - -#' Fit Gaussian Process -#' -#' @param parameter Parameters -#' @param inputs Inputs -#' @param x Response variable -#' @param n.sim Number of simulations -#' @return list -#' @importFrom stats dist -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -fit.gp <- function(parameter, - inputs, - x, - n.sim) { - tic <- proc.time() - p <- length(parameter) - input.matrix <- as.matrix(inputs[, parameter, drop = FALSE]) - colmin <- apply(input.matrix, 2, min) - colmax <- apply(input.matrix, 2, max) - colrange <- colmax - colmin - input.matrix <- sweep(input.matrix, 2, colmin, "-") - input.matrix <- sweep(input.matrix, 2, colrange, "/") - N <- nrow(input.matrix) - H <- cbind(1, input.matrix) - q <- ncol(H) - - hyperparams <- - estimate.hyperparams( - x = x, - input.matrix = input.matrix, - parameter = parameter, - n.sim = n.sim) - - delta.hat <- hyperparams[1:p] - nu.hat <- hyperparams[p + 1] - A <- exp(-(as.matrix(dist(t(t(input.matrix)/delta.hat), - upper = TRUE, diag = TRUE))^2)) - Astar <- A + nu.hat * diag(N) - Astarinv <- chol2inv(chol(Astar)) - - rm(Astar) - gc() - - AstarinvY <- Astarinv %*% x - tHAstarinv <- t(H) %*% Astarinv - tHAHinv <- solve(tHAstarinv %*% H) - betahat <- tHAHinv %*% (tHAstarinv %*% x) - Hbetahat <- H %*% betahat - resid <- x - Hbetahat - fitted <- Hbetahat + A %*% (Astarinv %*% resid) - - ##TODO: should this be used somewhere? - # AAstarinvH <- A %*% t(tHAstarinv) - - sigmasqhat <- - as.numeric(t(resid) %*% Astarinv %*% resid)/(N - q - 2) - - rm(A, Astarinv, AstarinvY, tHAstarinv, tHAHinv, - Hbetahat, resid, sigmasqhat) - gc() - - toc <- proc.time() - tic - time <- toc[3] - names(time) <- "Time to fit GP regression (seconds)" - - list(fitted = fitted, - time = time, - fit = NULL, - formula = NULL) -} - - -#' INLA Fitting -#' -#' @param parameter Parameter -#' @param inputs Inputs -#' @param x Response variable -#' @param k k -#' @param l l -#' @return list -#' @seealso \code{\link{evppi}} -#' @importFrom cli cli_alert_warning -#' @keywords internal -#' -make.proj <- function(parameter, - inputs, - x, k, l) { - tic <- proc.time() - scale <- 8 / (range(x)[2] - range(x)[1]) - scale.x <- scale * x - mean(scale * x) - - # generate basic function - bx <- bf(scale.x, case = "poly", 2) - - # principle fitted components - fit1 <- - pfc(scale(inputs[, parameter]), scale.x, bx, structure = "iso") - fit2 <- - pfc(scale(inputs[, parameter]), scale.x, bx, structure = "aniso") - fit3 <- - pfc(scale(inputs[, parameter]), scale.x, bx, structure = "unstr") - - struc <- - c("iso", "aniso", "unstr")[ - which(c(fit1$aic, fit2$aic, fit3$aic) == min(fit1$aic, fit2$aic, fit3$aic))] - AIC.deg <- array() - - for (i in 2:7) { - bx <- bf(scale.x, case = "poly", i) - fit <- - pfc(scale(inputs[, parameter]), scale.x, bx, structure = struc) - AIC.deg[i] <- fit$aic - } - - deg <- which(AIC.deg == min(AIC.deg, na.rm = TRUE)) - d <- min(dim(inputs[, parameter])[2], deg) - - ##TODO: where should this be used? - # by <- bf(scale.x, case = "poly", deg) - - # likelihood-based dimension reduction - comp.d <- - ldr( - scale(inputs[, parameter]), - scale.x, - bx, - structure = struc, - model = "pfc", - numdir = d, - numdir.test = TRUE) - - dim.d <- which(comp.d$aic == min(comp.d$aic)) - 1 - - comp <- - ldr( - scale(inputs[, parameter]), - scale.x, - bx, - structure = struc, - model = "pfc", - numdir = 2) - - toc <- proc.time() - tic - time <- toc[3] - - if (dim.d > 2){ - cur <- c("effects", "costs") - - cli::cli_alert_warning( - "The dimension of the sufficient reduction for the incremental - {.code cur[k]} column {.code l} is {.code dim.d}. - Dimensions greater than 2 imply that the EVPPI approximation using INLA may be inaccurate. - Full residual checking using {.fn diag.evppi} is required.") - } - names(time) <- "Time to fit find projections (seconds)" - - list(data = comp$R, - time = time, - dim = dim.d) -} - - -#' Mesh Plot -#' -#' Option of interactively saving the plot. -#' -#' @param mesh Mesh -#' @param data Data -#' @param plot Create plot? logical -#' -#' @importFrom utils select.list -#' @importFrom grDevices dev.off -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -plot.mesh <- function(mesh, data, plot) { - - if (!plot) return() - - cat("\n") - choice <- select.list(c("yes", "no"), - title = "Would you like to save the graph?", - graphics = FALSE) - if (choice == "yes") { - exts <- c("jpeg", "pdf", "bmp", "png", "tiff") - ext <- select.list(exts, - title = "Please select file extension", - graphics = FALSE) - name <- paste0(getwd(), "/mesh.", ext) - txt <- paste0(ext, "('", name, "')") - eval(parse(text = txt)) - txt <- paste0("Graph saved as: ", name) - cat(txt) - on.exit(dev.off()) - } - - cat("\n") - plot(mesh) - points(data, - col = "blue", - pch = 19, - cex = 0.8) -} - - -#' Make Mesh -#' -#' Fit using INLA methods. -#' -#' @param data Data -#' @param convex.inner convex.inner -#' @param convex.outer convex.outer -#' @param cutoff Cut-off value -#' @param max.edge Maximum edge -#' @return list -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -make.mesh <- function(data, - convex.inner, - convex.outer, - cutoff, - max.edge) { - tic <- proc.time() - inner <- suppressMessages({ - INLA::inla.nonconvex.hull(data, convex = convex.inner) - }) - outer <- INLA::inla.nonconvex.hull(data, convex = convex.outer) - mesh <- - INLA::inla.mesh.2d( - loc = data, - boundary = list(inner, outer), - max.edge = c(max.edge, max.edge), - cutoff = c(cutoff)) - toc <- proc.time() - tic - time <- toc[3] - names(time) <- "Time to fit determine the mesh (seconds)" - - list(mesh = mesh, - pts = data, - time = time) -} - - -#' Fit INLA -#' -#' @param parameter Parameters -#' @param inputs Inputs -#' @param x Response variable -#' @param mesh Mesh -#' @param data.scale data.scale -#' @param int.ord int.ord -#' @param convex.inner convex.inner -#' @param convex.outer convex.outer -#' @param cutoff Cut-off -#' @param max.edge Maximum edge -#' @param h.value h.value -#' @param family family -#' @return list -#' @importFrom stats as.formula -#' -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -fit.inla <- function(parameter, - inputs, - x, - mesh, - data.scale, - int.ord, - convex.inner, - convex.outer, - cutoff, - max.edge, - h.value, - family) { - # @importFrom INLA inla.stack.data - # @importFrom INLA inla.stack.A - # @importFrom INLA inla - - tic <- proc.time() - inputs.scale <- - scale(inputs, apply(inputs, 2, mean), apply(inputs, 2, sd)) - scale <- 8 / (range(x)[2] - range(x)[1]) - scale.x <- scale * x - mean(scale * x) - A <- - INLA::inla.spde.make.A(mesh = mesh, - loc = data.scale, - silent = 2L) - spde <- - INLA::inla.spde2.matern( - mesh = mesh, - alpha = 2) - stk.real <- - INLA::inla.stack( - tag = "est", - data = list(y = scale.x), - A = list(A, 1), - effects = list(s = 1:spde$n.spde, - data.frame( - b0 = 1, x = cbind(data.scale, inputs.scale) - )) - ) - data <- INLA::inla.stack.data(stk.real) - ctr.pred <- INLA::inla.stack.A(stk.real) - - inp <- names(stk.real$effects$data)[parameter + 4] - form <- paste(inp, "+", sep = "", collapse = "") - formula <- paste("y~0+(", - form, - "+0)+b0+f(s,model=spde)", - sep = "", - collapse = "") - if (int.ord[1] > 1) { - formula <- paste( - "y~0+(", - form, - "+0)^", - int.ord[1], - "+b0+f(s,model=spde)", - sep = "", - collapse = "") - } - Result <- suppressMessages({ - INLA::inla( - as.formula(formula), - data = data, - family = family, - control.predictor = list(A = ctr.pred, link = 1), - control.inla = list(h = h.value), - control.compute = list(config = TRUE) - ) - }) - fitted <- - (Result$summary.linear.predictor[seq_along(x), "mean"] + mean(scale * x)) / scale - fit <- Result - toc <- proc.time() - tic - time <- toc[3] - names(time) <- "Time to fit INLA/SPDE (seconds)" - - list( - fitted = fitted, - model = fit, - time = time, - formula = formula, - mesh = list(mesh = mesh, pts = data.scale)) -} - - -#' Compute EVPPI -#' @template args-he -#' @param fit.full fit.full -#' @return list -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -compute.evppi <- function(he, fit.full) { - EVPPI <- array() - tic <- proc.time() - for (i in seq_along(he$k)) { - NB.k <- -(he$k[i]*fit.full[[1]] - fit.full[[2]]) - EVPPI[i] <- (mean(apply(NB.k, 1, max, na.rm = TRUE)) - - max(apply(NB.k, 2, mean, na.rm = TRUE))) - } - toc <- proc.time() - tic - time <- toc[3] - names(time) <- "Time to compute the EVPPI (in seconds)" - - list(EVPPI = EVPPI, - time = time) -} - - -#' Prepare output -#' -#' @param parameters Parameters -#' @param inputs Inputs -#' @return name -#' @seealso \code{\link{evppi}} -#' @keywords internal -#' -prepare.output <- function(parameters, - inputs) { - - if (length(parameters) == 1) { - - name <- - if (is.numeric(parameters)) { - colnames(inputs)[parameters] - } else { - parameters} - - } else { - if (is.numeric(parameters)) { - n_params <- length(parameters) - end <- colnames(inputs)[parameters[n_params]] - name.mid <- - paste(colnames(inputs)[parameters[1:n_params - 1]], - ", ", - sep = "", - collapse = " ") - name <- paste(name.mid, "and ", end, sep = "", - collapse = " ") - } else { - n_params <- length(parameters) - end <- parameters[n_params] - name.mid <- paste(parameters[1:n_params - 1], ", ", - sep = "", - collapse = " ") - name <- paste(name.mid, "and ", end, - sep = "", - collapse = " ") - } - } - - name -} - diff --git a/R/evppi_plot_base.R b/R/evppi_plot_base.R index f406fa6f..08fb2f09 100644 --- a/R/evppi_plot_base.R +++ b/R/evppi_plot_base.R @@ -20,7 +20,6 @@ evppi_plot_base <- function(evppi_obj, pos_legend, col = NULL, annot = FALSE) { - legend_params <- evppi_legend_base(evppi_obj, pos_legend, col) diff --git a/R/evppi_base_params.R b/R/evppi_plot_helpers.R similarity index 100% rename from R/evppi_base_params.R rename to R/evppi_plot_helpers.R diff --git a/R/helper_base_params.R b/R/helper_base_params.R index 67c64a7b..a173b3d4 100644 --- a/R/helper_base_params.R +++ b/R/helper_base_params.R @@ -1,6 +1,7 @@ #' @importFrom grDevices colors #' @keywords dplot +#' helper_base_params <- function(he, graph_params) { diff --git a/R/helper_ggplot_params.R b/R/helper_ggplot_params.R index 1b6c517c..8ccc6d4c 100644 --- a/R/helper_ggplot_params.R +++ b/R/helper_ggplot_params.R @@ -35,7 +35,7 @@ helper_ggplot_params <- function(he, is_enough_types <- length(types) >= n_lines is_enough_colours <- length(cols) >= n_lines - is_enough_sizes <- length(sizes) >= n_lines + is_enough_sizes <- length(sizes) == n_lines if (!is_enough_types) { graph_params$line$type <- rep_len(types, n_lines)} diff --git a/R/ib_plot_base.R b/R/ib_plot_base.R index ffdf5692..bf96d1a3 100644 --- a/R/ib_plot_base.R +++ b/R/ib_plot_base.R @@ -52,9 +52,8 @@ ib_plot_base <- function(he, sep = "") } if (he$n_comparisons > 1) { - if (is.null(comparison)) { - comparison <- 1 - } + + comparison <- comparison %||% 1 # nbw <- sd(he$ib[w, , comparison])/1.5 @@ -65,9 +64,9 @@ ib_plot_base <- function(he, he$interventions[he$comp[comparison]], sep = "") } - if (is.null(xlim)) { - xlim <- range(d$x) - } + + xlim <- xlim %||% range(d$x) + plot( d$x, d$y, diff --git a/R/ib_plot_ggplot.R b/R/ib_plot_ggplot.R index a1db649c..4daa6c6a 100644 --- a/R/ib_plot_ggplot.R +++ b/R/ib_plot_ggplot.R @@ -27,9 +27,7 @@ ib_plot_ggplot <- function(he, n, xlim) { - if (is.null(comparison)) { - comparison <- 1 - } + comparison <- comparison %||% 1 if (max(he$k) < wtp) { wtp <- max(he$k) @@ -75,9 +73,9 @@ ib_plot_ggplot <- function(he, df <- data.frame(x = density$x, y = density$y) } - if (is.null(xlim)) { - xlim <- range(df$x) - } + + xlim <- xlim %||% range(df$x) + ib <- ggplot(df, aes(.data$x, .data$y)) + theme_bw() + diff --git a/R/inforank_params.R b/R/inforank_params.R index 3e42fec0..cdfa5fe4 100644 --- a/R/inforank_params.R +++ b/R/inforank_params.R @@ -24,6 +24,7 @@ inforank_params <- function(he, ##TODO: what to do for multiple ICER? if (is.null(wtp)) wtp <- he$k[min(which(he$k >= he$ICER[1]))] + if (!wtp %in% he$k) stop("wtp must be from values computed in call to bcea().", call. = FALSE) diff --git a/R/misc_helpers.R b/R/misc_helpers.R new file mode 100644 index 00000000..9d9f261a --- /dev/null +++ b/R/misc_helpers.R @@ -0,0 +1,5 @@ + +# +`%||%` <- function(x, y) { + if (is.null(x)) y else x +} diff --git a/R/mixedAn.default.R b/R/mixedAn.default.R index 5c63a3cd..268bde37 100644 --- a/R/mixedAn.default.R +++ b/R/mixedAn.default.R @@ -6,9 +6,7 @@ Ubar <- OL.star <- evi.star <- NULL - if (is.null(value)) { - value <- rep(1, he$n_comparators)/he$n_comparators - } + value <- value %||% rep(1, he$n_comparators)/he$n_comparators Ubar <- compute_Ubar(he, value) OL.star <- he$Ustar - Ubar diff --git a/R/multiplot.R b/R/multiplot.R index a0059268..08e79a4e 100644 --- a/R/multiplot.R +++ b/R/multiplot.R @@ -17,11 +17,10 @@ multiplot <- function(plotlist = NULL, n_plots <- length(plotlist) - if (is.null(layout_config)) { - layout_config <- matrix(seq(1, cols*ceiling(n_plots/cols)), - ncol = cols, - nrow = ceiling(n_plots/cols)) - } + layout_config <- + layout_config %||% matrix(seq(1, cols*ceiling(n_plots/cols)), + ncol = cols, + nrow = ceiling(n_plots/cols)) grid_params <- c(plotlist, list(layout_matrix = layout_config)) diff --git a/R/num_lines.R b/R/num_lines.R index b3cf6557..ab7ac92d 100644 --- a/R/num_lines.R +++ b/R/num_lines.R @@ -1,4 +1,5 @@ -# num_lines +# Number of lines +# which depends on the type of plot #' @title Get number of lines diff --git a/R/orthonorm.R b/R/orthonorm.R index 907266e6..e902f1c1 100644 --- a/R/orthonorm.R +++ b/R/orthonorm.R @@ -18,6 +18,7 @@ orthonorm <- function (u) { if (is.null(u)) return(NULL) + if (!(is.matrix(u))) u <- as.matrix(u) dd <- dim(u) diff --git a/R/pfc.R b/R/pfc.R index 1d3b3f28..02c7178d 100644 --- a/R/pfc.R +++ b/R/pfc.R @@ -21,417 +21,425 @@ # #' @importFrom stats cov #' -pfc <- - function(X, y, fy=NULL, numdir=NULL, structure=c("iso", "aniso", "unstr", "unstr2"), eps_aniso=1e-3, numdir.test=FALSE,...) - { - "%^%"<-function(M, pow) - { - if (prod(dim(M)==list(1,1))) return( as.matrix(M^pow) ) - eigenM = eigen(M) - return(eigenM$vectors%*%diag(c(eigenM$values)^pow)%*%t(eigenM$vectors)) - } - - Trace<-function(X) - { - if (!is.matrix(X)) stop("Argument to Trace is not a matrix in pfc") - return(sum(diag(X))) - } - if (is.null(fy)) {fy <- scale(y, TRUE, TRUE); numdir <- 1} - r <- dim(fy)[2] - X <- as.matrix(X) - op <- dim(X) - n <- op[1] - p <- op[2] - eff.numdir <- min(numdir, r, p) - - vnames <- dimnames(X)[[2]] - if (is.null(vnames)) vnames <- paste("X", 1:p, sep="") - if (p==1) return(onepfc(X=X, y=y, fy=fy, p, numdir.test)) - - Muhat <- apply(X, 2, mean) - Xc <- scale(X, TRUE, FALSE) - - P_F <- fy%*%solve(t(fy)%*%fy)%*%t(fy) - - Sigmahat <- cov(X) - Sigmahat_fit <- cov(P_F%*%X) - Sigmahat_res <- Sigmahat - Sigmahat_fit - - if (structure=="iso") +pfc <- function(X, + y, + fy = NULL, + numdir + = NULL, + structure = c("iso", "aniso", "unstr", "unstr2"), + eps_aniso = 1e-3, + numdir.test = FALSE, + ...) { + + "%^%"<-function(M, pow) + { + if (prod(dim(M)==list(1,1))) return( as.matrix(M^pow) ) + eigenM = eigen(M) + return(eigenM$vectors%*%diag(c(eigenM$values)^pow)%*%t(eigenM$vectors)) + } + + Trace<-function(X) + { + if (!is.matrix(X)) stop("Argument to Trace is not a matrix in pfc") + return(sum(diag(X))) + } + + if (is.null(fy)) {fy <- scale(y, TRUE, TRUE); numdir <- 1} + r <- dim(fy)[2] + X <- as.matrix(X) + op <- dim(X) + n <- op[1] + p <- op[2] + eff.numdir <- min(numdir, r, p) + + vnames <- dimnames(X)[[2]] + if (is.null(vnames)) vnames <- paste("X", 1:p, sep="") + if (p==1) return(onepfc(X=X, y=y, fy=fy, p, numdir.test)) + + Muhat <- apply(X, 2, mean) + Xc <- scale(X, TRUE, FALSE) + + P_F <- fy%*%solve(t(fy)%*%fy)%*%t(fy) + + Sigmahat <- cov(X) + Sigmahat_fit <- cov(P_F%*%X) + Sigmahat_res <- Sigmahat - Sigmahat_fit + + if (structure=="iso") + { + iso <- function(i) { - iso <- function(i) - { - ev <- eigen(Sigmahat) - ev.fit <- eigen(Sigmahat_fit) - all_evalues <-ev.fit$values - evalues <- all_evalues[1:i] - sigma2hat <- Re(sum(ev$values)/p) - - Gammahat <- Re(matrix(ev.fit$vectors[,1:i], ncol=i)) - dimnames(Gammahat) <- list(vnames, paste("Dir", 1:i, sep="")) - Betahat <-Re(t(Gammahat)%*%t(Xc)%*%fy%*%solve(t(fy)%*%fy)) - sigma2hat <- Re((sum(ev$values)-sum(evalues))/p) - Deltahat <- sigma2hat*diag(1, p) - dimnames(Deltahat) <- list(vnames, vnames) - - loglik <- - 0.5*n*p*(1+log(2*pi*sigma2hat)) - numpar <- p + (p-i)*i + i*dim(fy)[2] + 1 - aic <- -2*loglik + 2*numpar - bic <- -2*loglik + log(n)*numpar - - return(list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, - evalues=evalues, loglik=loglik, aic=aic, bic=bic, numpar=numpar)) - } + ev <- eigen(Sigmahat) + ev.fit <- eigen(Sigmahat_fit) + all_evalues <-ev.fit$values + evalues <- all_evalues[1:i] + sigma2hat <- Re(sum(ev$values)/p) - if (identical(numdir.test, FALSE)) - { - out <- iso(eff.numdir) - - ans <- list(R=X%*%orthonorm(out$Gammahat), Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, - Deltahat=out$Deltahat, loglik=out$loglik, aic=out$aic, bic=out$bic, numpar=out$numpar, - numdir=eff.numdir, evalues=out$evalues, structure="iso", y=y, fy=fy, - Xc=Xc, call=match.call(expand.dots=TRUE), numdir.test=numdir.test) - - class(ans) <- "pfc" - return(ans) - } + Gammahat <- Re(matrix(ev.fit$vectors[,1:i], ncol=i)) + dimnames(Gammahat) <- list(vnames, paste("Dir", 1:i, sep="")) + Betahat <-Re(t(Gammahat)%*%t(Xc)%*%fy%*%solve(t(fy)%*%fy)) + sigma2hat <- Re((sum(ev$values)-sum(evalues))/p) + Deltahat <- sigma2hat*diag(1, p) + dimnames(Deltahat) <- list(vnames, vnames) - if (identical(numdir.test, TRUE)) - { - aic <- bic <- numpar <- loglik <- vector(length=eff.numdir+1) - Betahat <- Deltahat <- Gammahat <-vector("list") - - # No fitting values (eff.numdir=0) - ev <- eigen(Sigmahat) - sigma2hat <- sum(ev$values)/p - loglik[1] <- - 0.5*n*p*(1+log(2*pi*sigma2hat)) - numpar[1] <- p + 1 - aic[1] <- -2*loglik[1] + 2*numpar[1] - bic[1] <- -2*loglik[1] + log(n)*numpar[1] - - for (i in 1:eff.numdir) - { - fit <- iso(i) - Betahat[[i]] <-fit$Betahat - Gammahat[[i]] <-fit$Gammahat - Deltahat[[i]] <- fit$Deltahat - loglik[i+1] <- fit$loglik - numpar[i+1] <- fit$numpar - aic[i+1] <- fit$aic - bic[i+1] <- fit$bic - } - ans <- list(R=X%*%orthonorm(Gammahat[[eff.numdir]]), Muhat=Muhat, Betahat=Betahat, Gammahat=Gammahat, - Deltahat=Deltahat, loglik=loglik, aic=aic, bic=bic, numpar=numpar, - numdir=eff.numdir, model="pfc", evalues=fit$evalues, structure="iso", - y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) - - class(ans)<- "pfc" - return(ans) - } + loglik <- - 0.5*n*p*(1+log(2*pi*sigma2hat)) + numpar <- p + (p-i)*i + i*dim(fy)[2] + 1 + aic <- -2*loglik + 2*numpar + bic <- -2*loglik + log(n)*numpar + + return(list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, + evalues=evalues, loglik=loglik, aic=aic, bic=bic, numpar=numpar)) } - if (structure=="aniso") + if (identical(numdir.test, FALSE)) { - aniso = function(X, y, fy, d, eps_aniso=1e-3, numdir.test) - { - vnames <- dimnames(X)[[2]] - if (is.null(vnames)) vnames <- paste("X", seq_len(ncol(X)), sep="") - op <- dim(X) - n <- op[1] - p <- op[2] - - # Initial Step - fit <- pfc(X=X, y=y, fy=fy, numdir=d, structure="iso", numdir.test=numdir.test) - - if (identical(numdir.test, FALSE)) - { - Betahatx <- fit$Betahat - Gammahatx <- fit$Gammahat - Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) - deltahat <- diag(cov(Xc)) - - repeat - { - Xnew = X%*%((1/sqrt(deltahat))*diag(p)) - fit <- pfc(X=Xnew, y=y, fy=fy, numdir=d, structure="iso", numdir.test=FALSE) - Betahatx <- fit$Betahat - Gammahatx <- (diag(p)*sqrt(deltahat))%*%fit$Gammahat - Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) - deltahat0 <- diag(t(Xc)%*%(Xc)/n) - if (sum(abs(deltahat-deltahat0)) < eps_aniso) break - deltahat <- deltahat0 - } - dimnames(Gammahatx) <- list(vnames, paste("Dir", 1:d, sep="")) - Deltahat <- deltahat*diag(p) - dimnames(Deltahat) <- list(vnames, vnames) - - loglik <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(deltahat)) - numpar <- p + d*(p-d) + ncol(fy)*d + p - aic <- -2*loglik + 2*numpar - bic <- -2*loglik + log(n)*numpar - - ans <- list(Betahat=Betahatx, Gammahat=orthonorm(Gammahatx), Deltahat=Deltahat, evalues=fit$evalues, - loglik=loglik, aic=aic, bic=bic, numpar=numpar, numdir.test=numdir.test) - - return(ans) - } - - Deltahat <- Betahat <- Gammahat <- vector("list") - aic <- bic <- numpar <- loglik <- vector(length=eff.numdir + 1) - - # No fitting values (eff.numdir=0) - ev <- eigen(Sigmahat) - loglik[1] <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(ev$values)) - numpar[1] <- p + p - aic[1] <- -2*loglik[1] + 2*numpar[1] - bic[1] <- -2*loglik[1] + log(n)*numpar[1] - - for (i in 1:eff.numdir) - { - Betahatx <- fit$Betahat[[i]] - Gammahatx <- fit$Gammahat[[i]] - Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) - deltahat <- diag(t(Xc)%*%(Xc)/n) - - repeat - { - Xnew = X%*%((1/sqrt(deltahat))*diag(p)) - fit2 <- pfc(X=Xnew, y=y, fy=fy, numdir=i, structure="iso", numdir.test=FALSE) - Betahatx <- fit2$Betahat - Gammahatx <- (diag(p)*sqrt(deltahat))%*%fit2$Gammahat - Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) - deltahat0 <- diag(t(Xc)%*%(Xc)/n) - if (sum(abs(deltahat-deltahat0)) < eps_aniso) break - deltahat <- deltahat0 - } - - Deltahat[[i]] <- deltahat*diag(p) - dimnames(Deltahat[[i]]) <- list(vnames, vnames) - - loglik[i+1] = - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(deltahat)) - numpar[i+1] <- p + (p-i)*i + i*dim(fy)[2] + p - aic[i+1] <- -2*loglik[i+1] + 2*numpar[i+1] - bic[i+1] <- -2*loglik[i+1] + log(n)*numpar[i+1] - Betahat[[i]] <- Betahatx - Gammahat[[i]] <- orthonorm(Gammahatx) - dimnames(Gammahat[[i]]) <- list(vnames, paste("Dir", 1:i, sep="")) - } - ans <- list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=fit2$evalues, - loglik=loglik, aic=aic, bic=bic, numpar=numpar, numdir.test=numdir.test) - - return(ans) - } - - fit <- aniso(X=X, y=y, fy=fy, d=eff.numdir, eps_aniso=eps_aniso, numdir.test=numdir.test) + out <- iso(eff.numdir) - ans <- list(Muhat=Muhat, Betahat=fit$Betahat, Gammahat=fit$Gammahat, Deltahat=fit$Deltahat, model="pfc", - loglik=fit$loglik, aic=fit$aic, bic=fit$bic, numpar=fit$numpar, numdir=eff.numdir, - evalues=fit$evalues, structure="aniso", Xc=Xc, y=y, fy=fy, call=match.call(), numdir.test=fit$numdir.test) - - if (numdir.test==FALSE) ans$R <- X%*%orthonorm(((fit$Deltahat)%^%(-1))%*%fit$Gammahat) else - ans$R <- X%*%orthonorm(((fit$Deltahat[[eff.numdir]])%^%(-1))%*%fit$Gammahat[[eff.numdir]]) - - class(ans)<- "pfc" + ans <- list(R=X%*%orthonorm(out$Gammahat), Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, + Deltahat=out$Deltahat, loglik=out$loglik, aic=out$aic, bic=out$bic, numpar=out$numpar, + numdir=eff.numdir, evalues=out$evalues, structure="iso", y=y, fy=fy, + Xc=Xc, call=match.call(expand.dots=TRUE), numdir.test=numdir.test) + class(ans) <- "pfc" return(ans) } - if (structure=="unstr") + if (identical(numdir.test, TRUE)) { - unstr<-function(i) - { - sqrt_Sigmahat_res <- Sigmahat_res%^%0.5 - Inv_Sqrt_Sigmahat_res <- solve(sqrt_Sigmahat_res) - lf_matrix <- Inv_Sqrt_Sigmahat_res%*%Sigmahat_fit%*%Inv_Sqrt_Sigmahat_res - all_evalues <- eigen(lf_matrix, symmetric=T)$values - evalues <- all_evalues[1:i] - - Vhat <- eigen(lf_matrix, symmetric=T)$vectors - Vhati <- matrix(Vhat[,1:i], ncol=i) - Gammahat <- (Sigmahat_res%^%0.5)%*%Vhati%*%solve((t(Vhati)%*%Sigmahat_res%*%Vhati)%^%0.5) - dimnames(Gammahat)<- list(vnames, paste("Dir", 1:i, sep="")) - - Khat<-diag(0, p) - if (i < min(ncol(fy),p)) {diag(Khat)[(i+1):min(ncol(fy), p )]<- all_evalues[(i+1):min(ncol(fy), p)]} - Deltahat <- sqrt_Sigmahat_res%*%Vhat%*%(diag(p)+Khat)%*%t(Vhat)%*%sqrt_Sigmahat_res - dimnames(Deltahat) <- list(vnames, vnames) - Betahat <- ((t(Vhati)%*%Sigmahat_res%*%Vhati)%^%0.5)%*%t(Vhati)%*%solve(Sigmahat_res%^%0.5)%*%t(Xc)%*%fy%*% solve(t(fy)%*%fy) - - temp0 <- -(n*p/2)*(1 + log(2*pi)) - temp1 <- -(n/2)*log(det(Sigmahat_res)) - temp2 <- 0 - - if (i < min(ncol(fy),p)) temp2 <- -(n/2)*sum(log(1 + all_evalues[(i+1):p])) - loglik <- temp0 + temp1 + temp2 - numpar <- p + (p-i)*i + i*ncol(fy) + p*(p+1)/2 - aic <- -2*loglik + 2*numpar - bic <- -2*loglik + log(n)*numpar - - return(list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=evalues, - loglik=loglik, aic=aic, bic=bic, numpar=numpar)) - } - - if (identical(numdir.test, FALSE)) - { - out <- unstr(eff.numdir) - - ans <- list(R=X%*%orthonorm(solve(out$Deltahat)%*%out$Gammahat), Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, - Deltahat=out$Deltahat, evalues=out$evalues, loglik=out$loglik, aic=out$aic, bic=out$bic, numpar=out$numpar, - numdir=eff.numdir, model="pfc", structure="unstr", y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) - - class(ans) <- "pfc" - return(ans) - } - aic <- bic <- numpar <- loglik <- vector(length=eff.numdir+1) - evalues <- vector(length=eff.numdir) - Betahat <- Deltahat <- Gammahat <-vector("list") + Betahat <- Deltahat <- Gammahat <-vector("list") - loglik[1] <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(det(Sigmahat)) - numpar[1] <- p + p*(p+1)/2 + # No fitting values (eff.numdir=0) + ev <- eigen(Sigmahat) + sigma2hat <- sum(ev$values)/p + loglik[1] <- - 0.5*n*p*(1+log(2*pi*sigma2hat)) + numpar[1] <- p + 1 aic[1] <- -2*loglik[1] + 2*numpar[1] bic[1] <- -2*loglik[1] + log(n)*numpar[1] - Deltahat[[1]] <- Sigmahat - dimnames(Deltahat[[1]]) <- list(vnames, vnames) for (i in 1:eff.numdir) { - fit <- unstr(i) - Betahat[[i]] <-fit$Betahat - Gammahat[[i]] <-fit$Gammahat + fit <- iso(i) + Betahat[[i]] <-fit$Betahat + Gammahat[[i]] <-fit$Gammahat Deltahat[[i]] <- fit$Deltahat - loglik[i+1] <- fit$loglik + loglik[i+1] <- fit$loglik numpar[i+1] <- fit$numpar aic[i+1] <- fit$aic - bic[i+1] <- fit$bic + bic[i+1] <- fit$bic } - ans <- list(R=X%*%orthonorm(solve(Deltahat[[eff.numdir]])%*%Gammahat[[eff.numdir]]), Muhat=Muhat, - Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=fit$evalues, loglik=loglik, - aic=aic, bic=bic, numpar=numpar, numdir=eff.numdir, model="pfc", structure="unstr", y=y, - fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) + ans <- list(R=X%*%orthonorm(Gammahat[[eff.numdir]]), Muhat=Muhat, Betahat=Betahat, Gammahat=Gammahat, + Deltahat=Deltahat, loglik=loglik, aic=aic, bic=bic, numpar=numpar, + numdir=eff.numdir, model="pfc", evalues=fit$evalues, structure="iso", + y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) class(ans)<- "pfc" return(ans) - } - else if (structure=="unstr2") + } + } + + if (structure=="aniso") + { + aniso = function(X, y, fy, d, eps_aniso=1e-3, numdir.test) { - unstr2 <- function(i) + vnames <- dimnames(X)[[2]] + if (is.null(vnames)) vnames <- paste("X", seq_len(ncol(X)), sep="") + op <- dim(X) + n <- op[1] + p <- op[2] + + # Initial Step + fit <- pfc(X=X, y=y, fy=fy, numdir=d, structure="iso", numdir.test=numdir.test) + + if (identical(numdir.test, FALSE)) { - objfun <- function(W) + Betahatx <- fit$Betahat + Gammahatx <- fit$Gammahat + Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) + deltahat <- diag(cov(Xc)) + + repeat { - Qt <- W$Qt - dc <- W$dim[1] - p <- ncol(Qt) - S <- W$Sigmas - U <- matrix(Qt[,1:dc], ncol=dc) - V <- matrix(Qt[,(dc+1):p], ncol=(p-dc)) - value <- -(n/2)*(p*log(2*pi)+p+log(det(t(V)%*%S$Sigmahat%*%V))+ log(det(t(U)%*%S$Sigmahat_res%*%U))) - - terme1 <- solve(t(U)%*%S$Sigmahat_res%*%U)%*%(t(U)%*%S$Sigmahat_res%*%V) - terme2 <- (t(U)%*%S$Sigmahat%*%V)%*%solve(t(V)%*%S$Sigmahat%*%V) - - gradient <- 2*(terme1 - terme2) - return(list(value=value, gradient=gradient)) + Xnew = X%*%((1/sqrt(deltahat))*diag(p)) + fit <- pfc(X=Xnew, y=y, fy=fy, numdir=d, structure="iso", numdir.test=FALSE) + Betahatx <- fit$Betahat + Gammahatx <- (diag(p)*sqrt(deltahat))%*%fit$Gammahat + Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) + deltahat0 <- diag(t(Xc)%*%(Xc)/n) + if (sum(abs(deltahat-deltahat0)) < eps_aniso) break + deltahat <- deltahat0 } - sigmas <- list(Sigmahat=Sigmahat, Sigmahat_fit=Sigmahat_fit, Sigmahat_res=Sigmahat_res, p=p, n=n) - - W <- list(Qt = svd(Sigmahat_fit)$u, dim=c(numdir, p), - Sigmas=list(Sigmahat=Sigmahat, Sigmahat_fit=Sigmahat_fit, - Sigmahat_res=Sigmahat_res, p=p, n=n)) - - objfun <- assign("objfun", objfun, envir=.BaseNamespaceEnv) - grassoptim <- GrassmannOptim(objfun, W,...) - - Gammahat <- matrix(grassoptim$Qt[,1:i], ncol=i, dimnames=list(vnames, paste("Dir", 1:i, sep=""))) - - Gammahat0 <- matrix(grassoptim$Qt[, (i+1):p], ncol=p-i, dimnames=list(vnames, paste("Dir", (i+1):p, sep=""))) - - Betahat <- t(Gammahat)%*%t(Xc)%*%fy%*%solve(t(fy)%*%fy) - Omegahat <- t(Gammahat)%*%Sigmahat_res%*%Gammahat - Omegahat0 <-t(Gammahat0)%*%Sigmahat%*%Gammahat0 - Deltahat <- Gammahat%*%Omegahat%*%t(Gammahat) + Gammahat0%*%Omegahat0%*%t(Gammahat0) + dimnames(Gammahatx) <- list(vnames, paste("Dir", 1:d, sep="")) + Deltahat <- deltahat*diag(p) dimnames(Deltahat) <- list(vnames, vnames) - temp0 <- -(n*p/2)*(1+log(2*pi)) - temp1 <- -(n/2)*log(det(t(Gammahat)%*%Sigmahat_res%*%Gammahat)) - temp2 <- -(n/2)*log(det(t(Gammahat0)%*%Sigmahat%*%Gammahat0)) - - loglik <- temp0 + temp1 + temp2 - numpar <- p + (p-i)*i + i*dim(fy)[2] + i*(i+1)/2 + (p-i)*(p-i+1)/2 + loglik <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(deltahat)) + numpar <- p + d*(p-d) + ncol(fy)*d + p aic <- -2*loglik + 2*numpar bic <- -2*loglik + log(n)*numpar - ev.fit <- eigen(Sigmahat_fit) - evalues <- ev.fit$values[1:i] - return(list(Betahat=Betahat, Gammahat=Gammahat, Gammahat0=Gammahat0, Omegahat=Omegahat, Omegahat0=Omegahat0, - Deltahat=Deltahat, evalues=evalues, loglik=loglik, aic=aic, bic=bic, numpar=numpar)) + ans <- list(Betahat=Betahatx, Gammahat=orthonorm(Gammahatx), Deltahat=Deltahat, evalues=fit$evalues, + loglik=loglik, aic=aic, bic=bic, numpar=numpar, numdir.test=numdir.test) + + return(ans) } - if (identical(numdir.test, FALSE)) + Deltahat <- Betahat <- Gammahat <- vector("list") + aic <- bic <- numpar <- loglik <- vector(length=eff.numdir + 1) + + # No fitting values (eff.numdir=0) + ev <- eigen(Sigmahat) + loglik[1] <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(ev$values)) + numpar[1] <- p + p + aic[1] <- -2*loglik[1] + 2*numpar[1] + bic[1] <- -2*loglik[1] + log(n)*numpar[1] + + for (i in 1:eff.numdir) { - out <- unstr2(numdir) + Betahatx <- fit$Betahat[[i]] + Gammahatx <- fit$Gammahat[[i]] + Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) + deltahat <- diag(t(Xc)%*%(Xc)/n) - ans <- list(R = X%*%out$Gammahat, Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, - Gammahat0=out$Gammahat0, Omegahat=out$Omegahat, Omegahat0=out$Omegahat0, - Deltahat=out$Deltahat, evalues=out$evalues, loglik=out$loglik, aic=out$aic, - bic=out$bic, numpar=out$numpar, numdir=numdir, model="pfc", structure="unstr2", - y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) + repeat + { + Xnew = X%*%((1/sqrt(deltahat))*diag(p)) + fit2 <- pfc(X=Xnew, y=y, fy=fy, numdir=i, structure="iso", numdir.test=FALSE) + Betahatx <- fit2$Betahat + Gammahatx <- (diag(p)*sqrt(deltahat))%*%fit2$Gammahat + Xc <- scale(X, TRUE, FALSE) - fy%*%t(Gammahatx%*%Betahatx) + deltahat0 <- diag(t(Xc)%*%(Xc)/n) + if (sum(abs(deltahat-deltahat0)) < eps_aniso) break + deltahat <- deltahat0 + } - class(ans) <- "pfc" - return(ans) + Deltahat[[i]] <- deltahat*diag(p) + dimnames(Deltahat[[i]]) <- list(vnames, vnames) + + loglik[i+1] = - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(prod(deltahat)) + numpar[i+1] <- p + (p-i)*i + i*dim(fy)[2] + p + aic[i+1] <- -2*loglik[i+1] + 2*numpar[i+1] + bic[i+1] <- -2*loglik[i+1] + log(n)*numpar[i+1] + Betahat[[i]] <- Betahatx + Gammahat[[i]] <- orthonorm(Gammahatx) + dimnames(Gammahat[[i]]) <- list(vnames, paste("Dir", 1:i, sep="")) } + ans <- list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=fit2$evalues, + loglik=loglik, aic=aic, bic=bic, numpar=numpar, numdir.test=numdir.test) - aic <- bic <- numpar <- loglik <- vector(length=numdir+1) - Betahat <- Deltahat <- Gammahat <- Gammahat0 <- Omegahat <- Omegahat0 <- vector("list") + return(ans) + } + + fit <- aniso(X=X, y=y, fy=fy, d=eff.numdir, eps_aniso=eps_aniso, numdir.test=numdir.test) + + ans <- list(Muhat=Muhat, Betahat=fit$Betahat, Gammahat=fit$Gammahat, Deltahat=fit$Deltahat, model="pfc", + loglik=fit$loglik, aic=fit$aic, bic=fit$bic, numpar=fit$numpar, numdir=eff.numdir, + evalues=fit$evalues, structure="aniso", Xc=Xc, y=y, fy=fy, call=match.call(), numdir.test=fit$numdir.test) + + if (numdir.test==FALSE) ans$R <- X%*%orthonorm(((fit$Deltahat)%^%(-1))%*%fit$Gammahat) else + ans$R <- X%*%orthonorm(((fit$Deltahat[[eff.numdir]])%^%(-1))%*%fit$Gammahat[[eff.numdir]]) + + class(ans)<- "pfc" + + return(ans) + } + + if (structure=="unstr") + { + unstr<-function(i) + { + sqrt_Sigmahat_res <- Sigmahat_res%^%0.5 + Inv_Sqrt_Sigmahat_res <- solve(sqrt_Sigmahat_res) + lf_matrix <- Inv_Sqrt_Sigmahat_res%*%Sigmahat_fit%*%Inv_Sqrt_Sigmahat_res + all_evalues <- eigen(lf_matrix, symmetric=T)$values + evalues <- all_evalues[1:i] - loglik[1] <- -(n*p/2)*(log(2*pi) + (1+log(Trace(Sigmahat)/p))) - numpar[1] <- p + p*(p+1)/2 - aic[1] <- -2*loglik[1] + 2*numpar[1] - bic[1] <- -2*loglik[1] + log(n)*numpar[1] + Vhat <- eigen(lf_matrix, symmetric=T)$vectors + Vhati <- matrix(Vhat[,1:i], ncol=i) + Gammahat <- (Sigmahat_res%^%0.5)%*%Vhati%*%solve((t(Vhati)%*%Sigmahat_res%*%Vhati)%^%0.5) + dimnames(Gammahat)<- list(vnames, paste("Dir", 1:i, sep="")) + + Khat<-diag(0, p) + if (i < min(ncol(fy),p)) {diag(Khat)[(i+1):min(ncol(fy), p )]<- all_evalues[(i+1):min(ncol(fy), p)]} + Deltahat <- sqrt_Sigmahat_res%*%Vhat%*%(diag(p)+Khat)%*%t(Vhat)%*%sqrt_Sigmahat_res + dimnames(Deltahat) <- list(vnames, vnames) + Betahat <- ((t(Vhati)%*%Sigmahat_res%*%Vhati)%^%0.5)%*%t(Vhati)%*%solve(Sigmahat_res%^%0.5)%*%t(Xc)%*%fy%*% solve(t(fy)%*%fy) - for(m in 1:numdir) + temp0 <- -(n*p/2)*(1 + log(2*pi)) + temp1 <- -(n/2)*log(det(Sigmahat_res)) + temp2 <- 0 + + if (i < min(ncol(fy),p)) temp2 <- -(n/2)*sum(log(1 + all_evalues[(i+1):p])) + loglik <- temp0 + temp1 + temp2 + numpar <- p + (p-i)*i + i*ncol(fy) + p*(p+1)/2 + aic <- -2*loglik + 2*numpar + bic <- -2*loglik + log(n)*numpar + + return(list(Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=evalues, + loglik=loglik, aic=aic, bic=bic, numpar=numpar)) + } + + if (identical(numdir.test, FALSE)) + { + out <- unstr(eff.numdir) + + ans <- list(R=X%*%orthonorm(solve(out$Deltahat)%*%out$Gammahat), Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, + Deltahat=out$Deltahat, evalues=out$evalues, loglik=out$loglik, aic=out$aic, bic=out$bic, numpar=out$numpar, + numdir=eff.numdir, model="pfc", structure="unstr", y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) + + class(ans) <- "pfc" + return(ans) + } + + aic <- bic <- numpar <- loglik <- vector(length=eff.numdir+1) + evalues <- vector(length=eff.numdir) + Betahat <- Deltahat <- Gammahat <-vector("list") + + loglik[1] <- - 0.5*n*p*(1+log(2*pi)) - 0.5*n*log(det(Sigmahat)) + numpar[1] <- p + p*(p+1)/2 + aic[1] <- -2*loglik[1] + 2*numpar[1] + bic[1] <- -2*loglik[1] + log(n)*numpar[1] + Deltahat[[1]] <- Sigmahat + dimnames(Deltahat[[1]]) <- list(vnames, vnames) + + for (i in 1:eff.numdir) + { + fit <- unstr(i) + Betahat[[i]] <-fit$Betahat + Gammahat[[i]] <-fit$Gammahat + Deltahat[[i]] <- fit$Deltahat + loglik[i+1] <- fit$loglik + numpar[i+1] <- fit$numpar + aic[i+1] <- fit$aic + bic[i+1] <- fit$bic + } + ans <- list(R=X%*%orthonorm(solve(Deltahat[[eff.numdir]])%*%Gammahat[[eff.numdir]]), Muhat=Muhat, + Betahat=Betahat, Gammahat=Gammahat, Deltahat=Deltahat, evalues=fit$evalues, loglik=loglik, + aic=aic, bic=bic, numpar=numpar, numdir=eff.numdir, model="pfc", structure="unstr", y=y, + fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) + + class(ans)<- "pfc" + return(ans) + } + else if (structure=="unstr2") + { + unstr2 <- function(i) + { + objfun <- function(W) { - fit <- unstr2(m) - Betahat[[m]] <-fit$Betahat - Gammahat[[m]] <-fit$Gammahat - Omegahat[[m]] <- fit$Omegahat - Omegahat0[[m]] <- fit$Omegahat0 - Deltahat[[m]] <- fit$Deltahat - loglik[m+1] <- fit$loglik - numpar[m+1] <- fit$numpar - aic[m+1] <- fit$aic - bic[m+1] <- fit$bic + Qt <- W$Qt + dc <- W$dim[1] + p <- ncol(Qt) + S <- W$Sigmas + U <- matrix(Qt[,1:dc], ncol=dc) + V <- matrix(Qt[,(dc+1):p], ncol=(p-dc)) + value <- -(n/2)*(p*log(2*pi)+p+log(det(t(V)%*%S$Sigmahat%*%V))+ log(det(t(U)%*%S$Sigmahat_res%*%U))) + + terme1 <- solve(t(U)%*%S$Sigmahat_res%*%U)%*%(t(U)%*%S$Sigmahat_res%*%V) + terme2 <- (t(U)%*%S$Sigmahat%*%V)%*%solve(t(V)%*%S$Sigmahat%*%V) + + gradient <- 2*(terme1 - terme2) + return(list(value=value, gradient=gradient)) } - ans <- list(R = X%*%Gammahat[[numdir]], evalues=fit$evalues, loglik =loglik, aic=aic, bic=bic, numdir=numdir, - numpar=numpar, Muhat=Muhat, Betahat=Betahat, Gammahat=Gammahat, - Gammahat0=Gammahat0, Omegahat=Omegahat, Omegahat0=Omegahat0, - Deltahat=Deltahat, model="pfc", structure="unstr2", + sigmas <- list(Sigmahat=Sigmahat, Sigmahat_fit=Sigmahat_fit, Sigmahat_res=Sigmahat_res, p=p, n=n) + + W <- list(Qt = svd(Sigmahat_fit)$u, dim=c(numdir, p), + Sigmas=list(Sigmahat=Sigmahat, Sigmahat_fit=Sigmahat_fit, + Sigmahat_res=Sigmahat_res, p=p, n=n)) + + objfun <- assign("objfun", objfun, envir=.BaseNamespaceEnv) + grassoptim <- GrassmannOptim(objfun, W,...) + + Gammahat <- matrix(grassoptim$Qt[,1:i], ncol=i, dimnames=list(vnames, paste("Dir", 1:i, sep=""))) + + Gammahat0 <- matrix(grassoptim$Qt[, (i+1):p], ncol=p-i, dimnames=list(vnames, paste("Dir", (i+1):p, sep=""))) + + Betahat <- t(Gammahat)%*%t(Xc)%*%fy%*%solve(t(fy)%*%fy) + Omegahat <- t(Gammahat)%*%Sigmahat_res%*%Gammahat + Omegahat0 <-t(Gammahat0)%*%Sigmahat%*%Gammahat0 + Deltahat <- Gammahat%*%Omegahat%*%t(Gammahat) + Gammahat0%*%Omegahat0%*%t(Gammahat0) + dimnames(Deltahat) <- list(vnames, vnames) + + temp0 <- -(n*p/2)*(1+log(2*pi)) + temp1 <- -(n/2)*log(det(t(Gammahat)%*%Sigmahat_res%*%Gammahat)) + temp2 <- -(n/2)*log(det(t(Gammahat0)%*%Sigmahat%*%Gammahat0)) + + loglik <- temp0 + temp1 + temp2 + numpar <- p + (p-i)*i + i*dim(fy)[2] + i*(i+1)/2 + (p-i)*(p-i+1)/2 + aic <- -2*loglik + 2*numpar + bic <- -2*loglik + log(n)*numpar + ev.fit <- eigen(Sigmahat_fit) + evalues <- ev.fit$values[1:i] + + return(list(Betahat=Betahat, Gammahat=Gammahat, Gammahat0=Gammahat0, Omegahat=Omegahat, Omegahat0=Omegahat0, + Deltahat=Deltahat, evalues=evalues, loglik=loglik, aic=aic, bic=bic, numpar=numpar)) + } + + if (identical(numdir.test, FALSE)) + { + out <- unstr2(numdir) + + ans <- list(R = X%*%out$Gammahat, Muhat=Muhat, Betahat=out$Betahat, Gammahat=out$Gammahat, + Gammahat0=out$Gammahat0, Omegahat=out$Omegahat, Omegahat0=out$Omegahat0, + Deltahat=out$Deltahat, evalues=out$evalues, loglik=out$loglik, aic=out$aic, + bic=out$bic, numpar=out$numpar, numdir=numdir, model="pfc", structure="unstr2", y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) - class(ans)<- "pfc" + class(ans) <- "pfc" return(ans) } + aic <- bic <- numpar <- loglik <- vector(length=numdir+1) + Betahat <- Deltahat <- Gammahat <- Gammahat0 <- Omegahat <- Omegahat0 <- vector("list") + + loglik[1] <- -(n*p/2)*(log(2*pi) + (1+log(Trace(Sigmahat)/p))) + numpar[1] <- p + p*(p+1)/2 + aic[1] <- -2*loglik[1] + 2*numpar[1] + bic[1] <- -2*loglik[1] + log(n)*numpar[1] + + for(m in 1:numdir) + { + fit <- unstr2(m) + Betahat[[m]] <-fit$Betahat + Gammahat[[m]] <-fit$Gammahat + Omegahat[[m]] <- fit$Omegahat + Omegahat0[[m]] <- fit$Omegahat0 + Deltahat[[m]] <- fit$Deltahat + loglik[m+1] <- fit$loglik + numpar[m+1] <- fit$numpar + aic[m+1] <- fit$aic + bic[m+1] <- fit$bic + } + ans <- list(R = X%*%Gammahat[[numdir]], evalues=fit$evalues, loglik =loglik, aic=aic, bic=bic, numdir=numdir, + numpar=numpar, Muhat=Muhat, Betahat=Betahat, Gammahat=Gammahat, + Gammahat0=Gammahat0, Omegahat=Omegahat, Omegahat0=Omegahat0, + Deltahat=Deltahat, model="pfc", structure="unstr2", + y=y, fy=fy, Xc=Xc, call=match.call(), numdir.test=numdir.test) + + class(ans)<- "pfc" + return(ans) } + +} -#' @importFrom stats lm +#' @importFrom stats lm as.formula #' onepfc <- function(X, y, fy, p, numdir.test) { # X is univariate predictor nobs <- length(X) r <- dim(fy)[2] - P_F <- fy%*%solve(t(fy)%*%fy)%*%t(fy) + P_F <- fy %*% solve(t(fy) %*% fy) %*% t(fy) Xc <- scale(X, TRUE, FALSE) - Sigmahat_fit <- (1/nobs)*t(Xc)%*%P_F%*%(Xc) + Sigmahat_fit <- (1/nobs)*t(Xc) %*% P_F %*% (Xc) ev.fit <- eigen(Sigmahat_fit) - temp.dat<-data.frame(cbind(X, fy)) - xnam<-paste("xx", 1:r, sep="") - names(temp.dat)<-c("yy", xnam) - fm.lm<- as.formula( paste("yy ~ ", paste(xnam, collapse= "+"))) + temp.dat <- data.frame(cbind(X, fy)) + xnam <- paste("xx", 1:r, sep="") + names(temp.dat) <- c("yy", xnam) + fm.lm <- as.formula( paste("yy ~ ", paste(xnam, collapse= "+"))) summary.fm <- summary(lm(fm.lm, data=temp.dat)) Betahat <- matrix(summary.fm$coefficients[2:(r+1),1], ncol=r) diff --git a/R/plot.bcea.R b/R/plot.bcea.R index e7dc27fc..6cf9b23d 100644 --- a/R/plot.bcea.R +++ b/R/plot.bcea.R @@ -35,6 +35,7 @@ #' \code{\link{ceac.plot}}, #' \code{\link{evi.plot}} #' @importFrom Rdpack reprompt +#' @importFrom purrr map_lgl #' #' @references #' diff --git a/R/plot.evppi.R b/R/plot.evppi.R index 9f70a3dd..053cea43 100644 --- a/R/plot.evppi.R +++ b/R/plot.evppi.R @@ -37,23 +37,25 @@ #' inp <- createInputs(vaccine_mat) #' #' # Compute the EVPPI using INLA/SPDE -#' x0 <- evppi(m, c("beta.1." , "beta.2."), input = inp$mat) +#' if (require("INLA")) { +#' x0 <- evppi(m, c("beta.1." , "beta.2."), input = inp$mat) +#' +#' plot(x0, pos = c(0,1)) #' -#' plot(x0, pos = c(0,1)) +#' x1 <- evppi(m, c(32,48,49), input = inp$mat) +#' plot(x1, pos = "topright") #' -#' x1 <- evppi(m, c(32,48,49), input = inp$mat) -#' plot(x1, pos = "topright") +#' plot(x0, col = c("black", "red"), pos = "topright") +#' plot(x0, col = c(2,3), pos = "bottomright") #' -#' plot(x0, col = c("black", "red"), pos = "topright") -#' plot(x0, col = c(2,3), pos = "bottomright") +#' plot(x0, pos = c(0,1), graph = "ggplot2") +#' plot(x1, pos = "top", graph = "ggplot2") #' -#' plot(x0, pos = c(0,1), graph = "ggplot2") -#' plot(x1, pos = "top", graph = "ggplot2") +#' plot(x0, col = c("black", "red"), pos = "right", graph = "ggplot2") +#' plot(x0, col = c(2,3), size = c(1,2), pos = "bottom", graph = "ggplot2") #' -#' plot(x0, col = c("black", "red"), pos = "right", graph = "ggplot2") -#' plot(x0, col = c(2,3), size = c(1,2), pos = "bottom", graph = "ggplot2") -#' -#' plot(x0, graph = "ggplot2", theme = ggplot2::theme_linedraw()) +#' plot(x0, graph = "ggplot2", theme = ggplot2::theme_linedraw()) +#' } #' #' if (FALSE) #' plot(x0, col = 3, pos = "topright") diff --git a/R/select_plot_type.R b/R/select_plot_type.R index 6951ad52..9ba48823 100644 --- a/R/select_plot_type.R +++ b/R/select_plot_type.R @@ -6,6 +6,7 @@ #' @param graph Type names; string #' @return Plot ID integer 1:base R; 2:ggplot2; 3:plotly #' @importFrom cli cli_alert_warning +#' @importFrom purrr map_lgl #' @keywords dplot internal #' select_plot_type <- function(graph) { @@ -15,7 +16,8 @@ select_plot_type <- function(graph) { graph_lup <- c(base = 1, ggplot2 = 2, plotly =3) graph_type <- graph_lup[graph] - is_req_pkgs <- map_lgl(c("ggplot2", "grid"), requireNamespace, quietly = TRUE) + is_req_pkgs <- purrr::map_lgl(c("ggplot2", "grid"), + requireNamespace, quietly = TRUE) if (graph_type == 2 && !all(is_req_pkgs)) { cli::cli_alert_warning( diff --git a/R/setComparisons.R b/R/setComparisons.R index 7f5812c3..aaa270d0 100644 --- a/R/setComparisons.R +++ b/R/setComparisons.R @@ -45,7 +45,12 @@ setComparisons <- function(he, comparison) { res$eib <- res$eib[, name_comp, drop = FALSE] res$ceac <- res$ceac[, name_comp, drop = FALSE] + ##TODO: is there a way not to recompute the whole thing? + res$best <- best_interv_given_k(res$eib, res$ref, res$comp) + res$kstar <- compute_kstar(res$k, res$best, res$ref) + ##TODO: currently compute _all_ interventions in compute_U() + ## change to this? # res$U <- res$U[, , name_comp, drop = FALSE] return(res) diff --git a/R/struct.psa.R b/R/struct.psa.R index 045941d2..a9029d83 100644 --- a/R/struct.psa.R +++ b/R/struct.psa.R @@ -115,12 +115,12 @@ struct.psa <- function(models, } dmin <- min(d) # minimum value to re-scale DICs - delta_dic <- abs(dmin-d) + delta_dic <- abs(dmin - d) # Only compute w using DIC if the user hasn't passed a suitable value - if(is.null(w)) { + if (is.null(w)) { w <- exp(-0.5*(delta_dic))/sum(exp(-0.5*(delta_dic))) # model weights (cfr BMHE) } - if ((!is.null(w) & length(w)!=n_models)) { + if ((!is.null(w) & length(w) != n_models)) { stop("If you are considering user-defined weights, you must pass a vector whose length is the same as the number of models to average!") } diff --git a/R/summary.bcea.R b/R/summary.bcea.R index 3df3e049..3eabb5e7 100644 --- a/R/summary.bcea.R +++ b/R/summary.bcea.R @@ -133,6 +133,7 @@ summary.bcea <- function(object, # cat(paste0(" : ", green(he$interventions[he$comp[i]]), "\n")) } } + cat("\n") if (length(he$kstar) == 0 && !is.na(he$step)) { cat( @@ -144,8 +145,12 @@ summary.bcea <- function(object, max(he$k), "] \n")) } + if (length(he$kstar) == 1 && !is.na(he$step)) { - kstar <- he$k[which(diff(he$best) == 1) + 1] + + ##TODO: why recalc when same as he$kstar? + kstar <- he$k[which(diff(he$best) != 0) + 1] + cat( paste0( "Optimal decision: choose ", @@ -158,6 +163,7 @@ summary.bcea <- function(object, kstar, "\n")) } + if (length(he$kstar) > 1 && !is.na(he$step)) { cat( paste0( @@ -183,6 +189,7 @@ summary.bcea <- function(object, he$kstar[length(he$kstar)], "\n")) } + cat("\n\n") cat(paste0("Analysis for willingness to pay parameter k = ", wtp, "\n")) cat("\n") diff --git a/R/zzz.R b/R/zzz.R index 49321e47..a7a0e5fd 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -10,6 +10,16 @@ ps.options(encoding = "CP1250") pdf.options(encoding = "CP1250") + Sys.setenv("_R_CHECK_LENGTH_1_CONDITION_" = "TRUE") + invisible() } +#' @title .onAttach +#' @description prints out a friendly reminder message to the user +#' @inheritParams base .onAttach +#' @return NULL +#' @noRd +.onAttach <- function(libname, pkgname) { + packageStartupMessage("The BCEA version loaded is: ", utils::packageVersion("BCEA")) +} diff --git a/README.md b/README.md index 3cc80b78..4e934bd2 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ > Perform Bayesian Cost-Effectiveness Analysis in R. -:rocket: **Version 2.4.2 out now!** [Check out the release notes here](https://github.com/n8thangreen/BCEA/releases/tag/v2.4.2). +:rocket: **Version 2.4.5 out now!** [Check out the release notes here](https://github.com/n8thangreen/BCEA/releases). ## Contents @@ -34,12 +34,12 @@ Main features of `BCEA` include: * EVPPI calculations and plots ## Installation -Install the released version from CRAN with +Install the [released version from CRAN](https://cran.r-project.org/package=BCEA) with ```r install.packages("BCEA") ``` -The development version can be installed using this GitHub repository. On Windows machines, you need to install a few dependencies, including [Rtools](https://cran.r-project.org/bin/windows/Rtools/) first, e.g. by running +The stable version (which can be updated more quickly) can be installed using this GitHub repository. On Windows machines, you need to install a few dependencies, including [Rtools](https://cran.r-project.org/bin/windows/Rtools/) first, e.g. by running ```r pkgs <- c("MASS", "Rtools", "remotes") @@ -49,13 +49,13 @@ install.packages(pkgs, repos=repos, dependencies = "Depends") before installing the package using `remotes`: ```r -remotes::install_github("giabaio/BCEA", ref="dev") +remotes::install_github("giabaio/BCEA") ``` Under Linux or MacOS, it is sufficient to install the package via `remotes`: ```r install.packages("remotes") -remotes::install_github("giabaio/BCEA", ref="dev") +remotes::install_github("giabaio/BCEA") ``` ## Articles @@ -72,14 +72,13 @@ Examples of using specific functions and their different arguments are given in ## Further details The `pkgdown` site is [here](https://n8thangreen.github.io/BCEA/). -More details on `BCEA` are available in our book [_Bayesian Cost-Effectiveness Analysis with the R Package BCEA_](https://gianluca.statistica.it/book/bcea/) (published in the UseR! Springer series). Also, details about the package, including some references and links to a pdf presentation and some posts on my own blog) are given [here](https://gianluca.statistica.it/software/bcea/). +More details on `BCEA` are available in our book [_Bayesian Cost-Effectiveness Analysis with the R Package BCEA_](https://gianluca.statistica.it/books/bcea/) (published in the UseR! Springer series). Also, details about the package, including some references and links to a pdf presentation and some posts on my own blog) are given [here](https://gianluca.statistica.it/software/bcea/). ## License [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) ## Contributing -Please submit contributions through `Pull Requests`, following the [contributing -guidelines](https://github.com/n8thangreen/BCEA/blob/dev/CONTRIBUTING.md). +Please submit contributions through `Pull Requests`, following the [contributing guidelines](https://github.com/n8thangreen/BCEA/blob/dev/CONTRIBUTING.md). To report issues and/or seek support, please file a new ticket in the [issue](https://github.com/n8thangreen/BCEA/issues) tracker. diff --git a/README_dev.md b/README_dev.md deleted file mode 100644 index 62facb3f..00000000 --- a/README_dev.md +++ /dev/null @@ -1,60 +0,0 @@ -# BCEA --- development version - - - -[![Build Status](https://app.travis-ci.com/n8thangreen/BCEA.svg?branch=dev)](https://app.travis-ci.com/n8thangreen/BCEA) -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/BCEA)](https://cran.r-project.org/package=BCEA) [![CRAN_Download_Badge](http://cranlogs.r-pkg.org/badges/BCEA)](https://cran.r-project.org/package=BCEA) -[![CRAN_Download_Badge](http://cranlogs.r-pkg.org/badges/grand-total/BCEA?color=orange)](https://cran.r-project.org/package=BCEA) -[![CodeFactor](https://www.codefactor.io/repository/github/n8thangreen/bcea/badge)](https://www.codefactor.io/repository/github/n8thangreen/bcea) - - -## Contents - -- [Overview](#introduction) -- [Features](#features) -- [Installation](#installation) -- [Further details](#further-details) - -## Overview - -Perform Bayesian Cost-Effectiveness Analysis in R. -Given the results of a Bayesian model (possibly based on MCMC) in the form of simulations from the posterior distributions of suitable variables of costs and clinical benefits for two or more interventions, produces a health economic evaluation. Compares one of the interventions (the "reference") to the others ("comparators"). - -## Features - -Main features of `BCEA` include: - -* Summary statistics and tables -* Cost-effectiveness analysis plots, such as CE planes and CEAC -* EVPPI calculations and plots - -This is the **development** version of BCEA (currently 2.4). It contains a major refactoring of the code to streamline the functions. - -## Installation -The development version can be installed using this GitHub repository. On Windows machines, you need to install a few dependencies, including [Rtools](https://cran.r-project.org/bin/windows/Rtools/) first, e.g. by running - -```r -pkgs <- c("MASS", "Rtools", "remotes") -repos <- c("https://cran.rstudio.com", "https://inla.r-inla-download.org/R/stable/") -install.packages(pkgs, repos=repos, dependencies = "Depends") -``` -before installing the package using `remotes`: - -```r -remotes::install_github("giabaio/BCEA", ref="dev") -``` -Under Linux or MacOS, it is sufficient to install the package via `remotes`: - -```r -install.packages("remotes") -remotes::install_github("giabaio/BCEA", ref="dev") -``` - -## Further details -The `pkgdown` site is [here](https://n8thangreen.github.io/BCEA/). -More details on `BCEA` are available in our book [_Bayesian Cost-Effectiveness Analysis with the R Package BCEA_](https://gianluca.statistica.it/book/bcea/) (published in the UseR! Springer series). Also, details about the package, including some references and links to a pdf presentation and some posts on my own blog) are given [here](https://gianluca.statistica.it/software/bcea/). - -## Licence -GPL-3 © [G Baio](https://github.com/giabaio/). - -Please note that this project is released with a [Contributor Code of Conduct](https://github.com/n8thangreen/BCEA/blob/dev/CONDUCT.md). By participating in this project you agree to abide by its terms. diff --git a/_pkgdown.yml b/_pkgdown.yml index e03de6e6..891e1399 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -9,6 +9,15 @@ template: bslib: font_scale: 1.0 +navbar: + structure: + left: [intro, reference, articles, tutorials, news, cheatsheet] + right: [search, github] + components: + cheatsheet: + text: "Cheatsheet" + href: "bcea_cheatsheet.pdf" + reference: - title: "Main" contents: diff --git a/cran-comments.md b/cran-comments.md index 01fa3693..1e88fbc5 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,19 +1,30 @@ -## Resubmission -This is a resubmission. In this version I have: -* Removed from the build large data files so that not more than 5 MB for a CRAN package. +# Resubmission (2.4.4) -* Reference about the method added in the Description with doi: prefix +* Suggested package wasn't used conditionally in unit test causing CRAN check error. `MCMCvis` is now moved to Required packages. +* Fixed R package CRAN webpage canonical form in README ## R CMD check results -There were no ERRORs or WARNINGs. + +0 errors | 0 warnings | 1 note * Any notes about using INLA has been condoned in previous versions, as they are only suggested. + +CHECK message: + + Suggests or Enhances not in mainstream repositories: + INLA + Availability using Additional_repositories specification: + INLA yes https://inla.r-inla-download.org/R/stable/ +Flavor: r-devel-linux-x86_64-debian-gcc, r-devel-windows-x86_64 +Check: package dependencies, Result: NOTE + Package suggested but not available for checking: 'INLA' ## Downstream dependencies -We checked 2 reverse dependencies, comparing R CMD check results across CRAN and dev versions of this package. +We checked 2 reverse dependencies (missingHE, heesim), comparing R CMD check results across CRAN and dev versions of this package. * We saw 0 new problems * We failed to check 0 packages + diff --git a/data/Smoking.RData b/data/Smoking.RData index c855d544..20dfc80f 100644 Binary files a/data/Smoking.RData and b/data/Smoking.RData differ diff --git a/data/Vaccine.RData b/data/Vaccine.RData index 86ce7677..c156a13b 100644 Binary files a/data/Vaccine.RData and b/data/Vaccine.RData differ diff --git a/data/datalist b/data/datalist index 31677fa0..c1af3b86 100644 --- a/data/datalist +++ b/data/datalist @@ -1,2 +1,2 @@ -Smoking: cost data eff life.years pi_post smoking treats -Vaccine: cost cost.time.off cost.trt2 eff N.resources QALYs.inf vaccine_mat cost.GP cost.time.vac cost.vac e.pts QALYs.adv QALYs.pne cost.hosp cost.travel c.pts N QALYs.death treats cost.otc cost.trt1 N.outcomes QALYs.hosp \ No newline at end of file +Smoking: cost data eff life.years pi_post smoking smoking_output treats +Vaccine: c.pts cost cost.GP cost.hosp cost.otc cost.time.off cost.time.vac cost.travel cost.trt1 cost.trt2 cost.vac e.pts eff N N.outcomes N.resources QALYs.adv QALYs.death QALYs.hosp QALYs.inf QALYs.pne treats vaccine_mat diff --git a/docs/404.html b/docs/404.html index f8bd2dba..77b316d0 100644 --- a/docs/404.html +++ b/docs/404.html @@ -6,33 +6,33 @@ Page not found (404) • BCEA - - - - - - - - - + + + + + + + + + - + - + - Skip to contents + Skip to contents
@@ -98,7 +101,7 @@ diff --git a/docs/CONDUCT.html b/docs/CONDUCT.html index 19005542..6d8e6cec 100644 --- a/docs/CONDUCT.html +++ b/docs/CONDUCT.html @@ -1,5 +1,5 @@ -Contributor Code of Conduct • BCEAContributor Code of Conduct • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.5
diff --git a/docs/CONTRIBUTING.html b/docs/CONTRIBUTING.html index 34ad2f1e..de8b1154 100644 --- a/docs/CONTRIBUTING.html +++ b/docs/CONTRIBUTING.html @@ -1,5 +1,5 @@ -Contributing to BCEA • BCEAContributing to BCEA • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.5 + + + + + +
+
+
+ + +

Please briefly describe your problem and what output you expect. If you have a general question, please don’t use this form. Instead, ask on https://stackoverflow.com/ or https://community.rstudio.com/.

+

Please include a minimal reproducible example (AKA a reprex). If you’ve never heard of a reprex before, start by reading https://www.tidyverse.org/help/#reprex.

+

Brief description of the problem

+
+# insert reprex here
+ + +
+ + +
+ + + + + + + diff --git a/docs/LICENSE.html b/docs/LICENSE.html index 4cbf2503..a5c26bde 100644 --- a/docs/LICENSE.html +++ b/docs/LICENSE.html @@ -1,5 +1,5 @@ -GNU General Public License • BCEAGNU General Public License • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.5 + + + + + +
+
+
+ +
+

Compute EVPPI

+
+ +
+

Usage

+
compute_evppi(he, fit.full)
+
+ +
+

Arguments

+
he
+

A bcea object containing the results of the Bayesian +modelling and the economic evaluation.

+ + +
fit.full
+

fit.full

+ +
+
+

Value

+ + +

list

+
+
+

See also

+ +
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/compute_kstar.html b/docs/reference/compute_kstar.html index 7e865c7c..917fcb56 100644 --- a/docs/reference/compute_kstar.html +++ b/docs/reference/compute_kstar.html @@ -1,5 +1,5 @@ -Compute k^* — compute_kstar • BCEACompute k^* — compute_kstar • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.5 + + + + + +
+
+
+ +
+

Q-Q Plot

+
+ +
+

Usage

+
evppi_qq_plot(evppi, he, interv)
+
+ + +
+ + +
+ + + + + + + diff --git a/docs/reference/evppi_residual_plot.html b/docs/reference/evppi_residual_plot.html new file mode 100644 index 00000000..ccc9e646 --- /dev/null +++ b/docs/reference/evppi_residual_plot.html @@ -0,0 +1,93 @@ + +Residual Plot — evppi_residual_plot • BCEA + Skip to contents + + +
+
+
+ +
+

Residual Plot

+
+ +
+

Usage

+
evppi_residual_plot(evppi, he, interv)
+
+ + +
+ + +
+ + + + + + + diff --git a/docs/reference/fit.gam.html b/docs/reference/fit.gam.html index 8e71f2e6..61053755 100644 --- a/docs/reference/fit.gam.html +++ b/docs/reference/fit.gam.html @@ -1,5 +1,5 @@ -Gaussian Additive Model Fitting — fit.gam • BCEAGaussian Additive Model Fitting — fit.gam • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.3 + + + + + +
+
+
+ +
+

Option of interactively saving the plot.

+
+ +
+

Usage

+
plot_mesh(mesh, data, plot, ...)
+
+ +
+

Arguments

+
mesh
+

Mesh

+ + +
data
+

Data

+ + +
plot
+

Create plot? logical

+ + +
...
+

Additional parameters

+ +
+
+

See also

+ +
+ +
+ + +
+ + + +
+ + + + + + + diff --git a/docs/reference/post.density.html b/docs/reference/post.density.html index 2337ef05..dc07d746 100644 --- a/docs/reference/post.density.html +++ b/docs/reference/post.density.html @@ -1,5 +1,5 @@ -Gaussian Process Fitting — post.density • BCEAGaussian Process Fitting — post.density • BCEA @@ -10,7 +10,7 @@ BCEA - 2.4.2 + 2.4.3