#' @title addnorm - adds a normal distribution to a histogram of a data set.
#'
#' @description addnorm - adds a normal distribution to a histogram of a data
#' set. This is generally to be used to illustrate whether log-transformation
#' normalizes a set of catch or cpue data.
#' @param inhist - is the output from a call to 'hist' (see examples)
#' @param xdata - is the data that is being plotted in the histogram.
#' @param inc - defaults to a value of 0.01; is the fine grain increment used to
#' define the normal curve. The histogram will be coarse grained relative to
#' this.
#' @return a list with a vector of 'x' values and a vector of 'y' values (to be
#' used to plot the fitted normal probability density function), and a vector
#' used two called 'stats' containing the mean and sandard deviation of the
#' input data
#' @export addnorm
#' @examples
#' x <- rnorm(1000,mean=5,sd=1)
#' dev.new(height=6,width=4,noRStudioGD = TRUE)
#' par(mfrow= c(1,1),mai=c(0.5,0.5,0.3,0.05))
#' par(cex=0.85, mgp=c(1.5,0.35,0), font.axis=7)
#' outH <- hist(x,breaks=25,col=3,main="")
#' nline <- addnorm(outH,x)
#' lines(nline$x,nline$y,lwd=3,col=2)
#' print(nline$stats)
addnorm <- function(inhist,xdata,inc=0.01) {
lower <- inhist$breaks[1]
upper <- tail(inhist$breaks,1)
cw <- inhist$breaks[2]-inhist$breaks[1]
x <- seq(lower,upper, inc) #+ (cw/2)
avCE <- mean(xdata,na.rm=TRUE)
sdCE <- sd(xdata,na.rm=TRUE)
N <- length(xdata)
ans <- list(x=x,y=(N*cw)*dnorm(x,avCE,sdCE),stats=c(avCE,sdCE,N))
return(ans)
} # end of addnorm
#' @title calcrmse calcuates the root mean square error for two vectors
#'
#' @description given a set of observations and predicted values, calcrmse
#' calculates the root mean square error sqrt(resid^2/n) to provide a
#' measure of relative fit. An assumption of normal errors is made so
#' if the observations and predicted values actually have a log-normal
#' distribution one should log-transform the input values.
#'
#' @param obs the observed data
#' @param pred the predicted values whose fit to the observations is to
#' measured
#'
#' @return a scalar value which is the rmse.
#' @export
#'
#' @examples
#' \dontrun{
#' x <- rep(1,10)
#' y <- c(1.109,1.210,0.947,0.933,0.832,0.864,0.633,0.820,1.004,1.049)
#' calcrmse(log(x),log(y)) # should be 0.1899495
#' }
calcrmse <- function(obs,pred){ # obs=fishery[,"CPUE"]; pred=fishery[,"PredCE"]
pickCE <- which(!is.na(obs))
resid <- obs[pickCE] - pred[pickCE]
rmse <- sqrt(sum(resid^2)/length(resid))
return(rmse)
}
#' @title checkdata determines which methods match the input data
#'
#' @description checkdata ensures the input fishery data contains the
#' data necessary for each of the analyses present and reports which
#' analyses are possible. It checks for simple errors
#'
#' @param indata this can be either the complete list of data objects of just
#' the matrix of fisheries data.
#' @param glob in case one inputs a data.frame of fishery data this allows one
#' to enter the globals list containing as a minimum the spsname and
#' resilience, so that a test can be made for the spm analysis
#'
#' @return a 3 x 2 matrix with vaiable and true or false for presence
#' @export
#'
#' @examples
#' \dontrun{
#' data(fishdat)
#' fish <- fishdat$fish
#' checkdata(fishdat)
#' checkdata(fish)
#' }
checkdata <- function(indata,glob=NA) { # indata=fish; glob=glb
rown <- c("catch-MSY","spm","aspm","catch-curves")
result <- as.data.frame(matrix(0,nrow=4,ncol=2,
dimnames=list(rown,c("Method","Possible"))))
result[,"Method"] <- rown
result[,"Possible"] <- FALSE
if (class(indata) == "data.frame") {
if ((all(c("year","catch") %in% colnames(indata))) &
(length(grep("cpue",colnames(indata))) > 0))
result["spm",] <- TRUE
if (is.list(glob)) {
if ((all(c("year","catch") %in% colnames(indata))) &
(glb$resilience %in% c("verylow","low","medium","high")))
result["catch-MSY",] <- TRUE
} else {
print("Only the fisheries data.frame checked for use with spm.")
}
return(result)
} else if (class(indata) == "list") {
fish <- indata$fish
glb <- indata$glb
glbnames <- c("maxage","M","Linf","K","t0","Waa","Wab","M50a","deltaM",
"steep","R0","sela50","deltaS","resilience","nages","ages",
"nyrs","spsname")
props <- indata$props
agedata <- indata$agedata
if ((all(c("year","catch") %in% colnames(fish))) &
(glb$resilience %in% c("verylow","low","medium","high")))
result["catch-MSY",] <- TRUE
if (all(c("year","catch","cpue") %in% colnames(fish)))
result["spm",] <- TRUE
if ((all(c("year","catch","cpue") %in% colnames(fish))) &
((all(glbnames %in% names(glb))) | (class(props) == "data.frame")))
result["aspm",] <- TRUE
if (class(agedata) != "NULL") result["catch-curves",] <- TRUE
return(result)
} else {
return(result)
}
} # end of checkspmdata
#' @title countgtzero used in apply to count the number of zeros in a vector
#'
#' @description countgtzero is designed to be used within the apply function to
#' count the number of zeros in a vector
#' @param invect vector of values
#' @return A single value of zero or the number of zeros
#' @export countgtzero
#' @examples
#' \dontrun{
#' x <- c(1,2,3,4,0,0,0,5,6,7,8) # 8 values > 0
#' countgtzero(x)
#' }
countgtzero <- function(invect) {
return(length(which(invect > 0)))
}
#' @title dataTemplate - generates a standard input datafile as a template
#'
#' @description dataTemplate - generates a standard input datafile to use as
#' a template. It is possible to define the number of blocks and then, once
#' the data file is created, go in and edit it appropriately to suit exactly
#' the data you have for your own species. Note that in the example commas
#' are used to separate individual data fields, this is required if you are
#' going to use the function 'readdata' to get your data into datalowSA,
#' which is recommended. But of course you are free to use it however you
#' wish! Each section in the data file, except the title, is identified
#' using CAPITAL letters, as in RESILIENCE, NYRS, YEARS, etc. These, again
#' are needed by 'readdata' to identify the various sections and proceed
#' to reading that data in appropriately. If you do not have data for a
#' section just delete it entirely, including the heading. DO NOT use a
#' section heading name anywhere else in the document as 'readdata' will
#' always try to use the first instance which could cause you trouble. The
#' example AGE data comes from Table 13.1 in Beverton and Holt (1957) and
#' consists of ageing data from North Sea plaice. Tw sexes have been given
#' here to illustrate the use of age by sex (female first row for a year,
#' males in the second row for a year), but zeros are given for males, The
#' numbers at length are pure fiction and again are merely here to
#' illustrate the format.
#'
#' @param filename - the name for the generated datafile, character string,
#' defaults to tmpdat.csv
#' @param title the first line of the data file; identify its contents
#' @return a standard definition data file ready to be read by readdata
#' and it also outputs the full address of the filename.
#' @export
#' @examples
#' \dontrun{
#' dataTemplate(filename="fishery1.csv")
#' ans <- readdata("fishery1.csv")
#' str(ans)
#' } # filename="C:/Users/had06a/Dropbox/Public/fishery1.csv"
dataTemplate <- function(filename="constants.csv",title="SpeciesName") {
cat(title, "\n\n",file=filename,append=FALSE)
cat("# The format of this file is such that the first line contains \n",
file=filename, append=TRUE)
cat("# the species name or title for the subsequent analyses \n",
file=filename,append=TRUE)
cat("# Each section contains a header fully in CAPITALS, which is used by \n",
file=filename, append=TRUE)
cat("# the function 'readdata' to identify the following data and read it \n",
file=filename, append=TRUE)
cat("# in appropriately. If you do not have data for a particular \n",
file=filename, append=TRUE)
cat("# section, just delete it and its header. \n",
file=filename, append=TRUE)
cat("# DO NOT confuse the software by having the same word used as \n",
file=filename, append=TRUE)
cat("# a section header occur more than once in the data file \n",
file=filename, append=TRUE)
cat("# Always save the file as a comma separated variable file as this \n",
file=filename, append=TRUE)
cat("# is the default format when using 'readdata'. \n\n",
file=filename, append=TRUE)
cat("RESILIENCE \n",file=filename, append=TRUE)
cat("verylow \n\n",file=filename, append=TRUE)
cat("NYRS \n",file=filename, append=TRUE)
cat("30, \n\n",file=filename, append=TRUE)
cat("YEARS,catch,cpue,SE \n",file=filename, append=TRUE)
cat("1986,3924.912,NA,NA \n",file=filename, append=TRUE)
cat("1987,5117.988,NA,NA \n",file=filename, append=TRUE)
cat("1988,4722.2,NA,NA \n",file=filename, append=TRUE)
cat("1989,1365.128,1.49,0.2 \n",file=filename, append=TRUE)
cat("1990,801.567,1.849,0.2 \n",file=filename, append=TRUE)
cat("1991,625.407,1.501,0.2 \n",file=filename, append=TRUE)
cat("1992,1108.241,2.997,0.2 \n",file=filename, append=TRUE)
cat("1993,964.409,1.415,0.2 \n",file=filename, append=TRUE)
cat("1994,800.618,1.139,0.2 \n",file=filename, append=TRUE)
cat("1995,962.399,0.7,0.2 \n",file=filename, append=TRUE)
cat("1996,1180.349,0.469,0.2 \n",file=filename, append=TRUE)
cat("1997,297.003,0.303,0.2 \n",file=filename, append=TRUE)
cat("1998,316.131,0.356,0.2 \n",file=filename, append=TRUE)
cat("1999,210.529,0.39,0.2 \n",file=filename, append=TRUE)
cat("2000,169.337,0.439,0.2 \n",file=filename, append=TRUE)
cat("2001,200.843,0.489,0.2 \n",file=filename, append=TRUE)
cat("2002,255.735,0.431,0.2 \n",file=filename, append=TRUE)
cat("2003,217.502,0.52,0.2 \n",file=filename, append=TRUE)
cat("2004,283.11,0.777,0.2 \n",file=filename, append=TRUE)
cat("2005,264.607,1.322,0.2 \n",file=filename, append=TRUE)
cat("2006,139.316,1.412,0.2 \n",file=filename, append=TRUE)
cat("2007,28.571,NA,NA \n",file=filename, append=TRUE)
cat("2008,3.331,NA,NA \n",file=filename, append=TRUE)
cat("2009,13.859,NA,NA \n",file=filename, append=TRUE)
cat("2010,21.44,NA,NA \n",file=filename, append=TRUE)
cat("2011,31.426,NA,NA \n",file=filename, append=TRUE)
cat("2012,17.253,NA,NA \n",file=filename, append=TRUE)
cat("2013,35.94,NA,NA \n",file=filename, append=TRUE)
cat("2014,22.087,NA,NA \n",file=filename, append=TRUE)
cat("2015,16.206,NA,NA \n\n\n",file=filename, append=TRUE)
cat("# Other variables can be added to the biology list below. The \n",
file=filename, append=TRUE)
cat("# standard format is to add the value followed by a comma, followed \n",
file=filename, append=TRUE)
cat("# by a descriptive name. This will then be able to be included in \n",
file=filename, append=TRUE)
cat("# the standard input function 'readdata' once we know what data you \n",
file=filename, append=TRUE)
cat("# may have that we did not think of. \n\n",file=filename, append=TRUE)
cat("BIOLOGY \n",file=filename, append=TRUE)
cat("80 , maxage the maximum age, usually a plus group \n",
file=filename, append=TRUE)
cat("0.036 , M natural mortality \n",file=filename, append=TRUE)
cat("39.6 , Linf vB asymptotic maximum length \n",file=filename, append=TRUE)
cat("0.06 , K vB Brody growth coefficient \n",file=filename, append=TRUE)
cat("-0.01 , t0 theoretical age at zero length\n",file=filename, append=TRUE)
cat("0.0365 , Waa weight at age a parameter \n",file=filename, append=TRUE)
cat("2.965 , Wab weight at age b parameter \n",file=filename, append=TRUE)
cat("31.0 , M50a age at 50% maturity \n",file=filename, append=TRUE)
cat("5 , deltaM diff between ages at 50 and 95% maturity \n",
file=filename, append=TRUE)
cat("31 , sela50 age at 50% selectivity \n",file=filename, append=TRUE)
cat("5 , deltaS diff between 50 and 95% selectivity \n",file=filename,
append=TRUE)
cat("0.6 , steep stock recruitment steepness \n",file=filename, append=TRUE)
cat("14 , Ln(R0) initial value but this will be estimated \n\n",
file=filename, append=TRUE)
cat("# The catch-at-age data derives from Table 13.1 in Beverton and \n",
file=filename, append=TRUE)
cat("# Holt, 1957; this version uses -99 to indicate missing date. \n\n",
file=filename, append=TRUE)
cat("AGE data \n",file=filename, append=TRUE)
cat("9 , years of numbers-at-age data \n",file=filename, append=TRUE)
cat("2 , sexes, females first=1 then males=2 \n",file=filename, append=TRUE)
cat("2, 3, 4, 5, 6, 7, 8, 9, 10, age classes \n",file=filename, append=TRUE)
cat("1929,1,328,2120,2783,1128,370,768,237,112,48 \n",file=filename, append=TRUE)
cat("1930,1,223,2246,1938,1620,302,106,181,58,18 \n",file=filename, append=TRUE)
cat("1931,1,95,2898,3017,1159,591,116,100,82,33 \n",file=filename, append=TRUE)
cat("1932,1,77,606,4385,1186,231,138,42,21,51 \n",file=filename, append=TRUE)
cat("1933,1,50,489,1121,4738,456,106,80,27,18 \n",file=filename, append=TRUE)
cat("1934,1,44,475,1666,1538,2510,160,50,43,14 \n",file=filename, append=TRUE)
cat("1935,1,131,1373,1595,1587,1326,883,144,30,28 \n",file=filename, append=TRUE)
cat("1936,1,38,691,2862,1094,864,382,436,27,15 \n",file=filename, append=TRUE)
cat("1937,1,138,1293,1804,1810,426,390,163,228,26 \n",file=filename, append=TRUE)
cat("1929,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1930,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1931,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1932,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1933,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1934,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1935,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1936,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n",file=filename, append=TRUE)
cat("1937,2,-99,-99,-99,-99,-99,-99,-99,-99,-99 \n\n\n",file=filename, append=TRUE)
cat("LENGTH data \n",file=filename, append=TRUE)
cat("2 , years of numbers-at-length data \n",file=filename, append=TRUE)
cat("2 , sexes, females=1 first then males=2 \n",file=filename, append=TRUE)
cat("15,18,21,24,27,30,33,36,39, length classes \n",file=filename, append=TRUE)
cat("1990,1,15,100,200,300,400,350,300,220,100 \n",file=filename, append=TRUE)
cat("1990,2,15,100,200,300,400,350,300,220,100 \n",file=filename, append=TRUE)
cat("1998,1,15,100,200,300,400,350,300,220,100 \n",file=filename, append=TRUE)
cat("1998,2,15,100,200,300,400,350,300,220,100 \n",file=filename, append=TRUE)
return(filename)
} # end of dataTemplate
#' @title facttonum converts a vector of numeric factors into numbers
#'
#' @description facttonum converts a vector of numeric factors into numbers.
#' If the factors are not numeric then the outcome will be a series of NA.
#' It is up to you to apply this function only to numeric factors. A warning
#' will be thrown if the resulting output vector contains NAs
#'
#' @param invect the vector of numeric factors to be converted back to numbers
#'
#' @return an output vector of numbers instead of the input factors
#' @export
#'
#' @examples
#' \dontrun{
#' DepCat <- as.factor(rep(seq(100,600,100),2)); DepCat
#' 5 * DepCat[3]
#' as.numeric(levels(DepCat)) # #only converts the levels not the replicates
#' DepCat <- facttonum(DepCat)
#' 5 * DepCat[3]
#' x <- factor(letters)
#' facttonum(x)
#' }
facttonum <- function(invect){
if (class(invect) == "factor") {
outvect <- suppressWarnings(as.numeric(levels(invect))[invect])
}
if (class(invect) == "numeric") outvect <- invect
if (any(is.na(outvect)))
warning("NAs produced, your input vector may have non-numbers present \n")
return(outvect)
} # end of facttonum
#' @title getLNCI gets the log-normal confidence intervals
#'
#' @description getLNCI takes the mean and the standard deviation and produces
#' the asymmetric log-normal confidence intervals around the mean values
#'
#' @param av the mean value or a vector of mean values
#' @param se the standard deviation
#' @param P the percent used for the CI, defaults to 95.
#'
#' @return a vector of three for a single input or a matrix of 3 columns for
#' input vectors
#' @export
#'
#' @examples
#' \dontrun{
#' av <- c(4.0,2.15)
#' se <- 0.33
#' getLNCI(av,se,P=95)
#' se <- c(0.33,0.4)
#' getLNCI(av,se)
#' }
getLNCI <- function(av,se,P=95){ # av=fissp[,"CPUE"]; se=rmse;P=0.95
Zmult <- -qnorm((1-(P/100))/2.0)
lower <- av * exp(-Zmult*se)
upper <- av * exp(Zmult*se)
if (length(av) > 1) {
result <- cbind(lower,av,upper)
colnames(result) <- c("lower","mean","upper")
} else {
result <- c(lower=lower,mean=av,upper=upper)
}
return(result)
} # end of getLNCI
#' @title getminy generates the lower bound for a plot
#'
#' @description getminy generates a lower bound for a plot where it is unknown
#' whether the minumum is less than zero of not. If less than 0 then
#' multiplying by the default mult of 1.05 works well but if the outcome if
#' > 0 then the multiplier needs to be adjusted appropriately so the minimum
#' is slightly lower than the minimum of the data
#'
#' @param x the vector of data to be tested for its minimum
#' @param mult the multiplier for both ends, defaults to 1.05 (=0.95 if >0)
#'
#' @return a suitable lower bound for a plot if required
#' @export
#'
#' @examples
#' vect <- rnorm(10,mean=0,sd=2)
#' sort(vect)
#' getminy(vect,mult=1.0)
getminy <- function(x,mult=1.05) {
ymin <- min(x,na.rm=TRUE)
if (ymin < 0) {
ymin <- ymin * mult
} else {
ymin <- ymin * (2 - mult)
}
return(ymin)
} # end of getminy
#' @title getmaxy generates the upper bound for a plot
#'
#' @description getmaxy generates an upper bound for a plot where it is unknown
#' whether the maximum is greater than zero of not. If > 0 then
#' multiplying by the default mult of 1.05 works well but if the outcome if
#' < 0 then the multiplier needs to be adjusted appropriately so the maximum
#' is slightly higher than the maximum of the data
#'
#' @param x the vector of data to be tested for its maximum
#' @param mult the multiplier for both ends, defaults to 1.05 (=0.95 if < 0)
#'
#' @return a suitable upper bound for a plot if required
#' @export
#'
#' @examples
#' \dontrun{
#' vect <- rnorm(10,mean=0,sd=2)
#' sort(vect,decreasing=TRUE)
#' getmaxy(vect,mult=1.0)
#' vect <- rnorm(10,mean = -5,sd = 1.5)
#' sort(vect,decreasing=TRUE)
#' getmaxy(vect,mult=1.0)
#' }
getmaxy <- function(x,mult=1.05) {
ymax <- max(x,na.rm=TRUE)
if (ymax > 0) {
ymax <- ymax * mult
} else {
ymax <- ymax * (2 - mult)
}
return(ymax)
} # end of getmaxy
#' @title getsingle extracts a single number from an input line of characters
#'
#' @description getsingle splits up a text line and translates the first non-
#' empty character string into a number.
#'
#' @param inline the line of text, usually taken after using readLines
#' @param sep the separator used to divide the numbers from descriptive text.
#' defaults to a comma.
#'
#' @return a single number
#' @export
#'
#' @examples
#' \dontrun{
#' x <- "12.3 , this is a number"
#' y <- "21.3 # 22.3 # here are two numbers"
#' getsingle(x)
#' getsingle(y,sep="#")
#' }
getsingle <- function(inline,sep=",") { # inline=dat[41]
tmp <- unlist(strsplit(inline,sep))
tmp <- gsub(" ","",tmp)
tmp <- tmp[nchar(tmp) > 0]
return(as.numeric(tmp[1]))
}
#' @title getvector extracts a vector of numbers from a line of characters
#'
#' @description getvector when reading in a csv file using readLines,
#' getvector extarcts a line of numbers from a specified line within
#' the readLine object.This function works out how many numbers there
#' are. If you wish to add a comment at the end of a vector of numbers
#' it must be separated from tehm by the separator. e.g. a comma
#' @param indat the readLines object
#' @param locate the line number from which to extract the numbers
#' @param sep the separator between numbers, defaults to ","
#'
#' @return a vector of numbers
#' @export
#'
#' @examples
#' \dontrun{
#' x <- "12.3, 15.1, 8.7,10.3, # this is a vector of numbers"
#' y <- "21.3 # 22.3 # 8.7 # 10.3 # here are four numbers"
#' getvector(x)
#' getvector(y,sep="#")
#' }
getvector <- function(indat,locate,sep=",") { # indat=dat; locate=pick+2;sep=","
vect <- indat[locate]
if (length(grep("\t",vect) > 0)) vect <- gsub("\t",sep,vect)
vect <- unlist(strsplit(vect,sep))
vect <- gsub(" ","",vect)
vect1 <- vect
vect <- suppressWarnings(as.numeric(vect))
vect <- vect[nchar(vect) > 0]
if (!any(vect1 == "NA")) {
vect <- vect[!is.na(vect)]
}
return(vect)
}
#' @title greplow - uses tolower in the search for the pattern
#'
#' @description greplow - a grep implementation that ignores the case of either
#' the search pattern or the object to be search. Both are converted to lower
#' case before using grep.
#' @param pattern - the text to search for in x
#' @param x - the vector or object within which to search for 'pattern' once
#' both have been converted to lowercase.
#'
#' @return the index location within x of 'pattern', if it is present, an empty
#' integer if not
#' @export greplow
#'
#' @examples
#' txt <- c("Long","Lat","LongE","LatE","Depth","Zone","Effort","Method")
#' greplow("zone",txt)
#' greplow("Zone",txt)
#' greplow("long",txt)
greplow <- function(pattern,x) {
return(grep(tolower(pattern),tolower(x)))
}
#' @title histyear plots a histogram of a variable for each year available
#'
#' @description histyear plots a histogram of a variable for each year
#' available.
#'
#' @param x the data.frame of data with at least a 'Year' and pickvar present
#' @param Lbound leftbound on all histograms, defaults to -3.5
#' @param Rbound right bound on all histograms, defaults to 12.25
#' @param inc the class width of the histogram, defaults to 0.25
#' @param new reuse a plotting device or output a new one
#' @param plots the number of rows and columns, defaults to c(11,3)
#' @param pickvar which variable in the input data.frame to plot each year
#' default = 'LnCE'
#' @param varlabel what label to use on x-axis, default = 'log(CPUE)'
#' @param col colour of each cell; defaults to 2 (red)
#' @param border colour of the border of each cell, default = 2 (red)
#' @param addnormal include the fitting of a normal distribution
#' @param right put the legend on the right hand side, default = TRUE
#' @param vline allows for the placing of a vertical line on each histogram
#' if set to "" then nothing is added. defaults to the average.
#'
#' @return a matrix of the year, mean value, stdev, and N number of
#' observations. It also plots a histogram for each year and fits a
#' normal distribution to each one.
#' @export
#'
#' @examples
#' \dontrun{
#' print("still to be developed")
#' }
histyear <- function(x,Lbound=-3.5,Rbound=12.25,inc=0.25,new=FALSE,
plots=c(11,3),pickvar="LnCE",varlabel="log(CPUE)",
col=2,border=2,addnormal=TRUE,right=TRUE,vline="mean") {
yrs <- sort(unique(x[,"Year"]))
nyr <- length(yrs)
columns <- c("Year","maxcount","Mean","StDev","N")
results <- matrix(0,nrow=nyr,ncol=length(columns),dimnames=list(yrs,columns))
rgeCE <- range(x[,pickvar],na.rm=TRUE)
avCE <- mean(x[,pickvar],na.rm=TRUE)
if (rgeCE[1] < Lbound) Lbound <- rgeCE[1] - inc
if (rgeCE[2] > Rbound) Rbound <- rgeCE[2] + inc
if (new) plotprep(width=6,height=9)
par(mfcol=plots,mai=c(0.225,0.25,0.025,0.05),oma=c(1.2,1.0,0.0,0.0))
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
for (yr in 1:nyr) {
pick <- which(x[,"Year"] == yrs[yr])
outh <- hist(x[pick,pickvar],breaks=seq(Lbound,Rbound,inc),col=col,
border=border,main="",xlab="",ylab="",panel.first=grid())
if (vline == "mean") abline(v=avCE,col=4,lwd=2)
if (is.numeric(vline)) abline(v=vline,col=4,lwd=2)
adjust <- 0
if (right) adjust <- 1
mtext(paste0(" ",yrs[yr]),side=3,outer=F,line=-1.5,font=7,cex=0.8,adj=adjust)
mtext(paste0(" ",length(pick)),side=3,outer=F,line=-2.5,font=7,cex=0.8,adj=adjust)
if (addnormal) {
pickmax <- which.max(outh$counts)
ans <- addnorm(outh,x[pick,pickvar])
lines(ans$x,ans$y,col=4,lwd=2)
results[yr,] <- c(yrs[yr],outh$mids[pickmax],ans$stats)
}
}
mtext("Frequency",side=2,outer=T,line=0.0,font=7,cex=1.0)
mtext(varlabel,side=1,outer=T,line=0.0,font=7,cex=1.0)
return(invisible(results))
} # end of histyear
#' @title incol is a utility to determine is a column is present in a matrix
#'
#' @description incol is a utility to determine whether a names columns is
#' present in a given matrix or data.frame.
#'
#' @param incol the name of the column; defaults to "year" as an example
#' @param inmat the matrix or data.frame within which to search for incol
#'
#' @return TRUE or FALSE
#' @export
#'
#' @examples
#' \dontrun{
#' test <- matrix(c(1,2,3,4),nrow=2,ncol=2,dimnames=list(1:2,c("year","Catch")))
#' print(test)
#' iscol("year",test)
#' iscol("Catch",test)
#' iscol("catch",test)
#' iscol("ages",test)
#' }
iscol <- function(incol="year",inmat) { # incol="ages"; inmat=dat
if (length(grep(incol,colnames(inmat))) < 1) return(FALSE)
else return(TRUE)
}
#' @title magnitude returns the magnitude of numbers
#'
#' @description magnitude is useful when using an
#' optimizer such as optim, which uses a parscale parameter.
#' magnitude can determine the respective parscale value for each
#' parameter value.
#'
#' @param x the vector of numbers (parameters) whose magnitudes are
#' needed
#'
#' @return a vector of magnitudes
#' @export
#'
#' @examples
#' \dontrun{
#' x <- c(0,0.03,0.3,3,30,300,3000)
#' magnitude(x)
#' }
magnitude <- function(x) {
return(10^(floor(log10(abs(x)))))
}
#' @title makedataset combines year, catch, and optionally cpue into a data set
#'
#' @description makedataset takes the components of a simplified data set and
#' combines them into a standard structure suitable for use with datalowSA.
#' It requires a vector called year, one called catch, optionally one called
#' cpue (all vectors need to be the same length with gaps filled with NA),
#' the species name, and its resilience. It returns a list containing
#' all the required parts.
#'
#' @param year a vector of years that contained fisheries data
#' @param catch a vector of catches in each year
#' @param cpue either a single NA (default) or a vector of cpue. Any gaps in the
#' cpue time-series should be filled with NA
#' @param species the species name in quotes
#' @param resilience either verylow, low, medium, or high
#'
#' @return a list of fish, glb, props, agedata, and lendata, the last 3 NULL
#' @export
#'
#' @examples
#' \dontrun{
#' year <- 1986:2016
#' catch <- c(112.9,206.3,95.7,183.1,147.4,198.9,102.1,235.5,247.8,426.8,448,577.4,
#' 558.5,427.9,509.3,502.4,429.6,360.2,306.2,195.7,210,287.3,214.2,260.6,
#' 272.2,356.9,345,282.7,285.1,237.8,233.3)
#' cpue <- c(1.2006,1.3547,1.0585,1.0846,0.9738,1.0437,0.7759,1.0532,1.284,1.3327,
#' 1.4014,NA,NA,1.142,0.9957,0.8818,0.7635,0.7668,0.7198,0.5997,0.6336,
#' 0.6936,0.8894,0.8644,0.8442,0.8427,0.8849,0.9964,0.9804,0.957,1.0629)
#' dat <- makedataset(year,catch,cpue,"dataspmbits","low")
#' dat
#' }
makedataset <- function(year, catch, cpue=NA, species="species",resilience="low") {
if (length(year) != length(catch))
stop("year vector different length to catch vector \n")
if (length(cpue) == 1) cpue <- rep(NA,length(year))
fish <- as.data.frame(cbind(year,catch,cpue),row.names=year)
glb <- list(spsname=species,resilience=resilience)
return(list(fish=fish,glb=glb,props=NULL,agedata=NULL,lendata=NULL))
} # end of makedataset
#' @title multinomLL is the negative log-likelihood for multinomial values
#'
#' @description multinomLL calculates the negative log-likelihood for the
#' multinomial distribution. Commonly used with length and age composition
#' data in assessment models. Zero observations in particular categories
#' do not matter as the equation is - Sum(Obs x log(expected-proportion)).
#' Thus zeros in the observations just contribute zero to the sum.
#'
#' @param counts the observed counts in each category (age or length) that is to
#' be fitted by the model
#' @param predprop the predicted proportions in each category. One needs a
#' separate function to calculate this.
#'
#' @return a single scalar value, the -veLL
#' @export
#'
#' @examples
#' \dontrun{
#' obs <- c(0,5,11,19,12,3,0,0)
#' expprop <- c(0.018,0.002,0.12,0.4,0.3,0.12,0.03,0.01)
#' multinomLL(obs,expprop) # should be 92.61393
#' }
multinomLL <- function(counts,predprop) {
LL <- -sum(counts*log(predprop)) #
return(LL)
}
#' @title outoptim provides a sucinct summary of the optim solution
#'
#' @description outoptim provides a more succinct printout of the optim solution
#' which only uses five lines on the screen rather than the standard print
#' of the output list structure
#'
#' @param inopt an object into which the output from optim is directed.
#'
#' @return nothing, but it does print the optim result to the console in five
#' lines without the list structure
#' @export
#'
#' @examples
#' \dontrun{
#' data(fishdat)
#' bestL <- optim(c(14.0,0.3),aspmLL,method="Nelder-Mead",infish=fishdat$fish,
#' inglb=fishdat$glb,inprops=fishdat$props,
#' control=list(maxit=1000,parscale=c(10,0.1)))
#' outoptim(bestL)
#' }
outoptim <- function(inopt){
cat("$par : ",inopt$par,"\n")
cat("$value : ",inopt$value,"\n")
cat("$counts : ",inopt$counts," iterations, gradient\n")
cat("$convergence : ",inopt$convergence,"\n")
cat("$message : ",inopt$message,"\n")
} # end of outoptim
#' @title parsyn types the standard syntax for the par command to the console
#'
#' @description parsyn types the standard syntax for the par command to the
#' console so it can be copied and pasted into your own code and modified.
#'
#' @return it writes two lines of R code to the console
#' @export
#'
#' @examples
#' parsyn()
parsyn <- function() {
cat("par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) \n")
cat("par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) \n")
}
#' @title penalty0 enables the adding of a large penalty as one approaches 0.0
#'
#' @description penalty0 allows for the option of adding a large penalty as
#' a parameter approaches 0.0 . See spmLL
#' for example code that contains such a parameter. For example, when
#' fitting an spm sometimes the optimal mathematical model fit can occur
#' by depressing the r value to 0 or even go negative. This is only
#' used internally to datalowSA and so is not formally exported.
#'
#' @param x the parameter value that potentially incurs a penalty
#'
#' @return a single value as a penalty to be added to a Log-Likelihood or SSQ
#'
#' @examples
#' datalowSA:::penalty0(0.5)
#' datalowSA:::penalty0(0.1)
#' datalowSA:::penalty0(0.01)
#' datalowSA:::penalty0(0.005)
penalty0 <- function(x){
ans <- 500*exp(-1000*x)
return(ans)
} # end of penalty0
#' @title penalty1 enables the adding of a large penalty as one approaches 1.0
#'
#' @description penalty1 allows for the option of adding a large penalty as
#' a parameter approaches 1.0 and moves to become larger than 1. See aspmLL
#' for example code that contains such a parameter. For example, when
#' fitting an ASPM sometimes the optimal mathematical model fit can occur
#' by depressing the R0 and having an initial depletion >>1. This is only
#' used internally to datalowSA and so is not formally exported.
#'
#' @param x the parameter value that potentially incurs a penalty
#'
#' @return a single value as a penalty to be added to a Log-Likelihood or SSQ
#'
#' @examples
#' datalowSA:::penalty1(0.5)
#' datalowSA:::penalty1(0.9)
#' datalowSA:::penalty1(0.99)
penalty1 <- function(x){
ans <- 50 * (abs((1 - abs(x) - 0.5))/0.5)^45
return(ans)
} # end of penalty1
#' @title plot1 a simple way to plot an xy line plot
#'
#' @description plot1 provides a quick way to plot out a single xy
#' line plot. It can be used with plotprep to generate a plot
#' outside of Rstudio or by itself to generate one within Rstudio.
#' It uses a standard par setup and permits custom labels, font,
#' and font size (cex). It checks the spread of y and if a ymax is
#' not given in the parameters finds the ymax and checks to see if
#' y goes negative in which case it uses getmin, so the
#' y-axis is set to 0 - ymax or ymin - ymax
#'
#' @param x The single vector of x data
#' @param y the single vector of y data. If more are required they can
#' be added spearately after calling plot1.
#' @param xlabel the label fot the x-axis, defaults to empty
#' @param ylabel the label fot the y-axis, defaults to empty
#' @param type the type of plot "l" is for line, the default, "p" is
#' points. If you want both plot a line and add points afterwards.
#' @param usefont which font to use, defaults to 7 which is Times bold
#' @param cex the size of the fonts used. defaults to 0.85
#' @param maxy defaults to 0, which does nothing. If a value is given
#' then this value is used rather than estimating from the input y
#' @param defpar if TRUE then plot1 will declare a par statement. If false it
#' will expect one outside the function. In this way plot1 can be
#' used when plotting multiple graphs, perhaps as mfrow=c(2,2)
#' @param inpch the input character type if using type="p", default=16
#' @param incol the colour to use for the line or points, default = black
#'
#' @return nothing but it does plot a graph and changes the par setting
#' @export
#'
#' @examples
#' \dontrun{
#' x <- rnorm(20,mean=5,sd=1)
#' plot1(x,x,xlabel="x-values",ylabel="yvalues")
#' }
plot1 <- function(x,y,xlabel="",ylabel="",type="l",usefont=7,cex=0.85,
maxy=0,defpar=TRUE,inpch=16,incol=1){
if (defpar) {
par(mfrow = c(1,1), mai = c(0.45,0.45,0.1,0.05),oma = c(0,0,0,0))
par(cex = cex, mgp = c(1.35, 0.35, 0), font.axis = usefont,
font = usefont, font.lab = usefont)
}
if (maxy > 0) ymax <- maxy else ymax <- getmaxy(y)
if (min(y,na.rm=TRUE) < 0.0) ymin <- getminy(y) else ymin <- 0.0
addline <- FALSE
plot(x,y,type=type,pch=inpch,lwd=2,col=incol,ylim=c(ymin,ymax),yaxs="i",
ylab=ylabel,xlab=xlabel,cex=cex,panel.first=grid())
} # end of plot1
#' @title plotfish plots the catch and optionally the cpue from fish
#'
#' @description plotfish uses the matrix of fishery data used in the
#' datalowSA standard data format. It requires the matrix or data.frame
#' to contain the columns 'year', 'catch', and optionally 'cpue'.
#'
#' @param fish the matrix or data.frame containing year, catch, and cpue.
#' @param glb the list of biologicals potentially containing the spsname
#' @param ce a logical parameter determining whether to plot the cpue or not.
#' the default = TRUE
#' @param title determines whether or not the spsname is printed at the top
#' of the plot. Default = TRUE but for a more formal publication it
#' might need to be set to FALSE, which also reallocates the room
#' given to the title to the plot.
#' @param fnt the font used in the plot and axes.
#' @param filename default is empty. If a filename is put here a .png file
#' with that name will be put into the working directory.
#' @param resol the resolution of the png file, defaults to 200 dpi
#'
#' @return prints the location of the png file produced to the console
#' @export
#'
#' @examples
#' \dontrun{
#' data(dataspm)
#' plotfish(dataspm$fish,glb,ce=TRUE,filename="")
#' }
plotfish <- function(fish,glb,ce=TRUE,title=TRUE,fnt=7,filename="",
resol=200) {
colnames(fish) <- tolower(colnames(fish))
rows <- 1
if (ce) rows <- 2
lenfile <- nchar(filename)
if (lenfile > 3) {
end <- substr(filename,(lenfile-3),lenfile)
if (end != ".png") filename <- paste0(filename,".png")
ifelse(rows == 1,hgt <- 3.0,hgt <- 3.5)
png(filename=filename,width=5.0,height=hgt,units="in",res=resol)
}
yrs <- fish[,"year"]
if (title) {
par(mfrow=c(rows,1),mai=c(0.5,0.45,0.1,0.05),oma=c(0.0,0,1.0,0.0))
} else {
par(mfrow=c(rows,1),mai=c(0.5,0.45,0.1,0.05),oma=c(0.0,0,0.0,0.0))
}
par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=fnt,font=fnt,font.lab=fnt)
ymax <- getmaxy(fish[,"catch"])
plot(yrs,fish[,"catch"],type="l",lwd=2,ylab="Catch",xlab="Year",
ylim=c(0,ymax),yaxs="i",panel.first = grid())
if (title)
mtext(glb$spsname,side=3,cex=1.0,line=0,font=fnt,outer=TRUE)
if (ce) {
ymax <- getmaxy(fish[,"cpue"])
plot(yrs,fish[,"cpue"],type="l",lwd=2,ylab="CPUE",xlab="Year",
ylim=c(0,ymax),yaxs="i",panel.first = grid())
}
if (lenfile > 0) {
outfile <- paste0(getwd(),"/",filename)
print(outfile)
dev.off()
}
} # end of plotfish
#' @title plotfishery plots the catch and optionally the cpue from fish
#'
#' @description plotfishery is now deprecated and only calls plotfish.
#'
#' @param fish the matrix or data.frame containing year, catch, and cpue.
#' @param glb the list of biologicals potentially containing the spsname
#' @param ce a logical parameter determining whether to plot the cpue
#' or not. the default = TRUE
#' @param title determines whether or not the spsname is printed at the top
#' of the plot. Default = TRUE but for a more formal publication it
#' might need to be set to FALSE, which also reallocates the room
#' given to the title to the plot.
#' @param fnt the font used in the plot and axes.
#' @param filename default is empty. If a filename is put here a .png file
#' with that name will be put into the working directory.
#' @param resol the resolution of the png file, defaults to 200 dpi
#'
#' @return prints the location of the png file produced to the console
#' @export
plotfishery <- function(fish,glb,ce=TRUE,title=TRUE,fnt=7,filename="",
resol=200) {
plotfish(fish,glb,ce,title,fnt,filename,resol)
}
#' @title plotprep: sets up a window and the par values for a single plot
#'
#' @description plotprep: sets up a window and the par values for a single plot.
#' it checks to see if a graphics device is open and opens a new one if not.
#' This is simply a utility function to save typing the standard syntax.
#' Some of the defaults can be changed. Typing the name without () will
#' provide a template for modification. If 'windows' is called repeatedly this
#' will generate a new active graphics device each time leaving the older ones
#' inactive but present. For quick exploratory plots this behaviour is not
#' wanted, hence the check if an active device exists already or not.
#' @param width defaults to 6 inches = 15.24cm - width of plot
#' @param height defaults to 3 inches = 7.62cm - height of plot
#' @param plots defaults to c(1,1), but arranges multiple plots. If used it may
#' be necessary to print out this code and adjust the mai and oma variables
#' @param usefont default is 7 (bold Times); 1 = sans serif, 2 = sans serif bold
#' @param cex default is 0.85, the size of font used for text within the plots
#' @param xmtext default is TRUE; if plots is not c(1,1) this alters the mai and
#' oma variables for the x-axis to allow for mtexting and avoid the x title
#' @param ymtext default is TRUE; if plots is not c(1,1) this alters the mai and
#' oma variables for the y-axis to allow for mtexting and avoid the y title
#' @param newdev reuse a previously defined graphics device or make a new one;
#' defaults to TRUE
#' @param rows defaults to TRUE, determines whether to use mfrow or mfcol
#' @param filename defaults to "" = do not save to a filename. If a filename is
#' @return Checks for and sets up a graphics device and sets the default plotting
#' par values. This changes the current plotting options!
#' @export plotprep
#' @examples
#' x <- rnorm(1000,mean=0,sd=1.0)
#' plotprep()
#' hist(x,breaks=30,main="",col=2)
plotprep <- function(width=6,height=3.6,plots=c(1,1),usefont=7,cex=0.85,
xmtext=TRUE,ymtext=TRUE,
newdev=TRUE,rows=TRUE,filename="") {
if ((names(dev.cur()) != "null device") & (newdev)) suppressWarnings(dev.off())
lenfile <- nchar(filename)
if (lenfile > 3) {
end <- substr(filename,(lenfile-3),lenfile)
if (end != ".png") filename <- paste0(filename,".png")
png(filename=filename,width=width,height=height,units="in",res=300)
} else {
if (names(dev.cur()) %in% c("null device","RStudioGD"))
dev.new(width=width,height=height,noRStudioGD = TRUE)
}
firstmai <- 0.45; secondmai <- 0.45; thirdmai <- 0.1
firstoma <- 0.0; secondoma <- 0.0; thirdoma <- 0.0
if (sum(plots) != 2) {
if (xmtext) {
firstmai <- 0.25
thirdmai <- 0.05
firstoma <- 1.0
thirdoma <- 0.1
}
if (ymtext) {
secondmai <- 0.25
secondoma <- 1.0
}
}
maival <- c(firstmai,secondmai,thirdmai,0.05)
omaval <- c(firstoma,secondoma,thirdoma,0.0)
if (rows) {
par(mfrow = plots,mai=maival,oma=omaval)
} else {
par(mfcol = plots,mai=maival,oma=omaval)
}
par(cex=cex, mgp=c(1.35,0.35,0), font.axis=usefont,font=usefont,font.lab=usefont)
if (lenfile > 0) cat("\n Remember to place 'graphics.off()' after the plot \n")
} # end of plotprep
#' @title quants used in apply to estimate quantiles across a vector
#'
#' @description quants used in 'apply' to estimate quantiles across a vector
#' @param invect vector of values
#' @return a vector of the c(0.025,0.05,0.1,0.5,0.9,0.95,0.975) quantiles
#' @export quants
#' @examples
#' \dontrun{
#' x <- rnorm(100,mean=5,sd=1)
#' quants(x)
#' }
quants <- function(invect) {
ans <- quantile(invect,probs = c(0.025,0.05,0.5,0.9,0.95,0.975),na.rm=T)
return(ans)
}
#' @title readdata reads in a standard format data file
#'
#' @description readdata reads in a standard format data file. An example of the
#' standard format is generated by using the function dataTemplate, which
#' generates an example datafile which can be used to demonstrate the
#' methods present in datalowSA, or can be used as a template to edit and
#' input a new or different dataset.
#'
#' @param filename the filename (including the full path if required) containing
#' the data in the standard format.
#' @param property does the data include a table of length-at-age, maturity-at-
#' age, weight-at-age, fecundity-at-age, and that kind of thing (see the
#' standard format for a list of data.)
#' @param verbose Default = TRUE, which prints out details as data is read in.
#'
#' @return a list of three objects, fish includes the Year, Catch, CPUE, and SE
#' of the CPUE, glb containing an array of biological properties that can
#' be global, and props containing age, laa, waa, maa, sela.
#' @export
#'
#' @examples
#' \dontrun{
#' dataTemplate(filename="test.csv", title="A test of the functions")
#' ans <- readdata("test.csv",property=FALSE)
#' str(ans)
#' } # filename="sardine.csv";
readdata <- function(filename,property=FALSE,verbose=TRUE) { # property=FALSE; filename=paste0(datadir,"eastdeepwatershark.csv")
dat <- readLines(filename)
spsname <- scan(file=filename,skip=0,what=character(),nlines = 1,quiet=TRUE)
spsname <- removeEmpty(gsub(",","",spsname))
pick <- grep("RESILIENCE",dat)[1]
resilience <- scan(file=filename,skip=pick,what=character(),nlines = 1,quiet=TRUE)
resilience <- removeEmpty(gsub(",","",resilience))
pick <- grep("NYRS",dat)[1] + 1
nyrs <- getsingle(dat[pick])
if (verbose) cat(spsname,"\n\n resilience = ",resilience,
" Num Years = ",nyrs,"\n\n")
# get the fishery data into as many columns as required
pick <- grep("YEARS",dat)[1]
columns <- tolower(removeEmpty(unlist(strsplit(dat[pick],","))))
numcols <- length(columns)
columns <- c("year",columns[2:numcols])
skips <- pick:(pick+nyrs-1)
fish <- as.data.frame(matrix(NA,nrow=nyrs,ncol=length(columns),
dimnames=list(1:nyrs,columns)))
for (i in 1:nyrs) { # i = 1
fish[i,] <- scan(file=filename,skip=skips[i],sep=",",nlines = 1,quiet=TRUE)[1:numcols]
}
if (verbose) {
cat("fish \n\n")
print(head(fish,10))
cat("\n")
}
glb=list(resilience=resilience,spsname=spsname)
pick <- grep("BIOLOGY",dat)[1] + 1
if (is.na(pick)) {
props <- NULL
} else {
maxage <- getsingle(dat[pick])
M <- getsingle(dat[pick+1])
Linf <- getsingle(dat[pick+2])
K <- getsingle(dat[pick+3])
t0 <- getsingle(dat[pick+4])
Waa <- getsingle(dat[pick+5])
Wab <- getsingle(dat[pick+6])
M50a <- getsingle(dat[pick+7])
deltaM <- getsingle(dat[pick+8])
sela50 <- getsingle(dat[pick+9])
deltaS <- getsingle(dat[pick+10])
steep <- getsingle(dat[pick+11])
R0 <- getsingle(dat[pick+12])
ages <- 0:maxage
nages <- length(ages)
glb <- list(maxage=maxage,M=M,
Linf=Linf, K=K, t0=t0,
Waa=Waa, Wab=Wab,
M50a=M50a, deltaM=deltaM,
steep=steep,R0=R0,
sela50=sela50, deltaS=deltaS,
resilience=resilience,
nages=nages,ages=ages,nyrs=nyrs,spsname=spsname)
columns <- c("age","laa","waa","maa","sela","feca")
props <- as.data.frame(matrix(NA,nrow=nages,ncol=length(columns),
dimnames=list(ages,columns)))
if (property) {
pick <- grep("PROPERTY",dat)[1] + 1
for (i in 1:nages) {
props[i,] <- getvector(dat,pick,sep=",")
pick <- pick + 1
}
} else {
# now calculate the properties
laa <- vB(c(Linf,K,t0),ages)
waa <- (Waa * laa ^ Wab)/1000
maa <- logist(M50a,deltaM,ages)
sela <- logist(sela50,deltaS,ages)
feca <- sela * waa
props <- as.data.frame(cbind(ages,laa,waa,maa,sela,feca))
}
}
if (verbose) {
cat("biology \n")
print(unlist(glb))
cat("\n properties \n")
print(head(props,10))
}
# dat <- readLines(filename)
pick <- grep("AGE",dat)[1] + 1
if (is.na(pick)) {
agedata <- NULL
} else {
yrage <- getsingle(dat[pick])
numsex <- getsingle(dat[(pick+1)])
ages <- getvector(dat,(pick+2),sep=",")
agemax <- max(ages)
nage <- length(ages)
naa <- matrix(0,nrow=(yrage*numsex),ncol=(nage+2))
colnames(naa) <- c("year","sex",ages)
pick <- pick+2 # i = 1
for (i in 1:(yrage*numsex)) naa[i,] <- getvector(dat,(pick+i),sep=",")
rownames(naa) <- naa[,1]
agedata <- list(yrage=yrage,ages=ages,agemax=agemax,nage=nage,naa=naa)
}
pick <- grep("LENGTH",dat)[1] + 1
if (is.na(pick)) {
lendata <- NULL
} else {
yrlen <- getsingle(dat[pick])
numsexl <- getsingle(dat[(pick+1)])
lengths <- getvector(dat,(pick+2),sep=",")
maxlen <- max(lengths)
nlength <- length(lengths)
if (verbose) {
cat("Number of years of Data: ",yrlen,"\n")
cat("Number of sexes with Data: ",numsexl,"\n")
cat("Length classes: ",lengths,"\n")
cat("Number of length classes ",nlength,"\n")
}
nal <- matrix(0,nrow=(yrlen*numsexl),ncol=(nlength+2))
colnames(nal) <- c("year","sex",lengths)
pick <- pick+2
for (i in 1:(yrlen*numsexl)) nal[i,] <- getvector(dat,(pick+i),sep=",")
rownames(nal) <- nal[,1]
lendata <- list(yrlen=yrlen,lengths=lengths,maxlen=maxlen,
nlength=nlength,nal=nal)
}
ans <- list(fish=fish,glb=glb,props=props,agedata=agedata,lendata=lendata)
return(ans)
} # end of readdata
#' @title removeEmpty removes empty strings from a vector of strings
#'
#' @description removeEmpty removes empty strings from a vector of strings.
#' Such spaces often created by spurious commas at the end of lines. It
#' also removes strings made up only of spaces and removes spaces from
#' inside of inidivdual chunks of text.
#'
#' @param invect a vector of input strings, possibly containing empty strings
#'
#' @return a possibly NULL vector of strings
#' @export
#'
#' @examples
#' x <- c("1","","2",""," ","3"," ","4","","a string","end")
#' x
#' length(x)
#' length(removeEmpty(x))
#' removeEmpty(x)
removeEmpty <- function(invect) {
tmp <- gsub(" ","",invect)
tmp <- tmp[nchar(tmp) > 0]
return(tmp)
}
#' @title toXL copies a data.frame or matrix to the clipboard
#'
#' @description toXL copies a data.frame or matrix to the clipboard
#' so one can then switch to Excel and just type <ctrl> + V to paste the
#' data in place
#'
#' @param x a vector or matrix
#' @param output - a boolean determining whether to print the object to the
#' screen as well as the clipboard; defaults to FALSE
#' @return Places the object 'x' into the clipboard ready for pasting
#' @export toXL
#' @examples
#' \dontrun{
#' datamatrix <- matrix(data=rnorm(100),nrow=10,ncol=10)
#' colnames(datamatrix) <- paste0("A",1:10)
#' rownames(datamatrix) <- paste0("B",1:10)
#' toXL(datamatrix,output=TRUE)
#' }
toXL <- function(x,output=FALSE) {
write.table(x,"clipboard",sep="\t")
if(output) print(x)
}
#' @title vB calculates the predicted von Bertalanffy length at age
#'
#' @description vB calculates length at age for the von Bertalanffy curve.
#'
#' @param par is a vector the first three cells of which are Linf, K, and t0
#' for the VB curve; the fourth parameter will be sigma, the standard
#' deviation of the normal likelihoods used with the residuals
#' @param ages is a vector of ages
#'
#' @return a vector of predicted lengths for the vector of ages in 'ages'
#' @export
#'
#' @examples
#' ages <- seq(0,20,1)
#' pars <- c(Linf=50,K=0.3,t0=-1.0,sigma=1.0) # Linf, K, t0, sigma
#' cbind(ages,vB(pars,ages))
vB <- function(par,ages) {
return(par[1] * (1 - exp(-par[2]*(ages-par[3]))))
}
#' @title which.closest find the number closest to a given value
#'
#' @description which.closest finds either the number in a vector which is
#' closest to the input value or its index value
#'
#' @param x the value to lookup
#' @param invect the vector in which to lookup the value x
#' @param index should the closest value be returned or its index; default=TRUE
#'
#' @return by default it returns the index in the vector of the value closest to
#' the input value
#' @export
#'
#' @examples
#' vals <- rnorm(100,mean=5,sd=2)
#' pick <- which.closest(5.0,vals,index=TRUE)
#' pick
#' vals[pick]
#' vals[c(pick-5):(pick+5)]
#' which.closest(5.0,vals,index=FALSE)
which.closest <- function(x,invect,index=T) {
pick <- which.min(abs(invect-x))
if (index) {
return(pick)
} else {
return(invect[pick])
}
} # end of which.closest
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.