Trends in all univariate local and regional proxy time series (that is, signal separated from stochastic noise) were recovered by following the modeling protocol for paleoecological time‑series (Simpson 2018). Due to the size of the dataset, the ‘bam’ implementation of the recommended generalized additive models (GAM) was applied using the package ‘mgcv’ 1.8-22 (Wood 2011) for the R 3.6-1 statistical programming language (R Core Team 2019). The auto-regressive option accounted for the temporal autocorrelation in the proxies time-series. Due to the multivariate nature of the grain size composition, a redundancy analysis (RDA) was performed with polynomials of time as predictors of isometric logarithmic ratio transformed values by using the R package ‘composition’ 1.40.1 (Boogaart et al. 2014) and ‘vegan’ 2.4.4 (Oksanen et al. 2017). The first derivative of the smoothing splines (GAM), the polynomials (RDA) and the linearly interpolated global proxies informed on the change of the climate and environmental parameters per 1000 years.

# Load required packages
library(palinsol)
library(mgcv)
library(itsadug)
library(compositions)
library(robCompositions)
library(vegan)

# Auxiliary functions
source("AuxiliaryFunctions/finitDeriv.R")
source("AuxiliaryFunctions/plotCovariate.R")
MaxTime <- 1366
TimeSeq <- seq(0, MaxTime+1, length = 2000+1)

# Add some time before 1.36 million years for PyRateContinuous
# to include the number of initial species at that time
# All covariate will be 0 before 1.36 million years
# and do not exert any influence on 
# speciation or extinction
AddTime <- seq(1367, 5000, length.out = 100)

# Glacial/interglacials background for plots
LR04 <- read.table("Data/LR04_MISboundaries.txt", header = TRUE, check.names = FALSE, skip = 1)
LR04 <- LR04[-c(5:9), ] # Omit MIS5
LR04 <- LR04[-c(2:3), ]
LR04 <- LR04[-c(76:nrow(LR04)), ]
LR04$Boundary <- as.character(LR04$Boundary)
LR04$Boundary <- as.numeric(unlist(lapply(strsplit(LR04$Boundary, "/"), function(x) x[2])))
LR04[, 2] <- LR04[, 2] / 1000

Medstack

Medstack δ18O planktonic isotope ratios in parts per thousand relative to VPDB

# Read data
Medstack <- read.table("Data/MEDSTACK_Wang2010.txt", sep = "\t", header = TRUE)

Medstack2 <- data.frame(Age = TimeSeq, 
                        d18O = approx(x = Medstack$Age, y = Medstack$d18Ostack, xout = TimeSeq)$y)
# Change per 1000 years:
Medstack2Change <- (Medstack2$d18O[1:(nrow(Medstack2)-1)] - Medstack2$d18O[2:nrow(Medstack2)]) / 1 
Medstack2AbsChange <- abs(Medstack2Change)
# Remove times that are only there for calculating the change in temperature:
Medstack2 <- Medstack2[-nrow(Medstack2), ]

plotCovariate(Medstack2$Age, Medstack2$d18O, LR04, 
              ylab = expression(paste(delta^{18}, O[planktonic], "(VPDB \u2030)")),
              col = "green4", ylim = c(2.8, -2.2))

# Write covariate for PyRateContinuous
WriteMed <- data.frame(time = c(Medstack2$Age, AddTime) / 100,
                       Medstack = c(scale(Medstack2$d18O), rep(0, length(AddTime))),
                       Medstackabschange = c(scale(Medstack2Change), rep(0, length(AddTime))))
