Comparison of standard residuals
residual.comparison <- multiplot(CDES.residuals.pri.90, CDES.residuals.sec.bp.90, cols=3)

Std.res.all <- bind_rows(CDES.stdres.pri, CDES.stdres.sec)
write.csv(Std.res.all, file=paste0(session, "_CITstdresiduals.csv"))
residual.comparison
## NULL
standard correcting bulk isotopes values
bulk.stds.acc <- read.csv("bulkisostds.csv")
bulk.stds <- data %>%
filter(D48excess.label!="D48 excess" & Donotuse != TRUE & Type=="standard" | Sample.ID == "NBS18") # make a subset of HG and EG with no analytical problems and no D48 excess
bulk.stds <- left_join(bulk.stds, select(bulk.stds.acc, Sample.ID, d18O.acc, d13C.acc), by = "Sample.ID")
## Warning: Column `Sample.ID` joining factors with different levels, coercing to
## character vector
d13C discrimination correction - First, need to pick which correction scheme you prefer, then update this with the correct data columns
# regression
m.d13C <- lm(bulk.stds$d13C.acc ~ bulk.stds$d13C)
disc.slopeC<-(coef(m.d13C)[[2]])
disc.interC<-(coef(m.d13C)[[1]])
R2 <- summary(m.d13C)$r.squared
disc.plot.C <- ggplot(bulk.stds, aes(x=d13C, y = d13C.acc)) +
stat_smooth(method = "lm") +
geom_point(size = 3, shape = 21) +
theme_bw()
disc.plot.C
## `geom_smooth()` using formula 'y ~ x'

# apply correction to data
data$d13C.disc.offset <- data$d13C - (disc.slopeC * data$d13C + disc.interC) #calculates offset values for all acid corrected values
data$d13C.disc <- data$d13C - data$d13C.disc.offset
18O discrimination correction - First, need to pick which correction scheme you prefer, then update this with the correct data columns
# regression
m.d18O <- lm(bulk.stds$d18O.acc ~ bulk.stds$d18O.VPDB.min)
disc.slopeO<-(coef(m.d18O)[[2]])
disc.interO<-(coef(m.d18O)[[1]])
R2 <- summary(m.d18O)$r.squared
disc.plot.O <- ggplot(bulk.stds, aes(x=d18O.VPDB.min, y = d18O.acc)) +
stat_smooth(method = "lm") +
geom_point(size = 3, shape = 21) +
theme_bw()
disc.plot.O
## `geom_smooth()` using formula 'y ~ x'

