Skip to content

Commit

Permalink
adding app and getting plots to pick up first term for year
Browse files Browse the repository at this point in the history
  • Loading branch information
quantifish committed Dec 4, 2023
1 parent 74dc975 commit 3d69d8d
Show file tree
Hide file tree
Showing 23 changed files with 190 additions and 136 deletions.
4 changes: 4 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
42 changes: 42 additions & 0 deletions R/app.R
Original file line number Diff line number Diff line change
@@ -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, ...)
}
61 changes: 12 additions & 49 deletions R/get-index.R
Original file line number Diff line number Diff line change
@@ -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)) %>%
Expand Down Expand Up @@ -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{[email protected]}
#'
#' @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)
Expand All @@ -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) %>%
Expand Down
50 changes: 42 additions & 8 deletions R/helper-functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -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{[email protected]}
#'
#' @export
#'
id_var_type <- function(fit, xfocus, hurdle = FALSE) {
Expand Down Expand Up @@ -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{[email protected]}
#'
#' @export
#'
geo_mean <- function(a) {
Expand All @@ -119,7 +155,6 @@ geo_mean <- function(a) {
#' Inverse logit
#'
#' @param a a vector.
#' @author Darcy Webber \email{[email protected]}
#' @export
#'
inv_logit <- function(a) {
Expand All @@ -130,7 +165,6 @@ inv_logit <- function(a) {
#' logit
#'
#' @param p a vector.
#' @author Darcy Webber \email{[email protected]}
#' @export
#'
logit <- function(p) {
Expand Down
4 changes: 2 additions & 2 deletions R/plot-compare.R
Original file line number Diff line number Diff line change
Expand Up @@ -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) {

Expand All @@ -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]
}
Expand Down
11 changes: 6 additions & 5 deletions R/plot-index.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,24 +9,25 @@
#' @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{[email protected]}
#'
#' @importFrom stats fitted
#' @import brms
#' @import ggplot2
#' @import dplyr
#' @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")
Expand Down
12 changes: 6 additions & 6 deletions R/plot-influ.R
Original file line number Diff line number Diff line change
Expand Up @@ -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{[email protected]}
#'
#' @seealso \code{\link{get_influ}}
#'
#' @examples
#' \dontrun{
#' data(epilepsy)
Expand All @@ -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]]
Expand Down
5 changes: 1 addition & 4 deletions R/plot-step.R
Original file line number Diff line number Diff line change
Expand Up @@ -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{[email protected]}
#'
#' @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)
Expand Down
Loading

0 comments on commit 3d69d8d

Please sign in to comment.