.Rhistory

# GR1 (G. radiolyrata)
GR1 <- subset(Johnson_data_long, Specimen_simple == "GR1")
GR1sg <- data.frame(
D = as.numeric(GR1$Height) * 1000,
d18Oc = as.numeric(GR1$d18O),
d18Oc_err = GR1$d18Oc_err,
D_err = GR1$D_err * 1000,
YEARMARKER = 0
)
x11(); ggplot(data = GR1sg) + geom_line(aes(x = D, y = d18Oc))
GR1sg$YEARMARKER[which(GR1sg$D  %in% c(28500, 36400, 43400, 48500))] <- 1 # Add manual year markers
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1",
input_from_file = FALSE, # Should input be read from a file?
object_name = GR1sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1"
)
GR1_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1/Age_model_results.csv")
View(GR1)
View(GR1_ShellChronage)
GR1_ShellChronage$Sample_nr2 <- GR1_ShellChronage$D / 1000
GR1_ShellChronage$Specimen <- "GR1"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(GR1_ShellChronage, Specimen, mean.day, se.day, Sample_nr2), by = c("Specimen", "Sample_nr2"))
GR1_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1/Age_model_results.csv")
GR1_ShellChronage$Height <- GR1_ShellChronage$D / 1000
GR1_ShellChronage$Specimen <- "GR1"
View(GR1_ShellChronage)
# Read in results of ShellChron for further processing
GR1_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1/Age_model_results.csv")
GR1_ShellChronage$Height <- GR1_ShellChronage$D / 1000
GR1_ShellChronage$Specimen_simple <- "GR1"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(GR1_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
View(Johnson_data_long)
Johnson_data <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson2022_isotope_data.csv", header = FALSE)
Johnson_data[Johnson_data == ""] <- NA
# Lengthen data and label specimens
Johnson_data_long <- data.frame(
Height = rep(NA, ncol(Johnson_data) / 3 * (nrow(Johnson_data) - 2)),
d13C = rep(NA, ncol(Johnson_data) / 3 * (nrow(Johnson_data) - 2)),
d18O = rep(NA, ncol(Johnson_data) / 3 * (nrow(Johnson_data) - 2)),
Specimen = rep(NA, ncol(Johnson_data) / 3 * (nrow(Johnson_data) - 2))
)
for(i in seq(1, ncol(Johnson_data) / 3, 1)){
Johnson_data_long[(1 + (i - 1) * (nrow(Johnson_data) - 2)):(i * (nrow(Johnson_data) - 2)), 1:3] <- Johnson_data[3:nrow(Johnson_data), (1 + (i - 1) * 3):(3 + (i - 1) * 3)]
Johnson_data_long$Specimen[(1 + (i - 1) * (nrow(Johnson_data) - 2)):(i * (nrow(Johnson_data) - 2))] <- Johnson_data[1, (1 + (i - 1) * 3)]
}
# Simplify specimen names
Johnson_data_long$Specimen_simple <- gsub(" .*", "", Johnson_data_long$Specimen)
# Remove empty rows
Johnson_data_long <- Johnson_data_long[complete.cases(Johnson_data_long), ]
# Add uncertainties on d18O and Height
Johnson_data_long$d18Oc_err <- 0.05 # per mille
Johnson_data_long$D_err <- 0.5 # mm
# Convert Height, d18O and d13C to numberic
Johnson_data_long$Height <- as.numeric(Johnson_data_long$Height)
Johnson_data_long$d13C <- as.numeric(Johnson_data_long$d13C)
Johnson_data_long$d18O <- as.numeric(Johnson_data_long$d18O)
# Add columns to store ShellChron results
Johnson_data_long$ShellChron_DOY <- NA # Create column to store Day of the Year results from ShellChron
Johnson_data_long$ShellChron_DOY_err <- NA # Create column to store uncertainties on Day of the Year results from ShellChron
GR1 <- subset(Johnson_data_long, Specimen_simple == "GR1")
GR1sg <- data.frame(
D = GR1$Height * 1000,
d18Oc = GR1$d18O,
d18Oc_err = GR1$d18Oc_err,
D_err = GR1$D_err * 1000,
YEARMARKER = 0
)
summary(GR1sq)
summary(GR1sg)
GR1sg$YEARMARKER[which(GR1sg$D  %in% c(28500, 36400, 43400, 48500))] <- 1 # Add manual year markers
# Read in results of ShellChron for further processing
GR1_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR1/Age_model_results.csv")
GR1_ShellChronage$Height <- GR1_ShellChronage$D / 1000
GR1_ShellChronage$Specimen_simple <- "GR1"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(GR1_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
View(joined)
View(Johnson_data_long)
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# GR2 (G. radiolyrata)
GR2 <- subset(Johnson_data_long, Specimen_simple == "GR2")
GR2sg <- data.frame(
D = GR2$Height * 1000,
d18Oc = GR2$d18O,
d18Oc_err = GR2$d18Oc_err,
D_err = GR2$D_err * 1000,
YEARMARKER = 0
)
x11(); ggplot(data = GR2sg) + geom_line(aes(x = D, y = d18Oc))
View(GR2sg)
GR2sg$YEARMARKER[which(GR2sg$D  %in% c(17900, 23900, 29900, 36200, 40600, 44700, 48400, 51900, 53800))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR2",
input_from_file = FALSE, # Should input be read from a file?
object_name = GR2sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR2"
)
# Read in results of ShellChron for further processing
GR2_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_GR2/Age_model_results.csv")
GR2_ShellChronage$Height <- GR2_ShellChronage$D / 1000
GR2_ShellChronage$Specimen_simple <- "GR2"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(GR2_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
View(Johnson_data_long)
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# devtools::install_github("nielsjdewinter/ShellChron")
require(tidyverse)
require(ShellChron)
Johnson_data <- read.csv("C:/Users/Niels/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson2022_isotope_data.csv", header = FALSE)
require(tidyverse)
require(ShellChron)
Johnson_data <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson2022_isotope_data.csv", header = FALSE)
Johnson_data[Johnson_data == ""] <- NA
# Pick up where we left off
Johnson_data_long <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
-------------------------------------------------------------------
# AO1 (G. radiolyrata)
AO1 <- subset(Johnson_data_long, Specimen_simple == "AO1")
View(Johnson_data_long)
AO1 <- subset(Johnson_data_long, Specimen_simple == "AO1")
# ------------------------------------------------------------------------------
# AO1 (G. radiolyrata)
AO1 <- subset(Johnson_data_long, Specimen_simple == "AO1")
AO1sg <- data.frame(
D = AO1$Height * 1000,
d18Oc = AO1$d18O,
d18Oc_err = AO1$d18Oc_err,
D_err = AO1$D_err * 1000,
YEARMARKER = 0
)
AO1sg$YEARMARKER[which(AO1sg$D  %in% c(16000, 41000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO1",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO1sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO1"
)
# Read in results of ShellChron for further processing
AO1_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO1/Age_model_results.csv")
AO1_ShellChronage$Height <- AO1_ShellChronage$D / 1000
AO1_ShellChronage$Specimen_simple <- "AO1"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AO1_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# AO2 (G. radiolyrata)
AO2 <- subset(Johnson_data_long, Specimen_simple == "AO2")
AO2sg <- data.frame(
D = AO2$Height * 1000,
d18Oc = AO2$d18O,
d18Oc_err = AO2$d18Oc_err,
D_err = AO2$D_err * 1000,
YEARMARKER = 0
)
AO2sg$YEARMARKER[which(AO2sg$D  %in% c(13100, 39000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO2",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO2sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO2"
)
# Read in results of ShellChron for further processing
AO2_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO2/Age_model_results.csv")
AO2_ShellChronage$Height <- AO2_ShellChronage$D / 1000
AO2_ShellChronage$Specimen_simple <- "AO2"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AO2_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# AO3 (G. radiolyrata)
AO3 <- subset(Johnson_data_long, Specimen_simple == "AO3")
AO3sg <- data.frame(
D = AO3$Height * 1000,
d18Oc = AO3$d18O,
d18Oc_err = AO3$d18Oc_err,
D_err = AO3$D_err * 1000,
YEARMARKER = 0
)
AO3sg$YEARMARKER[which(AO3sg$D  %in% c(13000, 42000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO3",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO3sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO3"
)
# ------------------------------------------------------------------------------
# AO4 (G. radiolyrata)
AO4 <- subset(Johnson_data_long, Specimen_simple == "AO4")
AO4sg <- data.frame(
D = AO4$Height * 1000,
d18Oc = AO4$d18O,
d18Oc_err = AO4$d18Oc_err,
D_err = AO4$D_err * 1000,
YEARMARKER = 0
)
AO4sg$YEARMARKER[which(AO4sg$D  %in% c(5000, 42000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO4",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO4sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO4"
)
# ------------------------------------------------------------------------------
# AO6 (G. radiolyrata)
AO6 <- subset(Johnson_data_long, Specimen_simple == "AO6")
AO6sg <- data.frame(
D = AO6$Height * 1000,
d18Oc = AO6$d18O,
d18Oc_err = AO6$d18Oc_err,
D_err = AO6$D_err * 1000,
YEARMARKER = 0
)
AO6sg$YEARMARKER[which(AO6sg$D  %in% c(7000, 34000, 49000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO6",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO6sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO6"
)
# Read in results of ShellChron for further processing
AO6_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO6/Age_model_results.csv")
AO6_ShellChronage$Height <- AO6_ShellChronage$D / 1000
AO6_ShellChronage$Specimen_simple <- "AO6"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AO6_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# AO7 (G. radiolyrata)
AO7 <- subset(Johnson_data_long, Specimen_simple == "AO7")
AO7sg <- data.frame(
D = AO7$Height * 1000,
d18Oc = AO7$d18O,
d18Oc_err = AO7$d18Oc_err,
D_err = AO7$D_err * 1000,
YEARMARKER = 0
)
AO7sg$YEARMARKER[which(AO7sg$D  %in% c(6000, 34000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO7",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO7sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO7"
)
# ------------------------------------------------------------------------------
# AO9 (G. radiolyrata)
AO9 <- subset(Johnson_data_long, Specimen_simple == "AO9")
AO9sg <- data.frame(
D = AO9$Height * 1000,
d18Oc = AO9$d18O,
d18Oc_err = AO9$d18Oc_err,
D_err = AO9$D_err * 1000,
YEARMARKER = 0
)
AO9sg$YEARMARKER[which(AO9sg$D  %in% c(9000, 42000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO9",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO9sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO9"
)
# Read in results of ShellChron for further processing
AO9_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO9/Age_model_results.csv")
AO9_ShellChronage$Height <- AO9_ShellChronage$D / 1000
AO9_ShellChronage$Specimen_simple <- "AO9"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AO9_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# AO10 (G. radiolyrata)
AO10 <- subset(Johnson_data_long, Specimen_simple == "AO10")
AO10sg <- data.frame(
D = AO10$Height * 1000,
d18Oc = AO10$d18O,
d18Oc_err = AO10$d18Oc_err,
D_err = AO10$D_err * 1000,
YEARMARKER = 0
)
AO10sg$YEARMARKER[which(AO10sg$D  %in% c(7000, 34000, 47000))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO10",
input_from_file = FALSE, # Should input be read from a file?
object_name = AO10sg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO10"
)
# Read in results of ShellChron for further processing
AO10_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AO10/Age_model_results.csv")
AO10_ShellChronage$Height <- AO10_ShellChronage$D / 1000
AO10_ShellChronage$Specimen_simple <- "AO10"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AO10_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# AI (A. islandica)
AI <- subset(Johnson_data_long, Specimen_simple == "AI")
AIsg <- data.frame(
D = AI$Height * 1000,
d18Oc = AI$d18O,
d18Oc_err = AI$d18Oc_err,
D_err = AI$D_err * 1000,
YEARMARKER = 0
)
AIsg$YEARMARKER[which(AIsg$D  %in% c(31200, 37700, 47300))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AI",
input_from_file = FALSE, # Should input be read from a file?
object_name = AIsg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AI"
)
# Read in results of ShellChron for further processing
AI_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_AI/Age_model_results.csv")
AI_ShellChronage$Height <- AI_ShellChronage$D / 1000
AI_ShellChronage$Specimen_simple <- "AI"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(AI_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertainty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")
# ------------------------------------------------------------------------------
# PR (P. rustica)
PR <- subset(Johnson_data_long, Specimen_simple == "PR")
PRsg <- data.frame(
D = PR$Height * 1000,
d18Oc = PR$d18O,
d18Oc_err = PR$d18Oc_err,
D_err = PR$D_err * 1000,
YEARMARKER = 0
)
PRsg$YEARMARKER[which(PRsg$D  %in% c(32300, 39800, 46500))] <- 1 # Add manual year markers
# Apply ShellChron
wrap_function(
path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_PR",
input_from_file = FALSE, # Should input be read from a file?
object_name = PRsg, # Name of object with input (only if input_from_file = FALSE)
transfer_function = "KimONeil97", # Set transfer function of the record, default is Kim and O'Neil 1997.
t_int = 1, # Set time interval in days
T_per = 365, # Set annual time period in days (default = 365)
d18Ow = 0, # Set d18Ow value or vector (default = constant year-round at 0 VSMOW). Alternative options are either one value (assumed constant year-round) or a vector with length T_per / t_int and interval t_int specifying d18Ow evolution through one year.
t_maxtemp = 182.5, # Define the day of the year at which temperature is heighest. Default = Assume that the day of maximum temperature is helfway through the year
SCEUApar = c(1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
sinfit = TRUE, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (TRUE/FALSE)
MC = 1000, # Number of MC simulations to include measurement error into error analysis. Default = 1000 (if MC = 0, error on D and d18O measurements not considered)
plot = TRUE, # Should intermediate plots be given to track progress? WARNING: plotting makes the script much slower, especially for long datasets.
plot_export = TRUE, # Should a plot of the results be saved as PDF?
export_raw = TRUE, # Should the results of all individual model runs be exported as CSV files?
export_path = "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_PR"
)
# Read in results of ShellChron for further processing
PR_ShellChronage <- read.csv("E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/ShellChron_PR/Age_model_results.csv")
PR_ShellChronage$Height <- PR_ShellChronage$D / 1000
PR_ShellChronage$Specimen_simple <- "PR"
# Add results to data sheet
joined <- left_join(Johnson_data_long, select(PR_ShellChronage, Specimen_simple, mean.day, se.day, Height), by = c("Specimen_simple", "Height"))
Johnson_data_long$ShellChron_DOY[which(!is.na(joined$mean.day))] <- joined$mean.day[which(!is.na(joined$mean.day))] %% 365 # Convert model result to day of year and add to column
Johnson_data_long$ShellChron_DOY_err[which(!is.na(joined$se.day))] <- joined$se.day[which(!is.na(joined$se.day))] # Add model uncertPRnty to column
# Save dated data intermediately
write.csv(Johnson_data_long, "E:/Dropbox/Research/postdoc/UNBIAS/PWP reconstructions/Combined_data/data_Johnson2022/Johnson_data_dated.csv")