Skip to content

Commit

Permalink
tests good again
Browse files Browse the repository at this point in the history
  • Loading branch information
quantifish committed Nov 28, 2023
1 parent 7b7147f commit 74dc975
Show file tree
Hide file tree
Showing 24 changed files with 92 additions and 68 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ License: MIT + file LICENSE
Encoding: UTF-8
Depends:
R (>= 3.5.0),
brms,
dplyr,
ggplot2,
patchwork
Expand All @@ -25,7 +26,6 @@ Imports:
readr,
gtable,
nlme,
brms,
rstan,
tidyr,
lubridate
Expand Down
4 changes: 4 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -49,16 +49,20 @@ importFrom(reshape2,melt)
importFrom(rstan,get_elapsed_time)
importFrom(rstan,get_num_upars)
importFrom(stats,as.formula)
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)
importFrom(stats,residuals)
importFrom(stats,setNames)
importFrom(stringi,stri_trim_right)
importFrom(stringr,str_detect)
importFrom(stringr,str_split)
Expand Down
11 changes: 9 additions & 2 deletions R/data.R
Original file line number Diff line number Diff line change
@@ -1,7 +1,14 @@
#' Simulated catch of lobsters per pot
#' Simulated CPUE data
#'
#' Simulated catch of lobsters per pot.
#'
#' @format A data.frame with four fields:
#' @format a \code{tibble} containing 5 fields including:
#' \describe{
#' \item{lobsters}{Age in years.}
#' \item{year}{Standard deviation of length at age.}
#' \item{month}{Standard deviation of length at age.}
#' \item{depth}{Standard deviation of length at age.}
#' \item{soak}{Standard deviation of length at age.}
#' }
#'
"lobsters_per_pot"
31 changes: 16 additions & 15 deletions R/get-index.R
Original file line number Diff line number Diff line change
@@ -1,23 +1,27 @@
#' Logit
#'
#' @return a \code{data.frame} or a \code{ggplot} object.
#'
#' @author Darcy Webber \email{[email protected]}
#'
#' @inheritParams stats::qlogis
#' @return the density.
#' @importFrom stats qlogis
#' @export
#'
logit <- qlogis
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
#'
#' @return a \code{data.frame} or a \code{ggplot} object.
#'
#' @author Darcy Webber \email{[email protected]}
#'
#' @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 <- plogis
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
Expand All @@ -26,9 +30,6 @@ logistic <- plogis
#' @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.
#'
#' @author Darcy Webber \email{[email protected]}
#'
#' @import brms
#' @import dplyr
#' @export
Expand Down Expand Up @@ -57,8 +58,8 @@ get_unstandarsied <- function(fit, year = "year", rescale = 1) {
gm <- geo_mean(unstd$cpue)

fout <- unstd %>%
mutate(Mean = cpue, Median = cpue) %>%
select(-cpue)
mutate(Mean = .data$cpue, Median = .data$cpue) %>%
select(-.data$cpue)

# Rescale the series
if (rescale == "raw") {
Expand Down
3 changes: 0 additions & 3 deletions R/get-influ.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,6 @@
#' @param group the variable to obtain
#' @param hurdle if a hurdle model then use the hurdle
#' @return a \code{data.frame}.
#'
#' @author Darcy Webber \email{[email protected]}
#'
#' @importFrom reshape2 melt
#' @importFrom readr parse_number
#' @importFrom stringr str_split str_detect
Expand Down
8 changes: 4 additions & 4 deletions R/get-tables.R
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,16 @@
#' @param criterion the criterion to use
#' @param sort to sort the table with the best models at the top
#' @param ... additional parameters passed to \code{add_criterion}
#'
#' @author Darcy Webber \email{[email protected]}
#'
#' @import brms
#' @import dplyr
#' @importFrom rstan get_elapsed_time get_num_upars
#' @importFrom lubridate seconds_to_period hour minute second
#' @importFrom stats setNames
#' @export
#'
table_criterion <- function(fits, criterion = c("loo", "loo_R2", "bayes_R2", "log_lik"), sort = TRUE, ...) {
table_criterion <- function(fits,
criterion = c("loo", "loo_R2", "bayes_R2", "log_lik"),
sort = TRUE, ...) {

n <- length(fits)

Expand Down
36 changes: 20 additions & 16 deletions R/influ.R
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,14 @@ Influence = proto()
#' }
#'
#' @name Influence$new
#' @param . Not sure how this works
#' @param model The model for which the influence of terms will be determined
#' @param data The data which was used to fit the model (required for some types of models that do not store the data)
#' @param response The response term in the model
#' @param focus The focus term for which influence is to be calculated.
#' @return A new Influence object
#'
Influence$new <- function(.,model,data=NULL,response=NULL,focus=NULL){
Influence$new <- function(., model, data = NULL, response = NULL, focus = NULL) {
instance = .$proto(
model = model,
data = data,
Expand All @@ -77,15 +78,15 @@ Influence$new <- function(.,model,data=NULL,response=NULL,focus=NULL){
#' @name Influence$init
#' @return None
#'
Influence$init <- function(.){
Influence$init <- function(.) {

# If no data was supplied....
if(is.null(.$data)){
if (is.null(.$data)) {
# See if the model has a data variable
.$data = .$model$data
if(is.null(.$data)){
if (is.null(.$data)) {
# See if it is available as a model attribute
.$data = attr(.$model,"data")
.$data = attr(.$model, "data")
}
if(is.null(.$data)){
stop("You need to provide the data that was fitted to for this type of model e. Influence$new(model,data,...)")
Expand Down Expand Up @@ -114,27 +115,30 @@ Influence$init <- function(.){
#' Extracts the coefficients
#'
#' @name coeffs
#' @param . The model to extract coefficients from
#' @param model The model to extract coefficients from
#' @param term The term in the model for which coefficients are extracted
#' @return A vector of coefficients
#' @importFrom stats coefficients
#' @export
#'
coeffs = function(.,model=.$model,term=.$focus){
coeffs = coefficients(model)
rows = substr(names(coeffs),1,nchar(term))==term
c(0,coeffs[rows])
#'
coeffs <- function(., model = .$model, term = .$focus) {
coeffs <- coefficients(model)
rows <- substr(names(coeffs), 1, nchar(term)) == term
c(0, coeffs[rows])
}

#' Obtain standard errors for the coefficients associated with a particular term in the model
#'
#' Extract standard errors using the method of Francis
#'
#' @name ses
#' @param . The model to extract standard errors for a coefficient from
#' @param model The model to extract standard errors for a coefficient from
#' @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') {
#The "cov.scaled" member of a glm object is the estimated covariance matrix of
Expand All @@ -147,10 +151,10 @@ ses = function(.,model=.$model, term=.$focus) {
V = model$var
}

rows = substr(row.names(V),1,nchar(term))==term
rows = substr(row.names(V), 1, nchar(term)) == term
V = V[rows,rows]
n = sum(rows)+1
Q = matrix(-1/n, nrow=n, ncol=n-1)
n = sum(rows) + 1
Q = matrix(-1/n, nrow = n, ncol = n - 1)
Q[-1,] = Q[-1,] + diag(rep(1,n-1))
V0 = (Q%*%V) %*% t(Q)
se = sqrt(diag(V0))
Expand All @@ -164,9 +168,9 @@ ses = function(.,model=.$model, term=.$focus) {
#' @param term The term in the model for which effects are extracted
#' @return A numeric vector of effects for the model term
#'
Influence$effects = function(.,model=.$model,term=.$focus){
Influence$effects = function(., model = .$model, term = .$focus){
coeffs = .$coeffs(model,term)
exp(coeffs-mean(coeffs))
exp(coeffs - mean(coeffs))
}

#' Perform necessary calculations
Expand Down
2 changes: 1 addition & 1 deletion R/plot-index.R
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ plot_index <- function(fit,
df$model <- factor(df$model, levels = c("Unstandardised", "Standardised"))

if (!show_unstandardised) {
df <- df %>% filter(model != "Unstandardised")
df <- df %>% filter(.data$model != "Unstandardised")
scale_col <- fill
scale_lin <- "solid"
} else {
Expand Down
Binary file modified data/lobsters_per_pot.rda
Binary file not shown.
2 changes: 2 additions & 0 deletions man/Influence-cash-new.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 2 additions & 0 deletions man/coeffs.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

3 changes: 0 additions & 3 deletions man/get_influ2.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

3 changes: 0 additions & 3 deletions man/get_unstandarsied.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

11 changes: 9 additions & 2 deletions man/lobsters_per_pot.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

15 changes: 11 additions & 4 deletions man/logistic.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/logit.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 2 additions & 0 deletions man/ses.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

3 changes: 0 additions & 3 deletions man/table_criterion.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

1 change: 0 additions & 1 deletion tests/testthat.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,2 @@
library(testthat)

test_check("influ2")
Binary file removed tests/testthat/Rplots.pdf
Binary file not shown.
Binary file modified tests/testthat/brm1.rds
Binary file not shown.
Binary file modified tests/testthat/brm2.rds
Binary file not shown.
17 changes: 9 additions & 8 deletions tests/testthat/test-cdi.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
library(influ2)
library(brms)

data(lobsters_per_pot)
Expand Down Expand Up @@ -87,23 +88,23 @@ test_that("this matches Nokome Bentley's influ package", {
filter(grepl("year", variable))

# Test Nokomes coeffs are the same as stats::coef function
plot(c1a, type = "b")
lines(c1b, col = 2)
# plot(c1a, type = "b")
# lines(c1b, col = 2)
expect_equal(c1a, c1b)

# Check that my get_coefs function is the same as the nlme::fixef function
plot(c2a$Estimate, type = "b")
lines(c2b$Estimate, col = 2)
# plot(c2a$Estimate, type = "b")
# lines(c2b$Estimate, col = 2)
expect_equal(c2a$Estimate, c2b$Estimate)

# Check that Nokomes coeffs are the same as mine. This is maximum likelihood vs Bayesian so tolerance needs to be a little higher.
plot(c1a, type = "b")
lines(c2a$Estimate, col = 2)
# plot(c1a, type = "b")
# lines(c2a$Estimate, col = 2)
expect_equal(as.numeric(c1a), c2a$Estimate, tolerance = 0.07)

# Check that stats::coef is the same as nlme::fixef. This is maximum likelihood vs Bayesian so tolerance needs to be a little higher.
plot(c1b, type = "b")
lines(c2b$Estimate, col = 2)
# plot(c1b, type = "b")
# lines(c2b$Estimate, col = 2)
expect_equal(as.numeric(c1b), c2b$Estimate, tolerance = 0.07)

# Check the influences are the same - I couldn't get this to work as Noko's code is broken
Expand Down
Loading

0 comments on commit 74dc975

Please sign in to comment.