From e46dc16b3694ebf5a114d11aba50d037c139a459 Mon Sep 17 00:00:00 2001
From: Ian Taylor <4992918+iantaylor-NOAA@users.noreply.github.com>
Date: Thu, 30 May 2024 09:34:26 -0700
Subject: [PATCH] run styler::style_pkg()

---
 R/AgeingError.R               |   4 +-
 R/CreateData.R                |   6 +-
 R/CreateSpecs.R               |   4 +-
 R/DoApplyAgeError.R           |   4 +-
 R/PlotOutputFn.R              | 495 +++++++++++++++++-----------------
 R/ProcessResults.R            |   4 +-
 R/RunFn.R                     | 349 ++++++++++++------------
 R/SimulatorFn.R               | 108 ++++----
 R/StepwiseFn.R                | 225 ++++++++--------
 R/ageing_comparison.R         |   9 +-
 R/estgrowth.vb.R              |  32 ++-
 R/plot_output.R               |   2 +-
 R/zzz.R                       |   1 -
 _pkgdown.yml                  |   2 +-
 tests/spelling.R              |   9 +-
 vignettes/getting_started.Rmd | 157 ++++++-----
 16 files changed, 735 insertions(+), 676 deletions(-)

diff --git a/R/AgeingError.R b/R/AgeingError.R
index dabdbdd..c1c9266 100644
--- a/R/AgeingError.R
+++ b/R/AgeingError.R
@@ -6,12 +6,12 @@
 #' @importFrom graphics hist lines mtext par
 #' @importFrom stats nlminb optim
 #' @importFrom utils read.table str write.csv
-#' 
+#'
 #' @references
 #' Punt, A.E., Smith, D.C., KrusicGolub, K., and Robertson, S. 2008.
 #' Quantifying age-reading error for use in fisheries stock assessments,
 #' with application to species in Australias southern and eastern scalefish and shark fishery.
 #' Canadian Journal of Fisheries and Aquatic Sciences 65: 1991-2005.
-#' 
+#'
 #' @keywords internal
 "_PACKAGE"
diff --git a/R/CreateData.R b/R/CreateData.R
index b2c2ed7..6c2a7cb 100644
--- a/R/CreateData.R
+++ b/R/CreateData.R
@@ -1,8 +1,8 @@
 #' Read the ageing error data
-#' 
+#'
 #' @param DataFile Filename for input data
 #' @param NDataSet Number of data sets within `DataFile`
-#' @param verbose Return messages to the console (in addition to any output to 
+#' @param verbose Return messages to the console (in addition to any output to
 #' `EchoFile`)
 #' @param EchoFile A file path to a file that will be created or appended to if
 #'   it already exists to store information about your data inputs. The default
@@ -339,4 +339,4 @@ CreateData <- function(DataFile = "data.dat",
     print(str(Outs))
   }
   return(Outs)
-}
\ No newline at end of file
+}
diff --git a/R/CreateSpecs.R b/R/CreateSpecs.R
index 6a7d5d1..be023dc 100644
--- a/R/CreateSpecs.R
+++ b/R/CreateSpecs.R
@@ -1,5 +1,5 @@
 #' Read the ageing error specifications
-#' 
+#'
 #' @param SpecsFile Filename for input specifications.
 #' @param DataSpecs The output from CreateData()
 #' @param verbose Return messages to the console (TRUE/FALSE)
@@ -254,4 +254,4 @@ CreateSpecs <- function(SpecsFile = "data.spc",
     print(str(Outs))
   }
   return(Outs)
-}
\ No newline at end of file
+}
diff --git a/R/DoApplyAgeError.R b/R/DoApplyAgeError.R
index fd51bfd..0f897eb 100644
--- a/R/DoApplyAgeError.R
+++ b/R/DoApplyAgeError.R
@@ -1,5 +1,5 @@
 #' Run the ageing error optimization routine
-#' 
+#'
 #' @param Species A string that will be used to create file names. Typically,
 #'   users will use the common name for the species of interest, especially if
 #'   you are saving files from multiple species in a single directory. Though,
@@ -296,4 +296,4 @@ DoApplyAgeError <- function(Species = "AgeingError",
   SaveAll$sdreport <- rep
   save(SaveAll, file = SaveFile)
   return(model)
-}
\ No newline at end of file
+}
diff --git a/R/PlotOutputFn.R b/R/PlotOutputFn.R
index f748727..b0b8a8e 100644
--- a/R/PlotOutputFn.R
+++ b/R/PlotOutputFn.R
@@ -10,7 +10,7 @@
 #'    shows the expected and standard deviation in age reads for that
 #'    reader/lab. This is displayed against a scatter plot of the read and
 #'    estimated ages for each otolith that was read by that reader/lab.
-#' 
+#'
 #' 1. Proportion-at-age histogram: The estimated proportion at age can be
 #'    plotted as a histogram and is displayed against the observed distribution
 #'    of read ages. This is useful to determine if hte estimated proportion at
