Export_results: Function to merge and export the results of the ShellChron...

In nielsjdewinter/ShellChron: Builds Chronologies from Oxygen Isotope Profiles in Shells

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

# 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[, -grep("window", colnames(resultarray[, , 3])), ] <- 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)

nielsjdewinter/ShellChron documentation built on Aug. 21, 2022, 4:31 p.m.