#' Error Check the Adult Index Data
#'
#' Check the adult sea lamprey trapping data (collected for estimation of
#' the Adult Index) for errors.
#' @param streamDat
#' A data frame of old and new stream mark-recapture estimates
#' used to estimate the lake-wide Adult Indices,
#' typically the output from \code{\link{AIprep}}. The data frame must
#' include: \code{year},
#' \code{lake}, lake-stream ID \code{lscode}
#' (see details), population estimate
#' \code{PEmr}, coefficient of variation \code{CVmr}
#' (100\% * sqrt(variance(PEmr)) / PEmr), \code{index}, a logical
#' identifying the index streams; \code{maintain} a logical identifying the
#' streams that will continue to have ongoing trapping even if not part of
#' the Adult Index; \code{indexContrib} a numeric, the stream population
#' estimate that will be used in the Adult Index (NA for new);
#' \code{indexContribCV} a numeric, the stream CV that will be used to
#' generate 95\% confidence intervals for the Adult Index (NA for new); and
#' \code{complete} a logical identifying streams and years for which the
#' Adult Index has already been estimated (FALSE for new).
#' @param csvDir
#' A character scalar identifying the path where the rtf file will be
#' stored, e.g., \code{csvDir = "C:\\temp\\mydir"}.
#' @param outFile
#' Name of the output rtf file, default NULL, in which case the file will be
#' named "YYYY Adult Index - error checking.doc" where YYYY
#' is the latest year represented in \code{streamDat}.
#' @param otherTabs
#' A list of other tables to be printed in error check report, default NULL.
#' The list names will be used as captions.
#' @details
#' Lake-stream IDs are combination of lake ID and stream ID
#' e.g., 1.064 = lake ID + (stream ID)/1000.
#' @return
#' An error checking document as an rtf file (with the file type *.doc,
#' so that MS Word will open it automatically).
#' @export
#'
AIcheck <- function(streamDat, csvDir, outFile=NULL, otherTabs=NULL) {
# library(GLFC)
# streamDat=stream1
# csvDir=DIRECTORY
# outFile=NULL
# otherTabs=othtabs
YEAR <- max(streamDat$year)
if (is.null(outFile)) {
outFile <- paste(YEAR, "Adult Index - error checking.doc")
}
# create a file with error checking results
doc <<- startrtf(file=outFile, dir=csvDir)#, width=11, height=8.5)
heading(paste("Error Checking the", YEAR, "Adult Sea Lamprey Data"))
heading(format(Sys.time(), "%Y %b %d %H:%M:%S"), 2)
# start by printing any "other tables" first
if (!is.null(otherTabs)) {
for(i in 1:length(otherTabs)) {
tabl(names(otherTabs)[i], TAB=otherTabs[[i]], row.names=FALSE)
}
}
# create a nice looking, formatted data frame for printing out errors ...
df.nice <- prettytable(
streamDat[, c("lake", "country", "strcode", "strname", "year", "trapcatch",
"PEmr", "CVmr")], 0, bigseps=c("", ",")[c(2, 2, 2, 2, 1, 2, 2, 2)])
with(streamDat, {
# show streams that are missing a current year population estimate
tab <- df.nice[year==YEAR & (is.na(PEmr) | is.na(CVmr)), ]
if (dim(tab)[1]>0) {
tabl("NOTE: Some streams were missing a ", YEAR,
" population estimate or CV.", TAB=tab, row.names=FALSE)
} else {
para("GREAT: All streams had ", YEAR, "population estimates and CVs.")
}
# check paired nature of mark-recapture estimates and trap catches
tab <- df.nice[((is.na(PEmr) & !is.na(CVmr)) |
(!is.na(PEmr) & is.na(CVmr)) | (is.na(trapcatch) & !is.na(PEmr))) &
lscode!=3.999, ]
if (dim(tab)[1]>0) {
tabl("ERROR: Something's missing ... either, trapcatch, PEmr, or CVmr.",
TAB=tab, row.names=FALSE)
} else para("OKAY: All mark-recapture estimates have trap catches and",
" CVs.")
# check values of trap catch and population estimate
tab <- df.nice[!is.na(PEmr) & !is.na(trapcatch) & PEmr<=trapcatch, ]
if (dim(tab)[1]>0) {
tabl(
"ERROR: Population estimates shouldn't be smaller than trap catches.",
TAB=tab, row.names=FALSE)
} else para("OKAY: Population estimates are greater than trap catches.")
sel1 <- year %in% (YEAR-(5:1))
nts <- tapply(!is.na(trapcatch[sel1]) & trapcatch[sel1]>0, lscode[sel1],
sum)
streamz1 <- names(nts)[nts>0]
sel2 <- year==YEAR
nts <- tapply(!is.na(trapcatch[sel2]) & trapcatch[sel2]>0, lscode[sel2],
sum)
streamz2 <- names(nts)[nts<1]
streamz <- as.numeric(intersect(streamz1, streamz2))
if (length(streamz)>0) {
tab <- df.nice[lscode %in% streamz &
((sel1 & !is.na(trapcatch) & trapcatch>0) | sel2), ]
tabl("Streams with no trap catch in ", YEAR,
", but at least one trap catch in previous five years.",
TAB=tab, row.names=FALSE)
} else para("OKAY: All streams with no trap catch in ", YEAR,
", also didn't have trap catch in the previous five years.")
tc <- ifelse(is.na(trapcatch), 0, trapcatch)
sel2 <- year==YEAR
tnow <- tapply(tc[sel2], lscode[sel2], mean, na.rm=TRUE)
sel1 <- year %in% (YEAR-(5:1)) & lscode %in% names(tnow)
tmin <- tapply(tc[sel1], lscode[sel1], min)
tmax <- tapply(tc[sel1], lscode[sel1], max)
streamz1 <- names(tnow)[tnow > 1.5*tmax]
streamz2 <- names(tnow)[tnow < tmin/2]
if (length(streamz1)>0) {
tab <- df.nice[lscode %in% streamz1 &
((sel1 & !is.na(trapcatch) & trapcatch>0) | sel2), ]
tabl("Streams with a ", YEAR,
" trap catch more than 50% greater than maximum from",
" the previous five years.",
TAB=tab, row.names=FALSE)
} else para("OKAY: No unusually large trap catches in ", YEAR, ".")
if (length(streamz2)>0) {
tab <- df.nice[lscode %in% streamz2 &
((sel1 & !is.na(trapcatch) & trapcatch>0) | sel2), ]
tabl("Streams with a ", YEAR,
" trap catch less than half of minimum from the previous five years.",
TAB=tab, row.names=FALSE)
} else para("OKAY: No unusually small trap catches in ", YEAR, ".")
})
streamDat$cle <- with(streamDat,
paste(casefold(substring(country, 1, 2), upper=TRUE),
Lakeabbs[lake], estr, sep=" - "))
sd <- with(streamDat, PEmr * CVmr/100)
streamDat$PEminusSD <- with(streamDat, PEmr - sd)
streamDat$PEminusSD[with(streamDat, !is.na(PEminusSD) & PEminusSD<0)] <- 0
streamDat$PEplusSD <- with(streamDat, PEmr + sd)
#### PLOTS (and tables) ####
# exploratory plots of the data
FIG.scatter <- function() {
with(streamDat, {
par(mfrow=c(3, 2), mar=c(4, 4, 2, 1), las=1, cex=1)
plotdf(streamDat[, c("lake", "country", "year")])
sel1 <- !is.na(trapcatch) & trapcatch>0 & !is.na(PEmr) & PEmr>0
plot(trapcatch[sel1], PEmr[sel1], log="xy",
xlab="Trap Catch", ylab="M-R Estimate")
sel1 <- !is.na(CVmr) & CVmr>0 & !is.na(PEmr) & PEmr>0
plot(PEmr[sel1], CVmr[sel1], log="xy", xlab="M-R Estimate", ylab="CV")
})
}
figu("Exploratory plots of ", YEAR, " adult sea lamprey data.",
FIG=FIG.scatter, newpage="port")
# plot trap efficiencies
FIG.trapEff <- function() {
sub <- streamDat[with(streamDat, !is.na(PEmr) & !is.na(trapcatch)), ]
with(sub, {
yrz <- min(year):max(year)
# use the latest observation for each stream
sub <- sub[with(sub, order(lscode, -year)), ]
suby <- sub[first(sub$lscode)==1, ]
suby <- suby[with(suby, order(lake, country, estr, strcode)), ]
sucle <- suby$cle
suls <- suby$lscode
tecol <- colr(suby$trapEff, "orange", "blue")
if(length(suls)>0) {
nrnc <- n2mfrow(length(suls))[2:1]
par(mfrow=nrnc, mar=c(0, 0, 0, 0), oma=c(4, 3, 1, 1), yaxs="i", cex=1.4)
for(i in seq(suls)) {
sel <- lscode==suls[i]
plot(1, 1, type="n", xlim=range(year), ylim=0:1, axes=FALSE,
xlab="", ylab="")
abline(h=seq(0.25, 0.75, 0.25), col="lightgray", lty=3)
abline(v=seq(1985, YEAR+2, 5), col="lightgray", lty=3)
lines(yrz, trapEff[sel][match(yrz, year[sel])],
type="o", pch=20, cex=1, col=tecol[i])
if (i <= nrnc[2]) {
axis(1, at=seq(1990, YEAR+2, 10), outer=TRUE, cex.axis=0.6, las=2,
tcl=-0.3)
}
if (i%%nrnc[2] == 1) {
axis(2, at=c(0.25, 0.75), labels=c("25", "75"), outer=TRUE,
cex.axis=0.6, las=1, tcl=-0.3)
}
mtext(paste(sucle[i], strname[sel][1], sep="\n"),
cex=1, adj=0.05, line=-2)
box(col="darkgray")
}
mtext("Year", side=1, outer=TRUE, cex=1.2, line=3)
mtext("Trap efficiency (%)", side=2, outer=TRUE, cex=1.2, line=2)
}
})
}
if(!is.null(streamDat$trapcatch)) {
# calculate trap efficiency
streamDat$trapEff <- with(streamDat, trapcatch/PEmr)
figu("Trap efficiency estimates (trap catch / mark-recapture PE), through ",
YEAR,
". Color is used to indicate streams with lowest (orange) and highest",
" (blue) trap efficiencies in the last year recorded for each stream.",
FIG=FIG.trapEff, newpage="port")
}
# plot stream PEs (long time series)
FIG.PElong <- function() {
sub <- streamDat[with(streamDat, !is.na(PEmr)), ]
with(sub, {
yrz <- min(year):max(year)
# use the latest observation for each stream
sub <- sub[with(sub, order(lscode, -year)), ]
suby <- sub[first(sub$lscode)==1, ]
suby <- suby[with(suby, order(lake, country, estr, strcode)), ]
sucle <- suby$cle
suls <- suby$lscode
pecol <- colr(suby$PEmr, "blue", "orange")
if(length(suls)>0) {
nrnc <- n2mfrow(length(suls))[2:1]
par(mfrow=nrnc, mar=c(0, 0, 0, 0), oma=c(4, 3, 1, 1), yaxs="i", cex=1.4)
for(i in seq(suls)) {
sel <- lscode==suls[i]
selsd <- sel & !is.na(PEplusSD)
plot(1, 1, type="n", xlim=range(year), ylim=1.1*range(0, PEmr)/1000,
axes=FALSE, xlab="", ylab="")
# ignore warnings from zero-length arrows
oldopt <- getOption("warn")
options(warn=-1)
arrows(year[selsd], PEminusSD[selsd]/1000,
year[selsd], PEplusSD[selsd]/1000,
col="gray", length=0.05, angle=90, code=3)
options(warn=oldopt)
lines(yrz, PEmr[sel][match(yrz, year[sel])]/1000,
type="o", pch=20, cex=1, col=pecol[i])
if (i <= nrnc[2]) {
axis(1, at=seq(1990, YEAR+2, 10), outer=TRUE, cex.axis=0.6, las=2,
tcl=-0.3)
}
if (i%%nrnc[2] == 1) {
axis(2, outer=TRUE, cex.axis=0.6, las=1, tcl=-0.3)
}
mtext(paste(sucle[i], strname[sel][1], sep="\n"),
cex=1, adj=0.05, line=-2)
box(col="darkgray")
}
mtext("Year", side=1, outer=TRUE, cex=1.2, line=3)
mtext("Adult PE (thousands)", side=2, outer=TRUE, cex=1.2, line=2)
}
})
}
figu("Adult mark-recapture estimates in streams,",
" plus or minus 1 standard deviation through ", YEAR,
". Color is used to indicate streams with lowest (blue) and highest",
" (orange) adult abundances in the last year recorded for each stream.",
FIG=FIG.PElong, newpage="port")
# line plots showing summary of stream PEs over past 5 years
FIG.strmest <- function() {
# selstreams <- with(streamDat, year > YEAR - 4.5 & !is.na(indexContrib))
selstreams <- with(streamDat, year > YEAR - 4.5 & !is.na(PEmr))
streams.pick <- sort(unique(streamDat$lscode[selstreams]))
sub <- streamDat[with(streamDat,
is.element(lscode, streams.pick) & year > YEAR - 4.5), ]
with(sub, {
suby <- sub[with(sub, year==YEAR), ]
suby <- suby[with(suby, order(lake, country, estr, strcode)), ]
sucle <- suby$cle
suls <- suby$lscode
if(length(suls)>0) {
nrnc <- n2mfrow(length(suls))
par(mfrow=nrnc, mar=c(0, 2, 2, 0), oma=c(3, 1, 1, 1), yaxs="i", cex=1.4)
for(i in seq(suls)) {
sel <- lscode==suls[i]
selsd <- sel & !is.na(PEplusSD)
mymax <- 1.1*max(c(PEplusSD[sel], PEmr[sel]), na.rm=TRUE)/1000
plot(year[sel], PEmr[sel]/1000, type="n", axes=FALSE,
xlim=range(year) + c(-1, 1)*0.3, ylim=c(0, mymax), xlab="", ylab="")
lines(year[sel], PEmr[sel]/1000, type="o", pch=16, cex=0.7)
# ignore warnings from zero-length arrows
oldopt <- getOption("warn")
options(warn=-1)
arrows(year[selsd], PEminusSD[selsd]/1000,
year[selsd], PEplusSD[selsd]/1000,
length=0.05, angle=90, code=3)
options(warn=oldopt)
if (i <= nrnc[2]) axis(1, outer=TRUE, cex.axis=0.6, tcl=-0.2)
axis(2, cex.axis=0.6, las=1, tcl=-0.2)
mtext(paste(sucle[i], strname[sel][1], sep="\n"), side=3,
cex=0.8, adj=0.1)
box(col="gray")
}
mtext("Year", side=1, outer=TRUE, cex=1.2, line=2)
mtext("Adult sea lampreys (thousands)", side=2, outer=TRUE, cex=1.2,
line=0)
}
})
}
figu("Recent adult sea lamprey mark-recapture estimates on index streams,",
" plus or minus 1 standard deviation, ", YEAR-4, "-", YEAR, ".",
FIG=FIG.strmest, newpage="port")
endrtf()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.