From 925b96a08fdadf7ef39fcde135718fbd8e70559e Mon Sep 17 00:00:00 2001 From: Niels J de Winter Date: Thu, 2 Jun 2022 14:44:22 +0200 Subject: [PATCH] Updated data_import with check of duplicate depth values. Modified Run_model to take modelled sinusoidal parameters of previous window as starting parameters for the next window --- .Rhistory | 1000 ++++++++++++++++++++++----------------------- R/Run_model.r | 77 ++-- R/Wrap_function.r | 2 +- R/data_import.r | 6 + 4 files changed, 547 insertions(+), 538 deletions(-) diff --git a/.Rhistory b/.Rhistory index 4bbc3e3..c17d3d9 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,512 +1,512 @@ -rownames(lmtable) <- c("intercept", "slope", "SD", "R2") -# Regress over timing -lm_ShellChron_LR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_windows_pred <- as.data.frame(predict(lm_ShellChron_LR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_windows_pred <- as.data.frame(predict(lm_ShellChron_HR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared), -) -rownames(lmtable) <- c("intercept", "slope", "SD", "R2") -lm_ShellChron_LR_windows$coefficients -c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared) -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared), -) -c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared) -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared) -) -# Regress over timing -lm_ShellChron_LR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_windows_pred <- as.data.frame(predict(lm_ShellChron_LR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_windows_pred <- as.data.frame(predict(lm_ShellChron_HR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared) -) -rownames(lmtable) <- c("intercept", "slope", "SD", "R2") -# Regress over timing -lm_ShellChron_LR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_windows_pred <- as.data.frame(predict(lm_ShellChron_LR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_windows_pred <- as.data.frame(predict(lm_ShellChron_HR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared) -) -rownames(lmtable_windows) <- c("intercept", "slope", "SD", "R2") -View(lmtable_windows) -# Regress over timing -lm_ShellChron_LR_samples <- lm(total/60 ~ samples, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_samples_pred <- as.data.frame(predict(lm_ShellChron_LR_samples, newdata = data.frame(samples = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_samples <- lm(total/60 ~ samples, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_samples_pred <- as.data.frame(predict(lm_ShellChron_HR_samples, newdata = data.frame(samples = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_samples <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_samples$coefficients, summary(lm_ShellChron_LR_samples)$sigma, summary(lm_ShellChron_LR_samples)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_samples$coefficients, summary(lm_ShellChron_HR_samples)$sigma, summary(lm_ShellChron_HR_samples)$r.squared) -) -rownames(lmtable_samples) <- c("intercept", "slope", "SD", "R2") -View(lmtable_samples) -rm(list=ls()) -setwd("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time") -# Calculate accuracy of Judd et al. (2018) model -t_real <- as.data.frame(read.csv("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_Realtime.csv", header = TRUE)) # Import real timing data -JD_real <- ((t_real + 0.25) * 365) %% 365 # Correct for phase lag in data and convert to julian day -JD_real_cum <- (t_real + 0.25) * 365 -# Create vector of datasets -recordnames <- c(paste("LR_", 3:12, sep = ""), paste("HR_", 3:12, sep = "")) -# Create dataframe to store timetrack results -timetracks <- data.frame(matrix(nrow = 20, ncol = 5)) -rownames(timetracks) <- recordnames -colnames(timetracks) <- paste("STEP", 1:5) -# Create dataframes to store differences between data and model results -Dt_Judd <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dt = rep(NA, length(which(!is.na(JD_real))) * 10) -) -Dt_Judd_cum <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dt = rep(NA, length(which(!is.na(JD_real))) * 10) -) -D18O_Judd <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dO = rep(NA, length(which(!is.na(JD_real))) * 10) -) -Dt_ShellChron <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dt = rep(NA, length(which(!is.na(JD_real))) * 10) -) -Dt_ShellChron_cum <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dt = rep(NA, length(which(!is.na(JD_real))) * 10) -) -D18O_ShellChron <- data.frame(record = rep(NA, length(which(!is.na(JD_real))) * 10), -dO = rep(NA, length(which(!is.na(JD_real))) * 10) -) -windows <- data.frame(matrix(nrow = 20, ncol = 2)) # Create matrix to store number of windows -for(i in 1:length(recordnames)){ -# Load timetrack data -timetrack <- read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/timetrack", recordnames[i], ".csv", sep = "")) -timetracks[i, ] <- timetrack$x -# Load input data -dat <- as.data.frame(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/Case1_", recordnames[i], ".csv", sep = ""), header = TRUE)) -# Load data Judd et al. model -Omod_Judd <- as.matrix(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/Judd_model/Omodsamp2.csv", sep = ""), header = FALSE)) # Load modelled d18O result -Omod_Judd <- as.vector(Omod_Judd[!is.na(Omod_Judd)]) -JDmod_Judd <- as.matrix(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/Judd_model/JDmodsampshift2.csv", sep = ""), header = FALSE)) # Load modelled time result -JDmod_Judd <- as.vector(JDmod_Judd[!is.na(JDmod_Judd)]) -# Create cumulative age from Judd et al. data -YM <- which(diff(JDmod_Judd) < -1) # Find year transitions, skip "years" in the model smaller than 1/2 year -x <- rep(0, length(JDmod_Judd)) -for(y in YM){ -x[(y + 1):length(x)] <- x[(y + 1):length(x)] + 1 -} -JDmod_Judd_cum <- JDmod_Judd + x * 365 -# Load ShellChron data -AMdat <- as.data.frame(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/Age_model_results.csv", sep = ""), header = TRUE)) # Load modelled time result -JDmod_ShellChron <- as.vector(AMdat$mean.day) -d18Odat <- as.data.frame(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/d18O_model_results.csv", sep = ""), header = TRUE)) # Load modelled d18O result -Omod_ShellChron <- as.vector(d18Odat$mean.d18O_mod) -JDraw <- as.data.frame(read.csv(paste("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_", recordnames[i], "/Day_of_year_raw.csv", sep = ""), header = TRUE)) # Load raw age model result -windows[i, ] <- c(ncol(JDraw) - 6, recordnames[i]) -# Calculate diffences between models and real data -if(grepl("LR", recordnames[i], fixed = TRUE)){ # Separate LR from HR datasets -# Judd et al., model -dt_Judd <- JDmod_Judd - JD_real$t_LR[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate time difference between model and data -dt_Judd[which(abs(dt_Judd) > 182.5)] <- 182.5 - dt_Judd[which(abs(dt_Judd) > 182.5)] # If difference is more than 1/2 year, take the other half -dt_Judd_cum <- JDmod_Judd_cum - JD_real_cum$t_LR[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate time difference between model and data -dO_Judd <- Omod_Judd - dat$d18Oc[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate d18Oc difference between model and data -# ShellChron -dt_ShellChron <- (JDmod_ShellChron %% 365) - JD_real$t_LR[1:length(JDmod_ShellChron)] # Calculate time difference between model and data -dt_ShellChron[which(abs(dt_ShellChron) > 182.5)] <- 182.5 - dt_ShellChron[which(abs(dt_ShellChron) > 182.5)] # If difference is more than 1/2 year, take the other half -dt_ShellChron_cum <- JDmod_ShellChron - JD_real_cum$t_LR[1:length(JDmod_ShellChron)] # Calculate time difference between model and data -dO_ShellChron <- Omod_ShellChron - dat$d18Oc # Calculate d18Oc difference between model and data -}else{ -dt_Judd <- JDmod_Judd - JD_real$t_HR[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate time difference between model and data -dt_Judd[which(abs(dt_Judd) > 182.5)] <- 182.5 - dt_Judd[which(abs(dt_Judd) > 182.5)] # If difference is more than 1/2 year, take the other half -dt_Judd_cum <- JDmod_Judd_cum - JD_real_cum$t_HR[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate time difference between model and data -dO_Judd <- Omod_Judd - dat$d18Oc[min(which(dat$YEARMARKER == 1)):(max(which(dat$YEARMARKER == 1)) - 1)] # Calculate d18Oc difference between model and data -# ShellChron -dt_ShellChron <- (JDmod_ShellChron %% 365) - JD_real$t_HR[1:length(JDmod_ShellChron)] # Calculate time difference between model and data -dt_ShellChron[which(abs(dt_ShellChron) > 182.5)] <- 182.5 - dt_ShellChron[which(abs(dt_ShellChron) > 182.5)] # If difference is more than 1/2 year, take the other half -dt_ShellChron_cum <- JDmod_ShellChron - JD_real_cum$t_HR[1:length(JDmod_ShellChron)] # Calculate time difference between model and data -dO_ShellChron <- Omod_ShellChron - dat$d18Oc # Calculate d18Oc difference between model and data -} -# Add results to dataframes -# Judd et al., model -Dt_Judd$dt[min(which(is.na(Dt_Judd$record))):(min(which(is.na(Dt_Judd$record))) + length(dt_Judd) - 1)] <- dt_Judd -Dt_Judd$record[min(which(is.na(Dt_Judd$record))):(min(which(is.na(Dt_Judd$record))) + length(dt_Judd) - 1)] <- rep(recordnames[i], length(dt_Judd)) -Dt_Judd_cum$dt[min(which(is.na(Dt_Judd_cum$record))):(min(which(is.na(Dt_Judd_cum$record))) + length(dt_Judd_cum) - 1)] <- dt_Judd_cum -Dt_Judd_cum$record[min(which(is.na(Dt_Judd_cum$record))):(min(which(is.na(Dt_Judd_cum$record))) + length(dt_Judd_cum) - 1)] <- rep(recordnames[i], length(dt_Judd_cum)) -D18O_Judd$dO[min(which(is.na(D18O_Judd$record))):(min(which(is.na(D18O_Judd$record))) + length(dO_Judd) - 1)] <- dO_Judd -D18O_Judd$record[min(which(is.na(D18O_Judd$record))):(min(which(is.na(D18O_Judd$record))) + length(dO_Judd) - 1)] <- rep(recordnames[i], length(dO_Judd)) -# ShellChron -Dt_ShellChron$dt[min(which(is.na(Dt_ShellChron$record))):(min(which(is.na(Dt_ShellChron$record))) + length(dt_ShellChron) - 1)] <- dt_ShellChron -Dt_ShellChron$record[min(which(is.na(Dt_ShellChron$record))):(min(which(is.na(Dt_ShellChron$record))) + length(dt_ShellChron) - 1)] <- rep(recordnames[i], length(dt_ShellChron)) -Dt_ShellChron_cum$dt[min(which(is.na(Dt_ShellChron_cum$record))):(min(which(is.na(Dt_ShellChron_cum$record))) + length(dt_ShellChron_cum) - 1)] <- dt_ShellChron_cum -Dt_ShellChron_cum$record[min(which(is.na(Dt_ShellChron_cum$record))):(min(which(is.na(Dt_ShellChron_cum$record))) + length(dt_ShellChron_cum) - 1)] <- rep(recordnames[i], length(dt_ShellChron_cum)) -D18O_ShellChron$dO[min(which(is.na(D18O_ShellChron$record))):(min(which(is.na(D18O_ShellChron$record))) + length(dO_ShellChron) - 1)] <- dO_ShellChron -D18O_ShellChron$record[min(which(is.na(D18O_ShellChron$record))):(min(which(is.na(D18O_ShellChron$record))) + length(dO_ShellChron) - 1)] <- rep(recordnames[i], length(dO_ShellChron)) -} -timetracks$total <- apply(timetracks, 1, sum) # Add column for total elapsed time -colnames(windows) <- c("windows", "record") # Rename column names for windows dataframe -Dt_Judd <- Dt_Judd[complete.cases(Dt_Judd), ] -# Summarize statistics -Dtstats_Judd <- Dt_Judd[complete.cases(Dt_Judd), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dt), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dt), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dt), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_mean = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "mean"), # Calculate means per sample +sd = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +N = sum(Nbin, na.rm = TRUE), # Calculate the number of modelled values, excluding NA's se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level +CL95 = qt(0.95, N) * se, # Calculate the 95% confidence level +param49_mean = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "mean"), # Calculate means per sample +param49_sd = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +param49_se = param49_sd / sqrt(N), +param49_CL95 = qt(0.95, N) * param49_se ) -Dt_cumstats_Judd <- Dt_Judd_cum[complete.cases(Dt_Judd_cum), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dt), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dt), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dt), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's +# write.csv(D47stats, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Aragonite_dataset_stats_grouped.csv") +# Summarize stats for Arctica islandica samples per specimen +# First propagate uncertainties on bins per specimen, keeping data from different machines separate +Aisstats_Machinebin <- dat[which(dat$D47_outlier != TRUE & dat$Analysis == "this study"), ] %>% # Summarize D47 statistics +group_by(Specimen, Machine) %>% +summarize( +sample = first(sample), +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_binmean = mean(D47, na.rm = TRUE), # Calculate means per sample +SD_bin = sd(D47, na.rm = TRUE), # Calculate stdevs per sample bin +SD_ext = first(D47_SD), # External standard deviation (reported or based on check STD) +binsd = max(SD_bin, SD_ext), # Find largest standard +Nbin = n(), # Calculate the number of modelled values, excluding NA's +binse = binsd / sqrt(Nbin), # Calculate the standard error +binCL95 = qt(0.95, Nbin) * binse, # Calculate the 95% confidence level +param49_binmean = mean(param49, na.rm = TRUE), # Propagate param49 mean per sample +param49_binsd = sd(param49, na.rm = TRUE), # Propagate param49 SD per sample +Machine = first(Machine) +) %>% +ungroup() +# write.csv(Aisstats_Machinebin, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_dataset_stats_by_Machine.csv") +# Now combine sample groups from different machines into one bin +Aisstats <- Aisstats_Machinebin %>% +group_by(Specimen) %>% +summarize( +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_mean = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "mean"), # Calculate means per sample +sd = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +N = sum(Nbin, na.rm = TRUE), # Calculate the number of modelled values, excluding NA's se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level +CL95 = qt(0.95, N) * se, # Calculate the 95% confidence level +param49_mean = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "mean"), # Calculate means per sample +param49_sd = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +param49_se = param49_sd / sqrt(N), +param49_CL95 = qt(0.95, N) * param49_se ) -D18Ostats_Judd <- D18O_Judd[complete.cases(D18O_Judd), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dO), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dO), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dO), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's -se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level -) -Dtstats_ShellChron <- Dt_ShellChron[complete.cases(Dt_ShellChron), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dt), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dt), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dt), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's -se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level +# write.csv(Aisstats, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_dataset_stats_grouped.csv") +# Update dat with new propagated SDs for regressions +newSDs <- D47stats[which(D47stats$Analysis == "this study"), which(colnames(D47stats) %in% c("ID", "sd"))] +dat$D47_SD[which(!is.na(match(dat$ID, newSDs$ID)))] <- newSDs$sd[match(dat$ID, newSDs$ID)[which(!is.na(match(dat$ID, newSDs$ID)))]] +# Monte Carlo sample from individual aliquot distributions for violin plots and polynomial regression including errors on D47 +Nsim <- 10 ^ 4 +violin_data <- data.frame(sample = rep(dat$sample, Nsim), +type = rep(dat$type, Nsim), +Analysis = rep(dat$Analysis, Nsim), +ID = rep(dat$ID, Nsim), +Specimen = rep(dat$Specimen, Nsim), +Temp = rep(dat$Temp, Nsim), +Temp_sampled = rnorm(Nsim * length(dat$Temp), dat$Temp, dat$Temp_SD), +D47 = rnorm(Nsim * length(dat$D47), dat$D47, dat$D47_SD), +param49 = rnorm(Nsim * length(dat$param49), dat$param49, dat$param49_sd), +outlier = rep(dat$D47_outlier, Nsim) ) -Dt_cumstats_ShellChron <- Dt_ShellChron_cum[complete.cases(Dt_ShellChron_cum), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dt), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dt), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dt), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's +violin_data$x <- 10 ^ 6 / (violin_data$Temp_sampled + 273.15) ^ 2 # Create 10^6/T^2 vector +Aisdata <- dat[which(dat$D47_outlier != TRUE & dat$Analysis == "this study"), ] # Isolate Arctica islandica data from this study +Ais_temp_aov <- aov(D47 ~ ID, data = Aisdata) # Conduct one-way ANOVA on temperature bins +Ais_spec1_aov <- aov(D47 ~ Specimen, data = Aisdata[which(Aisdata$Temp == 1.1), ]) # Conduct one-way ANOVA on Specimen bins within the 1 degree temperature treatment +Ais_spec18_aov <- aov(D47 ~ Specimen, data = Aisdata[which(Aisdata$Temp == 18.0), ]) # Conduct one-way ANOVA on Specimen bins within the 18 degree temperature treatment +write.csv(TukeyHSD(Ais_temp_aov)$ID, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_temp.csv") # Export summary of Tukey multiple pairwise-comparisons +write.csv(TukeyHSD(Ais_spec1_aov)$Specimen, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec1.csv") # Export summary of Tukey multiple pairwise-comparisons +write.csv(TukeyHSD(Ais_spec18_aov)$Specimen, "E:/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec18.csv") # Export summary of Tukey multiple pairwise-comparisons +require(seasonalclumped) +install.packages("seasonalclumped") +?carbmodel +require(seasonalclumped) +?carbmodel +require(seasonalclumped) +setwd("C:/Users/niels/Dropbox/Research/Proposals/VU application/interview 2/teaching session") +dat <- seasonalclumped::Case15 +View(dat) +require(ggplot2) # for plotting +require(ggpubr) # for violin plots +require(tidyverse) # for data treatment +require(gridExtra) # for combining plots +require(bfsl) # for York regressions +require(RColorBrewer) # to define color scales +require(ggrepel) # for plot labeling +source("C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Average_error_propagation.r") +# -------------------------Temperature offsets---------------------------------- +source("C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Temperature_offset.r") +dat <- as.data.frame(read.csv("C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Aragonite_compilation_4.csv", header = TRUE)) +dat$sample <- paste(dat$Temp, dat$Analysis, sep = "_") +# First propagate uncertainties on bins per temperature, keeping data from different machines separate +D47stats_Machinebin <- dat[which(dat$D47_outlier != TRUE),] %>% # Summarize D47 statistics +group_by(ID, Machine) %>% +summarize( +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_binmean = mean(D47, na.rm = TRUE), # Calculate means per sample +SD_bin = sd(D47, na.rm = TRUE), # Calculate stdevs per sample bin +SD_ext = first(D47_SD), # External standard deviation (reported or based on check STD) +binsd = if(is.na(SD_bin)){first(SD_ext)}else{max(SD_bin, SD_ext)}, # Find largest standard deviation +Nbin = n(), # Calculate the number of modelled values, excluding NA's +binse = binsd / sqrt(Nbin), # Calculate the standard error +binCL95 = qt(0.95, Nbin) * binse, # Calculate the 95% confidence level +param49_binmean = mean(param49, na.rm = TRUE), # Propagate param49 mean per sample +param49_binsd = sd(param49, na.rm = TRUE), # Propagate param49 SD per sample +Machine = first(Machine) +) %>% +ungroup() +# write.csv(D47stats_Machinebin, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Aragonite_dataset_stats_by_Machine.csv") +# Now combine sample groups from different machines into one bin +D47stats <- D47stats_Machinebin %>% +group_by(ID) %>% +summarize( +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_mean = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "mean"), # Calculate means per sample +sd = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +N = sum(Nbin, na.rm = TRUE), # Calculate the number of modelled values, excluding NA's se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level +CL95 = qt(0.95, N) * se, # Calculate the 95% confidence level +param49_mean = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "mean"), # Calculate means per sample +param49_sd = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +param49_se = param49_sd / sqrt(N), +param49_CL95 = qt(0.95, N) * param49_se ) -D18Ostats_ShellChron <- D18O_ShellChron[complete.cases(D18O_ShellChron), ] %>% # Summarize modelled time offset statistics -ggpubr::group_by(record) %>% -dplyr::summarize( -mean = mean(abs(dO), na.rm = TRUE), # Calculate means per sample -sd = sd(abs(dO), na.rm = TRUE), # Calculate stdevs per sample -N = dplyr::n_distinct(abs(dO), na.rm = TRUE), # Calculate the number of modelled values, excluding NA's +# write.csv(D47stats, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Aragonite_dataset_stats_grouped.csv") +# Summarize stats for Arctica islandica samples per specimen +# First propagate uncertainties on bins per specimen, keeping data from different machines separate +Aisstats_Machinebin <- dat[which(dat$D47_outlier != TRUE & dat$Analysis == "this study"), ] %>% # Summarize D47 statistics +group_by(Specimen, Machine) %>% +summarize( +sample = first(sample), +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_binmean = mean(D47, na.rm = TRUE), # Calculate means per sample +SD_bin = sd(D47, na.rm = TRUE), # Calculate stdevs per sample bin +SD_ext = first(D47_SD), # External standard deviation (reported or based on check STD) +binsd = max(SD_bin, SD_ext), # Find largest standard +Nbin = n(), # Calculate the number of modelled values, excluding NA's +binse = binsd / sqrt(Nbin), # Calculate the standard error +binCL95 = qt(0.95, Nbin) * binse, # Calculate the 95% confidence level +param49_binmean = mean(param49, na.rm = TRUE), # Propagate param49 mean per sample +param49_binsd = sd(param49, na.rm = TRUE), # Propagate param49 SD per sample +Machine = first(Machine) +) %>% +ungroup() +# write.csv(Aisstats_Machinebin, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_dataset_stats_by_Machine.csv") +# Now combine sample groups from different machines into one bin +Aisstats <- Aisstats_Machinebin %>% +group_by(Specimen) %>% +summarize( +Temp = mean(Temp, na.rm = TRUE), +Temp_SD = first(Temp_SD), +Analysis = first(Analysis), +type = first(type), +sample = first(sample), +D47_mean = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "mean"), # Calculate means per sample +sd = binmeans(x = D47_binmean, x_sd = binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +N = sum(Nbin, na.rm = TRUE), # Calculate the number of modelled values, excluding NA's se = sd / sqrt(N), # Calculate the standard error -CL95 = qt(0.95, N) * se # Calculate the 95% confidence level +CL95 = qt(0.95, N) * se, # Calculate the 95% confidence level +param49_mean = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "mean"), # Calculate means per sample +param49_sd = binmeans(x = param49_binmean, x_sd = param49_binsd, n = Nbin, output = "SD"), # Calculate stdevs per sample bin +param49_se = param49_sd / sqrt(N), +param49_CL95 = qt(0.95, N) * param49_se ) -# ------------------------------------------------------------------------------ -# Plot results of performance time test -# Create plot of accuracy of cumulative model vs. number of modelled years -# Add numeric years and HR/LR record to Dt dataframes -Dt_cumstats_Judd$years <- as.numeric(sapply(Dt_cumstats_Judd$record, substring, first = 4)) -Dt_cumstats_Judd$res <- sapply(Dt_cumstats_Judd$record, substring, first = 1, last = 2) -Dt_cumstats_ShellChron$years <- as.numeric(sapply(Dt_cumstats_ShellChron$record, substring, first = 4)) -Dt_cumstats_ShellChron$res <- sapply(Dt_cumstats_ShellChron$record, substring, first = 1, last = 2) -accplot <- ggplot(data = Dt_cumstats_Judd, aes(years, mean)) + -geom_point(aes(col = res)) + -geom_line(aes(col = res), linetype = "dashed", size = 1) + -geom_errorbar(aes(x = years, ymin = mean - CL95, ymax = mean + CL95, col = res), width = 0.3) + -geom_point(data = Dt_cumstats_ShellChron, aes(col = res)) + -geom_line(data = Dt_cumstats_ShellChron, aes(col = res), size = 1) + -geom_errorbar(data = Dt_cumstats_ShellChron, aes(x = years, ymin = mean - CL95, ymax = mean + CL95, col = res), width = 0.3) + -scale_y_log10("offset from actual age (d)", breaks = c(seq(1, 10, 1), seq(20, 100, 10), seq(200, 1000, 100))) + -scale_x_continuous("Length of record (yr)", breaks = seq(3, 12, 1)) + -scale_colour_brewer(palette = "Paired") -# Import timing for Judd et al., model -Timing_Judd <- as.data.frame(read.csv("Timing_Judd.csv", header = TRUE)) -# Plot duration of modelling against record -# Extract years and record types for timing -timetracks$years <- as.numeric(sapply(rownames(timetracks), substring, first = 4)) -timetracks$res <- sapply(rownames(timetracks), substring, first = 1, last = 2) -# Regress over timing -lm_ShellChron_LR <- lm(total/60 ~ years, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_pred <- as.data.frame(predict(lm_ShellChron_LR, newdata = data.frame(years = seq(3, 12 ,0.5)), interval = "confidence")) -lm_ShellChron_HR <- lm(total/60 ~ years, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_pred <- as.data.frame(predict(lm_ShellChron_HR, newdata = data.frame(years = seq(3, 12 ,0.5)), interval = "confidence")) -lm_Judd_LR <- lm(Judd/60 ~ years, data = Timing_Judd[Timing_Judd$res == "LR", ]) -lm_Judd_HR <- lm(Judd/60 ~ years, data = Timing_Judd[Timing_Judd$res == "HR", ]) -# Group regression stats -lmtable <- data.frame(ShellChron_LR = c(lm_ShellChron_LR$coefficients, summary(lm_ShellChron_LR)$sigma, summary(lm_ShellChron_LR)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR$coefficients, summary(lm_ShellChron_HR)$sigma, summary(lm_ShellChron_HR)$r.squared), -Judd_LR = c(lm_Judd_LR$coefficients, summary(lm_Judd_LR)$sigma, summary(lm_Judd_LR)$r.squared), -Judd_HR = c(lm_Judd_HR$coefficients, summary(lm_Judd_HR)$sigma, summary(lm_Judd_HR)$r.squared) +# write.csv(Aisstats, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_dataset_stats_grouped.csv") +# Update dat with new propagated SDs for regressions +newSDs <- D47stats[which(D47stats$Analysis == "this study"), which(colnames(D47stats) %in% c("ID", "sd"))] +dat$D47_SD[which(!is.na(match(dat$ID, newSDs$ID)))] <- newSDs$sd[match(dat$ID, newSDs$ID)[which(!is.na(match(dat$ID, newSDs$ID)))]] +# Monte Carlo sample from individual aliquot distributions for violin plots and polynomial regression including errors on D47 +Nsim <- 10 ^ 4 +violin_data <- data.frame(sample = rep(dat$sample, Nsim), +type = rep(dat$type, Nsim), +Analysis = rep(dat$Analysis, Nsim), +ID = rep(dat$ID, Nsim), +Specimen = rep(dat$Specimen, Nsim), +Temp = rep(dat$Temp, Nsim), +Temp_sampled = rnorm(Nsim * length(dat$Temp), dat$Temp, dat$Temp_SD), +D47 = rnorm(Nsim * length(dat$D47), dat$D47, dat$D47_SD), +param49 = rnorm(Nsim * length(dat$param49), dat$param49, dat$param49_sd), +outlier = rep(dat$D47_outlier, Nsim) ) -rownames(lmtable) <- c("intercept", "slope", "SD", "R2") -Timeplot <- ggplot(timetracks, aes(years, total / 60)) + -stat_smooth(data = timetracks[timetracks$res == "LR", ], method = lm, aes(col = res)) + -stat_smooth(data = timetracks[timetracks$res == "HR", ], method = lm, aes(col = res)) + -geom_point(aes(col = res)) + -stat_smooth(data = Timing_Judd[Timing_Judd$res == "LR", ], method = lm, aes(y = Judd / 60, col = res), linetype = "dashed") + -stat_smooth(data = Timing_Judd[Timing_Judd$res == "HR", ], method = lm, aes(y = Judd / 60, col = res), linetype = "dashed") + -geom_point(data = Timing_Judd, aes(years, Judd / 60, col = res), shape = 1) + -scale_y_continuous("Total time (minutes)", breaks = seq(0, 225, 25)) + -scale_x_continuous("Length of record (yr)", breaks = seq(3, 12, 1)) -write.csv(lmtable, "Timing_years_regression_stats.csv") -# Compare real model accuracy with accuracy output -# Load AMdat for the longest HR record -AMdat <- as.data.frame(read.csv("E:/Dropbox/Research/Manuscripts/[Review] GMD - Bivalve age model/tests/revision/Performance time/Case1_HR_12/Age_model_results.csv", header = TRUE)) # Load modelled time result -Acccomp <- data.frame(D = AMdat$D, model = AMdat$CL95.day, real = abs(Dt_ShellChron_cum$dt[which(Dt_ShellChron_cum$record == "HR_12")])) -Acccomp <- Acccomp[complete.cases(Acccomp),] -acc_comp_plot <- ggplot(data = Acccomp, aes(D, model)) + -geom_ribbon(aes(x = D, ymin = model, ymax = real), fill = "grey", alpha = 0.3) + -geom_line(aes(D, real, col = "real offset"), size = 0.5) + -geom_line(aes(col = "model uncertainty"), size = 0.5) + -scale_x_continuous("Depth along record (mm)", breaks = seq(0, 120, 10) * 1000, labels = seq(0, 120, 10)) + -scale_y_continuous("Uncertainty (95% CL; days)", breaks = seq(0, 100, 10), limits = c(0, 100)) + -scale_colour_brewer(palette = "Paired") + -theme(legend.position = "NONE") -acc_box <- ggplot(data = stack(Acccomp[, -1])) + -stat_boxplot(aes(x = ind, y= values), geom = "errorbar", width = 0.5) + -geom_boxplot(aes(x = ind, y= values, fill = ind), size = 0.5, coef = 1.5, outlier.colour = "black", outlier.size = 1) + -scale_fill_brewer(palette = "Paired") + -scale_y_continuous("", breaks = seq(0, 50, 10), limits = c(0, 50)) + -xlab("") + -theme(legend.position = "none") -combined_acc_comp_plots <- grid.arrange(acc_comp_plot, acc_box, ncol = 2, nrow = 1, widths = c(8, 1)) -acc_hist <- ggplot(data = Acccomp) + -geom_histogram(aes(real, fill = "real offset"), alpha = 0.5, bins = 100) + -geom_histogram(aes(model, fill = "model uncertainty"), alpha = 0.5, bins = 100) + -scale_x_continuous("Uncertainty (95% CL; days)", breaks = seq(0, 50, 5), limits = c(0, 50)) + -scale_y_continuous("Frequency", breaks = NULL) + -scale_fill_brewer(palette = "Paired") -# Export results -write.csv(Dt_cumstats_Judd, "Summary_dt_Judd.csv") -write.csv(Dt_cumstats_ShellChron, "Summary_dt_ShellChron.csv") -write.csv(D18Ostats_Judd, "Summary_dO_Judd.csv") -write.csv(D18Ostats_ShellChron, "Summary_dO_ShellChron.csv") -write.csv(timetracks, "Timing_results.csv") -timetracks$samples <- Dtstats_ShellChron$N[match(rownames(timetracks), Dtstats_ShellChron$record)] # Add nr of samples to timetrack -timetracks$windows <- as.numeric(windows$windows) -Windowplot <- ggplot(timetracks, aes(windows, total / 60)) + -stat_smooth(data = timetracks[timetracks$res == "LR", ], method = lm, aes(col = res)) + -stat_smooth(data = timetracks[timetracks$res == "HR", ], method = lm, aes(col = res)) + -geom_point(aes(col = res)) + -scale_y_continuous("Total time (minutes)", breaks = seq(0, 225, 25)) + -scale_x_continuous("# Windows in record", breaks = seq(0, 600, 50)) -# Regress over timing -lm_ShellChron_LR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_windows_pred <- as.data.frame(predict(lm_ShellChron_LR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_windows <- lm(total/60 ~ windows, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_windows_pred <- as.data.frame(predict(lm_ShellChron_HR_windows, newdata = data.frame(windows = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_windows <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_windows$coefficients, summary(lm_ShellChron_LR_windows)$sigma, summary(lm_ShellChron_LR_windows)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_windows$coefficients, summary(lm_ShellChron_HR_windows)$sigma, summary(lm_ShellChron_HR_windows)$r.squared) +violin_data$x <- 10 ^ 6 / (violin_data$Temp_sampled + 273.15) ^ 2 # Create 10^6/T^2 vector +# ------------------------Pairwise comparisons---------------------------------- +Aisdata <- dat[which(dat$D47_outlier != TRUE & dat$Analysis == "this study"), ] # Isolate Arctica islandica data from this study +Ais_temp_aov <- aov(D47 ~ ID, data = Aisdata) # Conduct one-way ANOVA on temperature bins +capture.output(summary(Ais_temp_aov), file = "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_temp_summary.txt") # Print summary of ANOVA +TukeyHSD(Ais_temp_aov) # Print results of Tukey multiple pairwise-comparisons (post-hoc Tukey's test) on temperature bins +# write.csv(TukeyHSD(Ais_temp_aov)$ID, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_temp.csv") # Export summary of Tukey multiple pairwise-comparisons +Ais_spec1_aov <- aov(D47 ~ Specimen, data = Aisdata[which(Aisdata$Temp == 1.1), ]) # Conduct one-way ANOVA on Specimen bins within the 1 degree temperature treatment +capture.output(summary(Ais_spec1_aov), file = "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec1_summary.txt") # Print summary of ANOVA +TukeyHSD(Ais_spec1_aov) # Print results of Tukey multiple pairwise-comparisons (post-hoc Tukey's test) on specimen bins +# write.csv(TukeyHSD(Ais_spec1_aov)$Specimen, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec1.csv") # Export summary of Tukey multiple pairwise-comparisons +Ais_spec18_aov <- aov(D47 ~ Specimen, data = Aisdata[which(Aisdata$Temp == 18.0), ]) # Conduct one-way ANOVA on Specimen bins within the 18 degree temperature treatment +capture.output(summary(Ais_spec18_aov), file = "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec18_summary.txt") # Print summary of ANOVA +TukeyHSD(Ais_spec18_aov) # Print results of Tukey multiple pairwise-comparisons (post-hoc Tukey's test) on specimen bins +# write.csv(TukeyHSD(Ais_spec18_aov)$Specimen, "C:/Users/Niels/Dropbox//Research//postdoc//UNBIAS//Clumped Temperature Calibration/Ais_Pairwise_comp_spec18.csv") # Export summary of Tukey multiple pairwise-comparisons +# ----------------------------Regressions--------------------------------------- +# Start with simple lm +D47m <- lm(D47 ~ I(10^6 / (Temp + 273.15) ^ 2), data = dat, subset = which(dat$D47_outlier == FALSE)) +newdat <- data.frame(Temp = seq(0, 1000, 0.1)) +D47m_pred <- predict.lm(D47m, newdata = newdat, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_result <- cbind(10^6 / (newdat + 273.15) ^2, D47m_pred$fit) +# Repeat excluding the high-T datapoints (Müller et al., 2017) +D47m_lowT <- lm(D47 ~ I(10^6 / (Temp + 273.15) ^ 2), data = dat, subset = which(dat$D47_outlier == FALSE & dat$Analysis != "Muller17")) +newdat_lowT <- data.frame(Temp = seq(0, 100, 0.1)) +D47m_lowT_pred <- predict.lm(D47m_lowT, newdata = newdat_lowT, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_lowT_result <- cbind(10^6 / (newdat_lowT + 273.15) ^2, D47m_lowT_pred$fit) +# Then try linear York regression on both +D47m_York <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE] + dat$Temp_SD[dat$D47_outlier == FALSE] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE] - dat$Temp_SD[dat$D47_outlier == FALSE] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE], +r = 0) +newdat_York <- data.frame(x = 10 ^6 / (seq(0, 1000, 0.1) + 273.15) ^ 2) +D47m_York_pred <- predict(D47m_York, newdata = newdat_York, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_York_result <- cbind(newdat_York, D47m_York_pred$fit) +D47m_lowT_York <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"] + dat$Temp_SD[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"] - dat$Temp_SD[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & dat$Analysis != "Muller17"], +r = 0) +newdat_lowT_York <- data.frame(x = 10 ^6 / (seq(0, 100, 0.1) + 273.15) ^ 2) +D47m_lowT_York_pred <- predict(D47m_lowT_York, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_lowT_York_result <- cbind(newdat_lowT_York, D47m_lowT_York_pred$fit) +# Also calculate linear York regression on A. islandica and all bivalve data +Ais_York <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")] + dat$Temp_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")] - dat$Temp_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica")], +r = 0) +Ais_York_pred <- predict(Ais_York, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Ais_York_result <- cbind(newdat_lowT_York, Ais_York_pred$fit) +Ais_York_MOTU <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")] + dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")] - dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica")], +r = 0) +Ais_York_MOTU_pred <- predict(Ais_York_MOTU, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Ais_York_MOTU_result <- cbind(newdat_lowT_York, Ais_York_MOTU_pred$fit) +Ais_York_PACMAN <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")] + dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")] - dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica")], +r = 0) +Ais_York_PACMAN_pred <- predict(Ais_York_PACMAN, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Ais_York_PACMAN_result <- cbind(newdat_lowT_York, Ais_York_PACMAN_pred$fit) +Mollusk_York <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + dat$Temp_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] - dat$Temp_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & (dat$Analysis == "this study" | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +r = 0) +Mollusk_York_pred <- predict(Mollusk_York, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Mollusk_York_result <- cbind(newdat_lowT_York, Mollusk_York_pred$fit) +Mollusk_York_MOTU <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] - dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "MOTU") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +r = 0) +Mollusk_York_MOTU_pred <- predict(Mollusk_York_MOTU, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Mollusk_York_MOTU_result <- cbind(newdat_lowT_York, Mollusk_York_MOTU_pred$fit) +Mollusk_York_PACMAN <- bfsl(x = 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2, +y = dat$D47[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +sd_x = abs(10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2 - 10^6 / (dat$Temp[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] - dat$Temp_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")] + 273.15) ^ 2) / 2, +sd_y = dat$D47_SD[dat$D47_outlier == FALSE & ((dat$Analysis == "this study" & dat$Machine == "PACMAN") | dat$ID == "A.islandica" | dat$Analysis == "Caldarescu21")], +r = 0) +Mollusk_York_PACMAN_pred <- predict(Mollusk_York_PACMAN, newdata = newdat_lowT_York, se.fit = TRUE, interval = "confidence", level = 0.95) +Mollusk_York_PACMAN_result <- cbind(newdat_lowT_York, Mollusk_York_PACMAN_pred$fit) +# Third order polynomial regression following Müller et al., 2019 and Jautzy et al., 2020 +D47m_poly <- lm(D47 ~ poly(I(10^6 / (Temp + 273.15) ^ 2), 3), data = dat, subset = which(dat$D47_outlier == FALSE)) +D47m_poly_pred <- predict.lm(D47m_poly, newdata = newdat, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_poly_result <- cbind(10^6 / (newdat + 273.15) ^2, D47m_poly_pred$fit) +# Polynomial regression with errors on D47 and Temp using MC simulation +D47m_poly_MC <- lm(D47 ~ poly(x, 3), data = violin_data, subset = which(violin_data$outlier == FALSE)) +D47m_poly_MC_pred <- predict.lm(D47m_poly_MC, newdata = newdat_York, se.fit = TRUE, interval = "confidence", level = 0.95) +D47m_poly_MC_result <- cbind(newdat_York, D47m_poly_MC_pred$fit) +# Readjust 95% CL calculations for actual degrees of freedom +D47m_poly_MC_result$lwr <- D47m_poly_MC_result$fit - (D47m_poly_MC_result$fit - D47m_poly_MC_result$lwr) * sqrt(D47m_poly_MC$df.residual) / sqrt(length(violin_data$D47) / Nsim - 4) +D47m_poly_MC_result$upr <- D47m_poly_MC_result$fit + (D47m_poly_MC_result$upr - D47m_poly_MC_result$fit) * sqrt(D47m_poly_MC$df.residual) / sqrt(length(violin_data$D47) / Nsim - 4) +# -----------------------Add Preexisting calibrations--------------------------- +# Add dummy data to plot Anderson calibration +Anderson21 <- data.frame(Temp = 10 ^ 6 / (seq(0, 1000, 0.1) + 273.15) ^ 2, +D47 = 0.0391 * 10 ^ 6 / (seq(0, 1000, 0.1) + 273.15) ^ 2 + 0.154) # Latest Anderson calibration +# Add dummy data to plot Meinicke calibration +MeinickeICDES <- data.frame(Temp = 10 ^ 6 / (seq(0, 1000, 0.1) + 273.15) ^ 2, +D47 = 0.0397 * 10 ^ 6 / (seq(0, 1000, 0.1) + 273.15) ^ 2 + 0.1518) # Recalculated Meinicke calibration +# Add theoretical calcite and aragonite calibration lines by Guo et al. 2009 recalculated to the I-CDES scale +Guo09 <- data.frame(Temp = 10 ^ 6 / (seq(0, 1000, 0.1) + 273.15) ^ 2, +D47_cc_original = -3.33040 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 5.54042 / (seq(0, 1000, 0.1) + 273.15), # Original formula published by Guo et al., 2009 +D47_ar_original = -3.43068 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 6.90300 / (seq(0, 1000, 0.1) + 273.15), # Original formula published by Guo et al., 2009 +D47_cc_CDES25 = -3.33040 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 5.54042 / (seq(0, 1000, 0.1) + 273.15) + 0.23252 + 0.268 - 0.232, # Corrected to CDES reference fram by updating the D47-D63 fractionation factor from 0.232 to 0.268 +D47_ar_CDES25 = -3.43068 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 6.90300 / (seq(0, 1000, 0.1) + 273.15) + 0.22893 + 0.268 - 0.232, # Corrected to CDES reference fram by updating the D47-D63 fractionation factor from 0.232 to 0.268 +D47_cc_CDES90 = -3.33040 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 5.54042 / (seq(0, 1000, 0.1) + 273.15) + 0.23252 + 0.268 - 0.232 - 0.088, # Bring to CDES90 by applying the 25-70 degrees acid fractionation factor by Petersen et al., 2019 +D47_ar_CDES90 = -3.43068 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 6.90300 / (seq(0, 1000, 0.1) + 273.15) + 0.22893 + 0.268 - 0.232 - 0.088, # Bring to CDES90 by applying the 25-70 degrees acid fractionation factor by Petersen et al., 2019 +D47_cc_CDES90_corr = (-3.33040 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 5.54042 / (seq(0, 1000, 0.1) + 273.15) + 0.23252 + 0.268 - 0.232 - 0.088) * 1.035, # Correct for D47-dependent D47-D63 fractionation factor of 35 ppm/per mille found by Guo et al. 2009 (and implemented in Jautzy et al., 2020) +D47_ar_CDES90_corr = (-3.43068 * 10 ^ 9 / (seq(0, 1000, 0.1) + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (seq(0, 1000, 0.1) + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (seq(0, 1000, 0.1) + 273.15) ^ 2 - 6.90300 / (seq(0, 1000, 0.1) + 273.15) + 0.22893 + 0.268 - 0.232 - 0.088) * 1.035 # Correct for D47-dependent D47-D63 fractionation factor of 35 ppm/per mille found by Guo et al. 2009 (and implemented in Jautzy et al., 2020) ) -rownames(lmtable_windows) <- c("intercept", "slope", "SD", "R2") -write.csv(lmtable_windows, "Timing_windows_regression_stats.csv") -Sampleplot <- ggplot(timetracks, aes(samples, total / 60)) + -stat_smooth(data = timetracks[timetracks$res == "LR", ], method = lm, aes(col = res)) + -stat_smooth(data = timetracks[timetracks$res == "HR", ], method = lm, aes(col = res)) + -geom_point(aes(col = res)) + -scale_y_continuous("Total time (minutes)", breaks = seq(0, 225, 25)) + -scale_x_continuous("# Samples in record", breaks = seq(0, 600, 50)) -# Regress over timing -lm_ShellChron_LR_samples <- lm(total/60 ~ samples, data = timetracks[timetracks$res == "LR", ]) -lm_ShellChron_LR_samples_pred <- as.data.frame(predict(lm_ShellChron_LR_samples, newdata = data.frame(samples = seq(0, 600, 50)), interval = "confidence")) -lm_ShellChron_HR_samples <- lm(total/60 ~ samples, data = timetracks[timetracks$res == "HR", ]) -lm_ShellChron_HR_samples_pred <- as.data.frame(predict(lm_ShellChron_HR_samples, newdata = data.frame(samples = seq(0, 600, 50)), interval = "confidence")) -# Group regression stats -lmtable_samples <- data.frame(ShellChron_LR = c(lm_ShellChron_LR_samples$coefficients, summary(lm_ShellChron_LR_samples)$sigma, summary(lm_ShellChron_LR_samples)$r.squared), -ShellChron_HR = c(lm_ShellChron_HR_samples$coefficients, summary(lm_ShellChron_HR_samples)$sigma, summary(lm_ShellChron_HR_samples)$r.squared) +# Values used to put Guo et al. 2009 data on the I-CDES scale +# ETH-1 - Iso A (Meckler et al., 2014) = 600 degr with D47 of 0.2052 +# ETH-2 - Iso B (Meckler et al., 2014) = 600 degr with D47 of 0.2085 +# ETH-3 - Iso C (Meckler et al., 2014) = ~20 degr with D47 of 0.6132 +ETH1_GuoCDES <- Guo09$D47_ar_CDES25[which(Guo09$Temp == 10 ^ 6 / (600 + 273.15) ^ 2)] +ETH3_GuoCDES <- Guo09$D47_ar_CDES25[which(Guo09$Temp == 10 ^ 6 / (20 + 273.15) ^ 2)] +ETH1_ICDES <- 0.2052 +ETH3_ICDES <- 0.6132 +# Place on I-CDES scale +Guo09$D47_cc_ICDES <- (Guo09$D47_cc_CDES25 - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES # Use linear correction through ETH1 and ETH3 values to transform to I-CDES scale +Guo09$D47_ar_ICDES <- (Guo09$D47_ar_CDES25 - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES # Use linear correction through ETH1 and ETH3 values to transform to I-CDES scale +Guo09$D47_cc <- Guo09$D47_cc_ICDES * 1.035 # Correct for D47-dependent D47-D63 fractionation factor of 35 ppm/per mille found by Guo et al. 2009 (and implemented in Jautzy et al., 2020) +Guo09$D47_ar <- Guo09$D47_ar_ICDES * 1.035 # Correct for D47-dependent D47-D63 fractionation factor of 35 ppm/per mille found by Guo et al. 2009 (and implemented in Jautzy et al., 2020) +View(D47m_York_result) +View(D47m_lowT_York_result) +D47m_lowT_York_result$Diff_D47m_York <- D47m_lowT_York_result$fit - D47m_York_result$fit[1:1001] +plot(D47m_lowT_York_result$Diff_D47m_York) +D47m_lowT_York_result$err <- D47m_lowT_York_result$fit - D47m_lowT_York_result$lwr +D47m_York_result$err <- D47m_York_result$fit - D47m_York_result$lwr +View(Temperature_offset) +View(Temperature_offset) +Temperature_offset(0.68755, 0.0052, 0, 0.0449, 0.0291) +Temperature_offset(0.68755, 0.005, 0, 0.0449, 0.0291) +Temperature_offset(0.6925, 0.005, 0, 0.0449, 0.0291) +Temperature_offset(0.6925, 0.005, 0, 0.0450, 0.0896) +Temperature_offset(0.4128477, 0.010, 0, 0.0450, 0.0896) +Temperature_offset(0.4128477, 0.00399, 0, 0.0450, 0.0896) +# -------------------------Regression residuals--------------------------------- +# Extract D47 values and their errors for temperatures of all samples +newdat_York <- data.frame(x = 10 ^6 / (dat$Temp + 273.15) ^ 2) +dat$D47res_York <- dat$D47 - predict(D47m_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +dat$D47res_lowT_York <- dat$D47 - predict(D47m_lowT_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +dat$D47res_poly <- dat$D47 - predict(D47m_poly_MC, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +dat$D47res_Anderson <- dat$D47 - (0.0391 * 10 ^ 6 / (dat$Temp + 273.15) ^ 2 + 0.154) +dat$D47res_Meinicke <- dat$D47 - (0.0397 * 10 ^ 6 / (dat$Temp + 273.15) ^ 2 + 0.1518) +dat$D47res_Guo_cc <- dat$D47 - (((-3.33040 * 10 ^ 9 / (dat$Temp + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (dat$Temp + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (dat$Temp + 273.15) ^ 2 - 5.54042 / (dat$Temp + 273.15) + 0.23252 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +dat$D47res_Guo_ar <- dat$D47 - (((-3.43068 * 10 ^ 9 / (dat$Temp + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (dat$Temp + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (dat$Temp + 273.15) ^ 2 - 6.90300 / (dat$Temp + 273.15) + 0.22893 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +newdat_York <- data.frame(x = 10 ^6 / (D47stats$Temp + 273.15) ^ 2) +D47stats$D47res_York <- D47stats$D47_mean - predict(D47m_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +D47stats$D47res_lowT_York <- D47stats$D47_mean - predict(D47m_lowT_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +D47stats$D47res_poly <- D47stats$D47_mean - predict(D47m_poly_MC, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +D47stats$D47res_Anderson <- D47stats$D47_mean - (0.0391 * 10 ^ 6 / (D47stats$Temp + 273.15) ^ 2 + 0.154) +D47stats$D47res_Meinicke <- D47stats$D47_mean - (0.0397 * 10 ^ 6 / (D47stats$Temp + 273.15) ^ 2 + 0.1518) +D47stats$D47res_Guo_cc <- D47stats$D47_mean - (((-3.33040 * 10 ^ 9 / (D47stats$Temp + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (D47stats$Temp + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (D47stats$Temp + 273.15) ^ 2 - 5.54042 / (D47stats$Temp + 273.15) + 0.23252 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +D47stats$D47res_Guo_ar <- D47stats$D47_mean - (((-3.43068 * 10 ^ 9 / (D47stats$Temp + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (D47stats$Temp + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (D47stats$Temp + 273.15) ^ 2 - 6.90300 / (D47stats$Temp + 273.15) + 0.22893 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +newdat_York <- data.frame(x = 10 ^6 / (Aisstats$Temp + 273.15) ^ 2) +Aisstats$D47res_York <- Aisstats$D47_mean - predict(D47m_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +Aisstats$D47res_lowT_York <- Aisstats$D47_mean - predict(D47m_lowT_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +Aisstats$D47res_poly <- Aisstats$D47_mean - predict(D47m_poly_MC, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +Aisstats$D47res_Anderson <- Aisstats$D47_mean - (0.0391 * 10 ^ 6 / (Aisstats$Temp + 273.15) ^ 2 + 0.154) +Aisstats$D47res_Meinicke <- Aisstats$D47_mean - (0.0397 * 10 ^ 6 / (Aisstats$Temp + 273.15) ^ 2 + 0.1518) +Aisstats$D47res_Guo_cc <- Aisstats$D47_mean - (((-3.33040 * 10 ^ 9 / (Aisstats$Temp + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (Aisstats$Temp + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (Aisstats$Temp + 273.15) ^ 2 - 5.54042 / (Aisstats$Temp + 273.15) + 0.23252 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +Aisstats$D47res_Guo_ar <- Aisstats$D47_mean - (((-3.43068 * 10 ^ 9 / (Aisstats$Temp + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (Aisstats$Temp + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (Aisstats$Temp + 273.15) ^ 2 - 6.90300 / (Aisstats$Temp + 273.15) + 0.22893 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +newdat_York <- data.frame(x = 10 ^6 / (violin_data$Temp + 273.15) ^ 2) +violin_data$D47res_York <- violin_data$D47 - predict(D47m_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +violin_data$D47res_lowT_York <- violin_data$D47 - predict(D47m_lowT_York, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +violin_data$D47res_poly <- violin_data$D47 - predict(D47m_poly_MC, newdata = newdat_York, se.fit = TRUE, interval = "none", level = 0.95)$fit +violin_data$D47res_Anderson <- violin_data$D47 - (0.0391 * 10 ^ 6 / (violin_data$Temp + 273.15) ^ 2 + 0.154) +violin_data$D47res_Meinicke <- violin_data$D47 - (0.0397 * 10 ^ 6 / (violin_data$Temp + 273.15) ^ 2 + 0.1518) +violin_data$D47res_Guo_cc <- violin_data$D47 - (((-3.33040 * 10 ^ 9 / (violin_data$Temp + 273.15) ^ 4 + 2.32415 * 10 ^ 7 / (violin_data$Temp + 273.15) ^ 3 - 2.91282 * 10 ^ 3 / (violin_data$Temp + 273.15) ^ 2 - 5.54042 / (violin_data$Temp + 273.15) + 0.23252 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +violin_data$D47res_Guo_ar <- violin_data$D47 - (((-3.43068 * 10 ^ 9 / (violin_data$Temp + 273.15) ^ 4 + 2.35766 * 10 ^ 7 / (violin_data$Temp + 273.15) ^ 3 - 8.06003 * 10 ^ 2 / (violin_data$Temp + 273.15) ^ 2 - 6.90300 / (violin_data$Temp + 273.15) + 0.22893 + 0.268 - 0.232) - ETH1_GuoCDES) * (ETH1_ICDES-ETH3_ICDES)/(ETH1_GuoCDES-ETH3_GuoCDES) + ETH1_ICDES) * 1.035 +# Calculate values for calibrations relative to regressions through our dataset +D47m_York_result_res <- D47m_York_result - D47m_York_result$fit +D47m_York_result_res$x <- D47m_York_result_res$x + D47m_York_result$fit +D47m_York_result_res$Anderson <- Anderson21$D47 - D47m_York_result$fit +D47m_York_result_res$Meinicke <- MeinickeICDES$D47 - D47m_York_result$fit +D47m_York_result_res$Guo_cc <- Guo09$D47_cc - D47m_York_result$fit +D47m_York_result_res$Guo_ar <- Guo09$D47_ar - D47m_York_result$fit +D47m_lowT_York_result_res <- D47m_lowT_York_result - D47m_lowT_York_result$fit +D47m_lowT_York_result_res$x <- D47m_lowT_York_result_res$x + D47m_lowT_York_result$fit +D47m_lowT_York_result_res$Anderson <- Anderson21$D47[1:1001] - D47m_lowT_York_result$fit +D47m_lowT_York_result_res$Meinicke <- MeinickeICDES$D47[1:1001] - D47m_lowT_York_result$fit +D47m_lowT_York_result_res$Guo_cc <- Guo09$D47_cc[1:1001] - D47m_lowT_York_result$fit +D47m_lowT_York_result_res$Guo_ar <- Guo09$D47_ar[1:1001] - D47m_lowT_York_result$fit +D47m_poly_MC_result_res <- D47m_poly_MC_result - D47m_poly_MC_result$fit +D47m_poly_MC_result_res$x <- D47m_poly_MC_result_res$x + D47m_poly_MC_result$fit +D47m_poly_MC_result_res$Anderson <- Anderson21$D47 - D47m_poly_MC_result$fit +D47m_poly_MC_result_res$Meinicke <- MeinickeICDES$D47 - D47m_poly_MC_result$fit +D47m_poly_MC_result_res$Guo_cc <- Guo09$D47_cc - D47m_poly_MC_result$fit +D47m_poly_MC_result_res$Guo_ar <- Guo09$D47_ar - D47m_poly_MC_result$fit +# Prepare summary of calibration offsets +calibration_offset <- data.frame(D47_offset = c(dat$D47res_Anderson[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))], +dat$D47res_Meinicke[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))], +dat$D47res_Guo_cc[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))], +dat$D47res_Guo_ar[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))]), +D47_SD = rep(dat$D47_SD[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))], 4), +Calibration = c(rep("Anderson et al., 2021", length(dat$D47res_Meinicke[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))])), +rep("Meinicke et al., 2020", length(dat$D47res_Meinicke[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))])), +rep("Guo et al., 2009 (calcite)", length(dat$D47res_Meinicke[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))])), +rep("Guo et al., 2009 (aragonite)", length(dat$D47res_Meinicke[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))]))) ) -rownames(lmtable_samples) <- c("intercept", "slope", "SD", "R2") -write.csv(lmtable_samples, "Timing_samples_regression_stats.csv") -Combined_plots <- ggpubr::ggarrange( # Combine plots -years_plot, -windows_plot, -samples_plot, -labels = c("A", "B", "C"), -ncol = 3, -nrow = 1) -Combined_plots <- ggpubr::ggarrange( # Combine plots -Timeplot, -Windowplot, -Sampleplot, -labels = c("A", "B", "C"), -ncol = 3, -nrow = 1) -x11(); plot(Combined_plots) -pdf("Performance time plots.pdf", width = 20, height = 5) -print(Combined_plots) -dev.off() -View(timetracks) -write.csv(timetracks, "Timing_results.csv") -devtools::document() -devtools::install() -devtools::document() -devtools::install() -require(rhub) -install.packages("rhub") -devtools::document() -devtools::check_rhub() -?check -devtools::check() -?rnorm -a_mean = 0.0391 -a_sd = 0.0004 -b_mean = 0.154 -b_sd = 0.004 -D_mean = 0.65 -D_sd = 0.03 -a = rnorm(10^3, a_mean, a_sd) -b = rnorm(10^3, b_mean, b_sd) -D = rnorm(10^3, D_mean, D_sd) -x11(); hist(a) -a = rnorm(10^6, a_mean, a_sd) -b = rnorm(10^6, b_mean, b_sd) -D = rnorm(10^6, D_mean, D_sd) -x11(); hist(a) -?hist -x11(); hist(a, breaks = 1000) -T = sqrt(a * 10^6 / (D - b)) -x11(); hist(T, breaks = 1000) -mean(T) -mean(T)-273.15 -sd(T) -sd(T)/SQRT(10) -sd(T)/sqrt(10) -D_mean = c(0.65, 0.63, 0.71) -D_sd = c(0.03, 0.032, 0.038) -D = matrix(ncol = 3, nrow = 10^6) -for(i in 1:length(D_mean)){ -D[, i] <- rnorm(10^6, D_mean[i], D_sd[i]) -} -View(D) -D = matrix(ncol = 3, nrow = 10^6) -T = matrix(ncol = 3, nrow = 10^6) -for(i in 1:length(D_mean)){ -D[, i] <- rnorm(10^6, D_mean[i], D_sd[i]) -T[, i] = sqrt(a * 10^6 / (D[, i] - b)) -} -View(`T`) -mean(T) -mean(T)-273.15 -sd(T) -sd(T)/sqrt(3) -?qt -qt(0.95, 3) -N = 10^6 -a_mean = 0.0391 -a_sd = 0.0004 -b_mean = 0.154 -b_sd = 0.004 -D_mean = c(0.65, 0.63, 0.71) -D_sd = c(0.03, 0.032, 0.038) -a = rnorm(N, a_mean, a_sd) -b = rnorm(N, b_mean, b_sd) -D = matrix(ncol = 3, nrow = N) -T = matrix(ncol = 3, nrow = N) -for(i in 1:length(D_mean)){ -D[, i] <- rnorm(N, D_mean[i], D_sd[i]) -T[, i] = sqrt(a * N / (D[, i] - b)) -} -T_mean = mean(T) -T_sd = sd(T) -T_se = T_sd / sqrt(length(D_mean)) -T_95CL <- qt(0.95, length(D_mean)) * T_se -N = 10^6 -a_mean = 0.0391 -a_sd = 0.000 -b_mean = 0.154 -b_sd = 0.00 -D_mean = c(0.65, 0.63, 0.71) -D_sd = c(0.03, 0.032, 0.038) -a = rnorm(N, a_mean, a_sd) -b = rnorm(N, b_mean, b_sd) -D = matrix(ncol = 3, nrow = N) -T = matrix(ncol = 3, nrow = N) -for(i in 1:length(D_mean)){ -D[, i] <- rnorm(N, D_mean[i], D_sd[i]) -T[, i] = sqrt(a * N / (D[, i] - b)) -} -T_mean = mean(T) -T_sd = sd(T) -T_se = T_sd / sqrt(length(D_mean)) -T_95CL <- qt(0.95, length(D_mean)) * T_se -devtools::check() -devtools::check_rhub() -require(rhub) -validate_email() -devtools::check_rhub() -devtools::release() -devtools::check_win_devel() -require(devtools) -?release_checks -release_checks() -devtools::release() -devtools::document() -install.packages("devtools") -install.packages("tidyverse") -install.packages("rtop") -install.packages("zoo") -install.packages("ggpubr") -install.packages(c("ggplot2", "tidyr", "scales", "dplyr", "magrittr")) -install.packages(c("knitr", "rmarkdown")) -devtools::document() -devtools:install() -devtools::install() -devtools::build() -devtools::check() -devtools::release() -require(ShellChron) -devtools::document() -devtools::build_manual() -devtools::build_manual() -devtools::document() -devtools::build_manual() -devtools::build_manual() -devtools::build_manual() -?build_manual -devtools::build_manual() -devtools::document() -devtools::build_manual() -devtools::document() -devtools::install() -?sinreg -require(ShellChron) -?sinreg -?run_model -require(ShellChron) -?run_model -?sinreg -devtools::build_manual() -devtools::document() -devtools::build_manual() -devtools::document() -devtools::build_manual() -devtools::check_win_devel() -require(devtools) -build_manual() -build_manual -?check -devtools::check(manual = TRUE) -devtools::check(manual = TRUE) -devtools::check() -devtools::document() -devtools::check() -devtools::document() -devtools::build_manual() -devtools::document() -devtools::build_manual() -devtools::check(manual = TRUE) -devtools::build_manual() -devtools::build_manual() -devtools::build_manual() -devtools::document() -devtools::check_man() -devtools::check() -devtools::release() -devtools::build_manual() -devtools::build() -devtools::document() -devtools::build_manual() +# Calculate the temperature equivalent of the calibration offset +mean_Temp <- mean(dat$Temp[which(dat$type == "bivalve" & (dat$Analysis == "this study" | dat$Analysis == "Bernasconi18"))]) +mean_D47 <- 0.0391 * 10 ^ 6 / (mean_Temp + 273.13) ^ 2 + 0.154 +# Statistics of difference with calibration +Calibration_offset_stats <- calibration_offset %>% # Summarize D47 offset statistics +group_by(Calibration) %>% +summarize( +N = n(), # Calculate the number of modelled values, excluding NA's +D47_offset_Average = binmeans(x = D47_offset, x_sd = D47_SD, output = "mean"), +D47_offset_SD = binmeans(x = D47_offset, x_sd = D47_SD, output = "SD"), +D47_offset_SE = D47_offset_SD / sqrt(N - 1), # Calculate the standard error +D47_offset_CL95 = qt(0.95, N - 1) * D47_offset_SE # Calculate the 95% confidence level +) %>% +ungroup() +Calibration_offset_stats$Temp_offset_Average <- sqrt(0.0391 * 10 ^ 6 / (mean_D47 - Calibration_offset_stats$D47_offset_Average - 0.154)) - 273.15 - mean_Temp +Calibration_offset_stats$Temp_offset_SD <- sqrt(0.0391 * 10 ^ 6 / (mean_D47 - Calibration_offset_stats$D47_offset_Average - Calibration_offset_stats$D47_offset_SD - 0.154)) - 273.15 - mean_Temp - Calibration_offset_stats$Temp_offset_Average +Calibration_offset_stats$Temp_offset_SE <- sqrt(0.0391 * 10 ^ 6 / (mean_D47 - Calibration_offset_stats$D47_offset_Average - Calibration_offset_stats$D47_offset_SE - 0.154)) - 273.15 - mean_Temp - Calibration_offset_stats$Temp_offset_Average +Calibration_offset_stats$Temp_offset_CL95 <- sqrt(0.0391 * 10 ^ 6 / (mean_D47 - Calibration_offset_stats$D47_offset_Average - Calibration_offset_stats$D47_offset_CL95 - 0.154)) - 273.15 - mean_Temp - Calibration_offset_stats$Temp_offset_Average +View(calibration_offset) +View(dat) +mean(dat$D47res_Anderson) +sd(dat$D47res_Anderson) +sd(dat$D47res_Anderson)/sqrt(length(dat$D47res_Anderson) - 1) +qt(0.05, length(dat$D47res_Anderson)) +qt(0.95, length(dat$D47res_Anderson)) +qt(0.95, length(dat$D47res_Anderson)) * sd(dat$D47res_Anderson)/sqrt(length(dat$D47res_Anderson) - 1) +mean(dat$D47res_Anderson[-which(dat$Analysis == "Muller17")]) +qt(0.95, length(dat$D47res_Anderson[-which(dat$Analysis == "Muller17")])) * sd(dat$D47res_Anderson[-which(dat$Analysis == "Muller17")])/sqrt(length(dat$D47res_Anderson[-which(dat$Analysis == "Muller17")]) - 1) +View(Calibration_offset_stats) +load("C:/Users/niels/Dropbox/Research/postdoc/UNBIAS/growth experiments/LAICPMS_Sr_spiking/Batch2/LA data/LA_combined_batch2.Rdata") +View(LA_combined) +Profile_plot_Sr_offset_all <- ggplot(LA_combined) + +geom_point(data = subset(LA_combined, !(Specimen_id %in% c(5, 7, 9))), +aes(Depth, +SrCa + as.numeric(Specimen_id) - 1, +col = Species), +alpha = 0.1, +size = 0.1) + +scale_y_continuous("[Sr]/[Ca] (mmol/mol)", +breaks = seq(0, 20, 1), +labels = seq(0, 20, 1), +limits = c(0, 20)) + +scale_x_continuous("Distance from ventral margin [mm]") + +ggtitle("Offset (+ 1) Sr/Ca curves") + +theme_bw() +require(tidyverse) +require(gridExtra) +Profile_plot_Sr_offset_all <- ggplot(LA_combined) + +geom_point(data = subset(LA_combined, !(Specimen_id %in% c(5, 7, 9))), +aes(Depth, +SrCa + as.numeric(Specimen_id) - 1, +col = Species), +alpha = 0.1, +size = 0.1) + +scale_y_continuous("[Sr]/[Ca] (mmol/mol)", +breaks = seq(0, 20, 1), +labels = seq(0, 20, 1), +limits = c(0, 20)) + +scale_x_continuous("Distance from ventral margin [mm]") + +ggtitle("Offset (+ 1) Sr/Ca curves") + +theme_bw() +x11(); plot(Profile_plot_Sr_offset_all) +Profile_plot_Sr_offset_all <- ggplot(LA_combined) + +geom_point(data = subset(LA_combined, !(Specimen_id %in% c(5, 7, 9))), +aes(Depth, +SrCa + as.numeric(Specimen_id) - 1, +col = Species), +alpha = 0.1, +size = 0.1) + +geom_smooth() + +scale_y_continuous("[Sr]/[Ca] (mmol/mol)", +breaks = seq(0, 20, 1), +labels = seq(0, 20, 1), +limits = c(0, 20)) + +scale_x_continuous("Distance from ventral margin [mm]") + +ggtitle("Offset (+ 1) Sr/Ca curves") + +theme_bw() +plot(Profile_plot_Sr_offset_all) +Profile_plot_Sr_offset_all <- ggplot(LA_combined) + +geom_point(data = subset(LA_combined, !(Specimen_id %in% c(5, 7, 9))), +aes(Depth, +SrCa + as.numeric(Specimen_id) - 1, +col = Species), +alpha = 0.1, +size = 0.1) + +geom_smooth(aes(x = Depth, +y = SrCa + as.numeric(Specimen_id) - 1, +col = Species)) + +scale_y_continuous("[Sr]/[Ca] (mmol/mol)", +breaks = seq(0, 20, 1), +labels = seq(0, 20, 1), +limits = c(0, 20)) + +scale_x_continuous("Distance from ventral margin [mm]") + +ggtitle("Offset (+ 1) Sr/Ca curves") + +theme_bw() +plot(Profile_plot_Sr_offset_all) +?geom_smooth +plot(LA_combined$Distance) +plot(LA_combined$Xpos, LA_combined$Ypos) +x11(); plot(LA_combined$Xpos, LA_combined$Ypos) diff --git a/R/Run_model.r b/R/Run_model.r index 1ee46d1..c7b50b5 100644 --- a/R/Run_model.r +++ b/R/Run_model.r @@ -76,7 +76,7 @@ run_model <- function(dat, # Core function to run the entire model on the data ( 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, # If errors = TRUE, give the number of iterations for Monte Carlo simulation used in error propagation (default = 1000, if MC = 0 no eror propagation is done) - plot = FALSE # Should the progress of model fitting be plotted? + plot = TRUE # Should the progress of model fitting be plotted? ){ d18Oc <- d18Oc_err <- Omod <- NULL # Predefine variables to circumvent global variable binding error @@ -170,49 +170,52 @@ run_model <- function(dat, # Core function to run the entire model on the data ( O_err <- rep(0, dynwindow$y) MC <- 0 } + if(i == 1){ # Estimate starting parameters for first data window + if(sinfit){ # If sinusoidal fitting is enabled + sinlist <- sinreg(Dsam, Osam) # Run sinusoidal regression to find initial parameter values + # Estimate starting parameters from regression results + O_av_start <- sinlist[[1]][1] # Export starting value for annual d18O average + O_amp_start <- sinlist[[1]][2] # Export starting value for d18O amplitude + O_pha_start <- sinlist[[1]][4] %% sinlist[[1]][3] # Estimate position (in depth of the first peak in d18O) + O_per_start <- sinlist[[1]][3] # Export starting value for period in distance domain + }else{ + O_av_start <- mean(Osam) # Estimate starting value for annual d18O average by mean of d18O in record + O_amp_start <- diff(range(Osam)) / 2 # Estimate starting value for d18O amplitude by half the difference between minimum and maximum d18Oc + O_per_start <- diff(range(Dsam)) # Estimate starting period as thickness of isolated year + O_pha_start <- 0.25 * O_per_start # Set starting phase to one quarter of a cycle if it cannot be estimated from sinusoidal regression + } - if(sinfit){ # If sinusoidal fitting is enabled - sinlist <- sinreg(Dsam, Osam) # Run sinusoidal regression to find initial parameter values - # Estimate starting parameters from regression results - O_av_start <- sinlist[[1]][1] # Export starting value for annual d18O average - O_amp_start <- sinlist[[1]][2] # Export starting value for d18O amplitude - O_pha_start <- sinlist[[1]][4] %% sinlist[[1]][3] # Estimate position (in depth of the first peak in d18O) - O_per_start <- sinlist[[1]][3] # Export starting value for period in distance domain - }else{ - O_av_start <- mean(Osam) # Estimate starting value for annual d18O average by mean of d18O in record - O_amp_start <- diff(range(Osam)) / 2 # Estimate starting value for d18O amplitude by half the difference between minimum and maximum d18Oc - O_per_start <- diff(range(Dsam)) # Estimate starting period as thickness of isolated year - O_pha_start <- 0.25 * O_per_start # Set starting phase to one quarter of a cycle if it cannot be estimated from sinusoidal regression - } + if(transfer_function == "KimONeil97"){ + T_av_start <- 18.03 * 1000 / (1000 * log((O_av_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 # Estimate mean temperature. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014) + T_amp_start <- 18.03 * 1000 / (1000 * log((O_av_start - O_amp_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 - T_av_start # Estimate temperature amplitude. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014) + }else if(transfer_function == "GrossmanKu86"){ + T_av_start <- 20.6 - 4.34 * (O_av_start - mean(d18Ow) - 0.2) # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999) + T_amp_start <- 20.6 - 4.34 * (O_av_start - O_amp_start - mean(d18Ow) - 0.2) - T_av_start # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999) + }else{ + print("ERROR: Supplied transfer function is not recognized") + } + + O_peak <- O_pha_start + Dsam[1] # Find position of d18O peak in distance domain + T_pha_start <- ((O_pha_start - 0.5 * O_per_start) %% O_per_start) / O_per_start * T_per # Estimate position of first peak in temperature (low in d18O) relative to annual cycle (days) + G_av_start <- O_per_start / G_per # Estimate average growth rate in distance/day + + # Collate starting parameters + par0 <- c( + T_amp = T_amp_start, + T_pha = T_pha_start, + T_av = T_av_start, - if(transfer_function == "KimONeil97"){ - T_av_start <- 18.03 * 1000 / (1000 * log((O_av_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 # Estimate mean temperature. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014) - T_amp_start <- 18.03 * 1000 / (1000 * log((O_av_start - O_amp_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 - T_av_start # Estimate temperature amplitude. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014) - }else if(transfer_function == "GrossmanKu86"){ - T_av_start <- 20.6 - 4.34 * (O_av_start - mean(d18Ow) - 0.2) # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999) - T_amp_start <- 20.6 - 4.34 * (O_av_start - O_amp_start - mean(d18Ow) - 0.2) - T_av_start # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999) + G_amp = G_av_start / 2, # Start by estimating growth rate changes by half the average + G_pha = T_pha_start, # Start by estimating that the peak in growth rate coincides with the peak in temperature + G_av = G_av_start, + G_skw = 50 # Start with a no skew + ) }else{ - print("ERROR: Supplied transfer function is not recognized") + par0 <- par1 # For remaining data windows, the starting parameters are based on the modelled parameters of the previous window } - O_peak <- O_pha_start + Dsam[1] # Find position of d18O peak in distance domain - T_pha_start <- ((O_pha_start - 0.5 * O_per_start) %% O_per_start) / O_per_start * T_per # Estimate position of first peak in temperature (low in d18O) relative to annual cycle (days) - G_av_start <- O_per_start / G_per # Estimate average growth rate in distance/day - years <- 3 # Set default number of years to 3 - # Collate starting parameters - par0 <- c( - T_amp = T_amp_start, - T_pha = T_pha_start, - T_av = T_av_start, - - G_amp = G_av_start / 2, # Start by estimating growth rate changes by half the average - G_pha = T_pha_start, # Start by estimating that the peak in growth rate coincides with the peak in temperature - G_av = G_av_start, - G_skw = 50 # Start with a no skew - ) - invisible(capture.output( # Suppress the details on converging SCEUA sceua_list <- rtop::sceua(growth_model, par0, diff --git a/R/Wrap_function.r b/R/Wrap_function.r index 8751526..f92fa24 100644 --- a/R/Wrap_function.r +++ b/R/Wrap_function.r @@ -96,7 +96,7 @@ wrap_function <- function(path = getwd(), # Wrapping function for the entire mod resultarray <- resultlist[[1]] parmat <- resultlist[[2]] - # STEP 3: Align model results to cumulative timescale + # STEP 3: Align model resultis s to cumulative timescale print("Calculating cumulative day of the year results...") suppressWarnings(resultarray[, , 3] <- cumulative_day(resultarray, TRUE, TRUE, export_path)) # Calculate cumulative day of the year for all model runs and replace matrix in result array diff --git a/R/data_import.r b/R/data_import.r index ee0ddf8..32f5224 100644 --- a/R/data_import.r +++ b/R/data_import.r @@ -50,6 +50,12 @@ data_import <- function(file_name){ return("ERROR: Input data does not match the default input data format") } + # Check for duplicate depth values + if(TRUE %in% duplicated(dat$D)){ + dat <- dat[-which(duplicated(dat$D) == TRUE), ] # Remove duplicated depth values + print("WARNING: Duplicated depth values were found and removed") + } + # Define sliding window based on indicated year markers YEARMARKER <- which(dat$YEARMARKER == 1) # Read position of yearmarkers in data. yearwindow <- diff(which(dat$YEARMARKER == 1)) # Calculate the number of datapoints in each year between consecutive year markers