write.table(WriteMed[, c(1, 2)], 
            file = "Results/Covariates/Medstack.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteMed[, c(1, 3)], 
            file = "Results/Covariates/Medstackabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

LR04

LR04 benthic δ18O stack isotope ratios in parts per thousand relative to VPDB

LR04Values <- read.table("Data/LR04stack.txt", header = TRUE, skip = 4, sep = "\t")
colnames(LR04Values) <- c("Time", "LR04", "SE")

LR04Values2 <- data.frame(time = TimeSeq,
                          LR04 = approx(x = LR04Values$Time, y = LR04Values$LR04, xout = TimeSeq)$y)
# Change per 1000 years:
LR04Change <- LR04Values2$LR04[1:(nrow(LR04Values2)-1)] - LR04Values2$LR04[2:nrow(LR04Values2)] 
LR04AbsChange <- abs(LR04Change)
# Remove times that are only there for calculating the change in temperature:
LR04Values2 <- LR04Values2[-nrow(LR04Values2), ]

plotCovariate(LR04Values2$time, LR04Values2$LR04, LR04, 
              ylab = expression(paste(delta^{18}, O[benthic], "(VPDB \u2030)")),
              col = "green4", ylim = c(5.1, 3.1))

# Write covariate for PyRateContinuous
WriteLR04 <- data.frame(time = c(LR04Values2$time, AddTime) / 100,
                        LR04 = c(scale(LR04Values2$LR04), rep(0, length(AddTime))),
                        LR04abschange = c(scale(LR04Change), rep(0, length(AddTime))))
write.table(WriteLR04[, c(1, 2)], 
            file = "Results/Covariates/LR04.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteLR04[, c(1, 3)], 
            file = "Results/Covariates/LR04abschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Insolation

Northern Hemisphere summer insolation at the latitude of Lake Ohrid (41° N)

data(BER78) 
SplineT <- seq(0, MaxTime, length = 2000)
Times <- -seq(-1, MaxTime + 1, length = 2000 + 2) * 1000
Obl <- function(t) {ber78 = ber78(t)}
Obls <- data.frame( t(sapply(Times, Obl)) )
Insolation41 <- rep(NA, length(Times))
for(i in 1:nrow(Obls)){
  Insolation41[i] <- Insol_d1d2(
    orbit = c(eps = as.numeric(Obls[i, 1]), 
              ecc = as.numeric(Obls[i, 2]), 
              varpi = as.numeric(Obls[i, 3])), 
    d1 = 171, d2 = 172, lat = 41*pi/180, S0 = 1368, avg = TRUE)
}
# Calculate change
H <- 1 # Change per kilo year
Insolation41Change <- (Insolation41[1:(length(Insolation41)-2)] - 
                         Insolation41[3:(length(Insolation41))]) / 2*H
Insolation41AbsChange <- abs(Insolation41Change)
# Remove times that are only there for calculating the change in insolation:
Insolation41 <- Insolation41[2:(length(Insolation41)-1)] 

plotCovariate(SplineT, Insolation41, LR04,
              ylab = expression(paste("Summer insolation ", "(W ", m^2,")")),
              col = "red")

# Write covariate for PyRateContinuous
WriteIns41 <- data.frame(time = c(SplineT, AddTime) / 100,
                         Ins41 = c(scale(Insolation41), rep(0, length(AddTime))),
                         Ins41abschange = c(scale(Insolation41Change), rep(0, length(AddTime))))
write.table(WriteIns41[, c(1, 2)], 
            file = "Results/Covariates/Ins41.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteIns41[, c(1, 3)], 
            file = "Results/Covariates/Ins41abschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Arboreal pollen

Percentages of arboreal pollen excluding Pinus pollen at Lake Ohrid

ArbPol <- read.table("Data/ArborealPollen.txt", header = TRUE, sep = "\t")
ArbPol2 <- ArbPol[nrow(ArbPol):1, ]
ArbPol2[, 1] <- -ArbPol2[, 1]
# Mark beginn for temporal autocorrelation
ArbPol2$Start <- c(TRUE, rep(FALSE, nrow(ArbPol2)-1))

# Fit GAM
APBam1 <- bam(ArborealPollen ~ s(Age, k = 500), data = ArbPol2)
Rho <- acf(resid_gam(APBam1), plot = FALSE)$acf[2]
APBam2 <- bam(ArborealPollen ~ s(Age, k = 500), data = ArbPol2, 
              rho = Rho, AR.start = ArbPol2$Start)
NewData <- data.frame(Age = -seq(0, MaxTime, length = 2000)) 
APPredict <- predict(APBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
APDeriv <- finitDeriv(Gam = APBam2, NewData, H = 0.001, Scale = FALSE)
APDeriv <- abs(APDeriv)

plotCovariate(-NewData$Age, APPredict, LR04,
              ylab = "Arboreal pollen (%)",
              col = "red", ylim = c(0, 100), Raw = ArbPol)

# Write covariate for PyRateContinuous
WriteArbPol <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          ArbPol = c(scale(APPredict), rep(0, length(AddTime))),
                          ArbPolabschange = c(scale(APDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteArbPol[, c(1, 2)], 
            file = "Results/Covariates/ArbPol.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteArbPol[, c(1, 3)], 
            file = "Results/Covariates/ArbPolabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Deciduous oaks

Percentage of pollen from deciduous oaks at Lake Ohrid

DecOaks <- read.table("Data/DeciduousOaks.txt", header = TRUE, sep = "\t")
DecOaks2 <- DecOaks[nrow(DecOaks):1, ]
DecOaks2[, 1] <- -DecOaks2[, 1]
# Mark beginn for temporal autocorrelation
DecOaks2$Start <- c(TRUE, rep(FALSE, nrow(DecOaks2)-1))

# Fit GAM
DecOaksBam1 <- bam(DecOak ~ s(Age, k = 500), data = DecOaks2)
Rho <- acf(resid_gam(DecOaksBam1), plot = FALSE)$acf[2]
DecOaksBam2 <- bam(DecOak ~ s(Age, k = 500), data = DecOaks2, 
                   rho = Rho, AR.start = DecOaks2$Start)
DecOaksPredict <- predict(DecOaksBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
DecOaksDeriv <- finitDeriv(Gam = DecOaksBam2, NewData, H = 0.001, Scale = FALSE)
DecOaksDeriv <- abs(DecOaksDeriv)

plotCovariate(-NewData$Age, DecOaksPredict, LR04,
              ylab = "Deciduous oaks (%)",
              col = "red", ylim = c(0, 60), Raw = DecOaks)

# Write covariate for PyRateContinuous
WriteDecOaks <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          DecOaks = c(scale(DecOaksPredict), rep(0, length(AddTime))),
                          DecOaksabschange = c(scale(DecOaksDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteDecOaks[, c(1, 2)], 
            file = "Results/Covariates/DecOaks.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteDecOaks[, c(1, 3)], 
            file = "Results/Covariates/DecOaksabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Lake Ohrid isotope

Lake Ohrid δ18Olakewater based on calcite and siderite δ18O information

Isotope <- read.table("Data/OxygenIsotope.txt", header = TRUE, sep = "\t")
Isotope2 <- Isotope[nrow(Isotope):1, ]
Isotope2$Age <- -Isotope2$Age
# Mark beginn for temporal autocorrelation
Isotope2$Start <- c(TRUE, rep(FALSE, nrow(Isotope2)-1))

# Fit GAM
IsotopeBam1 <- bam(Lakewater ~ s(Age, k = 500), data = Isotope2)
Rho <- acf(resid_gam(IsotopeBam1), plot = FALSE)$acf[2]
IsotopeBam2 <- bam(Lakewater ~ s(Age, k = 500), data = Isotope2, 
                   rho = Rho, AR.start = Isotope2$Start)
IsotopePredict <- predict(IsotopeBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
IsotopeDeriv <- finitDeriv(Gam = IsotopeBam2, NewData, H = 0.001, Scale = FALSE)
IsotopeDeriv <- abs(IsotopeDeriv)

plotCovariate(-NewData$Age, IsotopePredict, LR04,
              ylab = expression(paste(delta^{18}, "O"[lakewater], "(\u2030)")),
              col = "blue", ylim = c(-10, -2), 
              Raw = Isotope[, c("Age", "Lakewater")])

# Write covariate for PyRateContinuous
WriteIsotope <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          Isotope = c(scale(IsotopePredict), rep(0, length(AddTime))),
                          Isotopeabschange = c(scale(IsotopeDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteIsotope[, c(1, 2)], 
            file = "Results/Covariates/Isotope.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteIsotope[, c(1, 3)], 
            file = "Results/Covariates/Isotopeabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

TIC

Lake Ohrid total inorganic carbon (TIC) concentrations

Tic <- read.table("Data/TIC.txt", header = TRUE, sep = "\t")
Tic2 <- Tic[nrow(Tic):1, ]
Tic2[, 1] <- -Tic2[, 1]
# Mark beginn for temporal autocorrelation
Tic2$Start <- c(TRUE, rep(FALSE, nrow(Tic2)-1))

# Fit GAM
TicBam1 <- bam(TIC ~ s(Age, k = 500), data = Tic2)
Rho <- acf(resid_gam(TicBam1), plot = FALSE)$acf[2]
TicBam2 <- bam(TIC ~ s(Age, k = 500), data = Tic2, 
               rho = Rho, AR.start = Tic2$Start)
TicPredict <- predict(TicBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
TicDeriv <- finitDeriv(Gam = TicBam2, NewData, H = 0.001, Scale = FALSE)
TicDeriv <- abs(TicDeriv)

plotCovariate(-NewData$Age, TicPredict, LR04,
              ylab = "TIC (%wt)",
              col = "blue", ylim = c(0, 12), Raw = Tic)

# Write covariate for PyRateContinuous
WriteTic <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          TIC = c(scale(TicPredict), rep(0, length(AddTime))),
                          TICabschange = c(scale(TicDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteTic[, c(1, 2)], 
            file = "Results/Covariates/TIC.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteTic[, c(1, 3)], 
            file = "Results/Covariates/TICabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

TOC

Lake Ohrid total organic carbon (TOC) concentrations

Toc <- read.table("Data/TOC.txt", header = TRUE, sep = "\t")
Toc2 <- Toc[nrow(Toc):1, ]
Toc2[, 1] <- -Toc2[, 1]
# Mark beginn for temporal autocorrelation
Toc2$Start <- c(TRUE, rep(FALSE, nrow(Toc2)-1))

# Fit GAM
TocBam1 <- bam(TOC ~ s(Age, k = 500), data = Toc2)
Rho <- acf(resid_gam(TocBam1), plot = FALSE)$acf[2]
TocBam2 <- bam(TOC ~ s(Age, k = 500), data = Toc2, 
               rho = Rho, AR.start = Toc2$Start)
TocPredict <- predict(TocBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
TocDeriv <- finitDeriv(Gam = TocBam2, NewData, H = 0.001, Scale = FALSE)
TocDeriv <- abs(TocDeriv)

plotCovariate(-NewData$Age, TocPredict, LR04,
              ylab = "TOC (%wt)",
              col = "blue", ylim = c(0, 4), Raw = Toc)

# Write covariate for PyRateContinuous
WriteToc <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          TOC = c(scale(TocPredict), rep(0, length(AddTime))),
                          TOCabschange = c(scale(TocDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteToc[, c(1, 2)], 
            file = "Results/Covariates/TOC.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteToc[, c(1, 3)], 
            file = "Results/Covariates/TOCabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Potassium

Lake Ohrid potassium (K) counts from XRF scanning

K <- read.table("Data/Potassium.txt", header = TRUE, sep = "\t")
K2 <- K[nrow(K):1, ]
K2$Age <- -K2$Age
# Mark beginn for temporal autocorrelation
K2$Start <- c(TRUE, rep(FALSE, nrow(K2)-1))
# 11 point runing mean as in Wagner et al., 2019
K$K11pt <- rep(NA, nrow(K2))
K$K11pt[6:(nrow(K)-6)] <- sapply(6:(nrow(K)-6), function(x) mean(K[(x-5):(x+5), "K"]))
K$K11pt <- K$K11pt / 1000

KBam1 <- bam(K ~ s(Age, k = 500), data = K2, samfrac = 0.1) 
Rho <- acf(resid_gam(KBam1), plot = FALSE)$acf[2]
KBam2 <- bam(K ~ s(Age, k = 500), data = K2, samfrac = 0.1, 
             rho = Rho, AR.start = K2$Start)
KPredict <- predict.bam(KBam2, newdata = NewData, type = "response", se.fit = FALSE)
KPredict <- KPredict / 1000 # Convert to kilocounts

# 1st derivative for change in covariate
KDeriv <- finitDeriv(Gam = KBam2, NewData, H = 0.001, Scale = FALSE)
KDeriv <- abs(KDeriv)

plotCovariate(-NewData$Age, KPredict, LR04,
              ylab = expression(paste(K["11pt"], "(kcts)")),
              col = "blue", ylim = c(10, 0), 
              Raw = K[, c("Age", "K11pt")])

# Write covariate for PyRateContinuous
WriteK <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          K = c(scale(KPredict), rep(0, length(AddTime))),
                          Kabschange = c(scale(KDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteK[, c(1, 2)], 
            file = "Results/Covariates/K.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteK[, c(1, 3)], 
            file = "Results/Covariates/Kabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
rm("K", "K2", "KBam1", "KBam2") # Safe memory

Calcium

Lake Ohrid calcium (Ca) counts from XRF scanning

Ca <- read.table("Data/Calcium.txt", header = TRUE, sep = "\t")
Ca2 <- Ca[nrow(Ca):1, ]
Ca2$Age <- -Ca2$Age
# MarCa beginn for temporal autocorrelation
Ca2$Start <- c(TRUE, rep(FALSE, nrow(Ca2)-1))
# 11 point runing mean as in Wagner et al., 2019
Ca$Ca11pt <- rep(NA, nrow(Ca2))
Ca$Ca11pt[6:(nrow(Ca)-6)] <- sapply(6:(nrow(Ca)-6), function(x) mean(Ca[(x-5):(x+5), "Ca"]))
Ca$Ca11pt <- Ca$Ca11pt / 1000

CaBam1 <- bam(Ca ~ s(Age, k = 500), data = Ca2, samfrac = 0.1) 
Rho <- acf(resid_gam(CaBam1), plot = FALSE)$acf[2]
CaBam2 <- bam(Ca ~ s(Age, k = 500), data = Ca2, samfrac = 0.1, 
             rho = Rho, AR.start = Ca2$Start)
CaPredict <- predict.bam(CaBam2, newdata = NewData, type = "response", se.fit = FALSE)
CaPredict <- CaPredict / 1000 # Convert to Kilocounts

# 1st derivative for change in covariate
CaDeriv <- finitDeriv(Gam = CaBam2, NewData, H = 0.001, Scale = FALSE)
CaDeriv <- abs(CaDeriv)

plotCovariate(-NewData$Age, CaPredict, LR04,
              ylab = expression(paste(Ca["11pt"], "(kcts)")),
              col = "blue", ylim = c(0, 200), 
              Raw = Ca[, c("Age", "Ca11pt")])

# Write covariate for PyRateContinuous
WriteCa <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          Ca = c(scale(CaPredict), rep(0, length(AddTime))),
                          Caabschange = c(scale(CaDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteCa[, c(1, 2)], 
            file = "Results/Covariates/Ca.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteCa[, c(1, 3)], 
            file = "Results/Covariates/Caabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
rm("Ca", "Ca2", "CaBam1", "CaBam2") # Safe memory

Quartz

Lake Ohrid relative sedimentary quartz content

Quartz <- read.table("Data/Quartz.txt", header = TRUE, sep = "\t")
Quartz2 <- Quartz[nrow(Quartz):1, ]
Quartz2[, 1] <- -Quartz2[, 1]
# Mark beginn for temporal autocorrelation
Quartz2$Start <- c(TRUE, rep(FALSE, nrow(Quartz2)-1))

# Fit GAM
QuartzBam1 <- bam(Quartz ~ s(Age, k = 500), data = Quartz2)
Rho <- acf(resid_gam(QuartzBam1), plot = FALSE)$acf[2]
QuartzBam2 <- bam(Quartz ~ s(Age, k = 500), data = Quartz2, 
                  rho = Rho, AR.start = Quartz2$Start)
QuartzPredict <- predict(QuartzBam2, newdata = NewData, type = "response", se.fit = FALSE)

# 1st derivative for change in covariate
QuartzDeriv <- finitDeriv(Gam = QuartzBam2, NewData, H = 0.001, Scale = FALSE)
QuartzDeriv <- abs(QuartzDeriv)

plotCovariate(-NewData$Age, QuartzPredict, LR04,
              ylab = "Quartz (peak area)",
              col = "blue", ylim = c(0.4, 0), Raw = Quartz)

# Write covariate for PyRateContinuous
WriteQuartz <- data.frame(time = c(-NewData$Age, AddTime) / 100,
                          Quartz = c(scale(QuartzPredict), rep(0, length(AddTime))),
                          Quartzabschange = c(scale(QuartzDeriv[,2]), rep(0, length(AddTime))))
write.table(WriteQuartz[, c(1, 2)], 
            file = "Results/Covariates/Quartz.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteQuartz[, c(1, 3)], 
            file = "Results/Covariates/Quartzabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Grain size

Lake Ohrid grain size information (summarized by ordination)

Grainsize <- read.table("Data/Grainsize.txt", header = TRUE, sep = "\t")
# Combine all types of sand into one category
Grainsize$SAND <- Grainsize$VERY_COARSE_SAND + Grainsize$COARSE_SAND + Grainsize$MEDIUM_SAND + Grainsize$FINE_SAND + Grainsize$VERY_FINE_SAND
GrainsizeRobPca <- pcaCoDa(Grainsize[,c("SAND","VERY_COARSE_SILT", "COARSE_SILT", "MEDIUM_SILT", 
                                        "FINE_SILT", "VERY_FINE_SILT", "CLAY")]) 

# Exclude all samples with very little sand
GrainsizeRobPca <- GrainsizeRobPca$scores
Grainsize2 <- Grainsize[-which(GrainsizeRobPca[, 1] < -1.5),] 
GrainsizeRobPca <- pcaCoDa(Grainsize2[,c("SAND","VERY_COARSE_SILT", "COARSE_SILT", "MEDIUM_SILT", 
                                         "FINE_SILT", "VERY_FINE_SILT", "CLAY")]) 
# 1st axis explains 78% of the grainsize variance
GrainsizeRobPca$princompOutputClr$sdev^2 / sum(GrainsizeRobPca$princompOutputClr$sdev^2) 
##       Comp.1       Comp.2       Comp.3       Comp.4       Comp.5 
## 0.7838074299 0.1897189122 0.0203257678 0.0051105364 0.0009224254 
##       Comp.6 
## 0.0001149283
GrainsizeRobPca <- GrainsizeRobPca$scores

# Suggestion Gavin Simpson: RDA
GrainsizeIlr <- ilr(Grainsize2[,c("SAND","VERY_COARSE_SILT", "COARSE_SILT", "MEDIUM_SILT", 
                                  "FINE_SILT", "VERY_FINE_SILT", "CLAY")])
# RDA with polynomials of time as predictors of 
# isometric logarithmic ratio transformed grainsize values
Poly <- poly(Grainsize2$Age, 25) 
Poly2 <- as.data.frame(Poly) # Keep the poly class for the prediction
colnames(Poly2) <- paste0("V", 1:ncol(Poly2)) # rda fails without colnames
GrainsizeRda0 <- rda(GrainsizeIlr ~ 1, data = Poly2) # Model with intercept only
GrainsizeRda1 <- rda(GrainsizeIlr ~ ., data = Poly2) 
GrainsizeRdaSel <- ordiR2step(GrainsizeRda0, scope = formula(GrainsizeRda1), 
                              direction = "forward", Pin = 0.05,
                              R2permutations= 999, trace = FALSE)

NewData <- predict(Poly, SplineT)
NewData <- as.data.frame(NewData)
colnames(NewData) <- colnames(Poly2)
GrainsizePred <- -predict(GrainsizeRdaSel, newdata = NewData, type = "lc")[, 1]
# Change in grainsize
H <- 1/1000 # 1ka/1000
NewData <- predict(Poly, SplineT + H)
NewData <- as.data.frame(NewData)
colnames(NewData) <- colnames(Poly2)
PlusH <- -predict(GrainsizeRdaSel, newdata = NewData, type = "lc")[, 1]
NewData <- predict(Poly, SplineT - H)
NewData <- as.data.frame(NewData)
colnames(NewData) <- colnames(Poly2)
MinusH <- -predict(GrainsizeRdaSel, newdata = NewData, type = "lc")[, 1]
GrainsizeDeriv <- (PlusH - MinusH) / 2*H
GrainsizePred <- scale(GrainsizePred)

plotCovariate(SplineT, GrainsizePred, LR04,
              ylab = "Grain size",
              col = "blue", ylim = c(-2, 6), 
              Raw = cbind(Grainsize2$Age, scale(GrainsizeRobPca[, 1])))

WriteGrainsize <- data.frame(time = c(SplineT, AddTime) / 100,
                             Grainsize = c(GrainsizePred, rep(0, length(AddTime))),
                             Grainsizeabschange = c(scale(GrainsizeDeriv), rep(0, length(AddTime))))
write.table(WriteGrainsize[, c(1, 2)], 
            file = "Results/Covariates/Grainsize.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)
write.table(WriteGrainsize[, c(1, 3)], 
            file = "Results/Covariates/Grainsizeabschange.txt", 
            sep = "\t", row.names = FALSE, quote = FALSE)

Plattform and resources

Sys.time()
## [1] "2020-01-05 11:27:10 CET"
sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] vegan_2.5-5           lattice_0.20-38       permute_0.9-4        
##  [4] robCompositions_2.1.0 pls_2.6-0             e1071_1.7-2          
##  [7] data.table_1.12.2     ggplot2_3.2.1         compositions_1.40-2  
## [10] bayesm_3.1-3          energy_1.7-6          robustbase_0.93-5    
## [13] tensorA_0.36.1        itsadug_2.3           plotfunctions_1.3    
## [16] mgcv_1.8-28           nlme_3.1-141          palinsol_0.93        
## [19] gsl_2.1-6             knitr_1.25           
## 
## loaded via a namespace (and not attached):
##  [1] pbkrtest_0.4-7        RColorBrewer_1.1-2    prabclus_2.2-6       
##  [4] tools_3.6.2           R6_2.4.0              lazyeval_0.2.2       
##  [7] colorspace_1.4-1      trimcluster_0.1-2     nnet_7.3-12          
## [10] NADA_1.6-1            withr_2.1.2           sp_1.3-1             
## [13] tidyselect_0.2.5      GGally_1.3.1          compiler_3.6.2       
## [16] quantreg_5.51         SparseM_1.77          sROC_0.1-2           
## [19] diptest_0.75-7        scales_1.0.0          DEoptimR_1.0-8       
## [22] lmtest_0.9-37         mvtnorm_1.0-11        stringr_1.4.0        
## [25] digest_0.6.21         minqa_1.2.4           rmarkdown_1.15       
## [28] rrcov_1.4-7           pkgconfig_2.0.3       htmltools_0.4.0      
## [31] lme4_1.1-21           rlang_0.4.0           rstudioapi_0.10      
## [34] zCompositions_1.3.2-1 zoo_1.8-6             mclust_5.4.5         
## [37] dplyr_0.8.3           car_2.1-5             magrittr_1.5         
## [40] modeltools_0.2-21     Matrix_1.2-17         Rcpp_1.0.2           
## [43] munsell_0.5.0         stringi_1.4.3         yaml_2.2.0           
## [46] cvTools_0.3.2         MASS_7.3-51.4         flexmix_2.3-14       
## [49] plyr_1.8.4            grid_3.6.2            parallel_3.6.2       
## [52] crayon_1.3.4          splines_3.6.2         pillar_1.4.2         
## [55] ranger_0.11.2         boot_1.3-23           fpc_2.1-10           
## [58] codetools_0.2-16      stats4_3.6.2          glue_1.3.1           
## [61] evaluate_0.14         laeken_0.4.6          vcd_1.4-3            
## [64] nloptr_1.2.1          MatrixModels_0.4-1    VIM_4.8.0            
## [67] gtable_0.3.0          purrr_0.3.2           tidyr_0.8.3          
## [70] reshape_0.8.8         kernlab_0.9-27        assertthat_0.2.1     
## [73] xfun_0.10             class_7.3-15          survival_2.44-1.1    
## [76] pcaPP_1.9-73          truncnorm_1.0-8       tibble_2.1.3         
## [79] cluster_2.1.0