In EarthSystemDiagnostics/spectraluncertainty:

knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE, cache = TRUE,
  fig.width = 7, fig.height = 4.5, fig.pos = "H", dpi = 300, autodep = TRUE
)

# detach("package:spectraluncertainty", unload = TRUE)
# detach("package:sedproxy", unload = TRUE)
library(spectraluncertainty)
library(sedproxy)
library(tidyverse)

GetClimProxyCorr <- function(pfm){
  with(pfm$simulated.proxy, cor.test(simulated.proxy,
                                     clim.timepoints.ssr)
       )
}

GetReplicateCorr <- function(pfm){

  n.reps <- length(unique(pfm$everything$replicate))

  dat <- pfm$everything %>% 
    filter(stage == "simulated.proxy") %>% 
    select(replicate, timepoints, value) %>% 
    pivot_wider(., names_from = replicate, values_from = value)

  cor.mat <- cor(dat[, 2:(n.reps+1)])
  as.vector(cor.mat[upper.tri(cor.mat)])
}


spec.pars <- GetSpecPars("Mg_Ca", tau_b = 1000* 10/10, T = 1e04, delta_t = 100)
spec.obj <- do.call(ProxyErrorSpectrum, spec.pars)

CompareSimExpCorr <- function(spec.pars){

  spec.obj <- do.call(ProxyErrorSpectrum, spec.pars)
  sed.pars <- SedproxyParsFromSpecPars(spec.pars)
  sed.pars$n.replicates <- 2
  sed.pars$plot.sig.res <- spec.pars$delta_t

  PFM <- do.call(ClimToProxyClim, sed.pars)

  proxy.rep.corr <- mean(GetReplicateCorr(PFM))
  proxy.clim.corr <- as.numeric(GetClimProxyCorr(PFM)$estimate)

  exp.corr <- ExpectedCorrelation(spec.obj) %>% 
    filter(smoothed.resolution == spec.pars$delta_t) %>% 
    spread(correlation.type, rho)

  out <- c(proxy.rep.corr=proxy.rep.corr,
    exp.proxy.rep.corr=exp.corr$`proxy-proxy`,
    proxy.clim.corr=proxy.clim.corr,
    exp.proxy.clim.corr=exp.corr$`proxy-clim`)

  return(out)
}

spec.pars <- GetSpecPars("Mg_Ca", tau_b = 1000 * 10/5, T = 2e04, delta_t = 50,
                         sigma.meas = 0.25, sigma.ind = 2, sig.sq_a = 0)

spec.obj <- do.call(ProxyErrorSpectrum, spec.pars)


out <- t(replicate(10, {
  #suppressMessages(CompareSimExpCorr(spec.pars))
  CompareSimExpCorr(spec.pars)
  }))%>% 
  as_tibble() #%>% 
  #mutate(proxy.clim.corr.sq = proxy.clim.corr^2)

#boxplot(out)

out.long <- out %>% 
  select(proxy.rep.corr, proxy.clim.corr) %>% 
  gather() %>% 
  mutate(expected.value = ifelse(key == "proxy.rep.corr",
                                 unique(out$exp.proxy.rep.corr),
                                 unique(out$exp.proxy.clim.corr)))

exp.values <- out.long %>% 
  select(-value) %>% 
  distinct()

out.long %>% 
  ggplot(aes(x = key, y = value, colour = key)) +
  geom_violin() +
  geom_point(position = position_jitter(width = 0.125)) +
  geom_point(data = exp.values, aes(x = key, y = expected.value),
             size = 4)


out.long %>% 
  ggplot(aes(x = value, fill = key)) +
  geom_histogram(bins = 10) +
  facet_wrap(~key) +
  geom_vline(data = exp.values,
             aes(xintercept = expected.value), colour = "Red")


out  %>% 
  summarise(rep.corr.bias = mean(proxy.rep.corr) - mean(exp.proxy.rep.corr),
            rep.corr.sd = sd(proxy.rep.corr),
            clim.corr.bias = mean(proxy.clim.corr) - mean(exp.proxy.clim.corr),
            clim.corr.sd = sd(proxy.clim.corr))


t.test(out$proxy.rep.corr - (out$exp.proxy.rep.corr))
t.test(out$proxy.clim.corr - out$exp.proxy.clim.corr)


spec.pars <- GetSpecPars("Mg_Ca", tau_b = 1000 * 10/10, T = 1e04, delta_t = 100,
                         sigma.meas = 0.25, sigma.ind = 2, sig.sq_a = 0)

spec.obj <- do.call(ProxyErrorSpectrum, spec.pars)



ExpectedCorrelation(spec.obj) %>% 
  filter(n > 5) %>% 
  #gather(correlation, rho, -smoothed.resolution) %>% 
  ggplot(aes(x = smoothed.resolution, y = rho, colour = correlation.type)) +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = correlation.type),
              alpha = 0.25, colour = NA) +
  geom_line() +
  expand_limits(y = c(0, 1))

EarthSystemDiagnostics/spectraluncertainty documentation built on Jan. 8, 2020, 9:43 a.m.