In EarthSystemDiagnostics/spectraluncertainty:
Alternative names for psem
- SPR(UCE) Spectral Proxy Record UnCErtainty
- SUPR Spectral Uncertainty for Proxy Records
# Detect OS
get_os <- function() {
if (.Platform$OS.type == "windows") {
"win"
} else if (Sys.info()["sysname"] == "Darwin") {
"mac"
} else if (.Platform$OS.type == "unix") {
"unix"
} else {
stop("Unknown OS")
}
}
os <- get_os()
knitr::opts_chunk$set(
echo = TRUE, message = FALSE, warning = FALSE, cache = TRUE,
fig.width = 7, fig.height = 4.5, fig.pos = "H", dpi = 300,
dev = ifelse(os == "win", "pdf", "cairo_pdf"),
autodep = TRUE
)
# detach("package:psem", unload = TRUE)
# detach("package:sedproxy", unload = TRUE)
library(psem)
#library(spectraluncertainty)
library(sedproxy)
library(tidyverse)
library(corit)
# lookup tables and plotting keys
# Formats parameter values for tables
FormatParValues <- function(v, max.digits = 4, max.char = 6) {
v <- as.numeric(v)
v <- as.list(v)
v2 <- v
matches <- c(as.character(pi), as.character(pi / 2),
as.character(pi / 3), as.character(pi / 4),
as.character(2*pi), as.character(exp(1)),
as.character(2*pi/3))
replacements <- c("pi", "pi/2", "pi/3", "pi/4", "2*pi", "e", "2pi/3")
hc <- lapply(v, function(x) if(as.character(x) %in% matches){replacements[match(x, matches)]}else NA)
rat <- lapply(v2, function(x) as.character(MASS::as.fractions(x)))
sig <- lapply(v2, function(x) signif(x, max.digits))
nc.hc <- lapply(hc, nchar)
nc.rat <- lapply(rat, nchar)
nc.sig <- lapply(sig, nchar)
nc.hc[is.na(nc.hc)] <- Inf
nc.rat[is.na(nc.rat)] <- Inf
nc.sig[is.na(nc.sig)] <- Inf
v4 <- lapply(1:length(v), function(i) {
if (is.na(hc[[i]]) == FALSE & nc.hc[[i]] <= max.char) {
hc[[i]]
} else if (nc.rat[[i]] < max.char) {
rat[[i]]
} else {
sig[[i]]
}
})
return(as.character(v4))
}
Error types
set.seed(4)
n = 20
df <- tibble(
Age = seq(1, 10, length.out = n),
Climate = PaleoSpec::SimPowerlaw(1, n) + SSlogis(Age, 4, 5, 0.1) - 2,
`Independent errors` = Climate + rnorm(n, 0, 1.5),
`Constant errors` = Climate - 1.5,
`Correlated errors` = Climate + (Age - 5) * -0.58
)
df.long <- df %>%
gather(Noise.type, value, -Age, -Climate) %>%
mutate(Noise.type = factor(Noise.type,
levels = (c(
"Independent errors",
"Correlated errors",
"Constant errors"
)),
ordered = TRUE
))
txt.lbls <- tibble(
Age = -Inf, Climate = Inf, Noise.type = factor(levels(df.long$Noise.type),
ordered = T),
label = letters[1:3]
)
p.error.wide <- df.long %>%
ggplot(aes(x = Age, y = Climate)) +
geom_line(aes(linetype = "Climate")) +
geom_line(aes(y = value, colour = Noise.type, linetype = "Proxy")) +
geom_segment(aes(x = Age, xend = Age, y = Climate, yend = value, colour = Noise.type),
arrow = arrow(length = unit(0.03, "npc"), type = "closed")
) +
facet_wrap(~Noise.type, ncol = 3) +
theme_bw() +
theme(legend.position = "left", panel.grid.minor = element_blank()) +
scale_colour_manual(
values = wesanderson::wes_palette("GrandBudapest1", 4)[2:4],
guide = FALSE
) +
scale_linetype("", guide = FALSE) +
labs(x = "Age [ka]") +
geom_text(data = txt.lbls, aes(label = label), hjust = -1, vjust = 1.5,
fontface = 4)
p.error.wide
Simplified method description
spec.pars.smpl <- GetSpecPars("Mg_Ca",
T = 1e04 + 100, delta_t = 100,
tau_r = 100,
tau_b = 500, tau_s = 50, N = 1
)
# Fourier transform of bioturbation and sediment slice thickness filter
TF.f_bs <- function(nu, tau_b, tau_s) {
i <- 0 + 1i
sinc(nu * tau_s) * (1 + i * 2 * pi * nu * tau_b)^-1 * exp(i * 2 * pi * nu * tau_b)
}
# The redistributed stochastic climate variance
smpl.spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.smpl)
redist.stoch.clim.var <- smpl.spec.obj$proxy.error.spec$Aliasing.stochastic[smpl.spec.obj$proxy.error.spec$nu<0]
dat <- tibble(
nu = GetNu(spec.pars.smpl$T, spec.pars.smpl$delta_t),
Climate.spectrum = 0.1 * nu^-1
) %>%
filter(nu > 0) %>%
mutate(
# the transfer fucntion
tf.f_bs = TF.f_bs(nu, spec.pars.smpl$tau_b, spec.pars.smpl$tau_s),
# squared absolute value of the transfer function
sq.abs.tf.f_bs = abs(tf.f_bs)^2,
Filtered.climate.spectrum = Climate.spectrum * sq.abs.tf.f_bs,
# transfer function for the reference climate, sinc function
tf.ref = sinc(nu * spec.pars.smpl$tau_r),
sq.abs.tf.ref = abs(tf.ref)^2,
Reference.climate.spectrum = sq.abs.tf.ref * Climate.spectrum,
Error.spectrum = abs(tf.f_bs - tf.ref)^2 * Climate.spectrum,
Redist.stoch.clim.var = redist.stoch.clim.var
)
dat.long <- dat %>%
gather(Component, Spec, -nu) %>%
mutate(
line.type = ifelse(Component == "Error.spectrum", "Error", "Signal"),
Spec = Re(Spec)
)
clr.pal <- c((RColorBrewer::brewer.pal(9, "Set1"))[c(1,4,2,7, 9, 9)])
names(clr.pal) <- c(
"Climate.spectrum", "Reference.climate.spectrum", "Filtered.climate.spectrum",
"Error.spectrum", "Redist.as.white.noise", "Redist.stoch.clim.var"
)
clr.lbls <- c("Climate spectrum", "Reference climate\nspectrum",
"Climate spectrum\nafter bioturbation", "Spectrum of\nbioturbation error",
"Redistributed\nas white noise", "Redistributed\nclimate variance")
names(clr.lbls) <- names(clr.pal)
ltyp.pal <- c(1,1,1,2,1,1)
names(ltyp.pal) <- names(clr.pal)
p.a <- dat.long %>%
filter(Component %in% c(
"Climate.spectrum", "Error.spectrum",
"Filtered.climate.spectrum",
"Reference.climate.spectrum"
)) %>%
ggplot(aes(x = nu, y = Spec, colour = Component, linetype = Component)) +
geom_line() +
scale_y_continuous("Power spectral density", trans = "log10") +
scale_x_continuous("Frequency [1/years]", trans = "log10") +
scale_linetype_manual("", values = ltyp.pal, labels = clr.lbls, limits = names(clr.pal)[1:4]) +
scale_colour_manual("", values = clr.pal, labels = clr.lbls, limits = names(clr.pal)[1:4]) +
theme_bw() +
theme(panel.grid.minor = element_blank(), legend.key.height = unit(2, "line"))
p.b <- dat.long %>%
filter(Component %in% c(
"Bioturbation",
"Reference.climate.spectrum",
"tf.ref", "tf.f_bs", "sq.abs.tf.f_bs", "sq.abs.tf.ref",
"Error.spectrum"
) == FALSE) %>%
ggplot(aes(x = nu, y = Spec, colour = Component)) +
geom_line() +
geom_ribbon(
data = dat,
aes(
x = nu, ymin = Filtered.climate.spectrum,
ymax = Climate.spectrum, fill = "Redist.as.white.noise"
),
alpha = 0.5,
inherit.aes = FALSE
) +
scale_y_continuous("Power spectral density", trans = "log10") +
scale_x_continuous("Frequency [1/years]", trans = "log10") +
scale_fill_manual("", values = clr.pal, labels = clr.lbls) +
scale_colour_manual("", values = clr.pal, labels = clr.lbls) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key.height = unit(2, "line")) +
guides(colour = guide_legend(order = 1))
p.a.b <- egg::ggarrange(p.a, p.b, ncol = 1, draw = FALSE, labels = c("a", "b"),
label.args = list(gp = grid::gpar(font = 1, cex = 1.2)))
print(p.a.b)
Power spectrum of proxy error
# PSD Climate
example.lat <- 20
clim.spec.ex1 <- ModelSpectrum(
#freq = GetNu(100000, 100),
freq = NULL,
latitude = example.lat,
variable = "temperature", beta = 1
)
p.clim.spec.ex1 <- PlotModelSpectrum(clim.spec.ex1)
p.clim.spec.ex1
# Proxy seasonality
# climproxyrecords::shakun.metadata %>%
# filter(Core == "ODP 658C")
#
# ecusdata::shak.fraile.seasonality %>%
# filter(ID.no == "N31") %>%
# gather(Month, value, starts_with("v"))
seasonal.amp <- AmpFromLocation(
longitude = -18,
latitude = example.lat,
proxy.type = "degC",
depth.upr = 0, depth.lwr = -120
)
orbital.pars <- RelativeAmplitudeModulation(
latitude = example.lat,
maxTimeKYear = 23,
minTimeKYear = 1,
bPlot = FALSE
)
# Look for good example location
# ssa <- sapply(-90:90, function(l) RelativeAmplitudeModulation(latitude = l,
# maxTimeKYear = 23,
# minTimeKYear = 1,
# bPlot = FALSE)$sig.sq_a)
#
# ssc <- sapply(-90:90, function(l) AmpFromLocation(longitude = -28,
# latitude = l,
# proxy.type = "degC",
# depth.upr = 0, depth.lwr = -200)$sig.sq_c)
# plot(-90:90, ssa, type = "l")
# abline(v = 20)
# plot(-90:90, ssc, type = "l")
# abline(v = 20)
ex.sed.acc.rate <- 10
spec.pars.ex1 <- GetSpecPars("Mg_Ca",
T = 1e04 + 100,
delta_t = 100,
tau_r = 100,
sig.sq_a = orbital.pars$sig.sq_a,
sig.sq_c = round(seasonal.amp$sig.sq_c, 3),
tau_b = 1000 * 10 / ex.sed.acc.rate,
tau_s = 1000 * 1 / ex.sed.acc.rate,
N = 30,
tau_p = 7/12,
phi_c = 0, delta_phi_c = 2 * pi / 3,
phi_a = pi / 2,
sigma.cal = 0.25,
sigma.meas = 0.25,
sigma.ind = 1,
clim.spec.fun = "ModelSpectrum",
clim.spec.fun.args =
list(latitude = example.lat, beta = 1)
)
spec.obj.ex1 <- do.call(ProxyErrorSpectrum, spec.pars.ex1)
p.spec.1 <- spec.obj.ex1 %>%
PlotSpecError(.)
p.spec.1 +
coord_cartesian(ylim = c(2e-1, 1e04)) #+ geom_vline(xintercept = 1/c(23e03))
# FormatParValues(c(2*pi, 1/3, 0.1, pi, pi/2, exp(1), 2.12457896, 2*pi/3))
par.df <- as.data.frame(spec.pars.ex1) %>%
gather(parameter, value) %>%
left_join(parameter.symbol, .) %>%
left_join(., parameter.description) %>%
left_join(., parameter.source) %>%
filter(complete.cases(latex)) %>%
select(latex, value, description, source) %>%
mutate(value = FormatParValues(value, max.digits = 2, max.char = 6)) %>%
tbl_df()
names(par.df) <- c("Parameter", "Value", "Description", "Source")
# par.df %>%
# knitr::kable(.)
pander::pandoc.table(par.df)
# TEX for direct pasting into ms.
par.df$Description <- gsub("pi", "\\(\\pi\\)", x = par.df$Description, fixed = TRUE)
par.df$Source <- gsub("pi", "\\(\\pi\\)", x = par.df$Source, fixed = TRUE)
par.df$Value <- gsub("pi", "\\(\\pi\\)", x = par.df$Value, fixed = TRUE)
cat(apply(par.df[, 1:4], 1,
function(x) paste0(paste0(x, collapse = " & "), " \\\\")))
Timescale dependent error
var.obj.ex1 <- IntegrateErrorSpectra(spec.obj.ex1)
ev.a <- var.obj.ex1 %>%
PlotTSDVariance(., include.constant.errors = TRUE) +
theme(legend.position = "none")
ev.b <- var.obj.ex1 %>%
PlotTSDVariance(.)
ev.a.b <- egg::ggarrange(
ev.a, ev.b, ncol = 2, draw = FALSE,
labels = c("a", "b"),
label.args = list(gp = grid::gpar(font = 1, cex = 1.2)))
ev.a.b.2 <- grid::forceGrob(ev.a.b)
ev.a.b.2
Err diff 2 timeslices example spectrum
tslice.width.1 <- 1000
tslice.width.2 <- 1000
d.ts <- spec.pars.ex1$T - (tslice.width.1/2 + tslice.width.2/2)
err.tslice.diff <- ErrDiff2TimeSlices(spec.obj.ex1,
tau_1 = tslice.width.1,
tau_2 = tslice.width.2,
delta_ts = d.ts)
err.tslice.diff <- err.tslice.diff %>%
filter(component == "Total.error") %>%
select(-cov, -sigma.diff, -var.tau.12) %>%
mutate(naive.var.diff = var.tau.1 + var.tau.2)
err.tslice.diff %>%
rename(`Error\ntime-slice a` = var.tau.1,
`Error\ntime-slice b` = var.tau.2,
`Error\ntime-slice\ndifference` = var.diff,
`Naive error\ntime-slice\ndifference` = naive.var.diff) %>%
gather(s, var, -component) %>%
mutate(sigma = sqrt(var),
s = factor(s, ordered = TRUE, levels = s)) %>%
ggplot(aes(x = s, y = sigma)) +
#geom_point() +
geom_col(alpha = 0.75, width = 0.75) +
theme_bw() +
#labs(y = expression("Uncertainty" ~""%+-%""~"1"*sigma~" [°C]"), x = "") +
labs(y = expression("Error - 1"*sigma~" [°C]"), x = "")
GeoB_10054_4
# Core GeoB 10054-4
core.info <- replicatecoredata::core.info %>%
filter(Core == "GeoB 10054-4")
world <- map_data("world")
# world$long[world$long > 180] <- (360 - world$long[world$long > 180]) * -1
coremap <- ggplot(world, aes(x = long, y = lat, group = group)) +
geom_polygon(colour = "Darkgrey", fill = "Darkgrey", size = 0.5, show.legend = F) +
coord_map("gall", 110,
xlim = c(50, 180),
ylim = c(-50, 50)
) +
geom_point(
data = core.info, aes(Longitude, Latitude),
inherit.aes = FALSE, show.legend = T, colour = "Red",
size = 2, shape = 8
) +
theme_bw() +
theme(legend.position = "top") +
scale_x_continuous("Longitude", labels = function(x) paste0(x, "°")) +
scale_y_continuous("Latitude", labels = function(x) paste0(x, "°"))
coremap
# PSD Climate
clim.spec <- ModelSpectrum(
freq = NULL, latitude = core.info$Latitude,
variable = "d18O", beta = 1
)
p.clim.spec <- PlotModelSpectrum(clim.spec)
p.clim.spec
breit.lons <- unique(breitkreuz.amp$longitude)
nearest.lon <- breit.lons[which.min(abs(core.info$Longitude - breit.lons))]
breit.lats <- unique(breitkreuz.amp$latitude)
nearest.lat <- breit.lats[which.min(abs(core.info$Latitude - breit.lats))]
seasonal.amp <- AmpFromLocation(
longitude = core.info$Longitude,
latitude = core.info$Latitude,
proxy.type = "d18O",
depth.upr = 0, depth.lwr = -50
)
core.dat <- replicatecoredata::GeoB.vars %>%
filter(
variable == "d18O",
age.14C.cal <= 1e04,
core %in% c("GeoB 10054-4"),
taxon == "G. ruber",
n.forams.grp %in% c("30-30", "5-5")
)
record.names <- tibble(
n.forams.grp = c("30-30", "5-5"),
record = c("Rep2", "Rep1")
)
core.dat <- left_join(core.dat, record.names)
timescale.pars <- core.dat %>%
group_by(n.forams.grp) %>%
summarise(
T = diff(range(age.14C.cal)),
mu_delta_t = round(T / n()),
N = mean(n.forams)
) %>%
left_join(., record.names)
orbital.pars <- RelativeAmplitudeModulation(
latitude = core.info$Latitude,
maxTimeKYear = 23,
minTimeKYear = 1,
bPlot = FALSE
)
approx.s <- round(1000 * diff(range(core.dat$depth_cm)) /
diff(range(core.dat$age.14C.cal)))
## estimate of tau_b from replicated 14C measurements
tau_b.10054_4 <- 720
tau_b.10054_4.se <- ecustools::SDofSD(tau_b.10054_4, 10)
1 / 1000 * c(tau_b.10054_4 - tau_b.10054_4.se, tau_b.10054_4, tau_b.10054_4 + tau_b.10054_4.se) * 20
spec.pars.lst <- timescale.pars %>%
plyr::dlply(., c("record"), function(x) {
GetSpecPars(
proxy.type = "d18O",
delta_t = round(x$mu_delta_t),
tau_r = x$mu_delta_t,
T = x$T,
# Bioturbation
tau_b = tau_b.10054_4,
tau_s = 1000 * 1 / approx.s,
N = x$N,
sig.sq_c = seasonal.amp$sig.sq_c,
sigma.meas = 0.1,
sigma.ind = 0.15,
tau_p = 1,
phi_c = 0,
sig.sq_a = orbital.pars$sig.sq_a,
clim.spec.fun = "ModelSpectrum",
clim.spec.fun.args = list(
latitude = core.info$Latitude, beta = 1,
variable = "d18O"
)
)
})
spec.obj.lst <- plyr::llply(spec.pars.lst, function(x)
do.call(ProxyErrorSpectrum, x))
reps.par.tab <- t(plyr::ldply(spec.pars.lst, function(x) {
y <- as.data.frame(x)
y
}))
colnames(reps.par.tab) <- reps.par.tab[1, ]
geob10054.pars <- reps.par.tab %>%
data.frame(., stringsAsFactors = FALSE) %>%
# gather(parameter, value) %>%
tibble::rownames_to_column("parameter") %>%
left_join(., parameter.symbol) %>%
filter(complete.cases(latex)) %>%
select(latex, Rep1, Rep2) %>%
mutate(
Rep1 = FormatParValues(Rep1),
Rep2 = FormatParValues(Rep2)
) %>%
# mutate(value = ifelse(value == as.character(pi), "pi", value)) %>%
tbl_df()
pander::pandoc.table(geob10054.pars)
spec.obj.df <- plyr::ldply(spec.obj.lst, function(x) x[[2]], .id = "record")
plot.symbol <- parameter.symbol %>%
rename(varname = parameter) %>%
add_row(varname = "record", label = "R")
p <- spec.obj.df %>%
filter(nu > 0) %>%
gather(component, spec, -nu, -record) %>%
left_join(., timescale.pars) %>%
filter(spec > 1e-20) %>%
mutate(component = psem:::OrderStages(component)) %>%
ggplot(., aes(x = nu, y = spec, colour = component, linetype = component)) +
geom_line() +
scale_x_continuous(trans = "log10") +
scale_y_continuous(trans = "log10") +
scale_color_manual("", values = spec.colrs) +
scale_linetype_manual("", values = spec.linetypes) +
theme_bw() +
#theme(strip.background = element_blank(), strip.text.y = element_blank()) +
annotation_logticks(sides = "lb") +
labs(x = "Frequency [1/years]", y = "Power Spectral Density") +
facet_wrap(~record + N + mu_delta_t, labeller = function(labs) {
sedproxydocs::label_symbols(labs, label.key = plot.symbol)})
p
var.obj.lst <- plyr::llply(spec.obj.lst, IntegrateErrorSpectra)
err.obj.lst <- plyr::llply(var.obj.lst, function(x) {
GetProxyError(x, timescale = x[["spec.pars"]]$delta_t)
})
# Smooth higher res record to res of lower res record
var.obj.1b <- IntegrateErrorSpectra(spec.obj.lst[["Rep1"]],
tau_smooth = spec.pars.lst[["Rep2"]]$delta_t
)
err.obj.1b <- GetProxyError(var.obj.1b, timescale = spec.pars.lst[["Rep2"]]$delta_t)
err.obj.1b$record <- "Rep1"
err.obj.1.2 <- plyr::ldply(err.obj.lst, .id = "record")
var.obj.1.2 <- plyr::ldply(var.obj.lst, function(x) x[[2]], .id = "record")
TSD Variance
var.492 <- lapply(spec.obj.lst, function(x) IntegrateErrorSpectra(x, tau_smooth = 492))
err.492 <- lapply(var.492, function(x) GetProxyError(x, timescale = 492))
err.492.df <- plyr::ldply(err.492, .id = "record") %>%
group_by(record, smoothed.resolution) %>%
filter(component %in% c("Total.error", "Reference.climate") == FALSE) %>%
summarise(Total = sum(exl.f.zero^2)) %>%
mutate(smoothed.resolution = as.numeric(smoothed.resolution))
# small df for red dots highlighting total error at tau_smooth = 246
err.df <- err.obj.1.2 %>%
filter(record == "Rep2") %>%
bind_rows(., err.obj.1b) %>%
filter(component == "Total.error") %>%
mutate(x = 246, y = exl.f.zero^2)
subplt.lbls <-
tibble(record = factor(c("Rep1", "Rep2")), x = 10000, y = Inf, label = letters[1:2])
PlotTSDVariance(var.obj.1.2, units = "d18O", include.constant.errors = FALSE) +
facet_wrap(~record) +
#geom_point(data = err.df, aes(x = x, y = y), colour = "Red", inherit.aes = F) +
#geom_point(data = err.492.df, aes(x = smoothed.resolution, y = Total),
# colour = "Red", inherit.aes = F) +
geom_vline(xintercept = 492, linetype = 3) +
theme(panel.spacing = unit(1, "lines")) +
#annotate(geom = "text", x = 10000, y = Inf, label = "a") +
scale_x_continuous(trans = ecustools::reverselog_trans(),
breaks = c(100, 492, 1000, 10000),
labels = c(100, 492, 1000, 10000)/1000) +
labs(y = expression("Error variance [ ‰ d"^18*"O"^2 ~"]"),
x = "Timescale [kyr]") +
## Add subplot labels a, b, etc outside facet plotting areas
geom_text(data = subplt.lbls, aes(x = x, y = y, label = label),
inherit.aes = FALSE, hjust = c(5.5, 3), vjust = 1, size = 5) +
# makes points / text outside plotting area visible
coord_cartesian(clip = "off") +
# compensate for the thicker border due to the above
theme(panel.border = element_rect(size = 1/4), panel.spacing.x = unit(2, "lines"),
strip.background = element_blank(), strip.text = element_blank())
err.obj.df <- plyr::ldply(var.obj.lst, function(x) {
GetProxyError(x,
exclude = "Bioturbation", timescale = x[["spec.pars"]]$delta_t,
format = "short"
)
}, .id = "record") %>%
mutate(sigma.total = Total.inc.error)
core.dat.2 <- left_join(core.dat, err.obj.df) %>%
group_by(depth_cm) %>%
mutate(both.sampling.strat = ifelse(length(unique(n.forams.grp)) == 2,
"Same slice", "Extra slice"
))
ts_with_spec_error <- core.dat.2 %>%
# filter(both.sampling.strat == TRUE) %>%
ggplot(aes(
x = age.14C.cal / 1000, y = value,
fill = record,
colour = record,
group = n.forams.grp
)) +
geom_ribbon(aes(
ymax = value + sigma.total,
ymin = value - sigma.total,
colour = NULL
),
alpha = 0.5
) +
geom_point() +
geom_line() +
scale_y_reverse(expression(delta^"18" * O)) +
scale_x_continuous("Age [ka]", limits = c(0, 10), breaks = seq(0, 10, 2)) +
scale_colour_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) +
scale_fill_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) +
expand_limits(x = 0) +
theme_bw() +
theme(panel.grid.minor = element_blank())
ts_with_spec_error_facets <- ts_with_spec_error +
facet_grid(record ~ .)
#ts_with_spec_error_facets
Look only at samples from the same slices.
# core.dat.2 %>%
# # filter(both.sampling.strat == TRUE) %>%
# ggplot(aes(
# x = age.14C.cal / 1000, y = value,
# group = n.forams.grp,
# fill = record,
# colour = record
# )) +
# geom_ribbon(aes(
# ymax = value + sigma.total,
# ymin = value - sigma.total,
# colour = NULL
# ),
# alpha = 0.5
# ) +
# geom_point() +
# geom_line() +
# scale_y_reverse(expression(delta^"18" * O)) +
# scale_x_continuous("Age [ka]") +
# facet_grid(both.sampling.strat ~ .) +
# expand_limits(x = 0) +
# theme_bw() +
# theme(panel.grid.minor = element_blank())
#
# core.dat.2.coverage <- core.dat.2 %>%
# filter(both.sampling.strat == TRUE) %>%
# group_by(depth_cm) %>%
# mutate(n.samples = n()) %>%
# filter(n.samples == 2) %>%
# summarise(
# value.diff = mean(diff(value)),
# sigma.total = sqrt(sum(sigma.total^2)),
# covered = abs(value.diff) <= sigma.total
# )
#
# core.dat.2.coverage %>%
# ungroup() %>%
# summarise(coverage = sum(covered) / n())
Expected correlation
expected.correlations <- plyr::ldply(spec.obj.lst, ExpectedCorrelation,
.id = "record") %>%
left_join(., record.names) %>%
mutate(n.forams.grp = paste0("Expected ", record, "-",
ifelse(correlation.type == "proxy-proxy", record, "Climate")))
# # Interpolated and average expected correlations
#
# expected.correlations %>%
# group_by(correlation.type, record) %>%
# filter(complete.cases(smoothed.resolution, rho)) %>%
# do({
# smoothed.resolution <- seq(min(.$smoothed.resolution), max(.$smoothed.resolution), by = 1)
# rho <- approx(.$smoothed.resolution, y = .$rho, xout = smoothed.resolution)$y
# data.frame(smoothed.resolution, rho)
# }) %>%
# spread(record, rho) %>%
# filter(complete.cases(Rep1, Rep2)) %>%
# mutate(rho = (Rep1 + Rep2) / 2) %>%
# ggplot(aes(x = smoothed.resolution, y = rho, colour = correlation.type)) +
# geom_line()
#
time.series1 <- plyr::dlply(core.dat.2, "record",
function(x) zoo(x$value, order.by = x$age.14C.cal))
Cor <-
#expected.correlations %>%
# select(smoothed.resolution, record) %>%
# rename(timescale = smoothed.resolution) %>%
# filter(timescale <= 1500) %>%
# distinct() %>%
data.frame(timescale = seq(246,
1500,
by = 1)) %>%
group_by(timescale) %>%
do({
tcor <- corit::CorIrregTimser(
timser1 = time.series1[[1]],
timser2 = time.series1[[2]],
detr = TRUE, # remove linear trend time series
method = "InterpolationMethod",
appliedFilter = "runmean",
fc = 1 / (2*.$timescale), # cut-off frequency
#appliedFilter = "gauss",
#fc = 1 / (.$timescale), # cut-off frequency
dt = .$timescale / 10, # time step for the interpolation
int.method = "linear", # kind of interpolation
filt.output = TRUE)
# number of valid points in smoothed timeseries
# used for SE/CI of rho
data.frame(rho = tcor$cor,
n = sum(is.na(tcor$ft1)==FALSE))
}) %>%
mutate(correlation.type = "Observed Rep1-Rep2\n(with corit)")
rho.CI <- Cor %>%
select(-rho, -correlation.type) %>%
rename(smoothed.resolution = timescale) %>%
# join to add n to expected correlations
left_join(expected.correlations, .) %>%
group_by(record, correlation.type) %>%
do({
rhoCI(.$rho, n = .$n, alpha = 1-0.682689492137)
}) %>%
filter(complete.cases(n))
expected.correlations.ci <- expected.correlations %>%
left_join(., rho.CI) %>%
filter(
smoothed.resolution <= 1500)
expected.correlations.ci.mean <- expected.correlations.ci %>%
group_by(correlation.type, record) %>%
filter(complete.cases(smoothed.resolution, rho)) %>%
do({
smoothed.resolution <- seq(min(.$smoothed.resolution), max(.$smoothed.resolution), by = 1)
rho <- approx(.$smoothed.resolution, y = .$rho, xout = smoothed.resolution)$y
lwr <- approx(.$smoothed.resolution, y = .$lwr, xout = smoothed.resolution)$y
upr <- approx(.$smoothed.resolution, y = .$upr, xout = smoothed.resolution)$y
data.frame(smoothed.resolution, rho, lwr, upr)
}) %>%
group_by(correlation.type, smoothed.resolution) %>%
mutate(n = n()) %>%
filter(n == 2) %>%
summarise_if(is.numeric, mean) %>%
ungroup() %>%
mutate(correlation.type = ifelse(correlation.type == "proxy-proxy",
"Expected Proxy-Proxy",
"Expected Proxy-Climate"))
corr.clrs <- wesanderson::wes_palette("GrandBudapest1", 4)[c(2,4,3)]
names(corr.clrs) <- c(
"Expected Proxy-Proxy",
"Expected Proxy-Climate",
"Observed Rep1-Rep2\n(with corit)"
)
corr.ltyp <- c(1,2,1)
names(corr.ltyp) <- names(corr.clrs)
leg.order <- c("Expected Proxy-Climate",
"Observed Rep1-Rep2\n(with corit)",
"Expected Proxy-Proxy")
expected.correlations.ci.mean %>%
ggplot(aes(
x = smoothed.resolution, y = rho,
colour = correlation.type
)) +
geom_ribbon(aes(
ymax = upr, ymin = lwr,
fill = correlation.type
),
colour = NA, alpha = 0.5, show.legend = FALSE
) +
geom_line(aes()) +
geom_line(data = filter(Cor,
timescale < 1000),
aes(
x = timescale, y = rho
)) +
theme_bw() +
expand_limits(x = c(0, 1500), y = 0) +
scale_color_manual("", values = corr.clrs,
limits = leg.order) +
scale_fill_manual("1 standard error", values = corr.clrs,
limits = leg.order) +
scale_linetype_manual("", values = corr.ltyp,
limits = leg.order)+
scale_x_continuous("Timescale [years]") +
scale_y_continuous("Correlation coefficient") +
#guides(linetype = guide_legend(override.aes = list(linetype = c(2,2,1,1,1)))) +
theme(panel.grid.minor = element_blank(), legend.key.height = unit(1.5, "line")) +
expand_limits(y = c(0, 1))
# #corr.clrs <- RColorBrewer::brewer.pal(5, name = "Set2")
# #corr.clrs <- wesanderson::wes_palette("GrandBudapest1", 3)[c(1,2,1,2,3)]
# corr.clrs <- c(rep(RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)],2), "#5B1A18")
#
# names(corr.clrs) <- c(
# "Expected Rep1-Rep1",
# "Expected Rep2-Rep2",
# "Expected Rep1-Climate",
# "Expected Rep2-Climate",
# "Observed Rep1-Rep2\n(with corit)"
# )
#
# corr.ltyp <- c(1,1,2,2,1)
# names(corr.ltyp) <- c(
# "Expected Rep1-Rep1",
# "Expected Rep2-Rep2",
# "Expected Rep1-Climate",
# "Expected Rep2-Climate",
# "Observed Rep1-Rep2\n(with corit)"
# )
#
# leg.order <- c("Expected Rep1-Climate",
# "Expected Rep2-Climate",
# "Expected Rep1-Rep1",
# "Observed Rep1-Rep2\n(with corit)",
# "Expected Rep2-Rep2")
#
# expected.correlations.ci %>%
# ggplot(aes(
# x = smoothed.resolution, y = rho,
# colour = n.forams.grp
# )) +
# geom_ribbon(aes(
# ymax = upr, ymin = lwr,
# fill = n.forams.grp
# ),
# colour = NA, alpha = 0.5, show.legend = FALSE
# ) +
# geom_line(aes(linetype = n.forams.grp)) +
# geom_line(data = Cor, aes(
# x = timescale, y = rho
# )) +
# # geom_point(data = Cor, aes(
# # x = timescale, y = rho
# # )) +
# theme_bw() +
# expand_limits(x = c(0, 1500), y = 0) +
# scale_color_manual("", values = corr.clrs,
# limits = leg.order) +
# scale_fill_manual("1 standard error", values = corr.clrs,
# limits = leg.order) +
# scale_linetype_manual("", values = corr.ltyp,
# limits = leg.order)+
# scale_x_continuous("Timescale [years]") +
# scale_y_continuous("Correlation coefficient") +
# guides(linetype = guide_legend(override.aes = list(linetype = c(2,2,1,1,1)))) +
# theme(panel.grid = element_blank(), legend.key.height = unit(1.5, "line")) +
# expand_limits(y = c(0, 1))
expected.correlations.ci %>%
select(correlation.type, smoothed.resolution, rho, n, lwr, upr, width) %>%
group_by(correlation.type, smoothed.resolution) %>%
summarise_if(is.numeric, mean)
# SmoothIrregular2 <- function(df, time.var, value.var, tau_smooth) {
# mean.d_t <- floor(mean(diff(df[[time.var]])))
#
# strt.time <- mean.d_t * round((min(df[[time.var]]) - mean.d_t / 2) / mean.d_t)
# end.time <- mean.d_t * round((max(df[[time.var]]) + mean.d_t / 2) / mean.d_t)
#
# age <- seq(strt.time, end.time, by = tau_smooth)
#
# lst <- lapply(age, function(a) {
# df.sub <- subset(df, df[[time.var]] > a - tau_smooth / 2 & df[[time.var]] < a + tau_smooth / 2)
# list(
# mean.value = mean(df.sub[[value.var]]), n = nrow(df.sub),
# mean.d_t = tau_smooth / nrow(df.sub), tau_smooth = tau_smooth
# )
# })
#
# df.2 <- do.call(rbind.data.frame, lst)
# cbind(age, df.2)
# }
core.dat %>%
group_by(record) %>%
summarise(
max.dT = max(diff(age.14C.cal)),
mean.dT = mean(diff(age.14C.cal))
)
tau_smooth <- 2 * spec.pars.lst[["Rep2"]]$delta_t
core.dat.3 <- core.dat.2 %>%
group_by(n.forams.grp, record, n.forams) %>%
do({
# detrend?
# .$value <- residuals(lm(.$value~.$age.14C.cal))
SmoothIrregular(., "age.14C.cal", "value",
tau_smooth = tau_smooth, d_t = tau_smooth
)
})
core.dat.4 <- lapply(spec.pars.lst, function(i) {
core.dat.3 <- core.dat.3 %>%
filter(n.forams == i$N)
d_t <- core.dat.3 %>%
pull(mean.d_t)
pars.i <- i
pars.i$tau_smooth <- mean(core.dat.3$tau_smooth)
pars.i$tau_r <- mean(core.dat.3$tau_smooth)
spec.err.lst <- lapply(1:length(d_t), function(i) {
if (is.finite(d_t[i])) {
pars.i$delta_t <- d_t[i]
spec.obj <- do.call(ProxyErrorSpectrum, pars.i)
var.obj <- IntegrateErrorSpectra(spec.obj,
method = "sincfilter",
tau_smooth = pars.i$tau_smooth
)
err.obj <- GetProxyError(var.obj,
timescale = pars.i$tau_smooth,
format = "short",
exclude = c("Seasonal.unc.", "Bioturbation")
)
return(err.obj)
} else NA
})
err.df <- do.call(rbind.data.frame, spec.err.lst)
bind_cols(core.dat.3, err.df)
})
core.dat.4 <- do.call(rbind.data.frame, core.dat.4)
# # check that error is scaling with n
# core.dat.4 %>%
# ggplot(aes(x = n, y = Total.error^2, colour = record)) +
# geom_point(position = position_jitter(width = 0.01)) +
# expand_limits(y = 0) +
# scale_x_log10() +
# scale_y_log10() +
# geom_smooth(method = "lm") +
# geom_smooth(method = "lm", formula = y~offset(-x), linetype = 2,
# se = F)
Note: Bioturbation error is excluded as this is deterministic and should be the same for both of these replicate records.
p.rep.records.smoothed <- core.dat.4 %>%
ungroup() %>%
mutate(record = factor(record)) %>%
# filter(complete.cases(mean.value)) %>%
ggplot(aes(
x = age / 1000, y = mean.value,
colour = record
)) +
geom_ribbon(aes(
ymax = mean.value + Total.error,
ymin = mean.value - Total.error,
fill = record, colour = NULL
),
alpha = 0.5
) +
# geom_point(aes(shape = factor(n))) +
geom_line() +
geom_text(aes(label = n), show.legend = FALSE) +
# scale_y_reverse(expression(delta^"18"*O~"(centred)"))+
scale_y_reverse(expression(delta^"18" * O)) +
scale_x_continuous("Age [ka]", limits = c(0, 10), breaks = seq(0, 10, 2)) +
# scale_colour_manual("", values = wesanderson::wes_palette("Rushmore1", 5)[c(4, 5)]) +
# scale_fill_manual("", values = wesanderson::wes_palette("Rushmore1", 5)[c(4, 5)]) +
scale_colour_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) +
scale_fill_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) +
# scale_color_brewer("", palette = "Dark2") +
# scale_fill_brewer("", palette = "Dark2") +
# scale_fill_viridis_d("", option = "viridis") +
# scale_color_viridis_d("", option = "viridis") +
# facet_grid(record~.) +
#expand_limits(x = c(0, 10)) +
theme_bw() +
# labs(subtitle = paste0("Smoothed to ", tau_smooth, " year resolution")) +
labs(caption = expression(Delta[t] ~ "=" ~ tau[smooth] ~ "= 492")) +
theme(panel.grid.minor = element_blank(), legend.position = "none")
#p.rep.records.smoothed
#ggsave("smoothed_r1_r2-1.pdf", p.rep.records.smoothed, width = 7, height = 4.5)
# Get y axis range from plot 1
y.range <- ggplot_build(ts_with_spec_error)$layout$panel_scales_y[[1]]$range$range
p.rep.records.smoothed.2 <- p.rep.records.smoothed +
expand_limits(y = -y.range)
p.geo.comb <- egg::ggarrange(ts_with_spec_error, p.rep.records.smoothed.2,
ncol = 1,
labels = c(letters[1:2]),
label.args = list(gp = grid::gpar(font = 1, cex = 1.2))
)
ggsave("ori_smoothed_r1_r2-1.pdf", p.geo.comb, width = 6, height = 6)
equi <- plyr::ddply(core.dat.4, "n.forams.grp", function(x) {
dt = 10
T = diff(range(x$age))
n = floor(T / dt)
data.frame(
tpt = 1:n,
value = PaleoSpec::MakeEquidistant(x$age, x$mean.value, dt)[1:n],
error = PaleoSpec::MakeEquidistant(x$age, x$Total.inc.error, dt)[1:n]
)
}) %>%
tbl_df() %>%
filter(complete.cases(value, error))
equi.wide <- equi %>%
pivot_wider(names_from = n.forams.grp, values_from = c(value, error))
equi.wide <- equi.wide %>%
mutate(
value.diff = `value_30-30` - `value_5-5`,
err.tot = sqrt(`error_30-30`^2 + `error_5-5`^2),
covered = abs(value.diff) <= err.tot
) %>%
filter(complete.cases(value.diff))
equi.wide %>%
summarise(p.covered = sum(covered) / n())
equi.wide %>%
ggplot(aes(x = tpt, y = value.diff)) +
geom_line() +
geom_point() +
geom_line(aes(y = err.tot)) +
geom_line(aes(y = -err.tot))
core.dat.4.wide <- core.dat.4 %>%
ungroup() %>%
select(age, n.forams.grp, mean.value, Total.error) %>%
pivot_wider(
names_from = n.forams.grp,
values_from = c(mean.value, Total.error)
)
core.dat.4.wide <- core.dat.4.wide %>%
mutate(
value.diff = `mean.value_30-30` - `mean.value_5-5`,
err.tot = sqrt(`Total.error_30-30`^2 + `Total.error_5-5`^2),
covered = abs(value.diff) <= err.tot
) %>%
filter(complete.cases(value.diff))
core.dat.4.wide %>%
summarise(p.covered = sum(covered) / n())
core.dat.4.wide %>%
ggplot(aes(x = age, y = value.diff)) +
geom_line() +
geom_point() +
geom_line(aes(y = err.tot)) +
geom_line(aes(y = -err.tot))
Appendix
Bioturbation and sampling strategy
spec.pars.tau_b.1 <- GetSpecPars("Mg_Ca",
T = 1e04 + 100,
delta_t = 100,
tau_r = 100, sig.sq_a = 0,
tau_b = 500, tau_s = 50,
seas.amp = 8, tau_p = 4 / 12,
N = 30,
phi_c = pi / 2,
phi_a = pi / 2, sigma.cal = 0.3,
clim.spec.fun.args = list(
a = 0.1, b = 1
)
)
var.obj.4 <- crossing(
tau_b = 1e03 * 10 / c(5, 50, 500)
) %>%
group_by(tau_b) %>%
do({
spec.pars.tau_b.1$tau_b <- .$tau_b
spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.tau_b.1)
data.frame(as.list(IntegrateErrorSpectra(spec.obj)[["proxy.error.var"]]))
}) %>%
ungroup() %>%
filter(smoothed.resolution <= 1000 | smoothed.resolution == Inf)
p.tsdv.tau_b <- var.obj.4 %>%
PlotTSDVariance(., c(
"Calibration.unc.", "Seasonal.bias",
"Seasonal.bias.unc."
), include = FALSE) +
facet_wrap(~tau_b, labeller = function(labs) {
sedproxydocs::label_symbols(labs, label.key = plot.symbol)
}) +
expand_limits(x = 10^(log10(range(var.obj.4$smoothed.resolution)) + c(-0.1, 0.1))) +
labs(x = "Smoothed and interpreted resolution [Years]")
# Change tau_r for each time-scale
spec.pars.tau_b.1 <- GetSpecPars("Mg_Ca",
T = 1e04, delta_t = 100,
tau_r = 100, sig.sq_a = 0,
tau_b = 500, tau_s = 50,
seas.amp = 8, tau_p = 4 / 12,
N = 30,
phi_c = pi / 2,
phi_a = pi / 2, sigma.cal = 0.3,
clim.spec.fun.args = list(a = 0.1, b = 1)
)
RoundOdd <- function(x) {
x <- round(x)
il <- x %% 2 < 1
x[il] <- x[il] + 1
return(x)
}
var.obj.5 <- crossing(
delta_t = seq(100, 1000, 100),
tau_b = 1e03 * 10 / c(5, 50, 500)
) %>%
mutate(
T = 10000,
tau_r = delta_t,
T = delta_t * RoundOdd(T / delta_t)
) %>%
group_by(tau_b, delta_t) %>%
do({
# Replace tau_r and delta_t with smoothed resolution
spec.pars.tau_b.1$T <- .$T
spec.pars.tau_b.1$delta_t <- .$delta_t
spec.pars.tau_b.1$tau_r <- .$tau_r
spec.pars.tau_b.1$tau_b <- .$tau_b
spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.tau_b.1)
data.frame(as.list((IntegrateErrorSpectra(spec.obj,
method = "sincfilter",
tau_smooth = .$tau_r
)[["proxy.error.var"]])))
}) %>%
ungroup()
lab.df <- data.frame(
varname = c(
"clim.pow.0.5", "clim.pow.scal", "sed.acc.rate",
"bio.depth", "amp.seas.cycle", "amp.seas.cycle.2", "n.samples",
"resolution", "sigma.meas", "tau_b"
),
label = c("gamma", "beta", "SAR", "d", "Alpha", "Alpha[2]", "N", "R", "sigma", "tau[b]")
)
p.tsdv.tau_b.tau_r <- var.obj.5 %>%
filter(smoothed.resolution != Inf) %>%
select(
-delta_t, -Climate, -Total.error, -Calibration.unc., -Reference.climate,
-Seasonal.bias, -Seasonal.bias.unc.
) %>%
gather(component, value, -smoothed.resolution, -tau_b) %>%
ungroup() %>%
mutate(component = psem:::OrderStages(component, rev = TRUE)) %>%
ggplot(aes(x = smoothed.resolution, y = value, fill = component, colour = component)) +
# geom_line()+
scale_color_manual("", values = spec.colrs) +
geom_area(alpha = 0.75) +
scale_x_continuous(trans = ecustools::reverselog_trans(10)) +
scale_color_manual("", values = spec.colrs) +
scale_fill_manual("", values = spec.colrs) +
labs(y = expression("Error variance [°C"^2 ~ "]"), x = "Physical sampling and interpreted resolution [Years]") +
theme_bw() +
theme(panel.grid.minor = element_blank(), legend.position = "none") +
facet_wrap(~tau_b, labeller = function(labs) {
sedproxydocs::label_symbols(labs, label.key = lab.df)
}) +
expand_limits(x = 10^(log10(range(var.obj.5$smoothed.resolution)) + c(-0.1, 0.1)))
egg::ggarrange(p.tsdv.tau_b, p.tsdv.tau_b.tau_r, ncol = 1, labels = c(letters[1:2]),
label.args = list(gp = grid::gpar(font = 1, cex = 1.2)))
EarthSystemDiagnostics/spectraluncertainty documentation built on Jan. 8, 2020, 9:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Alternative names for psem
- SPR(UCE) Spectral Proxy Record UnCErtainty
- SUPR Spectral Uncertainty for Proxy Records
# Detect OS get_os <- function() { if (.Platform$OS.type == "windows") { "win" } else if (Sys.info()["sysname"] == "Darwin") { "mac" } else if (.Platform$OS.type == "unix") { "unix" } else { stop("Unknown OS") } } os <- get_os() knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, cache = TRUE, fig.width = 7, fig.height = 4.5, fig.pos = "H", dpi = 300, dev = ifelse(os == "win", "pdf", "cairo_pdf"), autodep = TRUE ) # detach("package:psem", unload = TRUE) # detach("package:sedproxy", unload = TRUE) library(psem) #library(spectraluncertainty) library(sedproxy) library(tidyverse) library(corit)
# lookup tables and plotting keys # Formats parameter values for tables FormatParValues <- function(v, max.digits = 4, max.char = 6) { v <- as.numeric(v) v <- as.list(v) v2 <- v matches <- c(as.character(pi), as.character(pi / 2), as.character(pi / 3), as.character(pi / 4), as.character(2*pi), as.character(exp(1)), as.character(2*pi/3)) replacements <- c("pi", "pi/2", "pi/3", "pi/4", "2*pi", "e", "2pi/3") hc <- lapply(v, function(x) if(as.character(x) %in% matches){replacements[match(x, matches)]}else NA) rat <- lapply(v2, function(x) as.character(MASS::as.fractions(x))) sig <- lapply(v2, function(x) signif(x, max.digits)) nc.hc <- lapply(hc, nchar) nc.rat <- lapply(rat, nchar) nc.sig <- lapply(sig, nchar) nc.hc[is.na(nc.hc)] <- Inf nc.rat[is.na(nc.rat)] <- Inf nc.sig[is.na(nc.sig)] <- Inf v4 <- lapply(1:length(v), function(i) { if (is.na(hc[[i]]) == FALSE & nc.hc[[i]] <= max.char) { hc[[i]] } else if (nc.rat[[i]] < max.char) { rat[[i]] } else { sig[[i]] } }) return(as.character(v4)) }
Error types
set.seed(4) n = 20 df <- tibble( Age = seq(1, 10, length.out = n), Climate = PaleoSpec::SimPowerlaw(1, n) + SSlogis(Age, 4, 5, 0.1) - 2, `Independent errors` = Climate + rnorm(n, 0, 1.5), `Constant errors` = Climate - 1.5, `Correlated errors` = Climate + (Age - 5) * -0.58 ) df.long <- df %>% gather(Noise.type, value, -Age, -Climate) %>% mutate(Noise.type = factor(Noise.type, levels = (c( "Independent errors", "Correlated errors", "Constant errors" )), ordered = TRUE ))
txt.lbls <- tibble( Age = -Inf, Climate = Inf, Noise.type = factor(levels(df.long$Noise.type), ordered = T), label = letters[1:3] ) p.error.wide <- df.long %>% ggplot(aes(x = Age, y = Climate)) + geom_line(aes(linetype = "Climate")) + geom_line(aes(y = value, colour = Noise.type, linetype = "Proxy")) + geom_segment(aes(x = Age, xend = Age, y = Climate, yend = value, colour = Noise.type), arrow = arrow(length = unit(0.03, "npc"), type = "closed") ) + facet_wrap(~Noise.type, ncol = 3) + theme_bw() + theme(legend.position = "left", panel.grid.minor = element_blank()) + scale_colour_manual( values = wesanderson::wes_palette("GrandBudapest1", 4)[2:4], guide = FALSE ) + scale_linetype("", guide = FALSE) + labs(x = "Age [ka]") + geom_text(data = txt.lbls, aes(label = label), hjust = -1, vjust = 1.5, fontface = 4) p.error.wide
Simplified method description
spec.pars.smpl <- GetSpecPars("Mg_Ca", T = 1e04 + 100, delta_t = 100, tau_r = 100, tau_b = 500, tau_s = 50, N = 1 ) # Fourier transform of bioturbation and sediment slice thickness filter TF.f_bs <- function(nu, tau_b, tau_s) { i <- 0 + 1i sinc(nu * tau_s) * (1 + i * 2 * pi * nu * tau_b)^-1 * exp(i * 2 * pi * nu * tau_b) } # The redistributed stochastic climate variance smpl.spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.smpl) redist.stoch.clim.var <- smpl.spec.obj$proxy.error.spec$Aliasing.stochastic[smpl.spec.obj$proxy.error.spec$nu<0] dat <- tibble( nu = GetNu(spec.pars.smpl$T, spec.pars.smpl$delta_t), Climate.spectrum = 0.1 * nu^-1 ) %>% filter(nu > 0) %>% mutate( # the transfer fucntion tf.f_bs = TF.f_bs(nu, spec.pars.smpl$tau_b, spec.pars.smpl$tau_s), # squared absolute value of the transfer function sq.abs.tf.f_bs = abs(tf.f_bs)^2, Filtered.climate.spectrum = Climate.spectrum * sq.abs.tf.f_bs, # transfer function for the reference climate, sinc function tf.ref = sinc(nu * spec.pars.smpl$tau_r), sq.abs.tf.ref = abs(tf.ref)^2, Reference.climate.spectrum = sq.abs.tf.ref * Climate.spectrum, Error.spectrum = abs(tf.f_bs - tf.ref)^2 * Climate.spectrum, Redist.stoch.clim.var = redist.stoch.clim.var ) dat.long <- dat %>% gather(Component, Spec, -nu) %>% mutate( line.type = ifelse(Component == "Error.spectrum", "Error", "Signal"), Spec = Re(Spec) ) clr.pal <- c((RColorBrewer::brewer.pal(9, "Set1"))[c(1,4,2,7, 9, 9)]) names(clr.pal) <- c( "Climate.spectrum", "Reference.climate.spectrum", "Filtered.climate.spectrum", "Error.spectrum", "Redist.as.white.noise", "Redist.stoch.clim.var" ) clr.lbls <- c("Climate spectrum", "Reference climate\nspectrum", "Climate spectrum\nafter bioturbation", "Spectrum of\nbioturbation error", "Redistributed\nas white noise", "Redistributed\nclimate variance") names(clr.lbls) <- names(clr.pal) ltyp.pal <- c(1,1,1,2,1,1) names(ltyp.pal) <- names(clr.pal) p.a <- dat.long %>% filter(Component %in% c( "Climate.spectrum", "Error.spectrum", "Filtered.climate.spectrum", "Reference.climate.spectrum" )) %>% ggplot(aes(x = nu, y = Spec, colour = Component, linetype = Component)) + geom_line() + scale_y_continuous("Power spectral density", trans = "log10") + scale_x_continuous("Frequency [1/years]", trans = "log10") + scale_linetype_manual("", values = ltyp.pal, labels = clr.lbls, limits = names(clr.pal)[1:4]) + scale_colour_manual("", values = clr.pal, labels = clr.lbls, limits = names(clr.pal)[1:4]) + theme_bw() + theme(panel.grid.minor = element_blank(), legend.key.height = unit(2, "line")) p.b <- dat.long %>% filter(Component %in% c( "Bioturbation", "Reference.climate.spectrum", "tf.ref", "tf.f_bs", "sq.abs.tf.f_bs", "sq.abs.tf.ref", "Error.spectrum" ) == FALSE) %>% ggplot(aes(x = nu, y = Spec, colour = Component)) + geom_line() + geom_ribbon( data = dat, aes( x = nu, ymin = Filtered.climate.spectrum, ymax = Climate.spectrum, fill = "Redist.as.white.noise" ), alpha = 0.5, inherit.aes = FALSE ) + scale_y_continuous("Power spectral density", trans = "log10") + scale_x_continuous("Frequency [1/years]", trans = "log10") + scale_fill_manual("", values = clr.pal, labels = clr.lbls) + scale_colour_manual("", values = clr.pal, labels = clr.lbls) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.key.height = unit(2, "line")) + guides(colour = guide_legend(order = 1))
p.a.b <- egg::ggarrange(p.a, p.b, ncol = 1, draw = FALSE, labels = c("a", "b"), label.args = list(gp = grid::gpar(font = 1, cex = 1.2))) print(p.a.b)
Power spectrum of proxy error
# PSD Climate example.lat <- 20 clim.spec.ex1 <- ModelSpectrum( #freq = GetNu(100000, 100), freq = NULL, latitude = example.lat, variable = "temperature", beta = 1 ) p.clim.spec.ex1 <- PlotModelSpectrum(clim.spec.ex1) p.clim.spec.ex1
# Proxy seasonality # climproxyrecords::shakun.metadata %>% # filter(Core == "ODP 658C") # # ecusdata::shak.fraile.seasonality %>% # filter(ID.no == "N31") %>% # gather(Month, value, starts_with("v")) seasonal.amp <- AmpFromLocation( longitude = -18, latitude = example.lat, proxy.type = "degC", depth.upr = 0, depth.lwr = -120 ) orbital.pars <- RelativeAmplitudeModulation( latitude = example.lat, maxTimeKYear = 23, minTimeKYear = 1, bPlot = FALSE ) # Look for good example location # ssa <- sapply(-90:90, function(l) RelativeAmplitudeModulation(latitude = l, # maxTimeKYear = 23, # minTimeKYear = 1, # bPlot = FALSE)$sig.sq_a) # # ssc <- sapply(-90:90, function(l) AmpFromLocation(longitude = -28, # latitude = l, # proxy.type = "degC", # depth.upr = 0, depth.lwr = -200)$sig.sq_c) # plot(-90:90, ssa, type = "l") # abline(v = 20) # plot(-90:90, ssc, type = "l") # abline(v = 20) ex.sed.acc.rate <- 10 spec.pars.ex1 <- GetSpecPars("Mg_Ca", T = 1e04 + 100, delta_t = 100, tau_r = 100, sig.sq_a = orbital.pars$sig.sq_a, sig.sq_c = round(seasonal.amp$sig.sq_c, 3), tau_b = 1000 * 10 / ex.sed.acc.rate, tau_s = 1000 * 1 / ex.sed.acc.rate, N = 30, tau_p = 7/12, phi_c = 0, delta_phi_c = 2 * pi / 3, phi_a = pi / 2, sigma.cal = 0.25, sigma.meas = 0.25, sigma.ind = 1, clim.spec.fun = "ModelSpectrum", clim.spec.fun.args = list(latitude = example.lat, beta = 1) ) spec.obj.ex1 <- do.call(ProxyErrorSpectrum, spec.pars.ex1)
p.spec.1 <- spec.obj.ex1 %>% PlotSpecError(.) p.spec.1 + coord_cartesian(ylim = c(2e-1, 1e04)) #+ geom_vline(xintercept = 1/c(23e03))
# FormatParValues(c(2*pi, 1/3, 0.1, pi, pi/2, exp(1), 2.12457896, 2*pi/3)) par.df <- as.data.frame(spec.pars.ex1) %>% gather(parameter, value) %>% left_join(parameter.symbol, .) %>% left_join(., parameter.description) %>% left_join(., parameter.source) %>% filter(complete.cases(latex)) %>% select(latex, value, description, source) %>% mutate(value = FormatParValues(value, max.digits = 2, max.char = 6)) %>% tbl_df() names(par.df) <- c("Parameter", "Value", "Description", "Source") # par.df %>% # knitr::kable(.) pander::pandoc.table(par.df)
# TEX for direct pasting into ms. par.df$Description <- gsub("pi", "\\(\\pi\\)", x = par.df$Description, fixed = TRUE) par.df$Source <- gsub("pi", "\\(\\pi\\)", x = par.df$Source, fixed = TRUE) par.df$Value <- gsub("pi", "\\(\\pi\\)", x = par.df$Value, fixed = TRUE) cat(apply(par.df[, 1:4], 1, function(x) paste0(paste0(x, collapse = " & "), " \\\\")))
Timescale dependent error
var.obj.ex1 <- IntegrateErrorSpectra(spec.obj.ex1) ev.a <- var.obj.ex1 %>% PlotTSDVariance(., include.constant.errors = TRUE) + theme(legend.position = "none") ev.b <- var.obj.ex1 %>% PlotTSDVariance(.) ev.a.b <- egg::ggarrange( ev.a, ev.b, ncol = 2, draw = FALSE, labels = c("a", "b"), label.args = list(gp = grid::gpar(font = 1, cex = 1.2))) ev.a.b.2 <- grid::forceGrob(ev.a.b) ev.a.b.2
Err diff 2 timeslices example spectrum
tslice.width.1 <- 1000 tslice.width.2 <- 1000 d.ts <- spec.pars.ex1$T - (tslice.width.1/2 + tslice.width.2/2) err.tslice.diff <- ErrDiff2TimeSlices(spec.obj.ex1, tau_1 = tslice.width.1, tau_2 = tslice.width.2, delta_ts = d.ts) err.tslice.diff <- err.tslice.diff %>% filter(component == "Total.error") %>% select(-cov, -sigma.diff, -var.tau.12) %>% mutate(naive.var.diff = var.tau.1 + var.tau.2) err.tslice.diff %>% rename(`Error\ntime-slice a` = var.tau.1, `Error\ntime-slice b` = var.tau.2, `Error\ntime-slice\ndifference` = var.diff, `Naive error\ntime-slice\ndifference` = naive.var.diff) %>% gather(s, var, -component) %>% mutate(sigma = sqrt(var), s = factor(s, ordered = TRUE, levels = s)) %>% ggplot(aes(x = s, y = sigma)) + #geom_point() + geom_col(alpha = 0.75, width = 0.75) + theme_bw() + #labs(y = expression("Uncertainty" ~""%+-%""~"1"*sigma~" [°C]"), x = "") + labs(y = expression("Error - 1"*sigma~" [°C]"), x = "")
GeoB_10054_4
# Core GeoB 10054-4 core.info <- replicatecoredata::core.info %>% filter(Core == "GeoB 10054-4") world <- map_data("world") # world$long[world$long > 180] <- (360 - world$long[world$long > 180]) * -1 coremap <- ggplot(world, aes(x = long, y = lat, group = group)) + geom_polygon(colour = "Darkgrey", fill = "Darkgrey", size = 0.5, show.legend = F) + coord_map("gall", 110, xlim = c(50, 180), ylim = c(-50, 50) ) + geom_point( data = core.info, aes(Longitude, Latitude), inherit.aes = FALSE, show.legend = T, colour = "Red", size = 2, shape = 8 ) + theme_bw() + theme(legend.position = "top") + scale_x_continuous("Longitude", labels = function(x) paste0(x, "°")) + scale_y_continuous("Latitude", labels = function(x) paste0(x, "°")) coremap
# PSD Climate clim.spec <- ModelSpectrum( freq = NULL, latitude = core.info$Latitude, variable = "d18O", beta = 1 ) p.clim.spec <- PlotModelSpectrum(clim.spec) p.clim.spec
breit.lons <- unique(breitkreuz.amp$longitude) nearest.lon <- breit.lons[which.min(abs(core.info$Longitude - breit.lons))] breit.lats <- unique(breitkreuz.amp$latitude) nearest.lat <- breit.lats[which.min(abs(core.info$Latitude - breit.lats))] seasonal.amp <- AmpFromLocation( longitude = core.info$Longitude, latitude = core.info$Latitude, proxy.type = "d18O", depth.upr = 0, depth.lwr = -50 )
core.dat <- replicatecoredata::GeoB.vars %>% filter( variable == "d18O", age.14C.cal <= 1e04, core %in% c("GeoB 10054-4"), taxon == "G. ruber", n.forams.grp %in% c("30-30", "5-5") ) record.names <- tibble( n.forams.grp = c("30-30", "5-5"), record = c("Rep2", "Rep1") ) core.dat <- left_join(core.dat, record.names) timescale.pars <- core.dat %>% group_by(n.forams.grp) %>% summarise( T = diff(range(age.14C.cal)), mu_delta_t = round(T / n()), N = mean(n.forams) ) %>% left_join(., record.names) orbital.pars <- RelativeAmplitudeModulation( latitude = core.info$Latitude, maxTimeKYear = 23, minTimeKYear = 1, bPlot = FALSE )
approx.s <- round(1000 * diff(range(core.dat$depth_cm)) / diff(range(core.dat$age.14C.cal))) ## estimate of tau_b from replicated 14C measurements tau_b.10054_4 <- 720 tau_b.10054_4.se <- ecustools::SDofSD(tau_b.10054_4, 10) 1 / 1000 * c(tau_b.10054_4 - tau_b.10054_4.se, tau_b.10054_4, tau_b.10054_4 + tau_b.10054_4.se) * 20 spec.pars.lst <- timescale.pars %>% plyr::dlply(., c("record"), function(x) { GetSpecPars( proxy.type = "d18O", delta_t = round(x$mu_delta_t), tau_r = x$mu_delta_t, T = x$T, # Bioturbation tau_b = tau_b.10054_4, tau_s = 1000 * 1 / approx.s, N = x$N, sig.sq_c = seasonal.amp$sig.sq_c, sigma.meas = 0.1, sigma.ind = 0.15, tau_p = 1, phi_c = 0, sig.sq_a = orbital.pars$sig.sq_a, clim.spec.fun = "ModelSpectrum", clim.spec.fun.args = list( latitude = core.info$Latitude, beta = 1, variable = "d18O" ) ) }) spec.obj.lst <- plyr::llply(spec.pars.lst, function(x) do.call(ProxyErrorSpectrum, x))
reps.par.tab <- t(plyr::ldply(spec.pars.lst, function(x) { y <- as.data.frame(x) y })) colnames(reps.par.tab) <- reps.par.tab[1, ] geob10054.pars <- reps.par.tab %>% data.frame(., stringsAsFactors = FALSE) %>% # gather(parameter, value) %>% tibble::rownames_to_column("parameter") %>% left_join(., parameter.symbol) %>% filter(complete.cases(latex)) %>% select(latex, Rep1, Rep2) %>% mutate( Rep1 = FormatParValues(Rep1), Rep2 = FormatParValues(Rep2) ) %>% # mutate(value = ifelse(value == as.character(pi), "pi", value)) %>% tbl_df() pander::pandoc.table(geob10054.pars)
spec.obj.df <- plyr::ldply(spec.obj.lst, function(x) x[[2]], .id = "record") plot.symbol <- parameter.symbol %>% rename(varname = parameter) %>% add_row(varname = "record", label = "R") p <- spec.obj.df %>% filter(nu > 0) %>% gather(component, spec, -nu, -record) %>% left_join(., timescale.pars) %>% filter(spec > 1e-20) %>% mutate(component = psem:::OrderStages(component)) %>% ggplot(., aes(x = nu, y = spec, colour = component, linetype = component)) + geom_line() + scale_x_continuous(trans = "log10") + scale_y_continuous(trans = "log10") + scale_color_manual("", values = spec.colrs) + scale_linetype_manual("", values = spec.linetypes) + theme_bw() + #theme(strip.background = element_blank(), strip.text.y = element_blank()) + annotation_logticks(sides = "lb") + labs(x = "Frequency [1/years]", y = "Power Spectral Density") + facet_wrap(~record + N + mu_delta_t, labeller = function(labs) { sedproxydocs::label_symbols(labs, label.key = plot.symbol)}) p
var.obj.lst <- plyr::llply(spec.obj.lst, IntegrateErrorSpectra) err.obj.lst <- plyr::llply(var.obj.lst, function(x) { GetProxyError(x, timescale = x[["spec.pars"]]$delta_t) }) # Smooth higher res record to res of lower res record var.obj.1b <- IntegrateErrorSpectra(spec.obj.lst[["Rep1"]], tau_smooth = spec.pars.lst[["Rep2"]]$delta_t ) err.obj.1b <- GetProxyError(var.obj.1b, timescale = spec.pars.lst[["Rep2"]]$delta_t) err.obj.1b$record <- "Rep1" err.obj.1.2 <- plyr::ldply(err.obj.lst, .id = "record") var.obj.1.2 <- plyr::ldply(var.obj.lst, function(x) x[[2]], .id = "record")
TSD Variance
var.492 <- lapply(spec.obj.lst, function(x) IntegrateErrorSpectra(x, tau_smooth = 492)) err.492 <- lapply(var.492, function(x) GetProxyError(x, timescale = 492)) err.492.df <- plyr::ldply(err.492, .id = "record") %>% group_by(record, smoothed.resolution) %>% filter(component %in% c("Total.error", "Reference.climate") == FALSE) %>% summarise(Total = sum(exl.f.zero^2)) %>% mutate(smoothed.resolution = as.numeric(smoothed.resolution))
# small df for red dots highlighting total error at tau_smooth = 246 err.df <- err.obj.1.2 %>% filter(record == "Rep2") %>% bind_rows(., err.obj.1b) %>% filter(component == "Total.error") %>% mutate(x = 246, y = exl.f.zero^2) subplt.lbls <- tibble(record = factor(c("Rep1", "Rep2")), x = 10000, y = Inf, label = letters[1:2]) PlotTSDVariance(var.obj.1.2, units = "d18O", include.constant.errors = FALSE) + facet_wrap(~record) + #geom_point(data = err.df, aes(x = x, y = y), colour = "Red", inherit.aes = F) + #geom_point(data = err.492.df, aes(x = smoothed.resolution, y = Total), # colour = "Red", inherit.aes = F) + geom_vline(xintercept = 492, linetype = 3) + theme(panel.spacing = unit(1, "lines")) + #annotate(geom = "text", x = 10000, y = Inf, label = "a") + scale_x_continuous(trans = ecustools::reverselog_trans(), breaks = c(100, 492, 1000, 10000), labels = c(100, 492, 1000, 10000)/1000) + labs(y = expression("Error variance [ ‰ d"^18*"O"^2 ~"]"), x = "Timescale [kyr]") + ## Add subplot labels a, b, etc outside facet plotting areas geom_text(data = subplt.lbls, aes(x = x, y = y, label = label), inherit.aes = FALSE, hjust = c(5.5, 3), vjust = 1, size = 5) + # makes points / text outside plotting area visible coord_cartesian(clip = "off") + # compensate for the thicker border due to the above theme(panel.border = element_rect(size = 1/4), panel.spacing.x = unit(2, "lines"), strip.background = element_blank(), strip.text = element_blank())
err.obj.df <- plyr::ldply(var.obj.lst, function(x) { GetProxyError(x, exclude = "Bioturbation", timescale = x[["spec.pars"]]$delta_t, format = "short" ) }, .id = "record") %>% mutate(sigma.total = Total.inc.error) core.dat.2 <- left_join(core.dat, err.obj.df) %>% group_by(depth_cm) %>% mutate(both.sampling.strat = ifelse(length(unique(n.forams.grp)) == 2, "Same slice", "Extra slice" )) ts_with_spec_error <- core.dat.2 %>% # filter(both.sampling.strat == TRUE) %>% ggplot(aes( x = age.14C.cal / 1000, y = value, fill = record, colour = record, group = n.forams.grp )) + geom_ribbon(aes( ymax = value + sigma.total, ymin = value - sigma.total, colour = NULL ), alpha = 0.5 ) + geom_point() + geom_line() + scale_y_reverse(expression(delta^"18" * O)) + scale_x_continuous("Age [ka]", limits = c(0, 10), breaks = seq(0, 10, 2)) + scale_colour_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) + scale_fill_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) + expand_limits(x = 0) + theme_bw() + theme(panel.grid.minor = element_blank()) ts_with_spec_error_facets <- ts_with_spec_error + facet_grid(record ~ .) #ts_with_spec_error_facets
Look only at samples from the same slices.
# core.dat.2 %>% # # filter(both.sampling.strat == TRUE) %>% # ggplot(aes( # x = age.14C.cal / 1000, y = value, # group = n.forams.grp, # fill = record, # colour = record # )) + # geom_ribbon(aes( # ymax = value + sigma.total, # ymin = value - sigma.total, # colour = NULL # ), # alpha = 0.5 # ) + # geom_point() + # geom_line() + # scale_y_reverse(expression(delta^"18" * O)) + # scale_x_continuous("Age [ka]") + # facet_grid(both.sampling.strat ~ .) + # expand_limits(x = 0) + # theme_bw() + # theme(panel.grid.minor = element_blank()) # # core.dat.2.coverage <- core.dat.2 %>% # filter(both.sampling.strat == TRUE) %>% # group_by(depth_cm) %>% # mutate(n.samples = n()) %>% # filter(n.samples == 2) %>% # summarise( # value.diff = mean(diff(value)), # sigma.total = sqrt(sum(sigma.total^2)), # covered = abs(value.diff) <= sigma.total # ) # # core.dat.2.coverage %>% # ungroup() %>% # summarise(coverage = sum(covered) / n())
Expected correlation
expected.correlations <- plyr::ldply(spec.obj.lst, ExpectedCorrelation, .id = "record") %>% left_join(., record.names) %>% mutate(n.forams.grp = paste0("Expected ", record, "-", ifelse(correlation.type == "proxy-proxy", record, "Climate")))
# # Interpolated and average expected correlations # # expected.correlations %>% # group_by(correlation.type, record) %>% # filter(complete.cases(smoothed.resolution, rho)) %>% # do({ # smoothed.resolution <- seq(min(.$smoothed.resolution), max(.$smoothed.resolution), by = 1) # rho <- approx(.$smoothed.resolution, y = .$rho, xout = smoothed.resolution)$y # data.frame(smoothed.resolution, rho) # }) %>% # spread(record, rho) %>% # filter(complete.cases(Rep1, Rep2)) %>% # mutate(rho = (Rep1 + Rep2) / 2) %>% # ggplot(aes(x = smoothed.resolution, y = rho, colour = correlation.type)) + # geom_line() #
time.series1 <- plyr::dlply(core.dat.2, "record", function(x) zoo(x$value, order.by = x$age.14C.cal)) Cor <- #expected.correlations %>% # select(smoothed.resolution, record) %>% # rename(timescale = smoothed.resolution) %>% # filter(timescale <= 1500) %>% # distinct() %>% data.frame(timescale = seq(246, 1500, by = 1)) %>% group_by(timescale) %>% do({ tcor <- corit::CorIrregTimser( timser1 = time.series1[[1]], timser2 = time.series1[[2]], detr = TRUE, # remove linear trend time series method = "InterpolationMethod", appliedFilter = "runmean", fc = 1 / (2*.$timescale), # cut-off frequency #appliedFilter = "gauss", #fc = 1 / (.$timescale), # cut-off frequency dt = .$timescale / 10, # time step for the interpolation int.method = "linear", # kind of interpolation filt.output = TRUE) # number of valid points in smoothed timeseries # used for SE/CI of rho data.frame(rho = tcor$cor, n = sum(is.na(tcor$ft1)==FALSE)) }) %>% mutate(correlation.type = "Observed Rep1-Rep2\n(with corit)") rho.CI <- Cor %>% select(-rho, -correlation.type) %>% rename(smoothed.resolution = timescale) %>% # join to add n to expected correlations left_join(expected.correlations, .) %>% group_by(record, correlation.type) %>% do({ rhoCI(.$rho, n = .$n, alpha = 1-0.682689492137) }) %>% filter(complete.cases(n)) expected.correlations.ci <- expected.correlations %>% left_join(., rho.CI) %>% filter( smoothed.resolution <= 1500)
expected.correlations.ci.mean <- expected.correlations.ci %>% group_by(correlation.type, record) %>% filter(complete.cases(smoothed.resolution, rho)) %>% do({ smoothed.resolution <- seq(min(.$smoothed.resolution), max(.$smoothed.resolution), by = 1) rho <- approx(.$smoothed.resolution, y = .$rho, xout = smoothed.resolution)$y lwr <- approx(.$smoothed.resolution, y = .$lwr, xout = smoothed.resolution)$y upr <- approx(.$smoothed.resolution, y = .$upr, xout = smoothed.resolution)$y data.frame(smoothed.resolution, rho, lwr, upr) }) %>% group_by(correlation.type, smoothed.resolution) %>% mutate(n = n()) %>% filter(n == 2) %>% summarise_if(is.numeric, mean) %>% ungroup() %>% mutate(correlation.type = ifelse(correlation.type == "proxy-proxy", "Expected Proxy-Proxy", "Expected Proxy-Climate"))
corr.clrs <- wesanderson::wes_palette("GrandBudapest1", 4)[c(2,4,3)] names(corr.clrs) <- c( "Expected Proxy-Proxy", "Expected Proxy-Climate", "Observed Rep1-Rep2\n(with corit)" ) corr.ltyp <- c(1,2,1) names(corr.ltyp) <- names(corr.clrs) leg.order <- c("Expected Proxy-Climate", "Observed Rep1-Rep2\n(with corit)", "Expected Proxy-Proxy") expected.correlations.ci.mean %>% ggplot(aes( x = smoothed.resolution, y = rho, colour = correlation.type )) + geom_ribbon(aes( ymax = upr, ymin = lwr, fill = correlation.type ), colour = NA, alpha = 0.5, show.legend = FALSE ) + geom_line(aes()) + geom_line(data = filter(Cor, timescale < 1000), aes( x = timescale, y = rho )) + theme_bw() + expand_limits(x = c(0, 1500), y = 0) + scale_color_manual("", values = corr.clrs, limits = leg.order) + scale_fill_manual("1 standard error", values = corr.clrs, limits = leg.order) + scale_linetype_manual("", values = corr.ltyp, limits = leg.order)+ scale_x_continuous("Timescale [years]") + scale_y_continuous("Correlation coefficient") + #guides(linetype = guide_legend(override.aes = list(linetype = c(2,2,1,1,1)))) + theme(panel.grid.minor = element_blank(), legend.key.height = unit(1.5, "line")) + expand_limits(y = c(0, 1))
# #corr.clrs <- RColorBrewer::brewer.pal(5, name = "Set2") # #corr.clrs <- wesanderson::wes_palette("GrandBudapest1", 3)[c(1,2,1,2,3)] # corr.clrs <- c(rep(RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)],2), "#5B1A18") # # names(corr.clrs) <- c( # "Expected Rep1-Rep1", # "Expected Rep2-Rep2", # "Expected Rep1-Climate", # "Expected Rep2-Climate", # "Observed Rep1-Rep2\n(with corit)" # ) # # corr.ltyp <- c(1,1,2,2,1) # names(corr.ltyp) <- c( # "Expected Rep1-Rep1", # "Expected Rep2-Rep2", # "Expected Rep1-Climate", # "Expected Rep2-Climate", # "Observed Rep1-Rep2\n(with corit)" # ) # # leg.order <- c("Expected Rep1-Climate", # "Expected Rep2-Climate", # "Expected Rep1-Rep1", # "Observed Rep1-Rep2\n(with corit)", # "Expected Rep2-Rep2") # # expected.correlations.ci %>% # ggplot(aes( # x = smoothed.resolution, y = rho, # colour = n.forams.grp # )) + # geom_ribbon(aes( # ymax = upr, ymin = lwr, # fill = n.forams.grp # ), # colour = NA, alpha = 0.5, show.legend = FALSE # ) + # geom_line(aes(linetype = n.forams.grp)) + # geom_line(data = Cor, aes( # x = timescale, y = rho # )) + # # geom_point(data = Cor, aes( # # x = timescale, y = rho # # )) + # theme_bw() + # expand_limits(x = c(0, 1500), y = 0) + # scale_color_manual("", values = corr.clrs, # limits = leg.order) + # scale_fill_manual("1 standard error", values = corr.clrs, # limits = leg.order) + # scale_linetype_manual("", values = corr.ltyp, # limits = leg.order)+ # scale_x_continuous("Timescale [years]") + # scale_y_continuous("Correlation coefficient") + # guides(linetype = guide_legend(override.aes = list(linetype = c(2,2,1,1,1)))) + # theme(panel.grid = element_blank(), legend.key.height = unit(1.5, "line")) + # expand_limits(y = c(0, 1))
expected.correlations.ci %>% select(correlation.type, smoothed.resolution, rho, n, lwr, upr, width) %>% group_by(correlation.type, smoothed.resolution) %>% summarise_if(is.numeric, mean)
# SmoothIrregular2 <- function(df, time.var, value.var, tau_smooth) { # mean.d_t <- floor(mean(diff(df[[time.var]]))) # # strt.time <- mean.d_t * round((min(df[[time.var]]) - mean.d_t / 2) / mean.d_t) # end.time <- mean.d_t * round((max(df[[time.var]]) + mean.d_t / 2) / mean.d_t) # # age <- seq(strt.time, end.time, by = tau_smooth) # # lst <- lapply(age, function(a) { # df.sub <- subset(df, df[[time.var]] > a - tau_smooth / 2 & df[[time.var]] < a + tau_smooth / 2) # list( # mean.value = mean(df.sub[[value.var]]), n = nrow(df.sub), # mean.d_t = tau_smooth / nrow(df.sub), tau_smooth = tau_smooth # ) # }) # # df.2 <- do.call(rbind.data.frame, lst) # cbind(age, df.2) # } core.dat %>% group_by(record) %>% summarise( max.dT = max(diff(age.14C.cal)), mean.dT = mean(diff(age.14C.cal)) ) tau_smooth <- 2 * spec.pars.lst[["Rep2"]]$delta_t core.dat.3 <- core.dat.2 %>% group_by(n.forams.grp, record, n.forams) %>% do({ # detrend? # .$value <- residuals(lm(.$value~.$age.14C.cal)) SmoothIrregular(., "age.14C.cal", "value", tau_smooth = tau_smooth, d_t = tau_smooth ) })
core.dat.4 <- lapply(spec.pars.lst, function(i) { core.dat.3 <- core.dat.3 %>% filter(n.forams == i$N) d_t <- core.dat.3 %>% pull(mean.d_t) pars.i <- i pars.i$tau_smooth <- mean(core.dat.3$tau_smooth) pars.i$tau_r <- mean(core.dat.3$tau_smooth) spec.err.lst <- lapply(1:length(d_t), function(i) { if (is.finite(d_t[i])) { pars.i$delta_t <- d_t[i] spec.obj <- do.call(ProxyErrorSpectrum, pars.i) var.obj <- IntegrateErrorSpectra(spec.obj, method = "sincfilter", tau_smooth = pars.i$tau_smooth ) err.obj <- GetProxyError(var.obj, timescale = pars.i$tau_smooth, format = "short", exclude = c("Seasonal.unc.", "Bioturbation") ) return(err.obj) } else NA }) err.df <- do.call(rbind.data.frame, spec.err.lst) bind_cols(core.dat.3, err.df) }) core.dat.4 <- do.call(rbind.data.frame, core.dat.4)
# # check that error is scaling with n # core.dat.4 %>% # ggplot(aes(x = n, y = Total.error^2, colour = record)) + # geom_point(position = position_jitter(width = 0.01)) + # expand_limits(y = 0) + # scale_x_log10() + # scale_y_log10() + # geom_smooth(method = "lm") + # geom_smooth(method = "lm", formula = y~offset(-x), linetype = 2, # se = F)
Note: Bioturbation error is excluded as this is deterministic and should be the same for both of these replicate records.
p.rep.records.smoothed <- core.dat.4 %>% ungroup() %>% mutate(record = factor(record)) %>% # filter(complete.cases(mean.value)) %>% ggplot(aes( x = age / 1000, y = mean.value, colour = record )) + geom_ribbon(aes( ymax = mean.value + Total.error, ymin = mean.value - Total.error, fill = record, colour = NULL ), alpha = 0.5 ) + # geom_point(aes(shape = factor(n))) + geom_line() + geom_text(aes(label = n), show.legend = FALSE) + # scale_y_reverse(expression(delta^"18"*O~"(centred)"))+ scale_y_reverse(expression(delta^"18" * O)) + scale_x_continuous("Age [ka]", limits = c(0, 10), breaks = seq(0, 10, 2)) + # scale_colour_manual("", values = wesanderson::wes_palette("Rushmore1", 5)[c(4, 5)]) + # scale_fill_manual("", values = wesanderson::wes_palette("Rushmore1", 5)[c(4, 5)]) + scale_colour_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) + scale_fill_manual("", values = RColorBrewer::brewer.pal(5, name = "Dark2")[c(2, 3)]) + # scale_color_brewer("", palette = "Dark2") + # scale_fill_brewer("", palette = "Dark2") + # scale_fill_viridis_d("", option = "viridis") + # scale_color_viridis_d("", option = "viridis") + # facet_grid(record~.) + #expand_limits(x = c(0, 10)) + theme_bw() + # labs(subtitle = paste0("Smoothed to ", tau_smooth, " year resolution")) + labs(caption = expression(Delta[t] ~ "=" ~ tau[smooth] ~ "= 492")) + theme(panel.grid.minor = element_blank(), legend.position = "none") #p.rep.records.smoothed #ggsave("smoothed_r1_r2-1.pdf", p.rep.records.smoothed, width = 7, height = 4.5)
# Get y axis range from plot 1 y.range <- ggplot_build(ts_with_spec_error)$layout$panel_scales_y[[1]]$range$range p.rep.records.smoothed.2 <- p.rep.records.smoothed + expand_limits(y = -y.range) p.geo.comb <- egg::ggarrange(ts_with_spec_error, p.rep.records.smoothed.2, ncol = 1, labels = c(letters[1:2]), label.args = list(gp = grid::gpar(font = 1, cex = 1.2)) ) ggsave("ori_smoothed_r1_r2-1.pdf", p.geo.comb, width = 6, height = 6)
equi <- plyr::ddply(core.dat.4, "n.forams.grp", function(x) { dt = 10 T = diff(range(x$age)) n = floor(T / dt) data.frame( tpt = 1:n, value = PaleoSpec::MakeEquidistant(x$age, x$mean.value, dt)[1:n], error = PaleoSpec::MakeEquidistant(x$age, x$Total.inc.error, dt)[1:n] ) }) %>% tbl_df() %>% filter(complete.cases(value, error))
equi.wide <- equi %>% pivot_wider(names_from = n.forams.grp, values_from = c(value, error)) equi.wide <- equi.wide %>% mutate( value.diff = `value_30-30` - `value_5-5`, err.tot = sqrt(`error_30-30`^2 + `error_5-5`^2), covered = abs(value.diff) <= err.tot ) %>% filter(complete.cases(value.diff)) equi.wide %>% summarise(p.covered = sum(covered) / n()) equi.wide %>% ggplot(aes(x = tpt, y = value.diff)) + geom_line() + geom_point() + geom_line(aes(y = err.tot)) + geom_line(aes(y = -err.tot))
core.dat.4.wide <- core.dat.4 %>% ungroup() %>% select(age, n.forams.grp, mean.value, Total.error) %>% pivot_wider( names_from = n.forams.grp, values_from = c(mean.value, Total.error) ) core.dat.4.wide <- core.dat.4.wide %>% mutate( value.diff = `mean.value_30-30` - `mean.value_5-5`, err.tot = sqrt(`Total.error_30-30`^2 + `Total.error_5-5`^2), covered = abs(value.diff) <= err.tot ) %>% filter(complete.cases(value.diff)) core.dat.4.wide %>% summarise(p.covered = sum(covered) / n()) core.dat.4.wide %>% ggplot(aes(x = age, y = value.diff)) + geom_line() + geom_point() + geom_line(aes(y = err.tot)) + geom_line(aes(y = -err.tot))
Appendix
Bioturbation and sampling strategy
spec.pars.tau_b.1 <- GetSpecPars("Mg_Ca", T = 1e04 + 100, delta_t = 100, tau_r = 100, sig.sq_a = 0, tau_b = 500, tau_s = 50, seas.amp = 8, tau_p = 4 / 12, N = 30, phi_c = pi / 2, phi_a = pi / 2, sigma.cal = 0.3, clim.spec.fun.args = list( a = 0.1, b = 1 ) ) var.obj.4 <- crossing( tau_b = 1e03 * 10 / c(5, 50, 500) ) %>% group_by(tau_b) %>% do({ spec.pars.tau_b.1$tau_b <- .$tau_b spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.tau_b.1) data.frame(as.list(IntegrateErrorSpectra(spec.obj)[["proxy.error.var"]])) }) %>% ungroup() %>% filter(smoothed.resolution <= 1000 | smoothed.resolution == Inf) p.tsdv.tau_b <- var.obj.4 %>% PlotTSDVariance(., c( "Calibration.unc.", "Seasonal.bias", "Seasonal.bias.unc." ), include = FALSE) + facet_wrap(~tau_b, labeller = function(labs) { sedproxydocs::label_symbols(labs, label.key = plot.symbol) }) + expand_limits(x = 10^(log10(range(var.obj.4$smoothed.resolution)) + c(-0.1, 0.1))) + labs(x = "Smoothed and interpreted resolution [Years]")
# Change tau_r for each time-scale spec.pars.tau_b.1 <- GetSpecPars("Mg_Ca", T = 1e04, delta_t = 100, tau_r = 100, sig.sq_a = 0, tau_b = 500, tau_s = 50, seas.amp = 8, tau_p = 4 / 12, N = 30, phi_c = pi / 2, phi_a = pi / 2, sigma.cal = 0.3, clim.spec.fun.args = list(a = 0.1, b = 1) ) RoundOdd <- function(x) { x <- round(x) il <- x %% 2 < 1 x[il] <- x[il] + 1 return(x) } var.obj.5 <- crossing( delta_t = seq(100, 1000, 100), tau_b = 1e03 * 10 / c(5, 50, 500) ) %>% mutate( T = 10000, tau_r = delta_t, T = delta_t * RoundOdd(T / delta_t) ) %>% group_by(tau_b, delta_t) %>% do({ # Replace tau_r and delta_t with smoothed resolution spec.pars.tau_b.1$T <- .$T spec.pars.tau_b.1$delta_t <- .$delta_t spec.pars.tau_b.1$tau_r <- .$tau_r spec.pars.tau_b.1$tau_b <- .$tau_b spec.obj <- do.call(ProxyErrorSpectrum, spec.pars.tau_b.1) data.frame(as.list((IntegrateErrorSpectra(spec.obj, method = "sincfilter", tau_smooth = .$tau_r )[["proxy.error.var"]]))) }) %>% ungroup() lab.df <- data.frame( varname = c( "clim.pow.0.5", "clim.pow.scal", "sed.acc.rate", "bio.depth", "amp.seas.cycle", "amp.seas.cycle.2", "n.samples", "resolution", "sigma.meas", "tau_b" ), label = c("gamma", "beta", "SAR", "d", "Alpha", "Alpha[2]", "N", "R", "sigma", "tau[b]") ) p.tsdv.tau_b.tau_r <- var.obj.5 %>% filter(smoothed.resolution != Inf) %>% select( -delta_t, -Climate, -Total.error, -Calibration.unc., -Reference.climate, -Seasonal.bias, -Seasonal.bias.unc. ) %>% gather(component, value, -smoothed.resolution, -tau_b) %>% ungroup() %>% mutate(component = psem:::OrderStages(component, rev = TRUE)) %>% ggplot(aes(x = smoothed.resolution, y = value, fill = component, colour = component)) + # geom_line()+ scale_color_manual("", values = spec.colrs) + geom_area(alpha = 0.75) + scale_x_continuous(trans = ecustools::reverselog_trans(10)) + scale_color_manual("", values = spec.colrs) + scale_fill_manual("", values = spec.colrs) + labs(y = expression("Error variance [°C"^2 ~ "]"), x = "Physical sampling and interpreted resolution [Years]") + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + facet_wrap(~tau_b, labeller = function(labs) { sedproxydocs::label_symbols(labs, label.key = lab.df) }) + expand_limits(x = 10^(log10(range(var.obj.5$smoothed.resolution)) + c(-0.1, 0.1)))
egg::ggarrange(p.tsdv.tau_b, p.tsdv.tau_b.tau_r, ncol = 1, labels = c(letters[1:2]), label.args = list(gp = grid::gpar(font = 1, cex = 1.2)))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.