# apply correction to data
data$d18O.disc.offset <- data$d18O.VPDB.min - (disc.slopeO * data$d18O.VPDB.min + disc.interO) #calculates offset values for all acid corrected values
data$d18O.disc <- data$d18O.VPDB.min - data$d18O.disc.offset
average reference data (standards, equilibrated gases), and plot D47,90 vs D48, 90 and compare to theoretical dual equilibrium line of Fiebig et al., 2019 (derived from Hill et al., 2014 Keq calculations for D63 and D64 of calcite)
data.avg <- data %>%
group_by(Sample.ID) %>%
filter(D48excess.label != "D48 excess") %>%
summarise(Project = Project[1],
N = n(),
d13C.sd = sqrt(calc_wgt_variance(d13C.disc, d13C.sterr)),
d13C.se = d13C.sd/sqrt(N),
d13C.avg = calc_wgt_means(d13C.disc, d13C.sterr),
d18O.sd = sqrt(calc_wgt_variance(d18O.disc, d18O.sterr)),
d18O.se = d18O.sd/sqrt(N),
d18O.avg = calc_wgt_means(d18O.disc, d18O.sterr),
D47.sd = sqrt(calc_wgt_variance(CDES.noAD.pri, D47.HG.prop.err)),
D47.se = D47.sd/sqrt(N),
D47.2se = 2*D47.se,
D47.90.noAD.avg = calc_wgt_means(CDES.noAD.pri, D47.HG.prop.err),
D47.90.avg = calc_wgt_means(CDES.D47stdcorr90.pri, D47.HG.prop.err),
D47.25D.avg = calc_wgt_means(CDES.D47stdcorrD.pri, D47.HG.prop.err),
D47.25P.avg = calc_wgt_means(CDES.D47stdcorrP.pri, D47.HG.prop.err),
D48.sd = sqrt(calc_wgt_variance(CDES.48noAD.pri, D48.sterr)),
D48.se = D48.sd/sqrt(N),
D48.2se = 2*D48.se,
D48.90.avg = calc_wgt_means(CDES.48noAD.pri, D48.sterr),
D48.sd.all = sqrt(calc_wgt_variance(CDES.48noAD.pri.all, D48.sterr)),
D48.se.all = D48.sd.all/sqrt(N),
D48.2se.all = 2*D48.se.all,
D48.90.avg.all = calc_wgt_means(CDES.48noAD.pri.all, D48.sterr),
T.B = convert_ARF.D47_to_ARF.Bonifacie.temp(D47.90.avg),
T.B.se = calc_D47.Temp_se(T.B, D47.90.avg, D47.se),
T.G = convert_ARF.D47_to_ARF.Ghosh.temp(D47.25D.avg),
T.G.se = calc_D47.Temp_se(T.B, D47.25D.avg, D47.se),
T.P = convert_ARF.D47_to_ARF.Peterson.temp(D47.25P.avg),
T.P.se = calc_D47.Temp_se(T.P, D47.25P.avg, D47.se),
d18Ow.B = calc_d18Ow(d18O.avg, T.B),
d18Ow.B.se = calc_d18Owse(d18O.avg, d18O.se, T.B, T.B.se),
d18Ow.G = calc_d18Ow(d18O.avg, T.G),
d18Ow.G.se = calc_d18Owse(d18O.avg, d18O.se, T.G, T.G.se),
d18Ow.P = calc_d18Ow(d18O.avg, T.P),
d18Ow.P.se = calc_d18Owse(d18O.avg, d18O.se, T.P, T.P.se))
## Warning in sqrt(calc_wgt_variance(d13C.disc, d13C.sterr)): NaNs produced
## Warning in sqrt(calc_wgt_variance(d18O.disc, d18O.sterr)): NaNs produced
## Warning in sqrt(calc_wgt_variance(CDES.noAD.pri, D47.HG.prop.err)): NaNs
## produced
pD47vD48 <- ggplot(data = data.frame(x=c(0.1, 0.7)), aes(x)) +
stat_function(fun = function(x) {0.1829*(x^2) + 0.2138*(x) - 0.0139}, colour = "light gray") +
stat_function(fun = function(x) {0.3821*(x^2) + 0.0044*(x) + 0.1198}) +
geom_errorbar(data = data.avg, aes(x=D47.90.noAD.avg, ymin=D48.90.avg - D48.2se, ymax=D48.90.avg + D48.2se, colour = Sample.ID)) +
geom_errorbarh(data = data.avg, aes(x=D47.90.noAD.avg, y=D48.90.avg, xmin=D47.90.noAD.avg - D47.2se, xmax = D47.90.noAD.avg + D47.2se , colour = Sample.ID)) +
geom_point(data = data.avg, aes(x=D47.90.noAD.avg, y=D48.90.avg, fill = Sample.ID, colour = NULL), size = 3, shape = 21) +
#scale_shape_manual(values = c(21,22,23,24,25)) +
theme_bw()
## Warning: Ignoring unknown aesthetics: x
pD47vD48.all <- ggplot(data = data.frame(x=c(0.1, 0.7)), aes(x)) +
stat_function(fun = function(x) {0.1829*(x^2) + 0.2138*(x) - 0.0139}, colour = "light gray") +
stat_function(fun = function(x) {0.3821*(x^2) + 0.0044*(x) + 0.1198}) +
geom_errorbar(data = data.avg, aes(x=D47.90.noAD.avg, ymin=D48.90.avg.all - D48.2se.all, ymax=D48.90.avg.all + D48.2se.all, colour = Sample.ID)) +
geom_errorbarh(data = data.avg, aes(x=D47.90.noAD.avg, y=D48.90.avg.all, xmin=D47.90.noAD.avg - D47.2se, xmax = D47.90.noAD.avg + D47.2se , colour = Sample.ID)) +
geom_point(data = data.avg, aes(x=D47.90.noAD.avg, y=D48.90.avg.all, fill = Sample.ID, colour = NULL), size = 3, shape = 21) +
#scale_shape_manual(values = c(21,22,23,24,25)) +
theme_bw()
## Warning: Ignoring unknown aesthetics: x
pD47vD48
## Warning: Removed 2 rows containing missing values (geom_errorbarh).

pD47vD48.all
## Warning: Removed 2 rows containing missing values (geom_errorbarh).

