RunFn | R Documentation |
Run the Punt et al. (2008) ADMB-based ageing error model from within R.
RunFn(
Data,
SigOpt,
KnotAges,
BiasOpt,
NDataSets = 1,
MinAge,
MaxAge,
RefAge,
MinusAge,
PlusAge,
MaxSd,
MaxExpectedAge,
SaveFile,
EffSampleSize = 0,
Intern = TRUE,
AdmbFile = NULL,
JustWrite = FALSE,
CallType = "system",
ExtraArgs = " -est",
verbose = TRUE
)
Data |
This is the data set with the first column being an integer
providing the number of otoliths that are included in the row and the
subsequent columns are the reader or lab estimated ag,e where each
reader/lab has a unique reading error and bias. The modeling framework
allows for, at most, 15 readers, i.e., 16 columns. There should not be any
identical rows in the data frame because otoliths that have the exact same
read from every reader/lab should be combined into a single row with the
count as the first column. If you failed to combine identical rows prior
to running the model, you will be alerted with an error and the |
SigOpt |
This a vector with one entry for each reader (i.e.,
|
KnotAges |
Ages associated with each knot. This is a necessary input
for |
BiasOpt |
A vector with one entry for each reader/lab specifying the
type of bias specific to each reader. Positive values lead to estimated
parameters and negative values are used for shared parameters between
readers, just like with
An example entry for the situation where you have seven readers and you
assume that the first reader is unbiased, readers 2-7 have a curvilinear
bias, reader 3 shares parameters with reader 2, reader 5 shares parameters
with reader 4, and reader 7 shares parameters with reader 6 would look
like |
NDataSets |
This is generally |
MinAge |
An integer, specifying the minimum possible "true" age. |
MaxAge |
An integer, specifying the maximum possible "true" age. |
RefAge |
An arbitrarily chosen age from which "true" age-composition fixed-effects are calculated as an offset. This has no effect on the answer but could potentially effect estimation speed. |
MinusAge |
The minimum age for which an age-specific age-composition is
estimated. Ages below
,
where beta is an estimated log-linear trend in the "true"
proportion-at-age. If |
PlusAge |
Identical to |
MaxSd |
An upper bound on possible values for the standard deviation of reading error. |
MaxExpectedAge |
Set to MaxAge. |
SaveFile |
Directory where |
EffSampleSize |
Indicating whether effective sample size should be
calculated. Missing values in the data matrix will cause this to be
ineffective, in which case this should be set to |
Intern |
A logical input that controls the amount of output displayed,
where |
AdmbFile |
An optional character entry that specifies the directory
from which |
JustWrite |
A logical input that allows just the data files to be written without running ADMB executable. |
CallType |
Either |
ExtraArgs |
A string of characters providing extra arguments passed to
ADMB. The default is |
verbose |
A logical input that controls the amount of feedback users
receive from the program. The default is to provide the most output as
possible with |
The premise of Punt et al. (2008) is to calculate the likelihood of model parameters given an observed data set of otolith age reads from multiple age readers. For each reader/lab, two parameters are defined, one for standard deviation and one for bias. The model calculates the expected age of each read and the standard deviation of a normally distributed reading error given the true age of an otolith. These relationships can be linear or curvilinear.
The true age is obviously an unobserved process and can be considered a
random effect. Thus, the software computes the likelihood while summing
across all possible discrete values for the true age of each otolith. This
true age requires a hyperdistribution that represents the prior probability
that an otolith is any given age. The hyperdistribution is controlled by a
set of hyperparameters and the parameters that govern the standard deviation
and bias of each age reader/lab. Specifically, one hyperparameter is
estimated for every age between and including the MinusAge
and PlusAge
.
Ages outside of this range have a prior proportion at age defined as a
loglinear deviation from the proportion at age for the extreme ages, i.e.,
MinusAge
and PlusAge
. The slope of these loglinear deviations thus
constitutes an additional 1 or 2 fixed effect parameters. The true
proportion at age is then calculated from these fixed effects and loglinear
slope parameters by normalizing the resulting distribution such that it sums
to one.
James T. Thorson, Ian J. Stewart, Andre E. Punt, Ian G. Taylor
StepwiseFn()
will run multiple models.
PlotOutputFn()
will help summarize the output.
example(SimulatorFn)
## Not run:
utils::write.csv(AgeReads,
file = file.path(getwd(), "Simulated_data_example.csv"))
## End(Not run)
##### Format data
Nreaders <- ncol(AgeReads)
# Change NA to -999 (which the Punt software considers missing data)
AgeReads <- ifelse(is.na(AgeReads), -999, AgeReads)
# Potentially eliminate rows that are only read once
# These rows have no information about reading error, but are potentially
# informative about latent age-structure. It is unknown whether eliminating
# these rows degrades estimation of error and bias, and is currently
# recommended to speed up computation
if (FALSE) {
KeepRow <- ifelse(
rowSums(ifelse(AgeReads == -999, 0, 1), na.rm = TRUE) <= 1,
FALSE, TRUE
)
AgeReads <- AgeReads[KeepRow, ]
}
# AgeReads2 is the correctly formatted data object
AgeReads2 <- rMx(c(1, AgeReads[1, ]))
# Combine duplicate rows
for (RowI in 2:nrow(AgeReads)) {
DupRow <- NA
for (PreviousRowJ in 1:nrow(AgeReads2)) {
if (all(
AgeReads[RowI,1:Nreaders] == AgeReads2[PreviousRowJ,1:Nreaders+1]
)) {
DupRow <- PreviousRowJ
}
}
if (is.na(DupRow)) {# Add new row to AgeReads2
AgeReads2 <- rbind(AgeReads2, c(1, AgeReads[RowI, ]))
}
if(!is.na(DupRow)){# Increment number of samples for previous duplicate
AgeReads2[DupRow,1] <- AgeReads2[DupRow,1] + 1
}
}
######## Determine settings for ADMB
# Define minimum and maximum ages for integral across unobserved ages
MinAge <- 1
MaxAge <- ceiling(max(AgeReads2[,-1])/10)*10
BiasOpt <- c(0, -1, 0, -3)
SigOpt <- c(1, -1, 6, -3)
# Necessary for SigOpt option 5 or 6
KnotAges <- list(NA, NA, c(1, 10, 20, MaxAge), NA)
##### Run the model (MAY TAKE 5-10 MINUTES)
## Not run:
fileloc <- file.path(tempdir(), "age")
dir.create(fileloc, showWarnings = FALSE)
RunFn(Data = AgeReads2, SigOpt = SigOpt, KnotAges = KnotAges,
BiasOpt = BiasOpt,
NDataSets = 1, MinAge = MinAge, MaxAge = MaxAge, RefAge = 10,
MinusAge = 1, PlusAge = 30, SaveFile = fileloc,
AdmbFile = file.path(system.file("executables",
package = "nwfscAgeingError"), .Platform$file.sep),
EffSampleSize = 0, Intern = FALSE, JustWrite = FALSE, CallType = "shell"
)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.