@@ -52,283 +52,282 @@ PlotOutputFn <- function(Data,
                          ReaderNames = NULL,
                          ...) {
   PlotType <- match.arg(PlotType)
-    # Interpret inputs
-    Nreaders <- ncol(Data) - 1
-    Ages <- Nages <- MaxAge + 1
+  # Interpret inputs
+  Nreaders <- ncol(Data) - 1
+  Ages <- Nages <- MaxAge + 1
 
-    # Reader names
-    if (is.null(ReaderNames)) {
-      ReaderNames <- paste("Reader", 1:Nreaders)
-    }
+  # Reader names
+  if (is.null(ReaderNames)) {
+    ReaderNames <- paste("Reader", 1:Nreaders)
+  }
 
-    # Read REP file
-    Rep <- scan(file.path(SaveFile, "agemat.rep"),
-      comment.char = "%", what = "character", quiet = TRUE
-    )
+  # Read REP file
+  Rep <- scan(file.path(SaveFile, "agemat.rep"),
+    comment.char = "%", what = "character", quiet = TRUE
+  )
 
-    # Read Misclassification rates
-    Grep <- grep("reader#", Rep)
-    # dimensions are Reader, TrueAge, EstAge
-    MisclassArray <- array(NA,
-      dim = c(Nreaders, Ages, Ages),
-      dimnames = list(
-        paste("Reader", 1:Nreaders),
-        paste("TrueAge", 0:MaxAge),
-        paste("EstAge", 0:MaxAge)
-      )
+  # Read Misclassification rates
+  Grep <- grep("reader#", Rep)
+  # dimensions are Reader, TrueAge, EstAge
+  MisclassArray <- array(NA,
+    dim = c(Nreaders, Ages, Ages),
+    dimnames = list(
+      paste("Reader", 1:Nreaders),
+      paste("TrueAge", 0:MaxAge),
+      paste("EstAge", 0:MaxAge)
     )
-    for (i in 1:Nreaders) {
-      MisclassArray[i, , ] <- matrix(as.numeric(Rep[Grep[i] + 1 + 1:(Ages^2)]),
-        ncol = Ages, byrow = TRUE
-      )
-    }
-
-    # Input estimated age-structure
-    Grep <- grep("age-structure", Rep)
-    AgeStruct <- matrix(as.numeric(Rep[Grep[1] + 7 + 1:(2 * Ages)]),
-      ncol = 2, byrow = TRUE
+  )
+  for (i in 1:Nreaders) {
+    MisclassArray[i, , ] <- matrix(as.numeric(Rep[Grep[i] + 1 + 1:(Ages^2)]),
+      ncol = Ages, byrow = TRUE
     )
+  }
 
-    # Input reader error and bias
-    Grep <- grep("age", Rep)[3]
-    Temp <- Rep[Grep + 1:(5 * Nages * Nreaders)]
-    ErrorAndBiasArray <- array(as.numeric(Temp),
-      dim = c(5, Nages, Nreaders),
-      dimnames = list(
-        c(
-          "Reader", "True_Age",
-          "CV", "SD", "Expected_age"
-        ),
-        paste("Age", 0:MaxAge),
-        paste("Reader", 1:Nreaders)
-      )
-    )
+  # Input estimated age-structure
+  Grep <- grep("age-structure", Rep)
+  AgeStruct <- matrix(as.numeric(Rep[Grep[1] + 7 + 1:(2 * Ages)]),
+    ncol = 2, byrow = TRUE
+  )
 
-    # Estimate unobserved age for each otolith
-    # This is done by assigning each otolith to the age which has
-    # maximum posterior probability (i.e. the conditional mode,
-    # as is typically done for random effects)
-    AgeProbs <- array(NA,
-      dim = c(nrow(Data), Ages),
-      dimnames = list(
-        paste("Otolith", 1:nrow(Data)),
-        paste("TrueAge", 0:MaxAge)
-      )
+  # Input reader error and bias
+  Grep <- grep("age", Rep)[3]
+  Temp <- Rep[Grep + 1:(5 * Nages * Nreaders)]
+  ErrorAndBiasArray <- array(as.numeric(Temp),
+    dim = c(5, Nages, Nreaders),
+    dimnames = list(
+      c(
+        "Reader", "True_Age",
+        "CV", "SD", "Expected_age"
+      ),
+      paste("Age", 0:MaxAge),
+      paste("Reader", 1:Nreaders)
     )
-    # Check AEP adjustment (lets discuss)
-    OtI <- AgeI <- ReadI <- 1
-    for (OtI in 1:nrow(Data)) {
-      for (AgeI in 1:Ages) {
-        # AgeProbs[OtI, AgeI] <- AgeStruct[AgeI, 2]
-        AgeProbs[OtI, AgeI] <- 1
-        for (ReadI in 1:Nreaders) {
-          if (Data[OtI, ReadI + 1] != -999) {
-            AgeRead <- Data[OtI, ReadI + 1]
-            # AgeProbs[OtI, AgeI] <- AgeStruct[AgeI, 2] *
-            #  (MisclassArray[ReadI, AgeI, AgeRead+1])^Data[OtI, 1]
+  )
 
-            # ANDRE: Ageread+1 because the first column is age0
-            AgeProbs[OtI, AgeI] <- AgeProbs[OtI, AgeI] *
-              (MisclassArray[ReadI, AgeI, AgeRead + 1])^Data[OtI, 1]
-          } # end check for value other than -999
-        } # end loop over readers
-      } # end loop over ages
-    } # end loop over rows of data
-
-    # Remove MaxAge before calculating "TrueAge"
-    # because the MaxAge is a plus-group, and ends up with
-    # maximum probability for most ages in the upper tail
-    # ANDRE - Found another issue - should be age-1 because the first age is zero
-    TrueAge <- apply(AgeProbs,
-      MARGIN = 1,
-      FUN = function(Vec) {
-        order(Vec[-length(Vec)],
-          decreasing = TRUE
-        )[1]
-      }
-    ) - 1
+  # Estimate unobserved age for each otolith
+  # This is done by assigning each otolith to the age which has
+  # maximum posterior probability (i.e. the conditional mode,
+  # as is typically done for random effects)
+  AgeProbs <- array(NA,
+    dim = c(nrow(Data), Ages),
+    dimnames = list(
+      paste("Otolith", 1:nrow(Data)),
+      paste("TrueAge", 0:MaxAge)
+    )
+  )
+  # Check AEP adjustment (lets discuss)
+  OtI <- AgeI <- ReadI <- 1
+  for (OtI in 1:nrow(Data)) {
+    for (AgeI in 1:Ages) {
+      # AgeProbs[OtI, AgeI] <- AgeStruct[AgeI, 2]
+      AgeProbs[OtI, AgeI] <- 1
+      for (ReadI in 1:Nreaders) {
+        if (Data[OtI, ReadI + 1] != -999) {
+          AgeRead <- Data[OtI, ReadI + 1]
+          # AgeProbs[OtI, AgeI] <- AgeStruct[AgeI, 2] *
+          #  (MisclassArray[ReadI, AgeI, AgeRead+1])^Data[OtI, 1]
 
-    DataExpanded <- Data[rep(1:nrow(Data), Data[, 1]), -1]
-    DataExpanded[DataExpanded == -999] <- NA
+          # ANDRE: Ageread+1 because the first column is age0
+          AgeProbs[OtI, AgeI] <- AgeProbs[OtI, AgeI] *
+            (MisclassArray[ReadI, AgeI, AgeRead + 1])^Data[OtI, 1]
+        } # end check for value other than -999
+      } # end loop over readers
+    } # end loop over ages
+  } # end loop over rows of data
 
-    ####################################################################
-    # Plot comparison of data for each pair of readers
-    if (1 %in% subplot) {
+  # Remove MaxAge before calculating "TrueAge"
+  # because the MaxAge is a plus-group, and ends up with
+  # maximum probability for most ages in the upper tail
+  # ANDRE - Found another issue - should be age-1 because the first age is zero
+  TrueAge <- apply(AgeProbs,
+    MARGIN = 1,
+    FUN = function(Vec) {
+      order(Vec[-length(Vec)],
+        decreasing = TRUE
+      )[1]
+    }
+  ) - 1
 
-      # make plots of input data for each reader pair
-      for (ireader in 1:(Nreaders - 1)) {
-        for (jreader in (ireader + 1):Nreaders) {
-          if (PlotType == "PDF") {
-            pdfname <- paste0(
-              ReaderNames[ireader],
-              " vs ", ReaderNames[jreader], ".pdf"
-            )
-            grDevices::pdf(pdfname,
-              width = 6, height = 6
-            )
-          }
+  DataExpanded <- Data[rep(1:nrow(Data), Data[, 1]), -1]
+  DataExpanded[DataExpanded == -999] <- NA
 
-          out <- ageing_comparison(
-            xvec = DataExpanded[, ireader],
-            yvec = DataExpanded[, jreader],
-            xlab = ReaderNames[ireader],
-            ylab = ReaderNames[jreader],
-            maxage = max(DataExpanded, na.rm = TRUE),
-            hist = TRUE,
-            png = (PlotType == "PNG"),
-            SaveFile = SaveFile,
-            filename = paste0(
-              ReaderNames[ireader],
-              " vs ", ReaderNames[jreader], ".png"
-            ),
-            verbose = FALSE,
-            ...
+  ####################################################################
+  # Plot comparison of data for each pair of readers
+  if (1 %in% subplot) {
+    # make plots of input data for each reader pair
+    for (ireader in 1:(Nreaders - 1)) {
+      for (jreader in (ireader + 1):Nreaders) {
+        if (PlotType == "PDF") {
+          pdfname <- paste0(
+            ReaderNames[ireader],
+            " vs ", ReaderNames[jreader], ".pdf"
           )
-          if (PlotType == "PDF") {
-            grDevices::dev.off()
-            if (is.null(out)) {
-              unlink(pdfname)
-            }
+          grDevices::pdf(pdfname,
+            width = 6, height = 6
+          )
+        }
+
+        out <- ageing_comparison(
+          xvec = DataExpanded[, ireader],
+          yvec = DataExpanded[, jreader],
+          xlab = ReaderNames[ireader],
+          ylab = ReaderNames[jreader],
+          maxage = max(DataExpanded, na.rm = TRUE),
+          hist = TRUE,
+          png = (PlotType == "PNG"),
+          SaveFile = SaveFile,
+          filename = paste0(
+            ReaderNames[ireader],
+            " vs ", ReaderNames[jreader], ".png"
+          ),
+          verbose = FALSE,
+          ...
+        )
+        if (PlotType == "PDF") {
+          grDevices::dev.off()
+          if (is.null(out)) {
+            unlink(pdfname)
           }
         }
       }
-    } # end check for whether subplot 1 was requested
+    }
+  } # end check for whether subplot 1 was requested
 
-    ####################################################################
-    # Plot estimated age structure
+  ####################################################################
+  # Plot estimated age structure
 
-    if (2 %in% subplot) {
-      if (PlotType == "PDF") {
-        grDevices::pdf(
-          file.path(SaveFile, "Estimated vs Observed Age Structure.pdf"),
-          width = 6, height = 6
-        )
-      }
-      if (PlotType == "PNG") {
-        grDevices::png(
-          file.path(SaveFile, "Estimated vs Observed Age Structure.png"),
-          width = 6, height = 6, units = "in", res = 200
-        )
-      }
-      graphics::par(
-        mar = c(3, 3, 2, 0), mgp = c(1.5, 0.25, 0),
-        tck = -0.02, oma = c(0, 0, 0, 0) + 0.1
+  if (2 %in% subplot) {
+    if (PlotType == "PDF") {
+      grDevices::pdf(
+        file.path(SaveFile, "Estimated vs Observed Age Structure.pdf"),
+        width = 6, height = 6
       )
-      plot(
-        x = AgeStruct[, 1], y = AgeStruct[, 2], type = "s", lwd = 2,
-        xlab = "Age", ylab = "Prop", main = "Estimated=Black, Observed=Red"
-      )
-      graphics::hist(as.matrix(DataExpanded),
-        add = TRUE, freq = FALSE, breaks = seq(0, MaxAge, by = 1),
-        col = grDevices::rgb(red = 1, green = 0, blue = 0, alpha = 0.30)
+    }
+    if (PlotType == "PNG") {
+      grDevices::png(
+        file.path(SaveFile, "Estimated vs Observed Age Structure.png"),
+        width = 6, height = 6, units = "in", res = 200
       )
-      grDevices::dev.off()
-    } # end check for whether subplot was requested
+    }
+    graphics::par(
+      mar = c(3, 3, 2, 0), mgp = c(1.5, 0.25, 0),
+      tck = -0.02, oma = c(0, 0, 0, 0) + 0.1
+    )
+    plot(
+      x = AgeStruct[, 1], y = AgeStruct[, 2], type = "s", lwd = 2,
+      xlab = "Age", ylab = "Prop", main = "Estimated=Black, Observed=Red"
+    )
+    graphics::hist(as.matrix(DataExpanded),
+      add = TRUE, freq = FALSE, breaks = seq(0, MaxAge, by = 1),
+      col = grDevices::rgb(red = 1, green = 0, blue = 0, alpha = 0.30)
+    )
+    grDevices::dev.off()
+  } # end check for whether subplot was requested
 
-    ####################################################################
-    # Plot true age against different age reads
+  ####################################################################
+  # Plot true age against different age reads
 
-    if (3 %in% subplot) {
-      Ncol <- ceiling(sqrt(Nreaders))
-      Nrow <- ceiling(Nreaders / Ncol)
-      if (PlotType == "PDF") {
-        grDevices::pdf(
-          file.path(SaveFile, "True vs Reads (by reader).pdf"),
-          width = Ncol * 3, height = Nrow * 3
-        )
-      }
-      if (PlotType == "PNG") {
-        grDevices::png(
-          file.path(SaveFile, "True vs Reads (by reader).png"),
-          width = Ncol * 3, height = Nrow * 3, units = "in", res = 200
-        )
-      }
-      graphics::par(
-        mfrow = c(Nrow, Ncol), mar = c(3, 3, 2, 0), mgp = c(1.5, 0.25, 0),
-        tck = -0.02, oma = c(0, 0, 5, 0) + 0.1
+  if (3 %in% subplot) {
+    Ncol <- ceiling(sqrt(Nreaders))
+    Nrow <- ceiling(Nreaders / Ncol)
+    if (PlotType == "PDF") {
+      grDevices::pdf(
+        file.path(SaveFile, "True vs Reads (by reader).pdf"),
+        width = Ncol * 3, height = Nrow * 3
       )
-      for (ReadI in 1:Nreaders) {
-        Main <- ReaderNames[ReadI]
-
-        # Add 0.5 to match convention in Punt model that otoliths are read
-        # half way through year
-        Temp <- cbind(TrueAge, Data[, ReadI + 1] + 0.5)
-        # Exclude rows with no read for this reader
-        Temp <- Temp[which(Data[, ReadI + 1] != -999), ]
-        plot(
-          x = Temp[, 1], y = Temp[, 2],
-          ylim = c(0, MaxAge), xlim = c(0, MaxAge),
-          col = grDevices::rgb(red = 0, green = 0, blue = 0, alpha = 0.2),
-          xlab = "Mode predicted age | parameters",
-          ylab = "Read age", lwd = 2, main = Main, pch = 21, cex = 0.2
-        )
-        graphics::lines(
-          x = c(0, MaxAge),
-          y = c(0, MaxAge),
-          lwd = 1, lty = "dashed"
-        )
-        graphics::lines(
-          x = ErrorAndBiasArray["True_Age", , ReadI],
-          y = ErrorAndBiasArray["Expected_age", , ReadI],
-          type = "l", col = "red", lwd = 1
-        )
-        graphics::lines(
-          x = ErrorAndBiasArray["True_Age", , ReadI],
-          y = ErrorAndBiasArray["SD", , ReadI],
-          type = "l", col = "blue", lwd = 1
-        )
-        graphics::lines(
-          x = ErrorAndBiasArray["True_Age", , ReadI],
-          y = ErrorAndBiasArray["Expected_age", , ReadI] +
-            2 * ErrorAndBiasArray["SD", , ReadI],
-          type = "l", col = "red", lwd = 1, lty = "dashed"
-        )
-        graphics::lines(
-          x = ErrorAndBiasArray["True_Age", , ReadI],
-          y = ErrorAndBiasArray["Expected_age", , ReadI] -
-            2 * ErrorAndBiasArray["SD", , ReadI],
-          type = "l", col = "red", lwd = 1, lty = "dashed"
-        )
-      }
-      graphics::mtext(
-        side = 3, outer = TRUE,
-        text = paste0(
-          "Reads(dot), Sd(blue), expected_read(red solid line),\n",
-          " and 95% CI for expected_read(red dotted line)"
-        ),
-        line = 1
+    }
+    if (PlotType == "PNG") {
+      grDevices::png(
+        file.path(SaveFile, "True vs Reads (by reader).png"),
+        width = Ncol * 3, height = Nrow * 3, units = "in", res = 200
       )
-      grDevices::dev.off()
-    } # end check for whether subplot was requested
-
-    ## AIC
-    Nll <- as.numeric(scan(file.path(SaveFile, "agemat.par"),
-      comment.char = "%", what = "character",
-      quiet = TRUE
-    )[11])
-    Df <- as.numeric(scan(file.path(SaveFile, "agemat.par"),
-      comment.char = "%", what = "character",
-      quiet = TRUE
-    )[6])
-    n <- sum(ifelse(Data[, -1] == -999, 0, 1))
-    Aic <- 2 * Nll + 2 * Df
-    Aicc <- Aic + 2 * Df * (Df + 1) / (n - Df - 1)
-    Bic <- 2 * Nll + Df * log(n)
-
-    # Write definitions to file
+    }
+    graphics::par(
+      mfrow = c(Nrow, Ncol), mar = c(3, 3, 2, 0), mgp = c(1.5, 0.25, 0),
+      tck = -0.02, oma = c(0, 0, 5, 0) + 0.1
+    )
     for (ReadI in 1:Nreaders) {
       Main <- ReaderNames[ReadI]
-      utils::write.csv(ErrorAndBiasArray[, , ReadI],
-        file = file.path(SaveFile, paste0("SS_format_", Main, ".csv"))
+
+      # Add 0.5 to match convention in Punt model that otoliths are read
+      # half way through year
+      Temp <- cbind(TrueAge, Data[, ReadI + 1] + 0.5)
+      # Exclude rows with no read for this reader
+      Temp <- Temp[which(Data[, ReadI + 1] != -999), ]
+      plot(
+        x = Temp[, 1], y = Temp[, 2],
+        ylim = c(0, MaxAge), xlim = c(0, MaxAge),
+        col = grDevices::rgb(red = 0, green = 0, blue = 0, alpha = 0.2),
+        xlab = "Mode predicted age | parameters",
+        ylab = "Read age", lwd = 2, main = Main, pch = 21, cex = 0.2
+      )
+      graphics::lines(
+        x = c(0, MaxAge),
+        y = c(0, MaxAge),
+        lwd = 1, lty = "dashed"
+      )
+      graphics::lines(
+        x = ErrorAndBiasArray["True_Age", , ReadI],
+        y = ErrorAndBiasArray["Expected_age", , ReadI],
+        type = "l", col = "red", lwd = 1
+      )
+      graphics::lines(
+        x = ErrorAndBiasArray["True_Age", , ReadI],
+        y = ErrorAndBiasArray["SD", , ReadI],
+        type = "l", col = "blue", lwd = 1
+      )
+      graphics::lines(
+        x = ErrorAndBiasArray["True_Age", , ReadI],
+        y = ErrorAndBiasArray["Expected_age", , ReadI] +
+          2 * ErrorAndBiasArray["SD", , ReadI],
+        type = "l", col = "red", lwd = 1, lty = "dashed"
+      )
+      graphics::lines(
+        x = ErrorAndBiasArray["True_Age", , ReadI],
+        y = ErrorAndBiasArray["Expected_age", , ReadI] -
+          2 * ErrorAndBiasArray["SD", , ReadI],
+        type = "l", col = "red", lwd = 1, lty = "dashed"
       )
     }
+    graphics::mtext(
+      side = 3, outer = TRUE,
+      text = paste0(
+        "Reads(dot), Sd(blue), expected_read(red solid line),\n",
+        " and 95% CI for expected_read(red dotted line)"
+      ),
+      line = 1
+    )
+    grDevices::dev.off()
+  } # end check for whether subplot was requested
 
+  ## AIC
+  Nll <- as.numeric(scan(file.path(SaveFile, "agemat.par"),
+    comment.char = "%", what = "character",
+    quiet = TRUE
+  )[11])
+  Df <- as.numeric(scan(file.path(SaveFile, "agemat.par"),
+    comment.char = "%", what = "character",
+    quiet = TRUE
+  )[6])
+  n <- sum(ifelse(Data[, -1] == -999, 0, 1))
+  Aic <- 2 * Nll + 2 * Df
+  Aicc <- Aic + 2 * Df * (Df + 1) / (n - Df - 1)
+  Bic <- 2 * Nll + Df * log(n)
 
-    # Return stuff
-    ModelSelection <- list("AIC" = Aic, "AICc" = Aicc, "BIC" = Bic)
-    Output <- list(
-      "ModelSelection" = ModelSelection,
-      "ErrorAndBiasArray" = ErrorAndBiasArray
+  # Write definitions to file
+  for (ReadI in 1:Nreaders) {
+    Main <- ReaderNames[ReadI]
+    utils::write.csv(ErrorAndBiasArray[, , ReadI],
+      file = file.path(SaveFile, paste0("SS_format_", Main, ".csv"))
     )
-    return(Output)
   }
+
+
+  # Return stuff
+  ModelSelection <- list("AIC" = Aic, "AICc" = Aicc, "BIC" = Bic)
+  Output <- list(
+    "ModelSelection" = ModelSelection,
+    "ErrorAndBiasArray" = ErrorAndBiasArray
+  )
+  return(Output)
+}
diff --git a/R/ProcessResults.R b/R/ProcessResults.R
index 8fa221d..0481a64 100644
--- a/R/ProcessResults.R
+++ b/R/ProcessResults.R
@@ -1,5 +1,5 @@
 #' Process results of the ageing error estimation
-#' 
+#'
 #' @inheritParams DoApplyAgeError
 #' @param CalcEff Calculate effective sample sizes (TRUE/FALSE)
 #' @export
@@ -460,4 +460,4 @@ ProcessResults <- function(Species = "AgeingError",
   }
   # print(str(Output))
   return(Output)
-}
\ No newline at end of file
+}
diff --git a/R/RunFn.R b/R/RunFn.R
index 7264331..13b833d 100644
--- a/R/RunFn.R
+++ b/R/RunFn.R
@@ -158,14 +158,15 @@
 #' example(SimulatorFn)
 #' \dontrun{
 #' utils::write.csv(AgeReads,
-#'   file = file.path(getwd(), "Simulated_data_example.csv"))
+#'   file = file.path(getwd(), "Simulated_data_example.csv")
+#' )
 #' }
 #'
 #' ##### 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
@@ -178,51 +179,52 @@
 #'   )
 #'   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]
+#'       AgeReads[RowI, 1:Nreaders] == AgeReads2[PreviousRowJ, 1:Nreaders + 1]
 #'     )) {
 #'       DupRow <- PreviousRowJ
 #'     }
 #'   }
-#'   if (is.na(DupRow)) {# Add new row to AgeReads2
+#'   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
+#'   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
+#' 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)
 #' \dontrun{
 #' fileloc <- file.path(tempdir(), "age")
 #' dir.create(fileloc, showWarnings = FALSE)
-#' RunFn(Data = AgeReads2, SigOpt = SigOpt, KnotAges = KnotAges,
+#' 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),
+#'     package = "nwfscAgeingError"
+#'   ), .Platform$file.sep),
 #'   EffSampleSize = 0, Intern = FALSE, JustWrite = FALSE, CallType = "shell"
 #' )
 #' }
-
 RunFn <- function(Data,
                   SigOpt,
                   KnotAges,
@@ -243,184 +245,183 @@ RunFn <- function(Data,
                   CallType = "system",
                   ExtraArgs = " -est",
                   verbose = TRUE) {
+  # add slash to end of directories so that nobody has to waste as much time
+  # debugging as Ian just did
+  SaveFile <- paste0(SaveFile, "/")
 
-    # add slash to end of directories so that nobody has to waste as much time
-    # debugging as Ian just did
-    SaveFile <- paste0(SaveFile, "/")
-
-    # Copy ADMB file
-    if (!is.null(AdmbFile)) {
-      AdmbFile <- paste0(AdmbFile, "/")
-      # Check for missing file before trying to copy
-      if (is.na(file.info(file.path(AdmbFile, "agemat.exe"))$size)) {
-        warning(
-          "executable 'agemat.exe' not found in\n",
-          AdmbFile
-        )
-      }
-      if (verbose) {
-        cat(
-          "copying 'agemat.exe' from\n", AdmbFile,
-          "\nto\n", SaveFile, "\n"
-        )
-      }
-      file.copy(
-        from = file.path(AdmbFile, "agemat.exe"),
-        to = file.path(SaveFile, "agemat.exe"), overwrite = TRUE
+  # Copy ADMB file
+  if (!is.null(AdmbFile)) {
+    AdmbFile <- paste0(AdmbFile, "/")
+    # Check for missing file before trying to copy
+    if (is.na(file.info(file.path(AdmbFile, "agemat.exe"))$size)) {
+      warning(
+        "executable 'agemat.exe' not found in\n",
+        AdmbFile
       )
     }
-    # Check for missing file
-    if (is.na(file.info(file.path(SaveFile, "agemat.exe"))$size)) {
-      stop(
-        "executable 'agemat.exe' not found in\n",
-        SaveFile
+    if (verbose) {
+      cat(
+        "copying 'agemat.exe' from\n", AdmbFile,
+        "\nto\n", SaveFile, "\n"
       )
     }
+    file.copy(
+      from = file.path(AdmbFile, "agemat.exe"),
+      to = file.path(SaveFile, "agemat.exe"), overwrite = TRUE
+    )
+  }
+  # Check for missing file
+  if (is.na(file.info(file.path(SaveFile, "agemat.exe"))$size)) {
+    stop(
+      "executable 'agemat.exe' not found in\n",
+      SaveFile
+    )
+  }
 
-    # Check for errors
-    Nreaders <- ncol(Data) - 1
-    for (ReaderI in 1:Nreaders) {
-      if ((SigOpt[ReaderI] == 5 | SigOpt[ReaderI] == 6) &
-        is.na(KnotAges[[ReaderI]][1])) {
-        stop("Must specify KnotAges for any reader with SigOpt 5 or 6")
-      }
+  # Check for errors
+  Nreaders <- ncol(Data) - 1
+  for (ReaderI in 1:Nreaders) {
+    if ((SigOpt[ReaderI] == 5 | SigOpt[ReaderI] == 6) &
+      is.na(KnotAges[[ReaderI]][1])) {
+      stop("Must specify KnotAges for any reader with SigOpt 5 or 6")
     }
+  }
 
-    # Check for specification errors
-    for (ReaderI in 1:Nreaders) {
-      if ((SigOpt[ReaderI] < 0 & SigOpt[ReaderI] <= (-ReaderI)) |
-        (BiasOpt[ReaderI] < 0 & BiasOpt[ReaderI] <= (-ReaderI))) {
-        stop("Mirrored readers must mirror a lower numbered reader")
-      }
+  # Check for specification errors
+  for (ReaderI in 1:Nreaders) {
+    if ((SigOpt[ReaderI] < 0 & SigOpt[ReaderI] <= (-ReaderI)) |
+      (BiasOpt[ReaderI] < 0 & BiasOpt[ReaderI] <= (-ReaderI))) {
+      stop("Mirrored readers must mirror a lower numbered reader")
     }
+  }
+
+  # Write DAT file
+  datfile <- file.path(SaveFile, "agemat.dat")
+  # simple functions to avoid repeated code
+  writeLine <- function(x, file = datfile, append = TRUE) {
+    write(x = x, file = file, append = append)
+  }
+  writeTable <- function(x, file = datfile, append = TRUE) {
+    utils::write.table(
+      x = x, file = file, append = append,
+      row.names = FALSE, col.names = FALSE
+    )
+  }
 
-    # Write DAT file
-    datfile <- file.path(SaveFile, "agemat.dat")
-    # simple functions to avoid repeated code
-    writeLine <- function(x, file = datfile, append = TRUE) {
-      write(x = x, file = file, append = append)
+  # first line creates new file (not appended to existing file)
+  writeLine(c("# Maximum number of readers", Nreaders), append = FALSE)
+  # subsequent lines use default append = TRUE
+  writeLine(c("# Number of data sets", NDataSets))
+  writeLine(c("# Number of points per data set", nrow(Data)))
+  writeLine(c("# Readers per data set", ncol(Data) - 1))
+  writeLine("# Which readers per data set")
+  writeTable(rMx(1:(ncol(Data) - 1)))
+  writeLine(c("# Minimum age", MinAge))
+  writeLine(c("# Maximum age", MaxAge))
+  writeLine(c("# Reference age", RefAge))
+  writeLine(c("# Minus groups", MinusAge))
+  writeLine(c("# Plus groups", PlusAge))
+  # write table of options for bias
+  writeLine("# Option for bias")
+  writeTable(rMx(BiasOpt))
+  # write table of options for SD
+  writeLine("# Option for standard deviation")
+  writeTable(rMx(SigOpt))
+  writeLine(c("# Option for effective sample size", EffSampleSize))
+  writeLine(c("# Use Par File (1=Yes)", 0))
+  # Write knots related to splines
+  writeLine("\n# Number and location of knots for any splines")
+  for (ReaderI in 1:Nreaders) {
+    if (SigOpt[ReaderI] == 5 | SigOpt[ReaderI] == 6) {
+      writeTable(rMx(c(length(KnotAges[[ReaderI]]), KnotAges[[ReaderI]])))
     }
-    writeTable <- function(x, file = datfile, append = TRUE) {
-      utils::write.table(
-        x = x, file = file, append = append,
-        row.names = FALSE, col.names = FALSE
-      )
+  }
+  # Write initial values
+  # Bias
+  writeLine("\n# Min, Max, Init, Phase for Bias")
+  for (BiasI in 1:Nreaders) {
+    # No bias
+    if (BiasOpt[BiasI] <= 0) {}
+    # Linear bias
+    if (BiasOpt[BiasI] == 1) {
+      writeTable(rMx(c(0.001, 3, 1, 2)))
     }
-
-    # first line creates new file (not appended to existing file)
-    writeLine(c("# Maximum number of readers", Nreaders), append = FALSE)
-    # subsequent lines use default append = TRUE
-    writeLine(c("# Number of data sets", NDataSets))
-    writeLine(c("# Number of points per data set", nrow(Data)))
-    writeLine(c("# Readers per data set", ncol(Data) - 1))
-    writeLine("# Which readers per data set")
-    writeTable(rMx(1:(ncol(Data) - 1)))
-    writeLine(c("# Minimum age", MinAge))
-    writeLine(c("# Maximum age", MaxAge))
-    writeLine(c("# Reference age", RefAge))
-    writeLine(c("# Minus groups", MinusAge))
-    writeLine(c("# Plus groups", PlusAge))
-    # write table of options for bias
-    writeLine("# Option for bias")
-    writeTable(rMx(BiasOpt))
-    # write table of options for SD
-    writeLine("# Option for standard deviation")
-    writeTable(rMx(SigOpt))
-    writeLine(c("# Option for effective sample size", EffSampleSize))
-    writeLine(c("# Use Par File (1=Yes)", 0))
-    # Write knots related to splines
-    writeLine("\n# Number and location of knots for any splines")
-    for (ReaderI in 1:Nreaders) {
-      if (SigOpt[ReaderI] == 5 | SigOpt[ReaderI] == 6) {
-        writeTable(rMx(c(length(KnotAges[[ReaderI]]), KnotAges[[ReaderI]])))
-      }
+    # Curvilinear bias = 0.5+Par1 + (Par3-Par1)/(1.0-mfexp(-Par2*(float(MaxAge)-1)))*(1.0-mfexp(-Par2*(float(Age1)-1)))
+    # Starting value must be non-zero
+    if (BiasOpt[BiasI] == 2) {
+      writeTable(rMx(c(0.001, 10, 1, 2)))
+      writeTable(rMx(c(-10, 1, 0.01, 2)))
+      writeTable(rMx(c(0.001, MaxAge * 2, MaxAge, 2)))
     }
-    # Write initial values
-    # Bias
-    writeLine("\n# Min, Max, Init, Phase for Bias")
-    for (BiasI in 1:Nreaders) {
-      # No bias
-      if (BiasOpt[BiasI] <= 0) {}
-      # Linear bias
-      if (BiasOpt[BiasI] == 1) {
-        writeTable(rMx(c(0.001, 3, 1, 2)))
-      }
-      # Curvilinear bias = 0.5+Par1 + (Par3-Par1)/(1.0-mfexp(-Par2*(float(MaxAge)-1)))*(1.0-mfexp(-Par2*(float(Age1)-1)))
-      # Starting value must be non-zero
-      if (BiasOpt[BiasI] == 2) {
-        writeTable(rMx(c(0.001, 10, 1, 2)))
-        writeTable(rMx(c(-10, 1, 0.01, 2)))
-        writeTable(rMx(c(0.001, MaxAge * 2, MaxAge, 2)))
-      }
+  }
+  # Sigma
+  writeLine("\n# Min, Max, Init, Phase for Sigma")
+  for (SigI in 1:Nreaders) {
+    # No error
+    if (SigOpt[SigI] <= 0) {}
+    # Linear CV
+    if (SigOpt[SigI] == 1) {
+      writeTable(rMx(c(0.001, 3, 0.1, 2)))
     }
-    # Sigma
-    writeLine("\n# Min, Max, Init, Phase for Sigma")
-    for (SigI in 1:Nreaders) {
-      # No error
-      if (SigOpt[SigI] <= 0) {}
-      # Linear CV
-      if (SigOpt[SigI] == 1) {
-        writeTable(rMx(c(0.001, 3, 0.1, 2)))
-      }
-      # Curvilinear SD = Par1 + (Par3-Par1)/(1.0-exp(-Par2*(100-1)))*(1.0-exp(-Par2*(1:100)))
-      # Starting value must be non-zero
-      if (SigOpt[SigI] == 2) {
-        writeTable(rMx(c(0.001, 100, 1, 2)))
-        writeTable(rMx(c(-10, 1, 0.01, 2)))
-        writeTable(rMx(c(0.001, 100, 10, 2)))
-      }
-      # Curvilinear CV
-      # Curvilinear CV = Par1 + (Par3-Par1)/(1.0-exp(-Par2*(100-1)))*(1.0-exp(-Par2*(1:100)))
-      # Starting value must be non-zero
-      if (SigOpt[SigI] == 3) {
-        writeTable(rMx(c(0.001, 3, 0.1, 2)))
-        writeTable(rMx(c(-10, 1, 0.01, 2)))
-        writeTable(rMx(c(0.001, 3, 0.1, 2)))
-      }
-      # Spline with estimated derivative at beginning and end
-      # (Params 1-N: knot parameters; N+1 and N+2: derivative at beginning and end)
-      if (SigOpt[SigI] == 5) {
-        for (ParI in 1:(2 + length(KnotAges[[SigI]]))) {
-          writeTable(rMx(c(-10.0, 10.0, 1.0, 1)))
-        }
-      }
-      # Spline with derivative at beginning and end fixed at zero
-      # (Params 1-N: knot parameters)
-      if (SigOpt[SigI] == 6) {
-        for (ParI in 1:length(KnotAges[[SigI]])) {
-          writeTable(rMx(c(-10.0, 10.0, 1.0, 1)))
-        }
-      }
+    # Curvilinear SD = Par1 + (Par3-Par1)/(1.0-exp(-Par2*(100-1)))*(1.0-exp(-Par2*(1:100)))
+    # Starting value must be non-zero
+    if (SigOpt[SigI] == 2) {
+      writeTable(rMx(c(0.001, 100, 1, 2)))
+      writeTable(rMx(c(-10, 1, 0.01, 2)))
+      writeTable(rMx(c(0.001, 100, 10, 2)))
+    }
+    # Curvilinear CV
+    # Curvilinear CV = Par1 + (Par3-Par1)/(1.0-exp(-Par2*(100-1)))*(1.0-exp(-Par2*(1:100)))
+    # Starting value must be non-zero
+    if (SigOpt[SigI] == 3) {
+      writeTable(rMx(c(0.001, 3, 0.1, 2)))
+      writeTable(rMx(c(-10, 1, 0.01, 2)))
+      writeTable(rMx(c(0.001, 3, 0.1, 2)))
     }
-    # Probs (i.e. age-composition probability relative to reference age)
-    writeLine("\n# Min, Max, Phase for Probs")
-    writeTable(rMx(c(-20, 20, 2)))
-    # Slopes
-    writeLine("\n# Min, Max, Init, Phase for slopes")
-    for (DataSetI in 1:NDataSets) {
-      if (MaxAge > PlusAge) {
-        writeTable(rMx(c(-10, 0, 0, 1)))
+    # Spline with estimated derivative at beginning and end
+    # (Params 1-N: knot parameters; N+1 and N+2: derivative at beginning and end)
+    if (SigOpt[SigI] == 5) {
+      for (ParI in 1:(2 + length(KnotAges[[SigI]]))) {
+        writeTable(rMx(c(-10.0, 10.0, 1.0, 1)))
       }
-      if (MinAge < MinusAge) {
-        writeTable(rMx(c(-10, 0, 0, 1)))
+    }
+    # Spline with derivative at beginning and end fixed at zero
+    # (Params 1-N: knot parameters)
+    if (SigOpt[SigI] == 6) {
+      for (ParI in 1:length(KnotAges[[SigI]])) {
+        writeTable(rMx(c(-10.0, 10.0, 1.0, 1)))
       }
     }
-    # Write dataset
-    writeLine("\n# Data set")
-    writeTable(Data)
-    writeLine(c("# Test number", 123456))
+  }
+  # Probs (i.e. age-composition probability relative to reference age)
+  writeLine("\n# Min, Max, Phase for Probs")
+  writeTable(rMx(c(-20, 20, 2)))
+  # Slopes
+  writeLine("\n# Min, Max, Init, Phase for slopes")
+  for (DataSetI in 1:NDataSets) {
+    if (MaxAge > PlusAge) {
+      writeTable(rMx(c(-10, 0, 0, 1)))
+    }
+    if (MinAge < MinusAge) {
+      writeTable(rMx(c(-10, 0, 0, 1)))
+    }
+  }
+  # Write dataset
+  writeLine("\n# Data set")
+  writeTable(Data)
+  writeLine(c("# Test number", 123456))
 
-    # Run ADMB file
-    if (JustWrite == FALSE) {
-      setwd(SaveFile)
-      if (CallType == "shell") {
-        # This may need to have the location pasted onto it depending upon file structure
-        Output <- shell(paste0("agemat.exe", ExtraArgs), intern = Intern)
-      }
-      if (CallType == "system") {
-        Output <- system(paste0("agemat.exe", ExtraArgs), intern = Intern)
-      }
-      # Admb = scan(paste(SaveFile,"agemat.par",sep=""),comment.char="#",quiet=TRUE)
+  # Run ADMB file
+  if (JustWrite == FALSE) {
+    setwd(SaveFile)
+    if (CallType == "shell") {
+      # This may need to have the location pasted onto it depending upon file structure
+      Output <- shell(paste0("agemat.exe", ExtraArgs), intern = Intern)
+    }
+    if (CallType == "system") {
+      Output <- system(paste0("agemat.exe", ExtraArgs), intern = Intern)
     }
-    # return(Output)
+    # Admb = scan(paste(SaveFile,"agemat.par",sep=""),comment.char="#",quiet=TRUE)
   }
+  # return(Output)
+}
diff --git a/R/SimulatorFn.R b/R/SimulatorFn.R
index 17d2faf..6f5ad56 100644
--- a/R/SimulatorFn.R
+++ b/R/SimulatorFn.R
@@ -35,70 +35,78 @@
 #' @examples
 #' # Parameters for generating data
 #' # This represents 2 unique readers
-#'   # Row 1 -- Otoliths read only once by reader 
-#'   # Row 2 -- Otoliths read twice by reader 1
-#'   # Row 2 -- Otoliths read only once by reader 2
-#'   # Row 4 -- Otoliths read twice by reader 2
-#'   # Row 5 -- Otoliths read once by reader 1 and once by reader 2
-#' ReadsMat <- structure(matrix(nrow = 5, ncol = 5,
-#'   c(rep(25, 5),
+#' # Row 1 -- Otoliths read only once by reader
+#' # Row 2 -- Otoliths read twice by reader 1
+#' # Row 2 -- Otoliths read only once by reader 2
+#' # Row 4 -- Otoliths read twice by reader 2
+#' # Row 5 -- Otoliths read once by reader 1 and once by reader 2
+#' ReadsMat <- structure(matrix(
+#'   nrow = 5, ncol = 5,
+#'   c(
+#'     rep(25, 5),
 #'     1, 1, 0, 0, 1,
 #'     0, 1, 0, 0, 0,
 #'     0, 0, 1, 1, 1,
-#'     0, 0, 0, 1, 0)
-#'   ), dimnames = list(
-#'     c("Reader1_Only", "Reader1_DoubleReads",
-#'       "Reader2_Only", "Reader2_DoubleReads",
-#'       "Reader1_&_Reader2"
-#'     ),
-#'    c("NumberOfReads",
+#'     0, 0, 0, 1, 0
+#'   )
+#' ), dimnames = list(
+#'   c(
+#'     "Reader1_Only", "Reader1_DoubleReads",
+#'     "Reader2_Only", "Reader2_DoubleReads",
+#'     "Reader1_&_Reader2"
+#'   ),
+#'   c(
+#'     "NumberOfReads",
 #'     "Reader1", "Reader1_DoubleReads",
 #'     "Reader2", "Reader2_DoubleReads"
-#'     ))
 #'   )
-#' 
+#' ))
+#'
 #' # Generate data
 #' set.seed(2)
-#' AgeReads <- SimulatorFn(Nreaders = 4, M = 0.2,
+#' AgeReads <- SimulatorFn(
+#'   Nreaders = 4, M = 0.2,
 #'   SelexForm = "Logistic",
 #'   SelexParams = c(5, 0.2), BiasParams = c(1, 1, 1.1, 1.1),
 #'   ErrorParams = c(0.2, 0.2, 0.2, 0.2), ReadsMat = ReadsMat,
-#'   RecCv = 0.6, RecAr1 = 0.8, Amax = 100)
-
-SimulatorFn <-function(Nreaders, M,
-  SelexForm, ErrorParams, BiasParams, SelexParams,
-  ReadsMat, RecCv = 0.6, RecAr1 = 0.8, Amax = 100) {
-    RecDev <- stats::rnorm(Amax, mean = 0, sd = RecCv)
-    for (i in 2:length(RecDev)) {
-      RecDev[i] <- RecDev[i] * sqrt(1 - RecAr1) + RecDev[i - 1] * sqrt(RecAr1)
-    }
-    AgeStruct <- 1 * exp(-M * 1:Amax) * exp(RecDev - RecCv^2 / 2)
-    if (SelexForm == "Logistic") {
-      SelexAtAge <- 1 / (1 + exp((SelexParams[1] - 1:Amax) * SelexParams[2]))
-    } else {
-      stop("Selex not implemented")
-    }
-    Ages <- sample(x = 1:Amax, size = sum(ReadsMat[, 1]), prob = AgeStruct * SelexAtAge, replace = TRUE)
+#'   RecCv = 0.6, RecAr1 = 0.8, Amax = 100
+#' )
+SimulatorFn <- function(
+    Nreaders, M,
+    SelexForm, ErrorParams, BiasParams, SelexParams,
+    ReadsMat, RecCv = 0.6, RecAr1 = 0.8, Amax = 100) {
+  RecDev <- stats::rnorm(Amax, mean = 0, sd = RecCv)
+  for (i in 2:length(RecDev)) {
+    RecDev[i] <- RecDev[i] * sqrt(1 - RecAr1) + RecDev[i - 1] * sqrt(RecAr1)
+  }
+  AgeStruct <- 1 * exp(-M * 1:Amax) * exp(RecDev - RecCv^2 / 2)
+  if (SelexForm == "Logistic") {
+    SelexAtAge <- 1 / (1 + exp((SelexParams[1] - 1:Amax) * SelexParams[2]))
+  } else {
+    stop("Selex not implemented")
+  }
+  Ages <- sample(x = 1:Amax, size = sum(ReadsMat[, 1]), prob = AgeStruct * SelexAtAge, replace = TRUE)
 
-    AgeReads <- vector()
+  AgeReads <- vector()
 
-    IndexI <- 0
-    for (RowI in 1:nrow(ReadsMat)) {
-      for (OtolithI in 1:ReadsMat[RowI, 1]) {
-        IndexI <- IndexI + 1
-        Row <- vector()
-        for (ReaderI in 1:Nreaders) {
-          if (ReadsMat[RowI, ReaderI + 1] == 0) {
-            Row <- c(Row, NA)
-          } else {
-            Row <- c(Row,
-              round(Ages[IndexI] * BiasParams[ReaderI] +
-                stats::rnorm(1, mean = 0, sd = Ages[IndexI] * ErrorParams[ReaderI])
-              ))
-          }
+  IndexI <- 0
+  for (RowI in 1:nrow(ReadsMat)) {
+    for (OtolithI in 1:ReadsMat[RowI, 1]) {
+      IndexI <- IndexI + 1
+      Row <- vector()
+      for (ReaderI in 1:Nreaders) {
+        if (ReadsMat[RowI, ReaderI + 1] == 0) {
+          Row <- c(Row, NA)
+        } else {
+          Row <- c(
+            Row,
+            round(Ages[IndexI] * BiasParams[ReaderI] +
+              stats::rnorm(1, mean = 0, sd = Ages[IndexI] * ErrorParams[ReaderI]))
+          )
         }
-        AgeReads <- rbind(AgeReads, Row)
       }
+      AgeReads <- rbind(AgeReads, Row)
     }
-    return(AgeReads)
   }
+  return(AgeReads)
+}
diff --git a/R/StepwiseFn.R b/R/StepwiseFn.R
index 0438069..677b151 100644
--- a/R/StepwiseFn.R
+++ b/R/StepwiseFn.R
@@ -57,19 +57,22 @@
 #' ##### Run the model (MAY TAKE 5-10 MINUTES)
 #' fileloc <- file.path(tempdir(), "age")
 #' dir.create(fileloc, showWarnings = FALSE)
-#' RunFn(Data = AgeReads2, SigOpt = SigOpt, KnotAges = KnotAges,
+#' 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),
+#'     package = "nwfscAgeingError"
+#'   ), .Platform$file.sep),
 #'   EffSampleSize = 0, Intern = FALSE, JustWrite = FALSE, CallType = "shell"
 #' )
 #' # Plot output
-#' PlotOutputFn(Data = AgeReads2, MaxAge = MaxAge,
+#' PlotOutputFn(
+#'   Data = AgeReads2, MaxAge = MaxAge,
 #'   SaveFile = fileloc, PlotType = "PDF"
 #' )
-#'}
+#' }
 #'
 #' ##### Stepwise selection
 #'
@@ -91,9 +94,13 @@
 #' # across for that reader
 #' SearchMat <- array(NA,
 #'   dim = c(Nreaders * 2 + 2, 7),
-#'   dimnames = list(c(paste("Error_Reader", 1:Nreaders),
-#'     paste("Bias_Reader", 1:Nreaders), "MinusAge", "PlusAge"),
-#'     paste("Option", 1:7))
+#'   dimnames = list(
+#'     c(
+#'       paste("Error_Reader", 1:Nreaders),
+#'       paste("Bias_Reader", 1:Nreaders), "MinusAge", "PlusAge"
+#'     ),
+#'     paste("Option", 1:7)
+#'   )
 #' )
 #' # Readers 1 and 3 search across options 1-3 for ERROR
 #' SearchMat[c(1, 3), 1:3] <- rep(1, 2) %o% c(1, 2, 3)
@@ -119,7 +126,7 @@
 #'   StartMinusAge + 4,
 #'   StartMinusAge + 10
 #' )
-#' SearchMat[9, 1:7] <- ifelse(SearchMat[9,1:7] < MinAge,
+#' SearchMat[9, 1:7] <- ifelse(SearchMat[9, 1:7] < MinAge,
 #'   NA, SearchMat[9, 1:7]
 #' )
 #' # PlusAge searches with a search kernal of -10,-4,-1,+0,+1,+4,+10
@@ -132,8 +139,9 @@
 #'   StartPlusAge + 4,
 #'   StartPlusAge + 10
 #' )
-#' SearchMat[10,1:7] <- ifelse(SearchMat[10, 1:7] > MaxAge,
-#'   NA, SearchMat[10, 1:7])
+#' SearchMat[10, 1:7] <- ifelse(SearchMat[10, 1:7] > MaxAge,
+#'   NA, SearchMat[10, 1:7]
+#' )
 #'
 #' # Run model selection
 #' # This outputs a series of files
@@ -150,7 +158,8 @@
 #' # WARNING: One run of this stepwise model building example can take
 #' # 8+ hours, and should be run overnight
 #' \dontrun{
-#' StepwiseFn(SearchMat = SearchMat, Data = AgeReads2,
+#' StepwiseFn(
+#'   SearchMat = SearchMat, Data = AgeReads2,
 #'   NDataSets = 1, MinAge = MinAge, MaxAge = MaxAge,
 #'   RefAge = 10, MaxSd = 40, MaxExpectedAge = MaxAge + 10,
 #'   SaveFile = fileloc, InformationCriterion = c("AIC", "AICc", "BIC")[3]
@@ -171,113 +180,113 @@ StepwiseFn <- function(SearchMat,
                        Intern = TRUE,
                        InformationCriterion = c("AIC", "AICc", "BIC"),
                        SelectAges = TRUE) {
-
   InformationCriterion <- match.arg(InformationCriterion)
-    # Define variables
-    Nages <- MaxAge + 1
-    Nreaders <- ncol(Data) - 1
-    ParamVecOpt <- SearchMat[, 1]
-    Stop <- FALSE
-    IcRecord <- NULL
-    StateRecord <- NULL
-    OuterIndex <- 0
-
-    # Continue searching until Stop==TRUE
-    while (Stop == FALSE) {
-
-      # Increment and intialize variables
-      OuterIndex <- OuterIndex + 1
-      Index <- 0
-      IcVec <- NULL
-      ParamMat <- NULL
-      Rep <- NULL
-      ParamVecOptPreviouslyEstimates <- FALSE
+  # Define variables
+  Nages <- MaxAge + 1
+  Nreaders <- ncol(Data) - 1
+  ParamVecOpt <- SearchMat[, 1]
+  Stop <- FALSE
+  IcRecord <- NULL
+  StateRecord <- NULL
+  OuterIndex <- 0
 
-      # Loop across all combinations of parameters obtained from one change in the current parameters
-      for (VarI in 1:nrow(SearchMat)) {
-        for (ValueI in 1:length(stats::na.omit(SearchMat[VarI, ]))) {
+  # Continue searching until Stop==TRUE
+  while (Stop == FALSE) {
+    # Increment and intialize variables
+    OuterIndex <- OuterIndex + 1
+    Index <- 0
+    IcVec <- NULL
+    ParamMat <- NULL
+    Rep <- NULL
+    ParamVecOptPreviouslyEstimates <- FALSE
 
-          # Update the current vector of parameters
-          ParamVecCurrent <- ParamVecOpt
-          ParamVecCurrent[VarI] <- stats::na.omit(SearchMat[VarI, ])[ValueI]
+    # Loop across all combinations of parameters obtained from one change in the current parameters
+    for (VarI in 1:nrow(SearchMat)) {
+      for (ValueI in 1:length(stats::na.omit(SearchMat[VarI, ]))) {
+        # Update the current vector of parameters
+        ParamVecCurrent <- ParamVecOpt
+        ParamVecCurrent[VarI] <- stats::na.omit(SearchMat[VarI, ])[ValueI]
 
-          # Decide if this ParamVecCurrent should be run
-          if (all(ParamVecCurrent == ParamVecOpt) & ParamVecOptPreviouslyEstimates == FALSE | !all(ParamVecCurrent == ParamVecOpt)) {
+        # Decide if this ParamVecCurrent should be run
+        if (all(ParamVecCurrent == ParamVecOpt) & ParamVecOptPreviouslyEstimates == FALSE | !all(ParamVecCurrent == ParamVecOpt)) {
+          # If running the current optimum, change so that it won't run again this loop
+          if (all(ParamVecCurrent == ParamVecOpt)) ParamVecOptPreviouslyEstimates <- TRUE
 
-            # If running the current optimum, change so that it won't run again this loop
-            if (all(ParamVecCurrent == ParamVecOpt)) ParamVecOptPreviouslyEstimates <- TRUE
+          # Make a new file for ADMB
+          if (!grepl(paste0(.Platform$file.sep, "$"), SaveFile)) {
+            SaveFile <- paste0(SaveFile, .Platform$file.sep)
+          }
+          RunFile <- paste(SaveFile, "Run\\", sep = "")
+          dir.create(RunFile, showWarnings = FALSE)
+          file.copy(from = paste(SaveFile, "agemat.exe", sep = ""), to = paste(RunFile, "agemat.exe", sep = ""))
 
-            # Make a new file for ADMB
-            if (!grepl(paste0(.Platform$file.sep, "$"), SaveFile)) {
-              SaveFile <- paste0(SaveFile, .Platform$file.sep)
-            }
-            RunFile <- paste(SaveFile, "Run\\", sep = "")
-            dir.create(RunFile, showWarnings = FALSE)
-            file.copy(from = paste(SaveFile, "agemat.exe", sep = ""), to = paste(RunFile, "agemat.exe", sep = ""))
+          # Increment Index
+          Index <- Index + 1
+          print(paste("Loop=", OuterIndex, " Run=", Index, " StartTime=", date(), sep = ""))
 
-            # Increment Index
-            Index <- Index + 1
-            print(paste("Loop=", OuterIndex, " Run=", Index, " StartTime=", date(), sep = ""))
+          # Configure and run model
+          SigOpt <- ParamVecCurrent[1:Nreaders]
+          BiasOpt <- ParamVecCurrent[Nreaders + 1:Nreaders]
+          MinusAge <- ParamVecCurrent[2 * Nreaders + 1]
+          PlusAge <- ParamVecCurrent[2 * Nreaders + 2]
+          RunFn(SigOpt = SigOpt, KnotAges = KnotAges, BiasOpt = BiasOpt, Data = Data, NDataSets = NDataSets, MinAge = MinAge, MaxAge = MaxAge, RefAge = RefAge, MinusAge = MinusAge, PlusAge = PlusAge, MaxSd = MaxSd, MaxExpectedAge = MaxExpectedAge, SaveFile = RunFile)
 
-            # Configure and run model
-            SigOpt <- ParamVecCurrent[1:Nreaders]
-            BiasOpt <- ParamVecCurrent[Nreaders + 1:Nreaders]
-            MinusAge <- ParamVecCurrent[2 * Nreaders + 1]
-            PlusAge <- ParamVecCurrent[2 * Nreaders + 2]
-            RunFn(SigOpt = SigOpt, KnotAges = KnotAges, BiasOpt = BiasOpt, Data = Data, NDataSets = NDataSets, MinAge = MinAge, MaxAge = MaxAge, RefAge = RefAge, MinusAge = MinusAge, PlusAge = PlusAge, MaxSd = MaxSd, MaxExpectedAge = MaxExpectedAge, SaveFile = RunFile)
+          # Compute information criteria
+          Df <- as.numeric(scan(paste(RunFile, "agemat.par", sep = ""), comment.char = "%", what = "character", quiet = TRUE)[6])
+          Nll <- as.numeric(scan(paste(RunFile, "agemat.par", sep = ""), comment.char = "%", what = "character", quiet = TRUE)[11])
+          n <- sum(ifelse(Data[, -1] == -999, 0, 1))
+          Aic <- 2 * Nll + 2 * Df
+          Aicc <- Aic + 2 * Df * (Df + 1) / (n - Df - 1)
+          Bic <- 2 * Nll + Df * log(n)
+          if (InformationCriterion == "AIC") IcVec <- c(IcVec, Aic)
+          if (InformationCriterion == "AICc") IcVec <- c(IcVec, Aicc)
+          if (InformationCriterion == "BIC") IcVec <- c(IcVec, Bic)
+          ParamMat <- rbind(ParamMat, ParamVecCurrent)
+          utils::write.table(cbind(IcVec, ParamMat),
+            paste0(SaveFile, "Stepwise - Model loop ", OuterIndex, ".txt"),
+            sep = "\t", row.names = FALSE
+          )
 
-            # Compute information criteria
-            Df <- as.numeric(scan(paste(RunFile, "agemat.par", sep = ""), comment.char = "%", what = "character", quiet = TRUE)[6])
-            Nll <- as.numeric(scan(paste(RunFile, "agemat.par", sep = ""), comment.char = "%", what = "character", quiet = TRUE)[11])
-            n <- sum(ifelse(Data[, -1] == -999, 0, 1))
-            Aic <- 2 * Nll + 2 * Df
-            Aicc <- Aic + 2 * Df * (Df + 1) / (n - Df - 1)
-            Bic <- 2 * Nll + Df * log(n)
-            if (InformationCriterion == "AIC") IcVec <- c(IcVec, Aic)
-            if (InformationCriterion == "AICc") IcVec <- c(IcVec, Aicc)
-            if (InformationCriterion == "BIC") IcVec <- c(IcVec, Bic)
-            ParamMat <- rbind(ParamMat, ParamVecCurrent)
-            utils::write.table(cbind(IcVec, ParamMat),
-              paste0(SaveFile, "Stepwise - Model loop ", OuterIndex, ".txt"),
-              sep = "\t", row.names = FALSE
-            )
-
-            # Input misclassification matrices
-            Rep[[Index]] <- readLines(paste(RunFile, "agemat.rep", sep = ""))
-          } # End if-statement for only running ParamVecOpt once per loop
-        }
-      } # End loop accross VarI and ValueI
+          # Input misclassification matrices
+          Rep[[Index]] <- readLines(paste(RunFile, "agemat.rep", sep = ""))
+        } # End if-statement for only running ParamVecOpt once per loop
+      }
+    } # End loop accross VarI and ValueI
 
-      IcRecord <- rbind(IcRecord, IcVec)
-      StateRecord <- rbind(StateRecord, ParamVecOpt)
-      utils::capture.output(list(IcRecord = IcRecord,
-        StateRecord = StateRecord),
-        file = paste0(SaveFile, "Stepwise - Record.txt"))
+    IcRecord <- rbind(IcRecord, IcVec)
+    StateRecord <- rbind(StateRecord, ParamVecOpt)
+    utils::capture.output(
+      list(
+        IcRecord = IcRecord,
+        StateRecord = StateRecord
+      ),
+      file = paste0(SaveFile, "Stepwise - Record.txt")
+    )
 
-      # Change current vector of optimum parameters
-      Max <- which.max(IcVec)
-      if (all(ParamMat[Max, ] == ParamVecOpt)) Stop <- TRUE
-      ParamVecOpt <- ParamMat[Max, ]
+    # Change current vector of optimum parameters
+    Max <- which.max(IcVec)
+    if (all(ParamMat[Max, ] == ParamVecOpt)) Stop <- TRUE
+    ParamVecOpt <- ParamMat[Max, ]
 
-      # Change boundaries for MinusAge parameter
-      CurrentMinusAge <- ParamVecOpt[length(ParamVecOpt) - 1]
-      if (SelectAges == TRUE) {
-        SearchMat[length(ParamVecOpt) - 1, 1:7] <- c(CurrentMinusAge, CurrentMinusAge - 10, CurrentMinusAge - 4, CurrentMinusAge - 1, CurrentMinusAge + 1, CurrentMinusAge + 4, CurrentMinusAge + 10)
-        SearchMat[length(ParamVecOpt) - 1, 1:7] <- ifelse(SearchMat[length(ParamVecOpt) - 1, 1:7] < MinAge, NA, SearchMat[length(ParamVecOpt) - 1, 1:7])
-      } else {
-        SearchMat[length(ParamVecOpt) - 1, 1:7] <- c(CurrentMinusAge, rep(NA, 6))
-      }
+    # Change boundaries for MinusAge parameter
+    CurrentMinusAge <- ParamVecOpt[length(ParamVecOpt) - 1]
+    if (SelectAges == TRUE) {
+      SearchMat[length(ParamVecOpt) - 1, 1:7] <- c(CurrentMinusAge, CurrentMinusAge - 10, CurrentMinusAge - 4, CurrentMinusAge - 1, CurrentMinusAge + 1, CurrentMinusAge + 4, CurrentMinusAge + 10)
+      SearchMat[length(ParamVecOpt) - 1, 1:7] <- ifelse(SearchMat[length(ParamVecOpt) - 1, 1:7] < MinAge, NA, SearchMat[length(ParamVecOpt) - 1, 1:7])
+    } else {
+      SearchMat[length(ParamVecOpt) - 1, 1:7] <- c(CurrentMinusAge, rep(NA, 6))
+    }
 
-      # Change boundaries for PlusAge parameter
-      CurrentPlusAge <- ParamVecOpt[length(ParamVecOpt)]
-      if (SelectAges == TRUE) {
-        SearchMat[length(ParamVecOpt), 1:7] <- c(CurrentPlusAge, CurrentPlusAge - 10, CurrentPlusAge - 4, CurrentPlusAge - 1, CurrentPlusAge + 1, CurrentPlusAge + 4, CurrentPlusAge + 10)
-        SearchMat[length(ParamVecOpt), 1:7] <- ifelse(SearchMat[length(ParamVecOpt), 1:7] > MaxAge, NA, SearchMat[length(ParamVecOpt), 1:7])
-      } else {
-        SearchMat[length(ParamVecOpt), 1:7] <- c(CurrentPlusAge, rep(NA, 6))
-      }
-      # Save image for each while loop
-      writeLines(Rep[[Max]], con = paste(SaveFile, "agemat.rep", sep = ""))
-      PlotOutputFn(Data = Data, MaxAge = MaxAge, SaveFile = SaveFile, PlotType = "JPG")
-    } # End while statement
-  } # End StepwiseFn
+    # Change boundaries for PlusAge parameter
+    CurrentPlusAge <- ParamVecOpt[length(ParamVecOpt)]
+    if (SelectAges == TRUE) {
+      SearchMat[length(ParamVecOpt), 1:7] <- c(CurrentPlusAge, CurrentPlusAge - 10, CurrentPlusAge - 4, CurrentPlusAge - 1, CurrentPlusAge + 1, CurrentPlusAge + 4, CurrentPlusAge + 10)
+      SearchMat[length(ParamVecOpt), 1:7] <- ifelse(SearchMat[length(ParamVecOpt), 1:7] > MaxAge, NA, SearchMat[length(ParamVecOpt), 1:7])
+    } else {
+      SearchMat[length(ParamVecOpt), 1:7] <- c(CurrentPlusAge, rep(NA, 6))
+    }
+    # Save image for each while loop
+    writeLines(Rep[[Max]], con = paste(SaveFile, "agemat.rep", sep = ""))
+    PlotOutputFn(Data = Data, MaxAge = MaxAge, SaveFile = SaveFile, PlotType = "JPG")
+  } # End while statement
+} # End StepwiseFn
diff --git a/R/ageing_comparison.R b/R/ageing_comparison.R
index 22b20d2..eb25df5 100644
--- a/R/ageing_comparison.R
+++ b/R/ageing_comparison.R
@@ -107,19 +107,22 @@ ageing_comparison <- function(xvec, yvec, scale.pts = 2,
   graphics::grid()
   # add lines at 0 if axes have padding around 0
   if (axs == "r") {
-    graphics::rect(xleft = 0, ybottom = 0,
+    graphics::rect(
+      xleft = 0, ybottom = 0,
       xright = graphics::par()$usr[2], ytop = graphics::par()$usr[4]
     )
   } else {
     graphics::box()
   }
   # add points
-  graphics::points(df1[, 1:2], col = col.pts,
+  graphics::points(df1[, 1:2],
+    col = col.pts,
     pch = 16, cex = scale.pts * sqrt(df1[, 3])
   )
   # add counts as text
   if (counts) {
-    graphics::text(df1[, 1:2], col = "white",
+    graphics::text(df1[, 1:2],
+      col = "white",
       lab = paste(df1[, 3]), cex = scale.pts / 3
     )
   }
diff --git a/R/estgrowth.vb.R b/R/estgrowth.vb.R
index 5adb8d6..99b474e 100644
--- a/R/estgrowth.vb.R
+++ b/R/estgrowth.vb.R
@@ -32,28 +32,34 @@
 #'
 #' @examples
 #' \dontrun{
-#' bio_dat <- data.frame(Age = rep(0:30, each = 20),
-#'   Length_cm = rnorm(n = 31 * 20, mean = 50, sd = 5))
+#' bio_dat <- data.frame(
+#'   Age = rep(0:30, each = 20),
+#'   Length_cm = rnorm(n = 31 * 20, mean = 50, sd = 5)
+#' )
 #' pars_in <- lapply(FUN = log, X = list(
 #'   "K" = 0.13,
 #'   "Linf" = 55,
 #'   "L0" = 5,
 #'   "CV0" = 0.1,
-#'   "CV1" = 0.1))
-#' solve <- optim(fn = estgrowth.vb, par = unlist(pars_in), hessian = FALSE,
+#'   "CV1" = 0.1
+#' ))
+#' solve <- optim(
+#'   fn = estgrowth.vb, par = unlist(pars_in), hessian = FALSE,
 #'   Ages = bio_dat[, "Age"],
-#'   Lengths = bio_dat[, "Length_cm"])
-#' predictions <- estgrowth.vb(Par = solve$par, ReturnType = "Pred",
+#'   Lengths = bio_dat[, "Length_cm"]
+#' )
+#' predictions <- estgrowth.vb(
+#'   Par = solve$par, ReturnType = "Pred",
 #'   sdFactor = 1,
 #'   Ages = bio_dat[, "Age"],
-#'   Lengths = bio_dat[, "Length_cm"])
+#'   Lengths = bio_dat[, "Length_cm"]
+#' )
 #' plot(bio_dat$Age, predictions[, "Lhat_pred"],
-#'   xlab = "Age (years)", ylab = "Predicted length (cm)")
+#'   xlab = "Age (years)", ylab = "Predicted length (cm)"
+#' )
 #' exp(solve$par)
 #' }
-
 estgrowth.vb <- function(Par, Ages, Lengths, ReturnType = c("NLL", "Pred"), sdFactor = 1) {
-
   #### Initialization
   if (is.null(names(Par))) {
     names(Par) <- c("K", "Linf", "L0", "CV0", "CV1")
@@ -96,11 +102,11 @@ estgrowth.vb <- function(Par, Ages, Lengths, ReturnType = c("NLL", "Pred"), sdFa
 
   #### Return
   if (ReturnType == "NLL") {
-    Return <- -1 * sum(dnorm(Lengths, mean=Lhat, sd=SD, log=TRUE), na.rm=TRUE)
+    Return <- -1 * sum(dnorm(Lengths, mean = Lhat, sd = SD, log = TRUE), na.rm = TRUE)
   }
   if (ReturnType == "Pred") {
-    Lhat_low   <- lhat.here(Par, Ages, variability = (-sdFactor * SD))
-    Lhat_high  <- lhat.here(Par, Ages, variability = (sdFactor * SD))
+    Lhat_low <- lhat.here(Par, Ages, variability = (-sdFactor * SD))
+    Lhat_high <- lhat.here(Par, Ages, variability = (sdFactor * SD))
     Return <- cbind(Lhat_low = Lhat_low, Lhat_pred = Lhat, Lhat_high = Lhat_high)
   }
   return(Return)
diff --git a/R/plot_output.R b/R/plot_output.R
index cb324d7..e63003f 100644
--- a/R/plot_output.R
+++ b/R/plot_output.R
@@ -332,4 +332,4 @@ plot_output <- function(Data,
     ErrorAndBiasArray = ErrorAndBiasArray
   )
   return(Output)
-}
\ No newline at end of file
+}
diff --git a/R/zzz.R b/R/zzz.R
index a65726d..16e5dd3 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,4 +1,3 @@
-
 .onUnload <- function(libpath) {
   library.dynam.unload("AgeingError", libpath)
 }
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 63f4f8e..a3a35e1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -1,4 +1,4 @@
-url: http://pfmc-assessments.github.io/nwfscAgeingError/
+url: http://pfmc-assessments.github.io/AgeingError/
 template:
   bootstrap: 5
 
diff --git a/tests/spelling.R b/tests/spelling.R
index 6713838..13f77d9 100644
--- a/tests/spelling.R
+++ b/tests/spelling.R
@@ -1,3 +1,6 @@
-if(requireNamespace('spelling', quietly = TRUE))
-  spelling::spell_check_test(vignettes = TRUE, error = FALSE,
-                             skip_on_cran = TRUE)
+if (requireNamespace("spelling", quietly = TRUE)) {
+  spelling::spell_check_test(
+    vignettes = TRUE, error = FALSE,
+    skip_on_cran = TRUE
+  )
+}
diff --git a/vignettes/getting_started.Rmd b/vignettes/getting_started.Rmd
index 4932794..bcf3592 100644
--- a/vignettes/getting_started.Rmd
+++ b/vignettes/getting_started.Rmd
@@ -363,8 +363,10 @@ the data file with a single data set so we set `NDataSet=1`.
 
 ```{r}
 data_dir <- system.file("extdata", package = "AgeingError")
-BG2022_dat <- AgeingError:::CreateData(file.path(data_dir, "BG2022.dat"), NDataSet=1, 
-                                       verbose=TRUE, EchoFile="BG2022echo.out")
+BG2022_dat <- AgeingError:::CreateData(file.path(data_dir, "BG2022.dat"),
+  NDataSet = 1,
+  verbose = TRUE, EchoFile = "BG2022echo.out"
+)
 ```
 When we set `verbose=TRUE` a summary of the loaded data is printed to the console.
 Check that the this matches with the input file, some things that are worth 
@@ -388,8 +390,10 @@ function is not exported either so again we need to call the package name and
 `:::` to access it. 
 
 ```{r}
-BG2022_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "BG2022.spc"), DataSpecs=BG2022_dat, 
-                                       verbose=TRUE)
+BG2022_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "BG2022.spc"),
+  DataSpecs = BG2022_dat,
+  verbose = TRUE
+)
 ```
 
 The function creates a nested list with one list element for each reader in the 
@@ -410,14 +414,15 @@ a large amount of output is printed to the screen.
 
 <!-- The output of this function isn't useful, set results="hide" -->
 ```{r, results="hide"}
-BG2022_mod <- AgeingError::DoApplyAgeError(Species = "BG2022",
-                                           DataSpecs = BG2022_dat,
-                                           ModelSpecsInp = BG2022_spc,
-                                           AprobWght = 1e-06,
-                                           SlopeWght = 0.01,
-                                           SaveDir = "Results",
-                                           verbose = FALSE)
-
+BG2022_mod <- AgeingError::DoApplyAgeError(
+  Species = "BG2022",
+  DataSpecs = BG2022_dat,
+  ModelSpecsInp = BG2022_spc,
+  AprobWght = 1e-06,
+  SlopeWght = 0.01,
+  SaveDir = "Results",
+  verbose = FALSE
+)
 ```
 
 The model is stored as a list in the object `BG2022_mod` and this
@@ -486,8 +491,9 @@ Load the data file.
 
 ```{r}
 BG2022final_dat <- AgeingError:::CreateData(file.path(data_dir, "BG2022_trim_8_1.dat"),
-                                            NDataSet=1, verbose=TRUE,
-                                            EchoFile="BG2022finalecho.out")
+  NDataSet = 1, verbose = TRUE,
+  EchoFile = "BG2022finalecho.out"
+)
 ```
 
 Note the number of lines of data has been reduced from 264 to 252. 
@@ -499,21 +505,24 @@ age-reading error (`SigmaOpt=7`) and fixed the first parameter at 0.2.
 
 ```{r}
 BG2022final_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "BG2022_trim_8_1.spc"),
-                                             DataSpecs=BG2022final_dat, 
-                                             verbose=TRUE)
+  DataSpecs = BG2022final_dat,
+  verbose = TRUE
+)
 ```
 
 We run the final model.
 
 <!-- The output of this function isn't useful, set results="hide" -->
 ```{r, results="hide"}
-BG2022final_mod <- AgeingError::DoApplyAgeError(Species = "BG2022_trim_8_1",
-                                           DataSpecs = BG2022final_dat,
-                                           ModelSpecsInp = BG2022final_spc,
-                                           AprobWght = 1e-06,
-                                           SlopeWght = 0.01,
-                                           SaveDir = "Results",
-                                           verbose = FALSE)
+BG2022final_mod <- AgeingError::DoApplyAgeError(
+  Species = "BG2022_trim_8_1",
+  DataSpecs = BG2022final_dat,
+  ModelSpecsInp = BG2022final_spc,
+  AprobWght = 1e-06,
+  SlopeWght = 0.01,
+  SaveDir = "Results",
+  verbose = FALSE
+)
 ```
 
 As with the earlier model we save the results.
@@ -534,8 +543,10 @@ readings from three separate readers who have each re-read their own reads (i.e.
 there are no inter-reader comparisons in the data). 
 
 ```{r}
-REB2022_dat <- AgeingError:::CreateData(file.path(data_dir, "REB2022.dat"), NDataSet=1, 
-                                       verbose=TRUE, EchoFile="REB2022echo.out")
+REB2022_dat <- AgeingError:::CreateData(file.path(data_dir, "REB2022.dat"),
+  NDataSet = 1,
+  verbose = TRUE, EchoFile = "REB2022echo.out"
+)
 ```
 
 We check the summary above against the information in the data file. There don't
@@ -546,7 +557,8 @@ Next we load the model specifications using the `CreateSpecs` function.
 
 ```{r}
 REB2022_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "REB2022.spc"),
-                                         DataSpecs=REB2022_dat, verbose=TRUE)
+  DataSpecs = REB2022_dat, verbose = TRUE
+)
 ```
 
 For Bight Redfish we fit a one parameter model assuming the expected age is a
@@ -559,14 +571,15 @@ We estimate ageing error for Bight Redfish.
 
 <!-- The output of this function isn't useful, set results="hide" -->
 ```{r, results="hide"}
-REB2022_mod <- AgeingError::DoApplyAgeError(Species = "REB2022",
-                                            DataSpecs = REB2022_dat,
-                                            ModelSpecsInp = REB2022_spc,
-                                            AprobWght = 1e-06,
-                                            SlopeWght = 0.01,
-                                            SaveDir = "Results",
-                                            verbose = FALSE)
-
+REB2022_mod <- AgeingError::DoApplyAgeError(
+  Species = "REB2022",
+  DataSpecs = REB2022_dat,
+  ModelSpecsInp = REB2022_spc,
+  AprobWght = 1e-06,
+  SlopeWght = 0.01,
+  SaveDir = "Results",
+  verbose = FALSE
+)
 ```
 
 We then save the Bight Redfish model output.
@@ -585,17 +598,23 @@ model (note we suppress the creation of output in this document).
 
 <!-- set results="hide" -->
 ```{r, results="hide"}
-REB2022_dat2 <- AgeingError:::CreateData(file.path(data_dir, "REB2022_trim.dat"), NDataSet=1, 
-                  verbose=FALSE, EchoFile="REB2022_trimecho.out")
+REB2022_dat2 <- AgeingError:::CreateData(file.path(data_dir, "REB2022_trim.dat"),
+  NDataSet = 1,
+  verbose = FALSE, EchoFile = "REB2022_trimecho.out"
+)
 REB2022_spc2 <- AgeingError:::CreateSpecs(file.path(data_dir, "REB2022_trim.spc"),
-                  DataSpecs=REB2022_dat2, verbose=FALSE)
-REB2022_mod2 <- AgeingError::DoApplyAgeError(Species = "REB2022_final",
-                  DataSpecs = REB2022_dat2, ModelSpecsInp = REB2022_spc2,
-                  AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results",
-                  verbose = FALSE)
-REB2022_out2 <- AgeingError:::ProcessResults(Species = "REB2022_final", 
-                  SaveDir = "Results", CalcEff = TRUE, verbose = FALSE)
-
+  DataSpecs = REB2022_dat2, verbose = FALSE
+)
+REB2022_mod2 <- AgeingError::DoApplyAgeError(
+  Species = "REB2022_final",
+  DataSpecs = REB2022_dat2, ModelSpecsInp = REB2022_spc2,
+  AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results",
+  verbose = FALSE
+)
+REB2022_out2 <- AgeingError:::ProcessResults(
+  Species = "REB2022_final",
+  SaveDir = "Results", CalcEff = TRUE, verbose = FALSE
+)
 ```
 
 
@@ -612,8 +631,10 @@ Load the School Whiting data file, make sure to set the number of data sets to
 three (`NDataSet=3`).
 
 ```{r}
-WHS2_dat <- AgeingError:::CreateData(file.path(data_dir, "WHS2.dat"), NDataSet=3, 
-                  verbose=TRUE, EchoFile="WHS2echo.out")
+WHS2_dat <- AgeingError:::CreateData(file.path(data_dir, "WHS2.dat"),
+  NDataSet = 3,
+  verbose = TRUE, EchoFile = "WHS2echo.out"
+)
 ```
 
 Make sure to check the summary above against the values in the data file.
@@ -622,7 +643,8 @@ Load the specifications for School Whiting.
 
 ```{r}
 WHS2_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "WHS2.spc"),
-                  DataSpecs=WHS2_dat, verbose=TRUE)
+  DataSpecs = WHS2_dat, verbose = TRUE
+)
 ```
 
 For School Whiting we fit a spline function (`SigmaOption=5`) with five knots
@@ -632,12 +654,16 @@ Fit the ageing error model and save the results
 
 <!-- Set results="hide" -->
 ```{r, results="hide"}
-WHS2_mod <- AgeingError::DoApplyAgeError(Species = "WHS2",
-                  DataSpecs = WHS2_dat, ModelSpecsInp = WHS2_spc,
-                  AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results",
-                  verbose = FALSE)
-WHS2_out <- AgeingError:::ProcessResults(Species = "WHS2", 
-                  SaveDir = "Results", CalcEff = TRUE, verbose = FALSE)
+WHS2_mod <- AgeingError::DoApplyAgeError(
+  Species = "WHS2",
+  DataSpecs = WHS2_dat, ModelSpecsInp = WHS2_spc,
+  AprobWght = 1e-06, SlopeWght = 0.01, SaveDir = "Results",
+  verbose = FALSE
+)
+WHS2_out <- AgeingError:::ProcessResults(
+  Species = "WHS2",
+  SaveDir = "Results", CalcEff = TRUE, verbose = FALSE
+)
 ```
 
 
@@ -653,8 +679,10 @@ appropriateness of this model should be discussed with the ageing technicians.
 Sablefish (_Anoplopoma fimbria_) provides an example of using the bias option.
 
 ```{r}
-Sable_dat <- AgeingError:::CreateData(file.path(data_dir, "Sable.dat"), NDataSet=1, 
-                                      verbose=TRUE, EchoFile="SableEcho.out")
+Sable_dat <- AgeingError:::CreateData(file.path(data_dir, "Sable.dat"),
+  NDataSet = 1,
+  verbose = TRUE, EchoFile = "SableEcho.out"
+)
 ```
 
 We check the summary above against the information in the data file. This data
@@ -664,7 +692,8 @@ We load the model specifications using the `CreateSpecs` function.
 
 ```{r}
 Sable_spc <- AgeingError:::CreateSpecs(file.path(data_dir, "Sable.spc"),
-                                         DataSpecs=Sable_dat, verbose=TRUE)
+  DataSpecs = Sable_dat, verbose = TRUE
+)
 ```
 
 The Sablefish example fits three separate series (1-2, 3-4 & 5-6). The model 
@@ -689,13 +718,15 @@ run it because this particular model is very slow.
 <!-- We don't run this code -->
 ```{r, eval=FALSE}
 ## run the Sablefish model
-Sable_mod <- AgeingError::DoApplyAgeError(Species = "Sable",
-                                          DataSpecs = Sable_dat,
-                                          ModelSpecsInp = Sable_spc,
-                                          AprobWght = 1e-06,
-                                          SlopeWght = 0.01,
-                                          SaveDir = "Results",
-                                          verbose = FALSE)
+Sable_mod <- AgeingError::DoApplyAgeError(
+  Species = "Sable",
+  DataSpecs = Sable_dat,
+  ModelSpecsInp = Sable_spc,
+  AprobWght = 1e-06,
+  SlopeWght = 0.01,
+  SaveDir = "Results",
+  verbose = FALSE
+)
 ## save the model results
 Sable_out <- AgeingError:::ProcessResults(Species = "Sable", SaveDir = "Results", CalcEff = TRUE, verbose = FALSE)
 ```