data.avg
## # A tibble: 24 x 36
## Sample.ID Project N d13C.sd d13C.se d13C.avg d18O.sd d18O.se
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3-021core WP2011 4 0.0757 0.0378 -5.96 0.148 0.0738
## 2 3-021rim WP2011 2 0.135 0.0956 -5.60 0.0170 0.0120
## 3 3A-097 WP2011 3 0.0273 0.0158 -6.75 0.160 0.0924
## 4 3B-021 WP2011 3 0.0394 0.0227 -3.68 0.401 0.231
## 5 3E-001cl… WP2011 3 0.0516 0.0298 -5.92 0.297 0.171
## 6 3F-019 WP2011 1 Inf Inf -6.95 Inf Inf
## 7 3H-014 WP2011 4 0.0801 0.0400 -5.94 0.168 0.0839
## 8 4-038 WP2011 3 0.0510 0.0294 -7.23 0.136 0.0784
## 9 6-003 WP2011 2 0.0182 0.0129 -6.85 0.00613 0.00434
## 10 6-042 WP2011 1 NaN NaN -6.27 NaN NaN
## # … with 14 more rows, and 28 more variables: d18O.avg <dbl>, D47.sd <dbl>,
## # D47.se <dbl>, D47.2se <dbl>, D47.90.noAD.avg <dbl>, D47.90.avg <dbl>,
## # D47.25D.avg <dbl>, D47.25P.avg <dbl>, D48.sd <dbl>, D48.se <dbl>,
## # D48.2se <dbl>, D48.90.avg <dbl>, D48.sd.all <dbl>, D48.se.all <dbl>,
## # D48.2se.all <dbl>, D48.90.avg.all <dbl>, T.B <dbl>, T.B.se <dbl>,
## # T.G <dbl>, T.G.se <dbl>, T.P <dbl>, T.P.se <dbl>, d18Ow.B <dbl>,
## # d18Ow.B.se <dbl>, d18Ow.G <dbl>, d18Ow.G.se <dbl>, d18Ow.P <dbl>,
## # d18Ow.P.se <dbl>
same plot of D47 v D48, but plotting all the individual values, instead of the averages
pD47vD48.all <- ggplot(data = data.frame(x=c(0.1, 0.7)), aes(x)) +
stat_function(fun = function(x) {0.1829*(x^2) + 0.2138*(x) - 0.0139}, colour = "light gray") + #D48 of C02, from Fiebig eq. for Wang et al., 2004
stat_function(fun = function(x) {0.3821*(x^2) + 0.0044*(x) + 0.1198}) + #from data in Hill et al 2014; regression from values in Table B2
geom_point(data = data, aes(x=CDES.noAD.pri, y=CDES.48noAD.pri, fill = Sample.ID, colour = NULL, shape = D48excess.label), size = 2) +
scale_shape_manual(values = c(21,23,24)) +
theme_bw()
ggplotly(pD47vD48.all)
write.csv(data, file=paste0(session, "_brandcorr_reprocessed.csv"))
Can subdivide the dataset into individual projects, each with all of the “all” data, to use for supplemental files for different projects. To use, you’ll need to manually update the separate project dataframe names and filtering criteria, and then the worksheet labels accordingly
Supp.table <- data %>%
filter(Project == "Marina" | Project == "all") %>%
select(Date, Sample.ID, Type, Donotuse, runinfo, spec.num, mass, num.acq, d45,d45.stdev,d45.sterr,d46,d46.stdev,d46.sterr,d47,d47.stdev,d47.sterr,d48,d48.stdev,d48.sterr,d49,d49.stdev,D47,D47.stdev,D47.sterr,D48,D48.stdev,D48.sterr,d13C,d13C.stdev,d13C.sterr,d18O.VPDB.min,d18O.stdev,d18O.sterr,D48.HGcorr,D48excess.label, HGcorr, D47.HG.prop.err,CDES.noAD.pri,CDES.D47stdcorr90.pri, CDES.ADP.pri,CDES.D47stdcorrP.pri,CDES.ghosh.temp.pri.D,CDES.peterson.temp.stdcorr.pri,CDES.Bonifacie.temp.pri)
samples <- data %>%
filter(Project == "Marina") %>%
select(Date, Sample.ID, Type, Donotuse, runinfo, spec.num, mass, num.acq, d45,d45.stdev,d45.sterr,d46,d46.stdev,d46.sterr,d47,d47.stdev,d47.sterr,d48,d48.stdev,d48.sterr,d49,d49.stdev,D47,D47.stdev,D47.sterr,D48,D48.stdev,D48.sterr,d13C,d13C.stdev,d13C.sterr,d18O.VPDB.min,d18O.stdev,d18O.sterr,D48.HGcorr,D48excess.label, HGcorr, D47.HG.prop.err,CDES.noAD.pri,CDES.D47stdcorr90.pri, CDES.ADP.pri,CDES.D47stdcorrP.pri,CDES.ghosh.temp.pri.D,CDES.peterson.temp.stdcorr.pri,CDES.Bonifacie.temp.pri)
Averages <- data.avg %>%
filter(Project == "all" | Project == "Marina" | Project == "int.standards")
add_ws_with_data <- function(wb, sheet, data) {
addWorksheet(wb, sheet)
writeData(wb, sheet=sheet, data)
return(wb)
}
wb <- createWorkbook("data")
wb <- add_ws_with_data(wb, "all data", data)
wb <- add_ws_with_data(wb, "all supp info", Supp.table)
wb <- add_ws_with_data(wb, "all replicates", samples)
wb <- add_ws_with_data(wb, "Averages", Averages)
saveWorkbook(wb, paste0(session, "_flatlist_project.xlsx"), overwrite = TRUE)