diff --git a/NAMESPACE b/NAMESPACE index ddd1076..201d55f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,12 +6,14 @@ export(geo_mean) export(get_bayes_R2) export(get_coefs) export(get_coefs_raw) +export(get_first_term) export(get_index) export(get_influ) export(get_influ2) export(get_marginal) export(get_unstandarsied) export(id_var_type) +export(influ_app) export(inv_logit) export(logistic) export(logit) @@ -34,7 +36,9 @@ import(dplyr) import(ggplot2) import(patchwork) import(proto) +import(shiny) importFrom(brms,add_criterion) +importFrom(brms,is.brmsfit) importFrom(brms,posterior_samples) importFrom(gtable,gtable_filter) importFrom(gtable,is.gtable) diff --git a/R/app.R b/R/app.R new file mode 100644 index 0000000..70a67ab --- /dev/null +++ b/R/app.R @@ -0,0 +1,42 @@ +#' Run shiny application +#' +#' @param fits a list of objects of class \code{brmsfit} that you want to compare. +#' @return a \code{shiny.appobj} object. +#' @import shiny +#' @export +#' +influ_app <- function(fits) { + + n <- length(fits) + flabs <- rep(NA, n) + + for (i in 1:n) { + str <- as.character(fits[[i]]$formula)[1] + left1 <- stringi::stri_trim_right(str = str, pattern = "[\u007E]", negate = FALSE) + flabs[i] <- substr(x = str, nchar(left1) + 2, nchar(str)) + } + + ui <- navbarPage( + title = "influ2", + tabPanel( + title = "Model comparisons", + tableOutput(outputId = "table_cf"), + selectInput(inputId = "pick_fits", label = "Select the fits", choices = flabs, selected = flabs, multiple = TRUE), + plotOutput(outputId = "plot_cf", click = "plot_click") + ) + ) + + server <- function(input, output, session) { + output$table_cf <- renderTable({ + get_bayes_R2(fits = fits) + }) + output$plot_cf <- renderPlot({ + ii <- which(input$pick_fits %in% flabs) + ifits <- fits[ii] + plot_compare(fits = ifits) + }) + } + + shinyApp(ui = ui, server = server) + # shinyApp(ui = ui, server = server, ...) +} diff --git a/R/get-index.R b/R/get-index.R index 632c886..ef87ea8 100644 --- a/R/get-index.R +++ b/R/get-index.R @@ -1,43 +1,21 @@ -#' Logit -#' -#' @inheritParams stats::qlogis -#' @return the density. -#' @importFrom stats qlogis -#' @export -#' -logit <- function(p, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) { - qlogis(p = p, location = location, scale = scale, lower.tail = lower.tail, log.p = log.p) -} - - -#' Logistic -#' -#' @param q vector of quantiles. -#' @inheritParams stats::plogis -#' @param log.p logical; if TRUE, probabilities p are given as log(p). -#' @return gives the distribution function. -#' @importFrom stats plogis -#' @export -#' -logistic <- function(q, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) { - plogis(q = q, location = location, scale = scale, lower.tail = lower.tail, log.p = log.p) -} - - #' Get the unstandardised indices #' #' @param fit An object of class \code{brmsfit}. #' @param year The year or time label (e.g. year, Year, fishing_year, etc). #' @param rescale How to re-scale the series. Choose from "raw" to retain the raw unstandardised series, or a number to re-scale by. #' @return a \code{data.frame} or a \code{ggplot} object. -#' @import brms +#' @importFrom brms is.brmsfit #' @import dplyr #' @export #' -get_unstandarsied <- function(fit, year = "year", rescale = 1) { +get_unstandarsied <- function(fit, year = NULL, rescale = 1) { if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.") + if (is.null(year)) { + year <- get_first_term(fit = fit) + } + if (fit$family$family %in% c("bernoulli", "negbinomial") | grepl("hurdle", fit$family$family)) { prop <- data.frame(y = fit$data[,1], Year = fit$data[,year]) %>% mutate(y = ifelse(.data$y > 0, 1, 0)) %>% @@ -84,20 +62,21 @@ get_unstandarsied <- function(fit, year = "year", rescale = 1) { #' @param do_plot Return a \code{ggplot} object instead of a \code{data.frame}. #' @param ... Additional parameters passed to \code{fitted}. #' @return a \code{data.frame} or a \code{ggplot} object. -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @importFrom stats fitted -#' @import brms +#' @importFrom brms is.brmsfit #' @import ggplot2 #' @import patchwork #' @import dplyr #' @export #' -get_index <- function(fit, year = "year", probs = c(0.025, 0.975), rescale = 1, do_plot = FALSE, ...) { +get_index <- function(fit, year = NULL, probs = c(0.025, 0.975), rescale = 1, do_plot = FALSE, ...) { if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.") + if (is.null(year)) { + year <- get_first_term(fit = fit) + } + # std <- get_coefs(fit = fit, var = year) yrs <- sort(unique(fit$data[,year])) n <- length(yrs) @@ -120,22 +99,6 @@ get_index <- function(fit, year = "year", probs = c(0.025, 0.975), rescale = 1, } newdata[,year] <- yrs newdata$pots <- 1 - - # fout1 <- fitted(object = fit, newdata = newdata, probs = c(probs[1], 0.5, probs[2]), re_formula = NA) - # newdata <- newdata[,1:5] - # newdata <- expand.grid(cpue = 1, period = unique(celr5$period), area2 = NA, vessel = NA, month = NA, "period:area2" = NA) - # names(newdata) <- names(fit5$data) - # head(fit$data) - # head(newdata) - # fout1 <- fitted(object = fit, newdata = newdata) - # fout1 <- fitted(object = fit4, newdata = newdata, probs = c(probs[1], 0.5, probs[2])) - # fout1 <- conditional_effects(x = fit, effects = "period")[[1]] - # conditions <- data.frame(area2 = c("916", "917", "933")) - # fout1 <- conditional_effects(x = fit, effects = "period", conditions = conditions)[[1]] - # fout1 <- posterior_epred(object = fit, newdata = newdata, probs = c(probs[1], 0.5, probs[2])) - # fout1 <- posterior_predict(object = fit, newdata = newdata) - # fout1 <- fitted(object = fit) - # pred1 <- predict(fit, newdata = newdata, re_formula = NULL, allow_new_levels = TRUE) # Get the predicted CPUE by year fout1 <- fitted(object = fit, newdata = newdata, probs = c(probs[1], 0.5, probs[2]), re_formula = NA) %>% diff --git a/R/helper-functions.R b/R/helper-functions.R index ee69f21..58b6f02 100644 --- a/R/helper-functions.R +++ b/R/helper-functions.R @@ -57,15 +57,54 @@ glm_step_plot <- function(data, mod_list, ibest = 5) { } +#' Get first model term +#' +#' @param fit An object of class \code{brmsfit}. +#' @return the first model term +#' @importFrom brms is.brmsfit +#' @export +#' +get_first_term <- function(fit) { + if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.") + f1 <- as.character(fit$formula)[1] + f2 <- gsub("~", "+", f1) + f3 <- str_split(f2, " \\+ ")[[1]] + return(f3[2]) +} + + +#' Logit +#' +#' @inheritParams stats::qlogis +#' @return the density. +#' @importFrom stats qlogis +#' @export +#' +logit <- function(p, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) { + qlogis(p = p, location = location, scale = scale, lower.tail = lower.tail, log.p = log.p) +} + + +#' Logistic +#' +#' @param q vector of quantiles. +#' @inheritParams stats::plogis +#' @param log.p logical; if TRUE, probabilities p are given as log(p). +#' @return gives the distribution function. +#' @importFrom stats plogis +#' @export +#' +logistic <- function(q, location = 0, scale = 1, lower.tail = TRUE, log.p = FALSE) { + plogis(q = q, location = location, scale = scale, lower.tail = lower.tail, log.p = log.p) +} + + #' Identify the variable type #' #' @param fit An object of class \code{brmsfit}. #' @param xfocus the x #' @param hurdle if hurdle or not #' @return The geometric mean of the vector. -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @export #' id_var_type <- function(fit, xfocus, hurdle = FALSE) { @@ -106,9 +145,6 @@ id_var_type <- function(fit, xfocus, hurdle = FALSE) { #' #' @param a a vector. #' @return The geometric mean of the vector. -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @export #' geo_mean <- function(a) { @@ -119,7 +155,6 @@ geo_mean <- function(a) { #' Inverse logit #' #' @param a a vector. -#' @author Darcy Webber \email{darcy@quantifish.co.nz} #' @export #' inv_logit <- function(a) { @@ -130,7 +165,6 @@ inv_logit <- function(a) { #' logit #' #' @param p a vector. -#' @author Darcy Webber \email{darcy@quantifish.co.nz} #' @export #' logit <- function(p) { diff --git a/R/plot-compare.R b/R/plot-compare.R index 8028867..16bcd18 100644 --- a/R/plot-compare.R +++ b/R/plot-compare.R @@ -16,7 +16,7 @@ #' @importFrom stringi stri_trim_right #' @export #' -plot_compare <- function(fits, labels = NULL, year = "year", +plot_compare <- function(fits, labels = NULL, year = NULL, probs = c(0.25, 0.75), show_probs = TRUE, rescale = "raw", rescale_series = NULL) { @@ -27,7 +27,7 @@ plot_compare <- function(fits, labels = NULL, year = "year", if (is.null(labels)) { str <- as.character(fits[[i]]$formula)[1] left1 <- stri_trim_right(str = str, pattern = "[\u007E]", negate = FALSE) - fout$Model <- substr(str, nchar(left1) + 2, nchar(str)) + fout$Model <- substr(x = str, nchar(left1) + 2, nchar(str)) } else { fout$Model <- labels[i] } diff --git a/R/plot-index.R b/R/plot-index.R index 6f75d8d..858da61 100644 --- a/R/plot-index.R +++ b/R/plot-index.R @@ -9,9 +9,6 @@ #' @param rescale the index of the series to rescale to. If set to NULL then no rescaling is done. #' @param show_unstandardised show the unstandardised series or not. #' @return a \code{ggplot} object. -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @importFrom stats fitted #' @import brms #' @import ggplot2 @@ -19,14 +16,18 @@ #' @export #' plot_index <- function(fit, - year = "Year", + year = NULL, fill = "purple", probs = c(0.25, 0.75), rescale = 1, show_unstandardised = TRUE) { if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.") - + + if (is.null(year)) { + year <- get_first_term(fit = fit) + } + # Get the standardised series fout <- get_index(fit = fit, year = year, probs = probs, rescale = rescale) %>% mutate(model = "Standardised") diff --git a/R/plot-influ.R b/R/plot-influ.R index e0873e4..2af940e 100644 --- a/R/plot-influ.R +++ b/R/plot-influ.R @@ -7,11 +7,7 @@ #' @param fill the colour to use in the plot. #' @param hurdle if a hurdle model then use the hurdle #' @return a \code{ggplot} object. -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @seealso \code{\link{get_influ}} -#' #' @examples #' \dontrun{ #' data(epilepsy) @@ -20,14 +16,18 @@ #' summary(fit1) #' plot_influ(fit = fit1, year = "Age") #' } -#' +#' @importFrom brms is.brmsfit #' @import ggplot2 #' @export #' -plot_influ <- function(fit, year = "fishing_year", fill = "purple", hurdle = FALSE) { +plot_influ <- function(fit, year = NULL, fill = "purple", hurdle = FALSE) { if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.") + if (is.null(year)) { + year <- get_first_term(fit = fit) + } + # Extract the models variable names # x1 <- gsub(paste0(as.character(fit$formula)[4], " ~ "), "", as.character(fit$formula)[1]) # x2 <- strsplit(x1, split = " + ", fixed = TRUE)[[1]] diff --git a/R/plot-step.R b/R/plot-step.R index 3551b33..17b7078 100644 --- a/R/plot-step.R +++ b/R/plot-step.R @@ -7,15 +7,12 @@ #' @param fill the colour of the credible interval ribbon #' @param probs the quantiles to plot #' @param show_probs plot the quantiles or not -#' -#' @author Darcy Webber \email{darcy@quantifish.co.nz} -#' #' @import brms #' @import ggplot2 #' @import dplyr #' @export #' -plot_step <- function(fits, year = "year", fill = "purple", +plot_step <- function(fits, year = NULL, fill = "purple", probs = c(0.25, 0.75), show_probs = TRUE) { m <- length(fits) diff --git a/_pkgdown.yml b/_pkgdown.yml index 53b7eb0..106bf45 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -10,6 +10,28 @@ navbar: right: [search, github] reference: +- title: influ2 functions + desc: Functions from the new influ2 R package. + contents: + - influ2 + - starts_with("get_") + - logit + - inv_logit + - logistic + - id_var_type + - geo_mean + +- title: influ2 plots + desc: Plots and tables from the new influ2 R package. + contents: + - starts_with("plot_") + - starts_with("table_") + +- title: Data + desc: Simulated data sets that are available. + contents: + - lobsters_per_pot + - title: influ desc: Functions from the original influ R package. contents: @@ -27,21 +49,4 @@ reference: - Influence$stepAndInfluPlot - coeffs - ses - -- title: influ2 - desc: Functions from the new influ2 R package. - contents: - - influ2 - - starts_with("get_") - - starts_with("plot_") - - starts_with("table_") - - logit - - inv_logit - - logistic - - id_var_type - - geo_mean - -- title: Data - desc: Simulated data sets that are available. - contents: - - lobsters_per_pot + \ No newline at end of file diff --git a/man/geo_mean.Rd b/man/geo_mean.Rd index 7a1282c..01bcfad 100644 --- a/man/geo_mean.Rd +++ b/man/geo_mean.Rd @@ -15,6 +15,3 @@ The geometric mean of the vector. \description{ Geometric mean } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/get_first_term.Rd b/man/get_first_term.Rd new file mode 100644 index 0000000..fc67b29 --- /dev/null +++ b/man/get_first_term.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helper-functions.R +\name{get_first_term} +\alias{get_first_term} +\title{Get first model term} +\usage{ +get_first_term(fit) +} +\arguments{ +\item{fit}{An object of class \code{brmsfit}.} +} +\value{ +the first model term +} +\description{ +Get first model term +} diff --git a/man/get_index.Rd b/man/get_index.Rd index 6ad760e..0e89ad1 100644 --- a/man/get_index.Rd +++ b/man/get_index.Rd @@ -6,7 +6,7 @@ \usage{ get_index( fit, - year = "year", + year = NULL, probs = c(0.025, 0.975), rescale = 1, do_plot = FALSE, @@ -32,6 +32,3 @@ a \code{data.frame} or a \code{ggplot} object. \description{ Get the standardised indices each year with associated uncertainty and return either as a table or a ggplot. } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/get_unstandarsied.Rd b/man/get_unstandarsied.Rd index 4030441..3a2b2ff 100644 --- a/man/get_unstandarsied.Rd +++ b/man/get_unstandarsied.Rd @@ -4,7 +4,7 @@ \alias{get_unstandarsied} \title{Get the unstandardised indices} \usage{ -get_unstandarsied(fit, year = "year", rescale = 1) +get_unstandarsied(fit, year = NULL, rescale = 1) } \arguments{ \item{fit}{An object of class \code{brmsfit}.} diff --git a/man/id_var_type.Rd b/man/id_var_type.Rd index 8a503cc..75ec724 100644 --- a/man/id_var_type.Rd +++ b/man/id_var_type.Rd @@ -19,6 +19,3 @@ The geometric mean of the vector. \description{ Identify the variable type } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/influ_app.Rd b/man/influ_app.Rd new file mode 100644 index 0000000..b42e17c --- /dev/null +++ b/man/influ_app.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/app.R +\name{influ_app} +\alias{influ_app} +\title{Run shiny application} +\usage{ +influ_app(fits) +} +\arguments{ +\item{fits}{a list of objects of class \code{brmsfit} that you want to compare.} +} +\value{ +a \code{shiny.appobj} object. +} +\description{ +Run shiny application +} diff --git a/man/inv_logit.Rd b/man/inv_logit.Rd index 869e858..a67392b 100644 --- a/man/inv_logit.Rd +++ b/man/inv_logit.Rd @@ -12,6 +12,3 @@ inv_logit(a) \description{ Inverse logit } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/logistic.Rd b/man/logistic.Rd index 7dbcda2..e0b3df8 100644 --- a/man/logistic.Rd +++ b/man/logistic.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/get-index.R +% Please edit documentation in R/helper-functions.R \name{logistic} \alias{logistic} \title{Logistic} diff --git a/man/logit.Rd b/man/logit.Rd index abcb40b..3a55a3a 100644 --- a/man/logit.Rd +++ b/man/logit.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/get-index.R, R/helper-functions.R +% Please edit documentation in R/helper-functions.R \name{logit} \alias{logit} \title{Logit} @@ -19,6 +19,3 @@ Logit logit } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/plot_compare.Rd b/man/plot_compare.Rd index fe97458..a02021c 100644 --- a/man/plot_compare.Rd +++ b/man/plot_compare.Rd @@ -7,7 +7,7 @@ plot_compare( fits, labels = NULL, - year = "year", + year = NULL, probs = c(0.25, 0.75), show_probs = TRUE, rescale = "raw", diff --git a/man/plot_index.Rd b/man/plot_index.Rd index c4f8331..8c6ef95 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -6,7 +6,7 @@ \usage{ plot_index( fit, - year = "Year", + year = NULL, fill = "purple", probs = c(0.25, 0.75), rescale = 1, @@ -32,6 +32,3 @@ a \code{ggplot} object. \description{ In this plot the unstandardised indices is the geometric mean of the data. } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/plot_influ.Rd b/man/plot_influ.Rd index b03716f..5dd766e 100644 --- a/man/plot_influ.Rd +++ b/man/plot_influ.Rd @@ -4,7 +4,7 @@ \alias{plot_influ} \title{Plot the influence metric for all variables in a model} \usage{ -plot_influ(fit, year = "fishing_year", fill = "purple", hurdle = FALSE) +plot_influ(fit, year = NULL, fill = "purple", hurdle = FALSE) } \arguments{ \item{fit}{An object of class \code{brmsfit}.} @@ -29,11 +29,7 @@ fit1 <- brm(count ~ Age + (1|patient), data = epilepsy, family = poisson(), iter summary(fit1) plot_influ(fit = fit1, year = "Age") } - } \seealso{ \code{\link{get_influ}} } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/man/plot_step.Rd b/man/plot_step.Rd index aa2c7ff..fc3be9f 100644 --- a/man/plot_step.Rd +++ b/man/plot_step.Rd @@ -6,7 +6,7 @@ \usage{ plot_step( fits, - year = "year", + year = NULL, fill = "purple", probs = c(0.25, 0.75), show_probs = TRUE @@ -26,6 +26,3 @@ plot_step( \description{ This requires that all steps be run using brms and then provided as a list of model fits. } -\author{ -Darcy Webber \email{darcy@quantifish.co.nz} -} diff --git a/vignettes/influ2.Rmd b/vignettes/influ2.Rmd index a107d80..4963e5e 100644 --- a/vignettes/influ2.Rmd +++ b/vignettes/influ2.Rmd @@ -56,7 +56,7 @@ theme_set(theme_bw()) options(mc.cores = 2) ``` -The simulated data set +The simulated data set: ```{r sim-data, echo=TRUE, message=FALSE} data(lobsters_per_pot) @@ -131,11 +131,11 @@ be modelled separately and then a list of models defined as below. ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "The new Bayesian step-plot from the influ2 package. Note that the new step-plot requires that all models/steps be run in brms before the function can be used."} fits <- list(fit1, fit2) -plot_step(fits = fits, year = "year") +plot_step(fits = fits) ``` ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE} -plot_index(fit2, year = "year") +plot_index(fit = fit2) ``` ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE} @@ -143,7 +143,7 @@ myInfl$stanPlot() ``` ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE} -plot_influ(fit = fit1, year = "year") +plot_influ(fit = fit1) ``` ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE} @@ -217,11 +217,11 @@ plot_bayesian_cdi(fit = fit2, xfocus = "month", yfocus = "year", ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "Comparison plot of all model runs."} fits <- list(fit1, fit2) -plot_compare(fits = fits, year = "year") +plot_compare(fits = fits) ``` ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "Unstandardised versus standardised series."} -plot_index(fit = fit1, year = "year") +plot_index(fit = fit1) ``` # Diagnostics @@ -312,8 +312,8 @@ can be produced using the `get_index` function. The `Estimate` and `Est.Error` columns are what should be used in stock assessment. ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE} -unstd <- get_unstandarsied(fit = fit2, year = "year") -cpue0 <- get_index(fit = fit2, year = "year") +unstd <- get_unstandarsied(fit = fit2) +cpue0 <- get_index(fit = fit2) cpue0 %>% mutate(Unstandardised = unstd$Median) %>% @@ -349,9 +349,8 @@ for (ii in 1:nrow(cpue1)) { cpue3$Qupper[ii] <- quantile(r1, probs = 0.975) } -df <- bind_rows(cpue1, cpue2, cpue3) - -ggplot(data = df, aes(x = Year, y = Median, colour = Type)) + +ggplot(data = bind_rows(cpue1, cpue2, cpue3), + aes(x = Year, y = Median, colour = Type)) + geom_pointrange(aes(ymin = Qlower, ymax = Qupper), position = position_dodge(width = 0.5)) + scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +