ShellChron (version 0.4.0)

export_results: Function to merge and export the results of the ShellChron model

Description

Takes the input data and model results and reformats them to tables of key parameters such as growth rate and shell age for each datapoint for easy plotting. This final function also combines uncertainties in the model result arising from uncertainties in input data (provided by the user) and uncertainties of the model (from overlapping modeling windows). Includes some optional plotting options.

Usage

export_results(
  path = getwd(),
  dat,
  resultarray,
  parmat,
  MC = 1000,
  dynwindow,
  plot = FALSE,
  plot_export = TRUE,
  export_raw = FALSE
)

Arguments

path

Path where result files are exported

dat

Matrix containing the input data

resultarray

Array containing the full results of the optimized growth model

parmat

Matrix listing all optimized growth rate and SST parameters used to model d18O in each data window

MC

Number of Monte Carlo simulations to apply for error propagation. Default = 1000

dynwindow

Information on the position and length of modeling windows

plot

Should an overview of the results of modeling be plotted? TRUE/FALSE

plot_export

Should the overview plot be exported as a PDF file? TRUE/FALSE

export_raw

Export tables containing all raw model results before being merged into tidy tables? TRUE/FALSE

Value

CSV tables of model results in the current working directory + optional plots in PDF format

References

package dependencies: tidyverse 1.3.0; ggpubr 0.4.0; magrittr function dependencies: sd_wt

Examples

# NOT RUN {
# Create dummy input data column by column
dat <- as.data.frame(seq(1000, 40000, 1000))
colnames(dat) <- "D"
dat$d18Oc <- sin((2 * pi * (seq(1, 40, 1) - 8 + 7 / 4)) / 7)
dat$YEARMARKER <- c(0, rep(c(0, 0, 0, 0, 0, 0, 1), 5), 0, 0, 0, 0)
dat$D_err <- rep(100, 40)
dat$d18Oc_err <- rep(0.1, 40)

testarray <- array(NA, dim = c(40, 36, 9)) # Create empty array
# with correct third dimension
windowfill <- seq(50, 500, 50) %% 365 # Create dummy simulation data 
# (ages) to copy through the array
for(i in 6:length(testarray[1, , 1])){
    testarray[, i, 3] <- c(windowfill, rep(NA, length(testarray[, 1, 3]) -
        length(windowfill)))
    windowfill <- c(NA, (windowfill + 51) %% 365)
}
# Add dummy /code{D} column.
testarray[, 1, 3] <- seq(1, length(testarray[, 1, 3]), 1)
# Add dummy YEARMARKER column
testarray[, 3, 3] <- c(0, rep(c(0, 0, 0, 0, 0, 0, 1), 5), 0, 0, 0, 0)
# Add dummy d18Oc column
testarray[, 2, 3] <- sin((2 * pi * (testarray[, 1, 3] - 8 + 7 / 4)) / 7)
# Create dummy seasonality data
seas <- as.data.frame(seq(1, 365, 1))
colnames(seas) <- "t"
seas$SST <- 15 + 10 * sin((2 * pi * (seq(1, 365, 1) - 182.5 +
    365 / 4)) / 365)
seas$GR <- 10 + 10 * sin((2 * pi * (seq(1, 365, 1) - 100 + 365 / 4)) / 365)
seas$d18O <- (exp((18.03 * 1000 / (seas$SST + 273.15) - 32.42) / 1000) - 1) *
    1000 + (0.97002 * 0 - 29.98)
# Apply dummy seasonality data to generate other tabs of testarray
testarray[, , 1] <- seas$d18O[match(testarray[, , 3], seas$t)] # d18O values
tab <- testarray[, , 1]
tab[which(!is.na(tab))] <- 0.1
testarray[, , 2] <- tab # dummy d18O residuals
testarray[, , 4] <- seas$GR[match(testarray[, , 3], seas$t)] # growth rates
testarray[, , 5] <- seas$SST[match(testarray[, , 3], seas$t)] # temperature
tab[which(!is.na(tab))] <- 0.1
testarray[, , 6] <- tab # dummy d18O SD
tab[which(!is.na(tab))] <- 20
testarray[, , 7] <- tab # dummy time SD
tab[which(!is.na(tab))] <- 3
testarray[, , 8] <- tab # dummy GR SD
tab[which(!is.na(tab))] <- 1
testarray[, , 9] <- tab # dummy temperature SD
darray <- array(rep(as.matrix(dat), 9), dim = c(40, 5, 9))
testarray[, 1:5, ] <- darray

# Create dummy dynwindow data
dynwindow <- as.data.frame(seq(1, 31, 1))
colnames(dynwindow) <- "x"
dynwindow$y <- rep(10, 31)

dimnames(testarray) <- list(
    paste("sample", 1:length(testarray[, 1, 3])),
    c(colnames(dat), paste("window", 1:length(dynwindow$x))),
    c("Modeled_d18O",
        "d18O_residuals",
        "Time_of_year",
        "Instantaneous_growth_rate",
        "Modeled temperature",
        "Modeled_d18O_SD",
        "Time_of_Year_SD",
        "Instantaneous_growth_rate_SD",
        "Modeled_temperature_SD")
)

# Set parameters
G_amp <- 20
G_per <- 365
G_pha <- 100
G_av <- 15
G_skw <- 70
T_amp <- 20
T_per <- 365
T_pha <- 150
T_av <- 15
pars <- c(T_amp, T_pha, T_av, G_amp, G_pha, G_av, G_skw)
parsSD <- c(3, 10, 3, 5, 10, 3, 5) # Artificial variability in parameters
parmat <- matrix(rnorm(length(pars) * length(dynwindow$x)), nrow =
    length(pars)) * parsSD + matrix(rep(pars, length(dynwindow$x)),
    nrow = length(pars))
rownames(parmat) <- c("T_amp", "T_pha", "T_av", "G_amp", "G_pha", "G_av",
    "G_skw")
# Run export function
test <- export_results(path = tempdir(),
    dat,
    testarray,
    parmat,
    MC = 1000,
    dynwindow,
    plot = FALSE,
    plot_export = FALSE,
    export_raw = FALSE)
# }