R/LRIBI.R
In droglenc/fishWiDNR: Support for Wisconsin DNR Fisheries Analyses
#' @title Computes the Wisconsin large rivers IBI
#'
#' @description Computes the Wisconsin large rivers IBI as described in Lyons \emph{et al.} (2001). Also provides a variety of summary information and interpretative plots.
#'
#' @param fmdb A data.frame of raw fish data from the Fisheries Management Database. This data.frame should be created with \code{\link{read.FMDB}} and filtered to include only the survey or visit number for which the IBI should be computed. See Details.
#' @param agg A data.frame of aggregate counts and weights from XXXX. This data.frame should be filtered to include only the survey or visit number for which the IBI should be computed. See Details.
#' @param effort A single numeric value that represents the effort (in miles) expended on the survey or visit represented in \code{fmdb} and \code{agg}.
#' @param region A single character string that indicates the ecoregion (must be either \code{"North"} or \code{"South"}) in which the survey was conducted. See Notes.
#' @param excludeSuckers A single logical that indicates that the two \dQuote{sucker} metrics should be excluded or not (DEFAULT) from the IBI calculation. See Notes.
#' @param round A logical that indicates whether the metric values should be rounded before comparing to Table 3 to define ratings. See Notes.
#' @param object,x An \code{LRIBI} object returned from \code{LRIBI}.
#' @param intermediate A logical that indicates whether intermediate results should be returned with \code{summary}. See Details.
#' @param show.data A logical that indicates whether summary data results should be returned with \code{summary}. See Details.
#' @param show.title A logical that indicates whether the title and IBI result should be returned with \code{summary}. See Details.
#' @param type A single character string that indicates the type of plot to construct (in \code{plot}) or problems to report (in \code{problemsLRIBI}). See details.
#' @param y Not used. Included here because required by main \code{plot}.
#' @param label.mar A single numeric that indicates the size of y-axis margin to be used for all \code{plot}s except for \code{type="dashboard"}. This margin holds the metric names and, thus, likely needs to be fairly large. Likely only need to change from default when the overall graphic size changes.
#' @param axis.mult A single numeric used when \code{type="dashboard"} to extend the y-axes so that all values or metric cutoff values are shown. Generally, should only be slightly greater than 1 (and likely does not need to change from the default).
#' @param rating.col Either a single or vector of length three that contains colors use to represent \dQuote{"Poor"}, \dQuote{"Fair"}, or \dQuote{"Good"} metric ratings. Use a vector of length three for different colors for the ratings or a single color to not distinguish among the ratings. See Examples.
#' @param xlab A single numeric used to label the x-axis of the plots.
#' @param main A single numeric used to produce a main title for the plots. If \code{NULL} (the DEFAULT) a title that includes the survey number, waterbody name, county, ecoregion, and year is produced. Use \code{main=""} to suppress the title.
#' @param \dots Optional additional arguments used in \code{plot}. Not yet implemented in \code{summary} or \code{coef} (needed here to match default function).
#'
#' @details This function derives the Wisconsin large rivers index of biotic integrity (IBI) value using the calculations described in Lyons \emph{et al.} (2001). This function requires data collected in a specific Wisconsin DNR survey and stored in two CSV files reported from WiDNR databases. The first is a "raw fish data" file from the Fisheries Management Database (hereafter, called the FMDB file) that is obtained from XXX FMDB DETAILS HERE XXX and can be read by \code{\link{read.FMDB}}. The second file is an "aggregate data" file (hereafter, called the aggregate file) that is obtained from XXX AGGREGATE FILES DETAILS HERE XXX. After reading both files should be filtered to include only data from the specific survey (could be a specific \dQuote{Survey.Seq.No} or \dQuote{Visit.Fish.Seq.No}) for which the IBI should be calculated. The main \code{LRIBI} function will automatically restrict the FMDB data.frame to only those surveys where the \dQuote{TARGET SPECIES} field is equal to \dQuote{ALL SPECIES}.
#'
#' Some discrepancies may exist between the FMDB and aggregate files. For example, the total counts or biomass of a species may differ between the files. Warning messages will be issued when these discrepancies are discovered. The main \code{LRIBI} function will continue with the counts from the FMDB file and the weights from the aggregate file when discrepancies occur. Help with determining the species for which the discrepancies occur is obtained by submitting the \code{LRIBI} object to \code{problemsLRIBI}.
#'
#' The overall IBI is returned by either \code{coef} (or \code{summary} with \code{show.metrics=FALSE}, \code{show.data=FALSE}, and \code{show.title=FALSE}). The intermediate metrics may be seen by including \code{show.metrics=TRUE} (DEFAULT) in \code{summary} or with \code{type="dashboard"}, \code{type="ratings"}, or \code{type="score"} in \code{plot}. The summarized data that led to the metrics may be seen by including \code{show.data=TRUE} (DEFAULT) in \code{summary} or with \code{type="counts"} or \code{type="weights"} in \code{plot}.
#'
#' @return The main \code{LRIBI} functions returns a list of class LRIBI with the following items:
#' \itemize{
#' \item survey The survey number given in \code{survey}.
#' \item region The ecoregion in which the survey was conducted (as determined from the county or as sent in \code{region}).
#' \item effort The number of miles samples (currently as given in \code{effort}).
#' \item short.title A character string title that includes the survey number, visit number, waterbody name, and year. Used with \code{plot}.
#' \item long.title A character string title that includes the survey number, visit number, waterbody name, county, ecoregion, and year. Used with \code{summary}
#' \item IBI The overall IBI score.
#' \item n.fish The total number (count) of fish sampled.
#' \item ttl.Biomass The total biomass (kg) of fish sampled.
#' \item sum.data A data.frame of species counts, biomasses, and DELT number, along with categories required for computing the IBI metrics.
#' \item IBI.data A data.frame of metric name, value, rating, and score for the 10 individual metrics in Lyons \emph{et al.} (2001) IBI.
#' \item noProblems A logical that is \code{TRUE} if there were no discrepancies between the FMDB and aggregate data files.
#' \item countsProblems A vector of species names for which there was a discrepancy in the counts between the FMDB and aggregate data files.
#' \item weightsProblems A vector of species names for which there was a discrepancy in the weights between the FMDB and aggregate data files.
#' \item fmdb.data A data.frame that contains the species name, total count, and total biomass for species in the FMDB data file.
#' \item agg.data A data.frame that contains the species name, total count, total biomass, and total DELT number for species in the aggregate data file.
#' }
#'
#' \code{summary} prints a variety of items (see Details) and returns the overall IBI value.
#'
#' \code{coef} returns the overall IBI value.
#'
#' \code{plot} produces a number of plots but does not return any value.
#'
#' \code{problemsLRIBI} prints a variety of items (see Details) but does not return any value.
#'
#' @note It is expected that the FMDB and aggregate files (see Details) will include the survey number in the file name and that number will be at the end of the file name (with the exception of the CSV extension). The filenames may have other text prior to the survey number. It is suggested to not change the filenames from those returned by the WiDNR database reports (e.g., \dQuote{fish_raw_data_Survey375049456.csv} and \dQuote{aggdata_Survey375049456.csv}). However, if the text prior to the survey number is changed from the default then that text should be entered into \code{fmdb.prefix} and \code{agg.prefix} as appropriate.
#'
#' The ecoregion can be automatically detected from the county in which the survey took place. However, there are a handful of counties that exist in both the \code{North} and \code{South} ecoregions. In these instances, the automatic detection of ecoregion will throw an error asking the user to supply the specific ecoregion in \code{region}. Alternatively the user may over-ride the automatically determined ecoregion by supplying a specific ecoregion in \code{region}.
#'
#' The values in Table 3 of Lyons \emph{et al.} (2001) imply that most metric values were rounded to whole numbers (\dQuote{WPUE} and \dQuote{\%DELT} were rounded to one decimal) for converting to IBI ratings. Metric values are rounded prior to comparing to Table 3 if \code{round=TRUE} (the DEFAULT) or are not rounded if \code{round=FALSE}.
#'
#' Lyons \emph{et al.} (2001) noted that the diversity of sucker species in tributaries to Lake Superior and some tributaries to Green Bay, Lake Michigan had very low diversities of sucker species. Thus, the two sucker metrics (\dQuote{Sucker species} and \dQuote{\% Round suckers (wt)}) should be excluded from the overall IBI calculation. The sum of the other eight metrics should be multipled by 1.25 (=10/8) to compensate for the use of fewer metrics and produce a result on a similar scale. By default the two sucker metrics are used in the overall IBI calculation, but the two metrics can be excluded (and the corrective multiplier applied) by using \code{excludeSuckers=TRUE}.
#'
#' @seealso \code{\link{LRIBI_AppendixA}}, \code{\link{LRIBI_Table3}}, and \code{\link{LRIBI_Ecoregions}}.
#'
#' @source Lyons, J., R.R. Piette, and K.W. Niermeyer. 2001. Development, validation, and application of a fish-based index of biotic integrity for Wisconsin's large warmwater rivers. Transactions of the American Fisheries Society 130:1077-1094.
#'
#' @author Derek H. Ogle (\email{derek@@derekogle.com}) and Joseph T. Mrnak
#'
#' @keywords manip plot
#'
#' @aliases LRIBI summary.LRIBI plot.LRIBI problemsLRIBI getLRIBISurveys
#'
#' @examples
#' ## Get location of external data files
#' dir <- dirname(system.file("extdata","aggdata_Survey375049456.csv",package="fishWiDNR"))
#'
#' ## Example of one survey
#' fmdb1 <- read.FMDB(paste0(dir,"/fish_raw_data_Survey375049456.csv"),expandCounts=TRUE)
#' agg1 <- read.csv(paste0(dir,"/aggdata_Survey375049456.csv"),stringsAsFactors=FALSE)
#' res <- LRIBI(fmdb1,agg1,effort=0.9)
#'
#' coef(res)
#' summary(res)
#' summary(res,show.data=FALSE)
#' summary(res,intermediate=FALSE)
#' summary(res,show.data=FALSE,intermediate=FALSE)
#' summary(res,show.data=FALSE,intermediate=FALSE,show.title=FALSE)
#'
#' plot(res)
#' plot(res,rating.col=NULL)
#' plot(res,rating.col="lightblue")
#' plot(res,rating.col=c("lightcyan","lightsalmon","lightblue"))
#' plot(res,type="counts",xlim=c(0,80))
#' plot(res,type="biomass",xlim=c(0,25))
#' plot(res,type="weights",xlim=c(0,25)) # same
#' plot(res,type="ratings")
#' plot(res,type="ratings",rating.col=NULL)
#' plot(res,type="ratings",rating.col="lightblue")
#' plot(res,type="ratings",rating.col=c("lightcyan","lightsalmon","lightblue"))
#'
#' problemsLRIBI(res)
#'
#' ## Another example of one survey
#' fmdb2 <- read.FMDB(paste0(dir,"/fish_raw_data_Survey375334018.csv"),expandCounts=TRUE)
#' agg2 <- read.csv(paste0(dir,"/aggdata_Survey375334018.csv"),stringsAsFactors=FALSE)
#'
#' res2 <- LRIBI(fmdb2,agg2,effort=0.9)
#' summary(res2)
#' plot(res2)
#' problemsLRIBI(res2)
#' problemsLRIBI(res2,type="count")
#' problemsLRIBI(res2,type="weight")
#'
#' # uses not rounded metric values
#' res2a <- LRIBI(fmdb2,agg2,effort=0.9,round=FALSE)
#' coef(res2); coef(res2a) # rounding can make a difference
#'
#' @rdname LRIBI
#' @export
LRIBI <- function(fmdb,agg,effort,
region=NULL,round=TRUE,excludeSuckers=FALSE){
## Load the internal lookup tables into this function's environment
## the data/get combination avoids the "no global binding" note at CHECK
appA <- get(utils::data("LRIBI_AppendixA",envir=environment()),envir=environment())
tab3 <- get(utils::data("LRIBI_Table3",envir=environment()),envir=environment())
ecoregs <- get(utils::data("LRIBI_Ecoregions",envir=environment()),envir=environment())
## Load FMDB and aggregate files (retruned are ready to be summarized)
fmdb <- iLRIBI_hndlFMDB(fmdb,appA)
agg <- iLRIBI_hndlAGG(agg,appA)
## Make sure efforts are reasonable
iCheckEffort(effort)
## Find county (make sure one county, capitalized properly, and not a factor)
cnty <- as.character(capFirst(unique(fmdb$County)))
## Find or handle the ecoregion
region <- iHndlRegion(region,cnty,ecoregs)
## Find the survey
survey <- iHndlSurvey(fmdb)
## Find the fish visit sequence number
visit <- iHndlVisit(fmdb)
## Make a title for the outputs
titles <- iMakeTitle(survey,visit,cnty,region,fmdb)
## Summarize and join two data types
joined <- iLRIBI_SumJoin(fmdb,agg,survey)
# Add grouping variables to the sum.data
joined$sum.data <- iLRIBI_AddGroupVars(joined$sum.data,appA)
## Compute individual metrics and lookup IBI scores
IBI.data <- iLRIBI_ComputeIBIMetrics(joined$sum.data,tab3,effort,
region,round,excludeSuckers)
## Return list of results in an LRIBI class
res <- list(survey=survey,region=region,effort=effort,
short.title=titles$short.title,
long.title=titles$long.title,
IBI=ifelse(excludeSuckers,sum(IBI.data$Score)*1.25,
sum(IBI.data$Score)),
n.fish=sum(joined$sum.data$Count),
ttl.Biomass=sum(joined$sum.data$Biomass),
sum.data=joined$sum.data,
IBI.data=IBI.data,
noProblems=is.null(c(joined$countsProblems,joined$weightsProblems)),
countsProblems=joined$countsProblems,
weightsProblems=joined$weightsProblems,
fmdb.data=joined$fmdb.data,
agg.data=joined$agg.data)
class(res) <- "LRIBI"
res
}
#' @rdname LRIBI
#' @export
summary.LRIBI <- function(object,intermediate=TRUE,show.data=TRUE,show.title=TRUE,...) {
if (show.title) {
cat(object$long.title,"\n")
cat("\nIBI =",object$IBI,"\n")
}
if (intermediate) {
if (show.title) cat("\n")
cat("IBI computed from these metric values, ratings and scores:\n")
print(object$IBI.data,row.names=FALSE)
}
if (show.data) {
if (show.title | intermediate) cat("\n")
cat("These observed data were used:\n")
print(object$sum.data[,c("Species1","Count","Biomass","DELTnum")],row.names=FALSE)
}
if (show.title | intermediate | show.data) invisible(object$IBI)
else return(object$IBI)
}
#' @rdname LRIBI
#' @export
coef.LRIBI <- function(object,...) object$IBI
#' @rdname LRIBI
#' @export
plot.LRIBI <- function(x,type=c("dashboard","counts","weights","biomass","ratings","scores"),
y=NULL,label.mar=9,axis.mult=1.25,
rating.col=FSA::col2rgbt(c("red","yellow","green"),1/2),
xlab=NULL,main=NULL,...) {
type <- match.arg(type)
switch(type,
dashboard= { iLRIBI_Dashboard(x,axis.mult,rating.col,main,...) },
counts= { iLRIBI_barData(x,type,label.mar,xlab,main,...) },
weights=,biomass= { iLRIBI_barData(x,type,label.mar,xlab,main,...) },
ratings=,scores= { iLRIBI_barMetrics(x,label.mar,rating.col,xlab,main,...) } )
}
#' @rdname LRIBI
#' @export
problemsLRIBI <- function(object,type=c("counts","weights")) {
type <- match.arg(type,several.ok=TRUE)
if (object$noProblems) message("No discrepancies between the FMDB and aggregate files\n were encountered for ",object$short.title)
else {
## Get fmdb and aggregate data together
tmp <- dplyr::left_join(object$fmdb.data,object$agg.data,by="Species1")[,-6]
names(tmp) <- sub(".x",".fmdb",names(tmp))
names(tmp) <- sub(".y",".agg",names(tmp))
## Address problems with counts
if ("counts" %in% type) {
if (is.null(object$countsProblems)) {
# don't send message if weights also in type
if (!"weights" %in% type) message("There was no discrepancies between the counts in the FMDB and\n aggregate files for ",object$short.title)
} else {
cat(object$long.title,"\n\n",sep="")
# create new variable that will have YES! to note discrepancy
tmp$PROBLEM <- ""
tmp$PROBLEM[tmp$Species1 %in% object$countsProblems] <- "YES!"
print(tmp[,-which(names(tmp) %in% c("Biomass.fmdb","Biomass.agg"))],
row.names=FALSE)
}
}
if ("weights" %in% type) {
if (is.null(object$weightsProblems)) {
# don't send message if counts also in type
if (!"counts" %in% type) message("There was no discrepancies between the weights in the FMDB and\n aggregate files for ",object$short.title)
} else {
cat(object$long.title,"\n\n",sep="")
# create new variable that will have YES! to note discrepancy
tmp$PROBLEM <- ""
tmp$PROBLEM[tmp$Species1 %in% object$weightsProblems] <- "YES!"
print(tmp[,-which(names(tmp) %in% c("Count.fmdb","Count.agg"))],
row.names=FALSE)
}
}
}
}
################################################################################
################################################################################
##### INTERNAL (not exported) COMPUTATION FUNCTIONS
##### Load FMDB data
iLRIBI_hndlFMDB <- function(fmdb,appA) {
if (!is.data.frame(fmdb)) stop("'fmdb' must be a data.frame from the FMDB.",
call.=FALSE)
# Reduce to only "all species" (Z100) surveys
fmdb <- fmdb[fmdb$Target.Species=="ALL SPECIES",,drop=TRUE]
# Convert g to kg
fmdb$Weight.KG <- fmdb$Weight.Grams/1000
# Make species a character variable (rather than factor)
fmdb$Species1 <- as.character(fmdb$Species1)
# Remove species not in Appendix A
fmdb <- fmdb[fmdb$Species1 %in% appA$common.name,,drop=TRUE]
# Group by species to make read for summarizing
fmdb <- dplyr::group_by_(fmdb,~Species1)
fmdb
}
##### Load aggregate data
iLRIBI_hndlAGG <- function(agg,appA) {
if (!is.data.frame(agg)) stop("'agg' must be an aggregate file data.frame.",
call.=FALSE)
# Convert species capitalization to match FMDB file (make sure it is character)
agg$Species1 <- as.character(FSA::capFirst(agg$Common.Name))
# Convert grams or lbs to kg
agg$Aggregate.KG <- ifelse(agg$Units=="GRAMS",agg$Aggregate.Weight/1000,
agg$Aggregate.Weight*454/1000)
# Remove species not in Appendix A
agg <- agg[agg$Species1 %in% appA$common.name,,drop=TRUE]
# Group by species to make read for summarizing
agg <- dplyr::group_by_(agg,~Species1)
agg
}
##### Join the FMDB and Aggregate files (and handle potential problems)
iLRIBI_SumJoin <- function(fmdb,agg,survey) {
## Initialize some variables
countsProblems <- weightsProblems <- NULL
## Count number of fish and sum weights in FMDB file
## dots needed to use with summarize_ for using dplyr in a function
dots <- list(~length(Species1),~sum(Weight.KG))
fmdb1 <- dplyr::summarize_(fmdb,.dots=stats::setNames(dots,
c("Count","Biomass")))
## Sum number of fish, aggregate weight, and number of DELT fish
## in aggregate file
dots <- list(~sum(Total.Count),~sum(Aggregate.KG),~sum(DELT..))
agg1 <- dplyr::summarize_(agg,.dots=stats::setNames(dots,
c("Count","Biomass","DELTnum")))
## Join the two summary date.frames by Species1
sum.data <- dplyr::full_join(fmdb1,agg1,by="Species1")
## Perform checks on counts
# problem if species in aggregate but not FMDB file
if (any(is.na(sum.data$Count.x))) {
warning("Some species in survey #",survey," had counts in the aggregate but not FMDB file.",call.=FALSE)
countsProblems <- c(countsProblems,
sum.data$Species1[is.na(sum.data$Count.x)])
}
# problem if counts between the two files do not match
tmp <- sum.data[!is.na(sum.data$Count.y),]
if (any(tmp$Count.x!=tmp$Count.y)) {
warning("Counts in FMDB and aggregate files for survey #",survey," do not match\n for some species! Continued with FMDB counts.",call.=FALSE)
countsProblems <- c(countsProblems,tmp$Species1[tmp$Count.x!=tmp$Count.y])
}
## Perform checks on weights
# problem if there are weights in both files and do they not match
tmp <- sum.data[!is.na(sum.data$Biomass.x) & !is.na(sum.data$Biomass.y),]
if (nrow(tmp)>0) {
if (nrow(tmp)==nrow(sum.data)) stop("There are no weights in either the FMDB or aggregate files.",call.=FALSE)
if (any(tmp$Biomass.x!=tmp$Biomass.y)) warning("Weights that exist in both FMDB and aggregate files for survey #",survey," do not match!\n Continued with aggregage weight.",call.=FALSE)
weightsProblems <- tmp$Species1[tmp$Biomass.x!=tmp$Biomass.y]
}
## Replace NA biomasses from FMDB with biomasses from aggregate file
tmp <- sum.data[is.na(sum.data$Biomass.x),]
tmp <- which(sum.data$Species1 %in% tmp$Species1)
sum.data$Biomass.x[tmp] <- sum.data$Biomass.y[tmp]
## Replace NAs in DELTnum with zeroes
sum.data$DELTnum[is.na(sum.data$DELTnum)] <- 0
## Rename the ".x" variables, Delete the ".y" variables
names(sum.data)[names(sum.data) %in% c("Count.x","Biomass.x")] <- c("Count","Biomass")
sum.data <- sum.data[,-which(names(sum.data) %in% c("Count.x","Count.y","Biomass.x","Biomass.y"))]
## Check if any weights are missing
if (any(is.na(sum.data$Biomass))) {
warning("Some species had counts but no weights recorded for survey #",survey,".",call.=FALSE)
weightsProblems <- unique(c(weightsProblems,
sum.data$Species1[is.na(sum.data$Biomass)]))
}
## Return result (use as.data.frame to strip off tibble attributes)
list(sum.data=as.data.frame(sum.data),
countsProblems=countsProblems,
weightsProblems=weightsProblems,
fmdb.data=as.data.frame(fmdb1),
agg.data=as.data.frame(agg1))
}
##### Add grouping variables required to create IBI metrics
iLRIBI_AddGroupVars <- function(sum.data,appA) {
## use mapvalues to lookup grouping variable levels by species common name
## used warn_missing=FALSE to suppress warnings when level did not exist
sum.data$habitat <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$habitat,
warn_missing=FALSE)
sum.data$habitat <- FSA::mapvalues(sum.data$habitat,
from="River-Large",to="River",warn_missing=FALSE)
sum.data$suckers <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$common.family,
warn_missing=FALSE)
sum.data$tolerance <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$tolerance,
warn_missing=FALSE)
sum.data$origin <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$origin,
warn_missing=FALSE)
sum.data$spawning <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$spawning,
warn_missing=FALSE)
sum.data$feeding <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$feeding,
warn_missing=FALSE)
sum.data$round.sucker <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$round.sucker,
warn_missing=FALSE)
sum.data
}
##### Compute the individual metrics
iLRIBI_ComputeIBIMetrics <- function(sum.data,tab3,effort,
region,round,excludeSuckers) {
## Individual metric summaries
metvals <- iLRIBI_SumValues(sum.data,effort)
## Individual metric ratings (Good, Fair, Poor)
metrats <- iLRIBI_Convert2Ratings(metvals,tab3,region,round)
## Create data.frame with metric names, values, and ratings
tmp <- data.frame(Metric=c("WPUE","# Native species","# Sucker species",
"# Intolerant species","# Riverine species","% DELT (n)",
"% Riverine (n)","% Lithophils (n)",
"% Insectivore (wt)","% Round suckers (wt)"),
Value=metvals,
Rating=metrats,stringsAsFactors=FALSE)
## Add new variable that has metric score (10,5,0)
tmp$Score <- FSA::fact2num(FSA::mapvalues(tmp$Rating,
from=c("Poor","Fair","Good"),to=c(0,5,10),
warn_missing=FALSE))
row.names(tmp) <- NULL
if (excludeSuckers) tmp <- tmp[!grepl("ucker",tmp$Metric),]
tmp
}
##### Extracts specific metrics for individual IBI metrics
iLRIBI_SumValues <- function(sum.data,effort) {
## Short helper function to make calculations above look cleaner
iPercTable <- function(formula,data,FUN) {
suppressWarnings(100*prop.table(FSA::sumTable(formula,data,FUN)))
}
## must exclude Tolerant species for WPUE calculation
c(WPUE=sum(sum.data$Biomass[sum.data$tolerance!="Tolerant"])/effort,
num.native=stats::xtabs(~origin,data=sum.data)[["Native"]],
num.sucker=stats::xtabs(~suckers,data=sum.data)[["Suckers"]],
num.intolerant=stats::xtabs(~tolerance,data=sum.data)[["Intolerant"]],
num.riverine=stats::xtabs(~habitat,data=sum.data)[["River"]],
perc.DELT=100*sum(sum.data$DELTNum)/sum(sum.data$Count),
perc.riverine=iPercTable(Count~habitat,data=sum.data,FUN=sum)[["River"]],
perc.lithophils=iPercTable(Count~spawning,data=sum.data,FUN=sum)[["Lithophil"]],
perc.insectivore=iPercTable(Biomass~feeding,data=sum.data,FUN=sum)[["Insectivore"]],
perc.rsuckers=iPercTable(Biomass~round.sucker,data=sum.data,FUN=sum)[["Yes"]])
}
##### Converts metric values into Poor, Fair, Good ratings
iLRIBI_Convert2Ratings <- function(metvals,tab3,region,round) {
## Internal function to make below cleaner
iTab3Lookup <- function(var,metvals,tab3,round,digits=0,less=TRUE) {
val <- metvals[[var]]
if (round) val <- round(val,digits)
if (less) {
ifelse(val<tab3$fair[tab3$metric==var],"Poor",
ifelse(val<=tab3$good[tab3$metric==var],"Fair","Good"))
} else {
ifelse(val>tab3$fair[tab3$metric==var],"Poor",
ifelse(val>=tab3$good[tab3$metric==var],"Fair","Good"))
}
} # end internal
## Lookup ratings for the specific ecoregion
tab3 <- tab3[tab3$ecoregion==region,,drop=TRUE]
c(WPUE_rate=iTab3Lookup("WPUE",metvals,tab3,round,digits=1),
num.native_rate=iTab3Lookup("num.native",metvals,tab3,round),
num.sucker_rate=iTab3Lookup("num.sucker",metvals,tab3,round),
num.intolerant_rate=iTab3Lookup("num.intolerant",metvals,tab3,round),
num.riverine_rate=iTab3Lookup("num.riverine",metvals,tab3,round),
perc.DELT_rate=iTab3Lookup("perc.DELT",metvals,tab3,round,digits=1,less=FALSE),
perc.riverine_rate=iTab3Lookup("perc.riverine",metvals,tab3,round),
perc.lithophils_rate=iTab3Lookup("perc.lithophils",metvals,tab3,round),
perc.insectivore_rate=iTab3Lookup("perc.insectivore",metvals,tab3,round),
perc.rsuckers_rate=iTab3Lookup("perc.rsuckers",metvals,tab3,round))
}
################################################################################
################################################################################
##### INTERNAL (not exported) PLOTTING FUNCTIONS
##### PLOT of data (WPUE, etc.)
iLRIBI_barData <- function(object,type,label.mar,xlab,main,...) {
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle1(object,main)
if (type=="counts") {
tmp <- object$sum.data$Count
if (is.null(xlab)) xlab="Observed Number"
} else {
tmp <- object$sum.data$Biomass
if (is.null(xlab)) xlab="Observed Weight (KG)"
}
graphics::par(mar=c(3,label.mar,0,1),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
bp <- graphics::barplot(tmp,names.arg=object$sum.data$Species1,
horiz=TRUE,xlab=xlab,...)
graphics::axis(2,at=bp,labels=NA)
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
}
##### Plot of metric ratings (Poor, Fair, Good)
iLRIBI_barMetrics <- function(object,label.mar,clrs,xlab,main,...) {
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle1(object,main)
clrs <- iHndlRatingColors(clrs)[factor(object$IBI.data$Rating,
levels=c("Poor","Fair","Good"))]
if (is.null(xlab)) xlab="Metric Rating (Score)"
graphics::par(mar=c(3,label.mar,0,2),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
bp <- graphics::barplot(object$IBI.data$Score,names.arg=object$IBI.data$Metric,
horiz=TRUE,xlab=xlab,xaxt="n",col=clrs,...)
graphics::axis(2,at=bp,labels=NA)
graphics::axis(1,at=c(0,5,10),labels=c("Poor (0)","Fair (5)","Good (10)"))
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
}
##### Dashboard plot of each metric value vs. its cutoff scores
iLRIBI_Dashboard <- function(object,axis.mult,clrs,main) {
## Internal function for individual bar plots
iLRIBI_plotBarIndiv <- function(object,tab3,axis.mult,metnum,clrs) {
df <- object$IBI.data
val <- df$Value[metnum]
# find the metric cutoffs
cuts <- tab3[tab3$ecoregion==object$region,,drop=TRUE]
cuts <- cuts[metnum,,drop=TRUE]
# adjust y-axis limits by type of metric
if (cuts$metric=="WPUE") {
# for WPUE be more than observed value of maximum cutoff
ymax <- max(axis.mult*val,axis.mult*max(cuts$good))
} else if (grepl("num",cuts$metric)) {
# for counts be more than observed, maximum cutoff+1, or 5 (helps with axis ticks)
ymax <- max(axis.mult*val,axis.mult*max(cuts$good+1),5)
} else {
# for percents be 100%
ymax <- 100
}
clrs <- clrs[factor(df$Rating[metnum],levels=c("Poor","Fair","Good"))]
op <- graphics::par(mar=c(1.2,3,1.2,0.7),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
graphics::barplot(val,ylab=df$Metric[metnum],ylim=c(0,ymax),col=clrs)
graphics::mtext(df$Rating[metnum],line=-0.5,cex=0.8)
graphics::abline(h=c(cuts$fair,cuts$good),col="gray50",lwd=3,lty=2)
graphics::par(op)
} # end internal
# get colors
clrs <- iHndlRatingColors(clrs)
# Make a layout for 9 plots (exclues perc.DELT ... see below)
tab3 <- get(utils::data("LRIBI_Table3",envir=environment()),envir=environment())
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle2(object,main)
# Plots individual metrics
iLRIBI_plotBarIndiv(object,tab3,axis.mult,1,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,2,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,3,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,4,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,5,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,7,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,8,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,9,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,10,clrs)
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
# did not plot DELT because nearly always zero???
# this sends a message reminding that it was not plotted.
message("The '% DELT (n)' metric of ",object$IBI.data$Rating[6],
" (",round(object$IBI.data$Value[6],2),"%) is not shown.")
}
##### Handle putting the main title on the plot for barplots
iHndlMainTitle1 <- function(object,main) {
if (is.null(main)) main <- object$short.title
if (main!="") {
graphics::layout(matrix(1:2,nrow=2),heights=c(0.5,9))
graphics::par(mar=c(0,0,0,0))
graphics::plot.new()
graphics::text(0.5,0.5,main,cex=0.85)
}
}
##### Handle putting the main title on the plot for dashboard
iHndlMainTitle2 <- function(object,main) {
if (is.null(main)) main <- object$short.title
if (main=="") graphics::layout(matrix(1:9,nrow=3,byrow=TRUE))
else {
graphics::layout(matrix(c(rep(1,3),2:10),nrow=4,byrow=TRUE),
heights=c(0.5,3,3,3))
graphics::par(mar=c(0,0,0,0))
graphics::plot.new()
graphics::text(0.5,0.5,main,cex=1.15)
}
}
##### Handle dealing with the colors used for the different ratings
iHndlRatingColors <- function(clrs) {
if (is.null(clrs)) clrs <- "gray70"
if (length(clrs)!=3) {
if (length(clrs)==1) clrs <- rep(clrs,3)
else stop("Must supply 1 or 3 colors to 'rating.col'.",call.=FALSE)
}
clrs
}
##### Make sure that effort is reasonable
iCheckEffort <- function(effort) {
if (length(effort)!=1L) stop("Only one 'effort' value can be used.",call.=FALSE)
if (effort<=0) stop("'effort' must be positive.",call.=FALSE)
if (effort>2) warning("Check that your 'effort' is in miles.",call.=FALSE)
}
##### Check or handle finding the ecoregion
iHndlRegion <- function(region,cnty,ecoregs) {
if (!is.null(region)) {
## user-supplied ecoregion ... make sure it is North or South
if (!region %in% c("North","South"))
stop("'region' must be 'North' or 'South'.",call.=FALSE)
} else {
## ecoregion not supplied by user ... find it from county
if (is.null(region)) region <- ecoregs$ecoregion[which(ecoregs$county==cnty)]
if (region=="Both") stop("This survey is in a county that is in both ecoregions.\n Please use 'region' to specify which ecoregions should be used.",call.=FALSE)
}
region
}
iHndlSurvey <- function(fmdb) {
survey <- unique(fmdb$Survey.Seq.No)
if (length(survey)>1L) stop("The IBI can only be computed for one survey number.",
call.=FALSE)
survey
}
iHndlVisit <- function(fmdb) {
visit <- unique(fmdb$Visit.Fish.Seq.No)
if (length(visit)>1L) warning("More than one 'Visit.Fish.Seq.No' is present in the 'fmdb' file.",call.=FALSE)
visit
}
#### Make a title for outputs
iMakeTitle <- function(survey,visit,cnty,region,fmdb) {
## both titles get survey and visit numbers
tmp <- paste0("Survey #",survey," (")
if (length(visit)>1L) tmp <- paste0(tmp,"multiple visits)")
else tmp <- paste0(tmp,"Visit #",visit,")")
## short title also gets watebordy name and survey year
stmp <- paste0(tmp," - ",as.character(capFirst(unique(fmdb$Waterbody.Name))),
", ",unique(fmdb$Survey.Year))
## long title also gets waterbody name, county, ecoregion, and survey year
ltmp <- paste0(tmp,"\n",as.character(capFirst(unique(fmdb$Waterbody.Name))),
", ",cnty," County, ",region," ecoregion, ",
unique(fmdb$Survey.Year))
list(short.title=stmp,long.title=ltmp)
}
droglenc/fishWiDNR documentation built on May 15, 2019, 2:51 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' @title Computes the Wisconsin large rivers IBI
#'
#' @description Computes the Wisconsin large rivers IBI as described in Lyons \emph{et al.} (2001). Also provides a variety of summary information and interpretative plots.
#'
#' @param fmdb A data.frame of raw fish data from the Fisheries Management Database. This data.frame should be created with \code{\link{read.FMDB}} and filtered to include only the survey or visit number for which the IBI should be computed. See Details.
#' @param agg A data.frame of aggregate counts and weights from XXXX. This data.frame should be filtered to include only the survey or visit number for which the IBI should be computed. See Details.
#' @param effort A single numeric value that represents the effort (in miles) expended on the survey or visit represented in \code{fmdb} and \code{agg}.
#' @param region A single character string that indicates the ecoregion (must be either \code{"North"} or \code{"South"}) in which the survey was conducted. See Notes.
#' @param excludeSuckers A single logical that indicates that the two \dQuote{sucker} metrics should be excluded or not (DEFAULT) from the IBI calculation. See Notes.
#' @param round A logical that indicates whether the metric values should be rounded before comparing to Table 3 to define ratings. See Notes.
#' @param object,x An \code{LRIBI} object returned from \code{LRIBI}.
#' @param intermediate A logical that indicates whether intermediate results should be returned with \code{summary}. See Details.
#' @param show.data A logical that indicates whether summary data results should be returned with \code{summary}. See Details.
#' @param show.title A logical that indicates whether the title and IBI result should be returned with \code{summary}. See Details.
#' @param type A single character string that indicates the type of plot to construct (in \code{plot}) or problems to report (in \code{problemsLRIBI}). See details.
#' @param y Not used. Included here because required by main \code{plot}.
#' @param label.mar A single numeric that indicates the size of y-axis margin to be used for all \code{plot}s except for \code{type="dashboard"}. This margin holds the metric names and, thus, likely needs to be fairly large. Likely only need to change from default when the overall graphic size changes.
#' @param axis.mult A single numeric used when \code{type="dashboard"} to extend the y-axes so that all values or metric cutoff values are shown. Generally, should only be slightly greater than 1 (and likely does not need to change from the default).
#' @param rating.col Either a single or vector of length three that contains colors use to represent \dQuote{"Poor"}, \dQuote{"Fair"}, or \dQuote{"Good"} metric ratings. Use a vector of length three for different colors for the ratings or a single color to not distinguish among the ratings. See Examples.
#' @param xlab A single numeric used to label the x-axis of the plots.
#' @param main A single numeric used to produce a main title for the plots. If \code{NULL} (the DEFAULT) a title that includes the survey number, waterbody name, county, ecoregion, and year is produced. Use \code{main=""} to suppress the title.
#' @param \dots Optional additional arguments used in \code{plot}. Not yet implemented in \code{summary} or \code{coef} (needed here to match default function).
#'
#' @details This function derives the Wisconsin large rivers index of biotic integrity (IBI) value using the calculations described in Lyons \emph{et al.} (2001). This function requires data collected in a specific Wisconsin DNR survey and stored in two CSV files reported from WiDNR databases. The first is a "raw fish data" file from the Fisheries Management Database (hereafter, called the FMDB file) that is obtained from XXX FMDB DETAILS HERE XXX and can be read by \code{\link{read.FMDB}}. The second file is an "aggregate data" file (hereafter, called the aggregate file) that is obtained from XXX AGGREGATE FILES DETAILS HERE XXX. After reading both files should be filtered to include only data from the specific survey (could be a specific \dQuote{Survey.Seq.No} or \dQuote{Visit.Fish.Seq.No}) for which the IBI should be calculated. The main \code{LRIBI} function will automatically restrict the FMDB data.frame to only those surveys where the \dQuote{TARGET SPECIES} field is equal to \dQuote{ALL SPECIES}.
#'
#' Some discrepancies may exist between the FMDB and aggregate files. For example, the total counts or biomass of a species may differ between the files. Warning messages will be issued when these discrepancies are discovered. The main \code{LRIBI} function will continue with the counts from the FMDB file and the weights from the aggregate file when discrepancies occur. Help with determining the species for which the discrepancies occur is obtained by submitting the \code{LRIBI} object to \code{problemsLRIBI}.
#'
#' The overall IBI is returned by either \code{coef} (or \code{summary} with \code{show.metrics=FALSE}, \code{show.data=FALSE}, and \code{show.title=FALSE}). The intermediate metrics may be seen by including \code{show.metrics=TRUE} (DEFAULT) in \code{summary} or with \code{type="dashboard"}, \code{type="ratings"}, or \code{type="score"} in \code{plot}. The summarized data that led to the metrics may be seen by including \code{show.data=TRUE} (DEFAULT) in \code{summary} or with \code{type="counts"} or \code{type="weights"} in \code{plot}.
#'
#' @return The main \code{LRIBI} functions returns a list of class LRIBI with the following items:
#' \itemize{
#' \item survey The survey number given in \code{survey}.
#' \item region The ecoregion in which the survey was conducted (as determined from the county or as sent in \code{region}).
#' \item effort The number of miles samples (currently as given in \code{effort}).
#' \item short.title A character string title that includes the survey number, visit number, waterbody name, and year. Used with \code{plot}.
#' \item long.title A character string title that includes the survey number, visit number, waterbody name, county, ecoregion, and year. Used with \code{summary}
#' \item IBI The overall IBI score.
#' \item n.fish The total number (count) of fish sampled.
#' \item ttl.Biomass The total biomass (kg) of fish sampled.
#' \item sum.data A data.frame of species counts, biomasses, and DELT number, along with categories required for computing the IBI metrics.
#' \item IBI.data A data.frame of metric name, value, rating, and score for the 10 individual metrics in Lyons \emph{et al.} (2001) IBI.
#' \item noProblems A logical that is \code{TRUE} if there were no discrepancies between the FMDB and aggregate data files.
#' \item countsProblems A vector of species names for which there was a discrepancy in the counts between the FMDB and aggregate data files.
#' \item weightsProblems A vector of species names for which there was a discrepancy in the weights between the FMDB and aggregate data files.
#' \item fmdb.data A data.frame that contains the species name, total count, and total biomass for species in the FMDB data file.
#' \item agg.data A data.frame that contains the species name, total count, total biomass, and total DELT number for species in the aggregate data file.
#' }
#'
#' \code{summary} prints a variety of items (see Details) and returns the overall IBI value.
#'
#' \code{coef} returns the overall IBI value.
#'
#' \code{plot} produces a number of plots but does not return any value.
#'
#' \code{problemsLRIBI} prints a variety of items (see Details) but does not return any value.
#'
#' @note It is expected that the FMDB and aggregate files (see Details) will include the survey number in the file name and that number will be at the end of the file name (with the exception of the CSV extension). The filenames may have other text prior to the survey number. It is suggested to not change the filenames from those returned by the WiDNR database reports (e.g., \dQuote{fish_raw_data_Survey375049456.csv} and \dQuote{aggdata_Survey375049456.csv}). However, if the text prior to the survey number is changed from the default then that text should be entered into \code{fmdb.prefix} and \code{agg.prefix} as appropriate.
#'
#' The ecoregion can be automatically detected from the county in which the survey took place. However, there are a handful of counties that exist in both the \code{North} and \code{South} ecoregions. In these instances, the automatic detection of ecoregion will throw an error asking the user to supply the specific ecoregion in \code{region}. Alternatively the user may over-ride the automatically determined ecoregion by supplying a specific ecoregion in \code{region}.
#'
#' The values in Table 3 of Lyons \emph{et al.} (2001) imply that most metric values were rounded to whole numbers (\dQuote{WPUE} and \dQuote{\%DELT} were rounded to one decimal) for converting to IBI ratings. Metric values are rounded prior to comparing to Table 3 if \code{round=TRUE} (the DEFAULT) or are not rounded if \code{round=FALSE}.
#'
#' Lyons \emph{et al.} (2001) noted that the diversity of sucker species in tributaries to Lake Superior and some tributaries to Green Bay, Lake Michigan had very low diversities of sucker species. Thus, the two sucker metrics (\dQuote{Sucker species} and \dQuote{\% Round suckers (wt)}) should be excluded from the overall IBI calculation. The sum of the other eight metrics should be multipled by 1.25 (=10/8) to compensate for the use of fewer metrics and produce a result on a similar scale. By default the two sucker metrics are used in the overall IBI calculation, but the two metrics can be excluded (and the corrective multiplier applied) by using \code{excludeSuckers=TRUE}.
#'
#' @seealso \code{\link{LRIBI_AppendixA}}, \code{\link{LRIBI_Table3}}, and \code{\link{LRIBI_Ecoregions}}.
#'
#' @source Lyons, J., R.R. Piette, and K.W. Niermeyer. 2001. Development, validation, and application of a fish-based index of biotic integrity for Wisconsin's large warmwater rivers. Transactions of the American Fisheries Society 130:1077-1094.
#'
#' @author Derek H. Ogle (\email{derek@@derekogle.com}) and Joseph T. Mrnak
#'
#' @keywords manip plot
#'
#' @aliases LRIBI summary.LRIBI plot.LRIBI problemsLRIBI getLRIBISurveys
#'
#' @examples
#' ## Get location of external data files
#' dir <- dirname(system.file("extdata","aggdata_Survey375049456.csv",package="fishWiDNR"))
#'
#' ## Example of one survey
#' fmdb1 <- read.FMDB(paste0(dir,"/fish_raw_data_Survey375049456.csv"),expandCounts=TRUE)
#' agg1 <- read.csv(paste0(dir,"/aggdata_Survey375049456.csv"),stringsAsFactors=FALSE)
#' res <- LRIBI(fmdb1,agg1,effort=0.9)
#'
#' coef(res)
#' summary(res)
#' summary(res,show.data=FALSE)
#' summary(res,intermediate=FALSE)
#' summary(res,show.data=FALSE,intermediate=FALSE)
#' summary(res,show.data=FALSE,intermediate=FALSE,show.title=FALSE)
#'
#' plot(res)
#' plot(res,rating.col=NULL)
#' plot(res,rating.col="lightblue")
#' plot(res,rating.col=c("lightcyan","lightsalmon","lightblue"))
#' plot(res,type="counts",xlim=c(0,80))
#' plot(res,type="biomass",xlim=c(0,25))
#' plot(res,type="weights",xlim=c(0,25)) # same
#' plot(res,type="ratings")
#' plot(res,type="ratings",rating.col=NULL)
#' plot(res,type="ratings",rating.col="lightblue")
#' plot(res,type="ratings",rating.col=c("lightcyan","lightsalmon","lightblue"))
#'
#' problemsLRIBI(res)
#'
#' ## Another example of one survey
#' fmdb2 <- read.FMDB(paste0(dir,"/fish_raw_data_Survey375334018.csv"),expandCounts=TRUE)
#' agg2 <- read.csv(paste0(dir,"/aggdata_Survey375334018.csv"),stringsAsFactors=FALSE)
#'
#' res2 <- LRIBI(fmdb2,agg2,effort=0.9)
#' summary(res2)
#' plot(res2)
#' problemsLRIBI(res2)
#' problemsLRIBI(res2,type="count")
#' problemsLRIBI(res2,type="weight")
#'
#' # uses not rounded metric values
#' res2a <- LRIBI(fmdb2,agg2,effort=0.9,round=FALSE)
#' coef(res2); coef(res2a) # rounding can make a difference
#'
#' @rdname LRIBI
#' @export
LRIBI <- function(fmdb,agg,effort,
region=NULL,round=TRUE,excludeSuckers=FALSE){
## Load the internal lookup tables into this function's environment
## the data/get combination avoids the "no global binding" note at CHECK
appA <- get(utils::data("LRIBI_AppendixA",envir=environment()),envir=environment())
tab3 <- get(utils::data("LRIBI_Table3",envir=environment()),envir=environment())
ecoregs <- get(utils::data("LRIBI_Ecoregions",envir=environment()),envir=environment())
## Load FMDB and aggregate files (retruned are ready to be summarized)
fmdb <- iLRIBI_hndlFMDB(fmdb,appA)
agg <- iLRIBI_hndlAGG(agg,appA)
## Make sure efforts are reasonable
iCheckEffort(effort)
## Find county (make sure one county, capitalized properly, and not a factor)
cnty <- as.character(capFirst(unique(fmdb$County)))
## Find or handle the ecoregion
region <- iHndlRegion(region,cnty,ecoregs)
## Find the survey
survey <- iHndlSurvey(fmdb)
## Find the fish visit sequence number
visit <- iHndlVisit(fmdb)
## Make a title for the outputs
titles <- iMakeTitle(survey,visit,cnty,region,fmdb)
## Summarize and join two data types
joined <- iLRIBI_SumJoin(fmdb,agg,survey)
# Add grouping variables to the sum.data
joined$sum.data <- iLRIBI_AddGroupVars(joined$sum.data,appA)
## Compute individual metrics and lookup IBI scores
IBI.data <- iLRIBI_ComputeIBIMetrics(joined$sum.data,tab3,effort,
region,round,excludeSuckers)
## Return list of results in an LRIBI class
res <- list(survey=survey,region=region,effort=effort,
short.title=titles$short.title,
long.title=titles$long.title,
IBI=ifelse(excludeSuckers,sum(IBI.data$Score)*1.25,
sum(IBI.data$Score)),
n.fish=sum(joined$sum.data$Count),
ttl.Biomass=sum(joined$sum.data$Biomass),
sum.data=joined$sum.data,
IBI.data=IBI.data,
noProblems=is.null(c(joined$countsProblems,joined$weightsProblems)),
countsProblems=joined$countsProblems,
weightsProblems=joined$weightsProblems,
fmdb.data=joined$fmdb.data,
agg.data=joined$agg.data)
class(res) <- "LRIBI"
res
}
#' @rdname LRIBI
#' @export
summary.LRIBI <- function(object,intermediate=TRUE,show.data=TRUE,show.title=TRUE,...) {
if (show.title) {
cat(object$long.title,"\n")
cat("\nIBI =",object$IBI,"\n")
}
if (intermediate) {
if (show.title) cat("\n")
cat("IBI computed from these metric values, ratings and scores:\n")
print(object$IBI.data,row.names=FALSE)
}
if (show.data) {
if (show.title | intermediate) cat("\n")
cat("These observed data were used:\n")
print(object$sum.data[,c("Species1","Count","Biomass","DELTnum")],row.names=FALSE)
}
if (show.title | intermediate | show.data) invisible(object$IBI)
else return(object$IBI)
}
#' @rdname LRIBI
#' @export
coef.LRIBI <- function(object,...) object$IBI
#' @rdname LRIBI
#' @export
plot.LRIBI <- function(x,type=c("dashboard","counts","weights","biomass","ratings","scores"),
y=NULL,label.mar=9,axis.mult=1.25,
rating.col=FSA::col2rgbt(c("red","yellow","green"),1/2),
xlab=NULL,main=NULL,...) {
type <- match.arg(type)
switch(type,
dashboard= { iLRIBI_Dashboard(x,axis.mult,rating.col,main,...) },
counts= { iLRIBI_barData(x,type,label.mar,xlab,main,...) },
weights=,biomass= { iLRIBI_barData(x,type,label.mar,xlab,main,...) },
ratings=,scores= { iLRIBI_barMetrics(x,label.mar,rating.col,xlab,main,...) } )
}
#' @rdname LRIBI
#' @export
problemsLRIBI <- function(object,type=c("counts","weights")) {
type <- match.arg(type,several.ok=TRUE)
if (object$noProblems) message("No discrepancies between the FMDB and aggregate files\n were encountered for ",object$short.title)
else {
## Get fmdb and aggregate data together
tmp <- dplyr::left_join(object$fmdb.data,object$agg.data,by="Species1")[,-6]
names(tmp) <- sub(".x",".fmdb",names(tmp))
names(tmp) <- sub(".y",".agg",names(tmp))
## Address problems with counts
if ("counts" %in% type) {
if (is.null(object$countsProblems)) {
# don't send message if weights also in type
if (!"weights" %in% type) message("There was no discrepancies between the counts in the FMDB and\n aggregate files for ",object$short.title)
} else {
cat(object$long.title,"\n\n",sep="")
# create new variable that will have YES! to note discrepancy
tmp$PROBLEM <- ""
tmp$PROBLEM[tmp$Species1 %in% object$countsProblems] <- "YES!"
print(tmp[,-which(names(tmp) %in% c("Biomass.fmdb","Biomass.agg"))],
row.names=FALSE)
}
}
if ("weights" %in% type) {
if (is.null(object$weightsProblems)) {
# don't send message if counts also in type
if (!"counts" %in% type) message("There was no discrepancies between the weights in the FMDB and\n aggregate files for ",object$short.title)
} else {
cat(object$long.title,"\n\n",sep="")
# create new variable that will have YES! to note discrepancy
tmp$PROBLEM <- ""
tmp$PROBLEM[tmp$Species1 %in% object$weightsProblems] <- "YES!"
print(tmp[,-which(names(tmp) %in% c("Count.fmdb","Count.agg"))],
row.names=FALSE)
}
}
}
}
################################################################################
################################################################################
##### INTERNAL (not exported) COMPUTATION FUNCTIONS
##### Load FMDB data
iLRIBI_hndlFMDB <- function(fmdb,appA) {
if (!is.data.frame(fmdb)) stop("'fmdb' must be a data.frame from the FMDB.",
call.=FALSE)
# Reduce to only "all species" (Z100) surveys
fmdb <- fmdb[fmdb$Target.Species=="ALL SPECIES",,drop=TRUE]
# Convert g to kg
fmdb$Weight.KG <- fmdb$Weight.Grams/1000
# Make species a character variable (rather than factor)
fmdb$Species1 <- as.character(fmdb$Species1)
# Remove species not in Appendix A
fmdb <- fmdb[fmdb$Species1 %in% appA$common.name,,drop=TRUE]
# Group by species to make read for summarizing
fmdb <- dplyr::group_by_(fmdb,~Species1)
fmdb
}
##### Load aggregate data
iLRIBI_hndlAGG <- function(agg,appA) {
if (!is.data.frame(agg)) stop("'agg' must be an aggregate file data.frame.",
call.=FALSE)
# Convert species capitalization to match FMDB file (make sure it is character)
agg$Species1 <- as.character(FSA::capFirst(agg$Common.Name))
# Convert grams or lbs to kg
agg$Aggregate.KG <- ifelse(agg$Units=="GRAMS",agg$Aggregate.Weight/1000,
agg$Aggregate.Weight*454/1000)
# Remove species not in Appendix A
agg <- agg[agg$Species1 %in% appA$common.name,,drop=TRUE]
# Group by species to make read for summarizing
agg <- dplyr::group_by_(agg,~Species1)
agg
}
##### Join the FMDB and Aggregate files (and handle potential problems)
iLRIBI_SumJoin <- function(fmdb,agg,survey) {
## Initialize some variables
countsProblems <- weightsProblems <- NULL
## Count number of fish and sum weights in FMDB file
## dots needed to use with summarize_ for using dplyr in a function
dots <- list(~length(Species1),~sum(Weight.KG))
fmdb1 <- dplyr::summarize_(fmdb,.dots=stats::setNames(dots,
c("Count","Biomass")))
## Sum number of fish, aggregate weight, and number of DELT fish
## in aggregate file
dots <- list(~sum(Total.Count),~sum(Aggregate.KG),~sum(DELT..))
agg1 <- dplyr::summarize_(agg,.dots=stats::setNames(dots,
c("Count","Biomass","DELTnum")))
## Join the two summary date.frames by Species1
sum.data <- dplyr::full_join(fmdb1,agg1,by="Species1")
## Perform checks on counts
# problem if species in aggregate but not FMDB file
if (any(is.na(sum.data$Count.x))) {
warning("Some species in survey #",survey," had counts in the aggregate but not FMDB file.",call.=FALSE)
countsProblems <- c(countsProblems,
sum.data$Species1[is.na(sum.data$Count.x)])
}
# problem if counts between the two files do not match
tmp <- sum.data[!is.na(sum.data$Count.y),]
if (any(tmp$Count.x!=tmp$Count.y)) {
warning("Counts in FMDB and aggregate files for survey #",survey," do not match\n for some species! Continued with FMDB counts.",call.=FALSE)
countsProblems <- c(countsProblems,tmp$Species1[tmp$Count.x!=tmp$Count.y])
}
## Perform checks on weights
# problem if there are weights in both files and do they not match
tmp <- sum.data[!is.na(sum.data$Biomass.x) & !is.na(sum.data$Biomass.y),]
if (nrow(tmp)>0) {
if (nrow(tmp)==nrow(sum.data)) stop("There are no weights in either the FMDB or aggregate files.",call.=FALSE)
if (any(tmp$Biomass.x!=tmp$Biomass.y)) warning("Weights that exist in both FMDB and aggregate files for survey #",survey," do not match!\n Continued with aggregage weight.",call.=FALSE)
weightsProblems <- tmp$Species1[tmp$Biomass.x!=tmp$Biomass.y]
}
## Replace NA biomasses from FMDB with biomasses from aggregate file
tmp <- sum.data[is.na(sum.data$Biomass.x),]
tmp <- which(sum.data$Species1 %in% tmp$Species1)
sum.data$Biomass.x[tmp] <- sum.data$Biomass.y[tmp]
## Replace NAs in DELTnum with zeroes
sum.data$DELTnum[is.na(sum.data$DELTnum)] <- 0
## Rename the ".x" variables, Delete the ".y" variables
names(sum.data)[names(sum.data) %in% c("Count.x","Biomass.x")] <- c("Count","Biomass")
sum.data <- sum.data[,-which(names(sum.data) %in% c("Count.x","Count.y","Biomass.x","Biomass.y"))]
## Check if any weights are missing
if (any(is.na(sum.data$Biomass))) {
warning("Some species had counts but no weights recorded for survey #",survey,".",call.=FALSE)
weightsProblems <- unique(c(weightsProblems,
sum.data$Species1[is.na(sum.data$Biomass)]))
}
## Return result (use as.data.frame to strip off tibble attributes)
list(sum.data=as.data.frame(sum.data),
countsProblems=countsProblems,
weightsProblems=weightsProblems,
fmdb.data=as.data.frame(fmdb1),
agg.data=as.data.frame(agg1))
}
##### Add grouping variables required to create IBI metrics
iLRIBI_AddGroupVars <- function(sum.data,appA) {
## use mapvalues to lookup grouping variable levels by species common name
## used warn_missing=FALSE to suppress warnings when level did not exist
sum.data$habitat <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$habitat,
warn_missing=FALSE)
sum.data$habitat <- FSA::mapvalues(sum.data$habitat,
from="River-Large",to="River",warn_missing=FALSE)
sum.data$suckers <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$common.family,
warn_missing=FALSE)
sum.data$tolerance <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$tolerance,
warn_missing=FALSE)
sum.data$origin <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$origin,
warn_missing=FALSE)
sum.data$spawning <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$spawning,
warn_missing=FALSE)
sum.data$feeding <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$feeding,
warn_missing=FALSE)
sum.data$round.sucker <- FSA::mapvalues(sum.data$Species1,
from=appA$common.name,to=appA$round.sucker,
warn_missing=FALSE)
sum.data
}
##### Compute the individual metrics
iLRIBI_ComputeIBIMetrics <- function(sum.data,tab3,effort,
region,round,excludeSuckers) {
## Individual metric summaries
metvals <- iLRIBI_SumValues(sum.data,effort)
## Individual metric ratings (Good, Fair, Poor)
metrats <- iLRIBI_Convert2Ratings(metvals,tab3,region,round)
## Create data.frame with metric names, values, and ratings
tmp <- data.frame(Metric=c("WPUE","# Native species","# Sucker species",
"# Intolerant species","# Riverine species","% DELT (n)",
"% Riverine (n)","% Lithophils (n)",
"% Insectivore (wt)","% Round suckers (wt)"),
Value=metvals,
Rating=metrats,stringsAsFactors=FALSE)
## Add new variable that has metric score (10,5,0)
tmp$Score <- FSA::fact2num(FSA::mapvalues(tmp$Rating,
from=c("Poor","Fair","Good"),to=c(0,5,10),
warn_missing=FALSE))
row.names(tmp) <- NULL
if (excludeSuckers) tmp <- tmp[!grepl("ucker",tmp$Metric),]
tmp
}
##### Extracts specific metrics for individual IBI metrics
iLRIBI_SumValues <- function(sum.data,effort) {
## Short helper function to make calculations above look cleaner
iPercTable <- function(formula,data,FUN) {
suppressWarnings(100*prop.table(FSA::sumTable(formula,data,FUN)))
}
## must exclude Tolerant species for WPUE calculation
c(WPUE=sum(sum.data$Biomass[sum.data$tolerance!="Tolerant"])/effort,
num.native=stats::xtabs(~origin,data=sum.data)[["Native"]],
num.sucker=stats::xtabs(~suckers,data=sum.data)[["Suckers"]],
num.intolerant=stats::xtabs(~tolerance,data=sum.data)[["Intolerant"]],
num.riverine=stats::xtabs(~habitat,data=sum.data)[["River"]],
perc.DELT=100*sum(sum.data$DELTNum)/sum(sum.data$Count),
perc.riverine=iPercTable(Count~habitat,data=sum.data,FUN=sum)[["River"]],
perc.lithophils=iPercTable(Count~spawning,data=sum.data,FUN=sum)[["Lithophil"]],
perc.insectivore=iPercTable(Biomass~feeding,data=sum.data,FUN=sum)[["Insectivore"]],
perc.rsuckers=iPercTable(Biomass~round.sucker,data=sum.data,FUN=sum)[["Yes"]])
}
##### Converts metric values into Poor, Fair, Good ratings
iLRIBI_Convert2Ratings <- function(metvals,tab3,region,round) {
## Internal function to make below cleaner
iTab3Lookup <- function(var,metvals,tab3,round,digits=0,less=TRUE) {
val <- metvals[[var]]
if (round) val <- round(val,digits)
if (less) {
ifelse(val<tab3$fair[tab3$metric==var],"Poor",
ifelse(val<=tab3$good[tab3$metric==var],"Fair","Good"))
} else {
ifelse(val>tab3$fair[tab3$metric==var],"Poor",
ifelse(val>=tab3$good[tab3$metric==var],"Fair","Good"))
}
} # end internal
## Lookup ratings for the specific ecoregion
tab3 <- tab3[tab3$ecoregion==region,,drop=TRUE]
c(WPUE_rate=iTab3Lookup("WPUE",metvals,tab3,round,digits=1),
num.native_rate=iTab3Lookup("num.native",metvals,tab3,round),
num.sucker_rate=iTab3Lookup("num.sucker",metvals,tab3,round),
num.intolerant_rate=iTab3Lookup("num.intolerant",metvals,tab3,round),
num.riverine_rate=iTab3Lookup("num.riverine",metvals,tab3,round),
perc.DELT_rate=iTab3Lookup("perc.DELT",metvals,tab3,round,digits=1,less=FALSE),
perc.riverine_rate=iTab3Lookup("perc.riverine",metvals,tab3,round),
perc.lithophils_rate=iTab3Lookup("perc.lithophils",metvals,tab3,round),
perc.insectivore_rate=iTab3Lookup("perc.insectivore",metvals,tab3,round),
perc.rsuckers_rate=iTab3Lookup("perc.rsuckers",metvals,tab3,round))
}
################################################################################
################################################################################
##### INTERNAL (not exported) PLOTTING FUNCTIONS
##### PLOT of data (WPUE, etc.)
iLRIBI_barData <- function(object,type,label.mar,xlab,main,...) {
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle1(object,main)
if (type=="counts") {
tmp <- object$sum.data$Count
if (is.null(xlab)) xlab="Observed Number"
} else {
tmp <- object$sum.data$Biomass
if (is.null(xlab)) xlab="Observed Weight (KG)"
}
graphics::par(mar=c(3,label.mar,0,1),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
bp <- graphics::barplot(tmp,names.arg=object$sum.data$Species1,
horiz=TRUE,xlab=xlab,...)
graphics::axis(2,at=bp,labels=NA)
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
}
##### Plot of metric ratings (Poor, Fair, Good)
iLRIBI_barMetrics <- function(object,label.mar,clrs,xlab,main,...) {
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle1(object,main)
clrs <- iHndlRatingColors(clrs)[factor(object$IBI.data$Rating,
levels=c("Poor","Fair","Good"))]
if (is.null(xlab)) xlab="Metric Rating (Score)"
graphics::par(mar=c(3,label.mar,0,2),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
bp <- graphics::barplot(object$IBI.data$Score,names.arg=object$IBI.data$Metric,
horiz=TRUE,xlab=xlab,xaxt="n",col=clrs,...)
graphics::axis(2,at=bp,labels=NA)
graphics::axis(1,at=c(0,5,10),labels=c("Poor (0)","Fair (5)","Good (10)"))
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
}
##### Dashboard plot of each metric value vs. its cutoff scores
iLRIBI_Dashboard <- function(object,axis.mult,clrs,main) {
## Internal function for individual bar plots
iLRIBI_plotBarIndiv <- function(object,tab3,axis.mult,metnum,clrs) {
df <- object$IBI.data
val <- df$Value[metnum]
# find the metric cutoffs
cuts <- tab3[tab3$ecoregion==object$region,,drop=TRUE]
cuts <- cuts[metnum,,drop=TRUE]
# adjust y-axis limits by type of metric
if (cuts$metric=="WPUE") {
# for WPUE be more than observed value of maximum cutoff
ymax <- max(axis.mult*val,axis.mult*max(cuts$good))
} else if (grepl("num",cuts$metric)) {
# for counts be more than observed, maximum cutoff+1, or 5 (helps with axis ticks)
ymax <- max(axis.mult*val,axis.mult*max(cuts$good+1),5)
} else {
# for percents be 100%
ymax <- 100
}
clrs <- clrs[factor(df$Rating[metnum],levels=c("Poor","Fair","Good"))]
op <- graphics::par(mar=c(1.2,3,1.2,0.7),mgp=c(1.9,0.5,0),tcl=-0.2,las=1)
graphics::barplot(val,ylab=df$Metric[metnum],ylim=c(0,ymax),col=clrs)
graphics::mtext(df$Rating[metnum],line=-0.5,cex=0.8)
graphics::abline(h=c(cuts$fair,cuts$good),col="gray50",lwd=3,lty=2)
graphics::par(op)
} # end internal
# get colors
clrs <- iHndlRatingColors(clrs)
# Make a layout for 9 plots (exclues perc.DELT ... see below)
tab3 <- get(utils::data("LRIBI_Table3",envir=environment()),envir=environment())
op <- graphics::par(no.readonly=TRUE)
iHndlMainTitle2(object,main)
# Plots individual metrics
iLRIBI_plotBarIndiv(object,tab3,axis.mult,1,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,2,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,3,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,4,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,5,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,7,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,8,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,9,clrs)
iLRIBI_plotBarIndiv(object,tab3,axis.mult,10,clrs)
# return graphing parameters and layout to original
graphics::par(op); graphics::layout(1)
# did not plot DELT because nearly always zero???
# this sends a message reminding that it was not plotted.
message("The '% DELT (n)' metric of ",object$IBI.data$Rating[6],
" (",round(object$IBI.data$Value[6],2),"%) is not shown.")
}
##### Handle putting the main title on the plot for barplots
iHndlMainTitle1 <- function(object,main) {
if (is.null(main)) main <- object$short.title
if (main!="") {
graphics::layout(matrix(1:2,nrow=2),heights=c(0.5,9))
graphics::par(mar=c(0,0,0,0))
graphics::plot.new()
graphics::text(0.5,0.5,main,cex=0.85)
}
}
##### Handle putting the main title on the plot for dashboard
iHndlMainTitle2 <- function(object,main) {
if (is.null(main)) main <- object$short.title
if (main=="") graphics::layout(matrix(1:9,nrow=3,byrow=TRUE))
else {
graphics::layout(matrix(c(rep(1,3),2:10),nrow=4,byrow=TRUE),
heights=c(0.5,3,3,3))
graphics::par(mar=c(0,0,0,0))
graphics::plot.new()
graphics::text(0.5,0.5,main,cex=1.15)
}
}
##### Handle dealing with the colors used for the different ratings
iHndlRatingColors <- function(clrs) {
if (is.null(clrs)) clrs <- "gray70"
if (length(clrs)!=3) {
if (length(clrs)==1) clrs <- rep(clrs,3)
else stop("Must supply 1 or 3 colors to 'rating.col'.",call.=FALSE)
}
clrs
}
##### Make sure that effort is reasonable
iCheckEffort <- function(effort) {
if (length(effort)!=1L) stop("Only one 'effort' value can be used.",call.=FALSE)
if (effort<=0) stop("'effort' must be positive.",call.=FALSE)
if (effort>2) warning("Check that your 'effort' is in miles.",call.=FALSE)
}
##### Check or handle finding the ecoregion
iHndlRegion <- function(region,cnty,ecoregs) {
if (!is.null(region)) {
## user-supplied ecoregion ... make sure it is North or South
if (!region %in% c("North","South"))
stop("'region' must be 'North' or 'South'.",call.=FALSE)
} else {
## ecoregion not supplied by user ... find it from county
if (is.null(region)) region <- ecoregs$ecoregion[which(ecoregs$county==cnty)]
if (region=="Both") stop("This survey is in a county that is in both ecoregions.\n Please use 'region' to specify which ecoregions should be used.",call.=FALSE)
}
region
}
iHndlSurvey <- function(fmdb) {
survey <- unique(fmdb$Survey.Seq.No)
if (length(survey)>1L) stop("The IBI can only be computed for one survey number.",
call.=FALSE)
survey
}
iHndlVisit <- function(fmdb) {
visit <- unique(fmdb$Visit.Fish.Seq.No)
if (length(visit)>1L) warning("More than one 'Visit.Fish.Seq.No' is present in the 'fmdb' file.",call.=FALSE)
visit
}
#### Make a title for outputs
iMakeTitle <- function(survey,visit,cnty,region,fmdb) {
## both titles get survey and visit numbers
tmp <- paste0("Survey #",survey," (")
if (length(visit)>1L) tmp <- paste0(tmp,"multiple visits)")
else tmp <- paste0(tmp,"Visit #",visit,")")
## short title also gets watebordy name and survey year
stmp <- paste0(tmp," - ",as.character(capFirst(unique(fmdb$Waterbody.Name))),
", ",unique(fmdb$Survey.Year))
## long title also gets waterbody name, county, ecoregion, and survey year
ltmp <- paste0(tmp,"\n",as.character(capFirst(unique(fmdb$Waterbody.Name))),
", ",cnty," County, ",region," ecoregion, ",
unique(fmdb$Survey.Year))
list(short.title=stmp,long.title=ltmp)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.