adding all new stuff - more to do

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: influ2
 Title: Influence plots
-Version: 0.0.0
+Version: 1.0.0
 Authors@R: 
     person(given = "Darcy",
            family = "Webber",
@@ -38,4 +38,4 @@ Suggests:
     bayesplot
 LazyData: false
 VignetteBuilder: knitr
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2

NAMESPACE

@@ -14,9 +14,6 @@ export(get_marginal)
 export(get_unstandarsied)
 export(id_var_type)
 export(influ_app)
-export(inv_logit)
-export(logistic)
-export(logit)
 export(plot_bayesian_cdi)
 export(plot_bayesian_cdi2)
 export(plot_bubble)
@@ -57,11 +54,9 @@ importFrom(stats,coefficients)
 importFrom(stats,fitted)
 importFrom(stats,median)
 importFrom(stats,model.matrix)
-importFrom(stats,plogis)
 importFrom(stats,poly)
 importFrom(stats,ppoints)
 importFrom(stats,qgamma)
-importFrom(stats,qlogis)
 importFrom(stats,qnorm)
 importFrom(stats,qqplot)
 importFrom(stats,quantile)

R/get-coefs.R

@@ -3,9 +3,6 @@
 #' @param fit An object of class \code{brmsfit}.
 #' @param var The variable to obtain.
 #' @return A \code{data.frame}.
-#' 
-#' @author Darcy Webber \email{[email protected]}
-#' 
 #' @importFrom reshape2 melt
 #' @importFrom readr parse_number
 #' @importFrom nlme fixef ranef
@@ -124,23 +121,21 @@ get_coefs_raw <- function(fit, var = "area") {
 #' @param hurdle If a hurdle model then use the hurdle.
 #' @param transform if the coefficients should be transformed using the link function.
 #' @return A \code{data.frame}.
-#' 
 #' @importFrom reshape2 melt
 #' @importFrom readr parse_number
 #' @importFrom nlme fixef ranef
 #' @importFrom brms posterior_samples
 #' @import dplyr
 #' @export
 #' 
-get_coefs <- function(fit, var = "area", 
+get_coefs <- function(fit, 
+                      var = "area", 
                       normalise = TRUE, 
                       hurdle = FALSE, 
                       transform = FALSE) {
 
   if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.")
-
   type <- id_var_type(fit = fit, xfocus = var, hurdle = hurdle)
-
   is_poly <- FALSE
 
   if (any(grepl("\\(1 \\|", var))) {
@@ -156,8 +151,7 @@ get_coefs <- function(fit, var = "area",
       mutate(variable = gsub(".*\\[|\\]", "", variable)) %>%
       mutate(variable = gsub(",Intercept", "", variable))
     if (!str_detect(var, ":")) {
-      ps <- ps %>%
-        filter(!str_detect(variable, ":"))
+      ps <- ps %>% filter(!str_detect(variable, ":"))
     }
   } else {
     # Population-level effects
@@ -182,11 +176,8 @@ get_coefs <- function(fit, var = "area",
       ps <- ps %>%
         filter(!str_detect(variable, ":"))
     }
-      # mutate(variable = paste0(var, gregexpr("[[:digit:]]+", variable)))
-      # mutate(variable = paste0(var, parse_number(as.character(variable))))
-    # unique(ps$variable)
-    # tail(ps)
   }
+  # unique(ps$variable)
 
   # Check to see if this is a polynomial
   if (any(grepl("poly", ps$variable))) {
@@ -203,8 +194,7 @@ get_coefs <- function(fit, var = "area",
         filter(grepl("hu", variable)) %>%
         mutate(variable = gsub("hu_", "", variable))
     } else {
-      ps <- ps %>%
-        filter(!grepl("hu", variable))
+      ps <- ps %>% filter(!grepl("hu", variable))
     }
   }
 
@@ -232,7 +222,12 @@ get_coefs <- function(fit, var = "area",
                       value = 0)
     coefs <- rbind(ps0, ps)
   } else {
-    coefs <- ps
+    data <- fit$data
+    data[,var] <- paste0(var, data[,var])
+    ps0 <- data.frame(iteration = 1:max(ps$iteration),
+                      variable = unique(data[,var])[!unique(data[,var]) %in% unique(ps$variable)],
+                      value = 0)
+    coefs <- rbind(ps0, ps)
   }
 
   # Arrange by vessel coefficient if vessel chosen

R/get-influ.R

@@ -18,8 +18,7 @@ get_influ2 <- function(fit,
   if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.")
 
   # Model data
-  data <- fit$data %>%
-    mutate(id = 1:n())
+  data <- fit$data %>% mutate(id = 1:n())
   y <- names(data)[1]
 
   # Identify the type of variable we are dealing with
@@ -59,10 +58,8 @@ get_influ2 <- function(fit,
     sdata <- standata(fit)
     X <- sdata$Xs
     Z <- sdata$Zs_1_1
-    s_coefs <- coefs %>%
-      filter(!str_detect(.data$variable, "bs_|sds_"))
-    coefs <- coefs %>%
-      filter(str_detect(.data$variable, "bs_"))
+    s_coefs <- coefs %>% filter(!str_detect(.data$variable, "bs_|sds_"))
+    coefs <- coefs %>% filter(str_detect(.data$variable, "bs_"))
   }
 
   # Do the matrix multiplication
@@ -82,7 +79,7 @@ get_influ2 <- function(fit,
   }
 
   influ_var <- melt(Xbeta) %>%
-    rename(iteration = .data$Var1, id = .data$Var2) %>%
+    rename(iteration = Var1, id = Var2) %>%
     left_join(data, by = "id")
   influ_rho <- influ_var %>%
     group_by(.data$iteration) %>%
@@ -153,7 +150,7 @@ get_influ <- function(fit, group = c("fishing_year", "area"), hurdle = FALSE) {
   }
 
   influ_var <- melt(Xbeta) %>%
-    rename(iteration = .data$Var1, id = .data$Var2) %>%
+    rename(iteration = Var1, id = .data$Var2) %>%
     left_join(data, by = "id")
   influ_rho <- influ_var %>%
     group_by(.data$iteration) %>%

R/helper-functions.R

@@ -73,32 +73,6 @@ get_first_term <- function(fit) {
 }
 
 
-#' 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}.
@@ -150,23 +124,3 @@ id_var_type <- function(fit, xfocus, hurdle = FALSE) {
 geo_mean <- function(a) {
   prod(a)^(1.0 / length(a))
 }
-
-
-#' Inverse logit
-#' 
-#' @param a a vector.
-#' @export
-#' 
-inv_logit <- function(a) {
-  exp(a) / (exp(a) + 1.0)
-}
-
-
-#' logit
-#' 
-#' @param p a vector.
-#' @export
-#' 
-logit <- function(p) {
-  log(p / (1.0 - p))
-}

R/influ.R

@@ -27,7 +27,6 @@
 #' This document provides minimal documentation and the best way to learn how to use this package is through the vignette:
 #'  "An example of how to use the 'influ' package"
 #'
-#' @docType package
 #' @name influ-package
 #' @aliases influ
 #' @title Influence in linear models 

R/influ2.R

@@ -22,7 +22,6 @@
 #' understanding fisheries catch-per-unit-effort standardizations. \emph{ICES Journal of Marine Science}. 
 #' \code{doi:10.1093/icesjms/fsr174}
 #' 
-#' @docType package
 #' @name influ2
 #' 
 NULL

R/plot-cdi.R

@@ -23,11 +23,7 @@
 #'   the same for all and raw, but the legend will change from numbers of records to a proportion.
 #' @param ... Further arguments passed to nothing.
 #' @return a \code{ggplot} object.
-#' 
-#' @author Darcy Webber \email{[email protected]}
-#' 
 #' @seealso \code{\link{get_coefs}}, \code{\link{get_influ}}, \code{\link{plot_bubble}}
-#' 
 #' @importFrom gtable is.gtable gtable_filter
 #' @importFrom stats poly
 #' @importFrom tidyselect all_of
@@ -49,7 +45,6 @@ plot_bayesian_cdi <- function(fit,
                               sum_by = "row", ...) {
 
   if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.")
-
   if (is.null(xlab)) xlab <- xfocus
   if (is.null(ylab)) ylab <- yfocus
   y_coefs <- "Conditional effect"
@@ -58,10 +53,15 @@ plot_bayesian_cdi <- function(fit,
   type <- id_var_type(fit = fit, xfocus = xfocus, hurdle = hurdle)
 
   # Posterior samples of coefficients
+  # if (type %in% c("random_effect")) {
   if (type %in% c("fixed_effect", "random_effect")) {
     coefs <- get_coefs(fit = fit, var = xfocus, hurdle = hurdle)
   } else {
-    coefs <- get_marginal(fit = fit, var = xfocus) # this would plot the marginal/conditional effect, but if it is a hurdle model it ignores the hurdle bit
+    # this would plot the marginal/conditional effect, but if it is a hurdle model it ignores the hurdle bit
+    coefs <- get_marginal(fit = fit, var = xfocus)# %>%
+      # mutate(value = log(value))
+    # library(marginaleffects)
+    # pred <- predictions(model = fit)
     # coefs <- get_coefs_raw(fit = fit, var = xfocus)
   }
 
@@ -73,9 +73,7 @@ plot_bayesian_cdi <- function(fit,
 
   # Model data
   if (is.numeric(coefs$variable)) {
-    data <- fit$data %>%
-      select(all_of(c(yfocus, xfocus)))
-
+    data <- fit$data %>% select(all_of(c(yfocus, xfocus)))
     length.out <- 15
     dmin <- min(data[,xfocus])
     dmax <- max(data[,xfocus])

_pkgdown.yml

@@ -15,11 +15,9 @@ reference:
   contents:
   - influ2
   - starts_with("get_")
-  - logit
-  - inv_logit
-  - logistic
   - id_var_type
   - geo_mean
+  - influ_app
 
 - title: influ2 plots
   desc: Plots and tables from the new influ2 R package.

tests/testthat/test-cdi.R

@@ -76,11 +76,8 @@ test_that("this matches Nokome Bentley's influ package", {
   c2a <- get_coefs(fit = brm1, var = "year", normalise = FALSE, transform = FALSE) %>%
     filter(variable != "year2000") %>%
     group_by(variable) %>%
-    summarise(Estimate = mean(value), 
-              Est.Error = sd(value), 
-              Q5 = quantile(value, probs = 0.05), 
-              Q50 = quantile(value, probs = 0.5), 
-              Q95 = quantile(value, probs = 0.95))
+    summarise(Estimate = mean(value), Est.Error = sd(value), 
+              Q5 = quantile(value, probs = 0.05), Q50 = quantile(value, probs = 0.5), Q95 = quantile(value, probs = 0.95))
 
   c2b <- fixef(object = brm1, probs = c(0.05, 0.95)) %>%
     data.frame() %>%