more work on vig and tests

.gitignore

@@ -1,3 +1,5 @@
+data-raw/
+
 # History files
 .Rhistory
 .Rapp.history

DESCRIPTION

@@ -34,11 +34,7 @@ Suggests:
     tidyverse,
     devtools,
     rmarkdown,
-    roxygen2,
-    testthat,
-    usethis,
     viridis,
-    proto,
     bayesplot
 LazyData: false
 VignetteBuilder: knitr

R/influ.R

@@ -134,14 +134,14 @@ coeffs = function(.,model=.$model,term=.$focus){
 #' @param term The term in the model for which coefficients SEs are extracted
 #' @export
 #' 
-ses = function(.,model=.$model,term=.$focus){
+ses = function(.,model=.$model, term=.$focus) {
   type = class(model)[1]
-  if(type=='glm'|type=='negbin'){
+  if (type == 'glm' | type == 'negbin') {
     #The "cov.scaled" member of a glm object is the estimated covariance matrix of 
     #the estimated coefficients scaled by dispersion"    
     V = summary(model)$cov.scaled
   }
-  else if(type=='survreg'){
+  else if (type == 'survreg') {
     #The member "var" for a survreg object is the "the variance-covariance matrix for the parameters, 
     #including the log(scale) parameter(s)" 
     V = model$var
@@ -179,18 +179,18 @@ Influence$effects = function(.,model=.$model,term=.$focus){
 Influence$calc <- function(.){
   #Get observed values
   observed = .$model$model[,.$response]
-  if(class(observed)[1]=='Surv') observed = as.numeric(observed[,1])
-  logged = substr(.$response,1,4)=='log('
+  if (class(observed)[1] == 'Surv') observed = as.numeric(observed[,1])
+  logged = substr(.$response,1,4) == 'log('
   #Create a data frame that is used to store various values for each level of focal term
   .$indices = data.frame(level=levels(.$model$model[,.$focus]))
   #Add an unstandardised index
-  if(logged | sum(observed<=0)==0){
-    if(logged) log_observed = observed else log_observed = log(observed)
+  if (logged | sum(observed<=0)==0){
+    if (logged) log_observed = observed else log_observed = log(observed)
     # Calculate geometic mean
     .$indices = merge(.$indices,aggregate(list(unstan=log_observed),list(level=.$model$model[,.$focus]),mean))
     # Turn into a relative index
     .$indices$unstan = with(.$indices,exp(unstan-mean(unstan)))
-  }else {
+  } else {
      # There are zeros (or even negative numbers) in the data so we cant calculate a geometric mean.
      # Use arithmetic mean instead
     .$indices = merge(.$indices,aggregate(list(unstan=observed),list(level=.$model$model[,.$focus]),mean))
@@ -209,8 +209,8 @@ Influence$calc <- function(.){
   .$summary = NULL
 
   #TODO calculate influence statistics in this loop too
-  for(termCount in 0:length(.$terms)){
-    if(termCount>0){
+  for (termCount in 0:length(.$terms)) {
+    if (termCount > 0) {
       term = .$terms[termCount]
       #Update both formula and model
       newFormula = formula(.$model)

R/plot-bubble.R

@@ -5,7 +5,6 @@
 #' @param data A \code{data.frame} containing at least two columns.
 #' @param xvar The name of the temporal column (e.g., year).
 #' @param yvar The names of the columns in the \code{data.frame} to plot.
-#' @param ... Further arguments passed to nothing.
 #' @return a ggplot object
 #' 
 #' @author Darcy Webber \email{[email protected]}

R/plot-cdi.R

@@ -4,11 +4,13 @@
 #' (bottom-left panel), and the influence statistic (bottom-right panel).
 #' 
 #' @param fit An object of class \code{brmsfit}.
-#' @param xfocus The column name of the variable to be plotted on the x axis. This column name must match one of the
-#'   column names in the \code{data.frame} that was passed to \code{brm} as the \code{data} argument.
-#' @param yfocus The column name of the variable to be plotted on the y axis. This column name must match one of the
-#'   column names in the \code{data.frame} that was passed to \code{brm} as the \code{data} argument. This is generally the
-#'   temporal variable in a generalised linear model (e.g. year).
+#' @param xfocus The column name of the variable to be plotted on the x axis. 
+#'   This column name must match one of the column names in the 
+#'   \code{data.frame} that was passed to \code{brm} as the \code{data} argument.
+#' @param yfocus The column name of the variable to be plotted on the y axis. 
+#'   This column name must match one of the column names in the 
+#'   \code{data.frame} that was passed to \code{brm} as the \code{data} argument. 
+#'   This is generally the temporal variable in a generalised linear model (e.g. year).
 #' @param hurdle If a hurdle model then use the hurdle.
 #' @param sort_coefs Should the coefficients be sorted from highest to lowest.
 #' @param axis.text.x.bl Include the x axis labels on the bottom-left (bl) bubble plot panel.
@@ -152,7 +154,7 @@ plot_bayesian_cdi <- function(fit,
     p1 <- p1 +
       # geom_point() +
       geom_violin(colour = colour, fill = colour, alpha = 0.5, draw_quantiles = 0.5, scale = "width") +
-      geom_hline(yintercept = 1, linetype = "dashed") +
+      geom_hline(yintercept = 0, linetype = "dashed") +
       scale_x_discrete(position = "top")# +
       # scale_x_discrete(position = "top", breaks = midpoints, minor_breaks = NULL, expand = expansion(mult = 0.05)) +
       # coord_cartesian(xlim = c(midpoints[1], midpoints[length(midpoints)]))

R/plot-index.R

@@ -7,6 +7,7 @@
 #' @param fill the fill colour for the percentiles.
 #' @param probs The percentiles to be computed by the \code{quantile} function.
 #' @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]}

R/plot-influ.R

@@ -29,24 +29,28 @@ plot_influ <- function(fit, year = "fishing_year", fill = "purple", hurdle = FAL
   if (!is.brmsfit(fit)) stop("fit is not an object of class brmsfit.")
 
   # 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]]
-  x <- x2[x2 != year]
-  x <- gsub("\\(1 \\| ", "", x)
-  x <- gsub("\\)", "", x)
-
+  # x1 <- gsub(paste0(as.character(fit$formula)[4], " ~ "), "", as.character(fit$formula)[1])
+  # x2 <- strsplit(x1, split = " + ", fixed = TRUE)[[1]]
+  # x <- x2[x2 != year]
+  # x <- gsub("\\(1 \\| ", "", x)
+  # x <- gsub("\\)", "", x)
+  # x <- gsub("poly\\(", "", x)
+  x <- names(fit$data)
+  x <- x[-1]
+  x <- x[x != year]
+  x <- x[!grepl("poly", x)]
+
   df <- NULL
   for (i in 1:length(x)) {
     inf1 <- get_influ2(fit = fit, group = c(year, x[i]), hurdle = hurdle) %>% 
       mutate(variable = x[i])
-
     df <- rbind(df, inf1)
   }
 
   # Order the factor levels
   df$variable <- factor(df$variable, levels = x)
 
-  p <- ggplot(data = df, aes_string(x = year)) +
+  p <- ggplot(data = df, aes(x = .data[[year]])) +
     geom_hline(yintercept = 1, linetype = "dashed") +
     geom_violin(aes(y = exp(.data$delta)), fill = fill, colour = fill, alpha = 0.5, draw_quantiles = 0.5, scale = "width") +
     facet_wrap(variable ~ ., ncol = 1, strip.position = "top") +

_pkgdown.yml

@@ -1,4 +1,47 @@
 url: http://www.quantifish.co.nz/influ2/
+
 template:
   bootstrap: 5
+  bootswatch: cerulean
+
+navbar:
+  structure:
+    left:  [intro, articles, tutorials, reference, news]
+    right: [search, github]
+
+reference:
+- title: influ
+  desc: Functions from the original influ R package.
+  contents:
+  - influ
+  - Influence
+  - Influence$new
+  - Influence$init
+  - Influence$calc
+  - Influence$effects
+  - Influence$stanPlot
+  - Influence$cdiPlot
+  - Influence$influPlot
+  - Influence$stepPlot
+  - Influence$cdiPlotAll
+  - 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

tests/testthat/test-cdi.R

@@ -1,12 +1,10 @@
 library(brms)
 
-options(mc.cores = parallel::detectCores())
-
 data(lobsters_per_pot)
 
-brm1 <- brm(lobsters ~ year + month, data = lobsters_per_pot, family = poisson)
+brm1 <- brm(lobsters ~ year + month, data = lobsters_per_pot, family = poisson, chains = 1, file = "brm1", file_refit = "never")
 glm1 <- glm(lobsters ~ year + month, data = lobsters_per_pot, family = poisson)
-brm2 <- brm(lobsters ~ year + (1 | month), data = lobsters_per_pot, family = negbinomial())
+brm2 <- brm(lobsters ~ year + (1 | month), data = lobsters_per_pot, family = negbinomial(), chains = 1, file = "brm2", file_refit = "never")
 
 
 context("CDI plot")

vignettes/hurdle.Rmd

@@ -57,7 +57,7 @@ head(sim_data)
 ```
 
 ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "Distribution of data by year and group."}
-plot_bubble(df = sim_data, group = c("year", "group"), fill = "green")
+plot_bubble(df = sim_data, group = c("year", "group"))
 ```
 
 ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE}
@@ -66,16 +66,14 @@ m1 <- brm(bf(y ~ year + group, hu ~ 1),
           file = "m1", file_refit = "never")
 ```
 
-```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "The new Bayesian CDI plot from the influ2 package."}
-# plot_bayesian_cdi(fit = m1, xfocus = "group", yfocus = "year", xlab = "Group")
-plot_index(m1, year = "year")
-```
-
 ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE}
 plot_index(m1, year = "year")
 ```
 
 ```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE}
 # plot_influ(m1, year = "year")
-plot_index(m1, year = "year")
+```
+
+```{r echo=TRUE, fig.height=6, fig.width=6, message=FALSE, fig.cap = "The new Bayesian CDI plot from the influ2 package."}
+# plot_bayesian_cdi(fit = m1, xfocus = "group", yfocus = "year", xlab = "Group")
 ```