Nothing
## ----setup, warning=FALSE, message=FALSE, echo=FALSE--------------------------
#checking if required packages are installed, and installing them if not
list.of.packages <- c("ggplot2", "cowplot", "knitr", "viridis", "tidyr", "formatR", "grid", "devtools", "magrittr", "kableExtra", "viridis")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, dep=TRUE)
#install virtualPollen if not installed
if(!("virtualPollen" %in% installed.packages())){
library(devtools)
install_github("blasbenito/virtualPollen")
}
# source("ecological_memory_functions.R")
library(virtualPollen) #to simulate pollen curves
library(ggplot2) #plotting library
library(cowplot) #plotting library
library(knitr) #report generation in R
library(viridis) #pretty plotting colors
library(grid)
library(tidyr)
library(formatR)
library(kableExtra) #to fit tables to pdf page size
library(magrittr) #kableExtra requires pipes
options(scipen=999)
#trying to line-wrap code in pdf output
#from https://github.com/yihui/knitr-examples/blob/master/077-wrap-output.Rmd
knitr::opts_chunk$set(echo = TRUE, fig.pos= "h")
opts_chunk$set(tidy.opts=list(width.cutoff=80), tidy=FALSE)
## ---- size="small"------------------------------------------------------------
#sets a state for the generator of pseudo-random numbers
set.seed(1)
#defines the variable "time"
time <- 1:10000
#samples (-1, 0, 1) with replacement
moves <- sample(x=c(-1, 0, 1), size=length(time), replace=TRUE)
#computes the cumulative sum of moves
random.walk <- cumsum(moves)
## ---- echo=FALSE , fig.height=2.5, fig.width=9, fig.cap="Random walk (a) and its temporal autocorrelation (b)."----
p1 <- ggplot(data=data.frame(Time=time, Value=random.walk), aes(x=Time, y=Value)) +
geom_line(color="gray40") +
ggtitle("Random walk") +
xlab("Time (years)")
p2 <- ggplot(data=acfToDf(random.walk, 5000, 50), aes(x=lag, y=acf)) +
geom_hline(aes(yintercept = 0)) +
geom_hline(aes(yintercept = ci.max), color="red", linetype="dashed") +
geom_hline(aes(yintercept = ci.min), color="red", linetype="dashed") +
geom_segment(mapping = aes(xend = lag, yend = 0)) +
ggtitle("Temporal autocorrelation") +
xlab("Lag (years)") +
ylab("Pearson correlation")
plot_grid(p1, p2, labels = c("a", "b"), align = "h")
## ---- size="small"------------------------------------------------------------
#setting a fixed seed for the generator of pseudo-random numbers
set.seed(1)
#generating 5 random numbers in the range [0, 1]
a <- runif(5)
#applying a convolution filter of length 2 and value 1
b <- filter(a, filter=c(1,1), method="convolution", circular=TRUE)
## ---- echo=FALSE , fig.height=2.5, fig.width=9, fig.cap="Original sequence (dashed line) and filtered sequence with filter (solid line)."----
ggplot(data=data.frame(Time=rep(1:length(a), 2), Value=c(a, b), Legend=c(rep("Original (a)", length(a)), rep("Filtered (b)", length(b)))), aes(x=Time, y=Value, group=Legend)) +
geom_line(aes(linetype=Legend), color="gray40", size=1) +
theme(legend.position="right") +
ggtitle("Effect of a convolution filter")
## ---- results="asis", echo=FALSE----------------------------------------------
temp.table <- data.frame(row=1:5, a=round(a, 2), b=round(b, 2))
temp.table$operation <- c("b1 = a1 x f2 + a2 x f1", "b2 = a2 x f2 + a3 x f1", "b3 = a3 x f2 + a4 x f1" , "b4 = a4 x f2 + a5 x f1", "b5 = a5 x f2 + a6 x f1")
kable(temp.table, caption = "Original sequence (a), filtered sequence (b), and filtering operations. Numbers beside letters a and b represent row numbers, while f1 and f2 represent the values of the convolution filter (both equal to 1 in this case).", booktabs = T) %>% kable_styling(latex_options = c("hold_position", "striped"))
## ---- size="small"------------------------------------------------------------
moves.10 <- filter(moves, filter=rep(1, 10), circular=TRUE)
moves.100 <- filter(moves, filter=rep(1, 100), circular=TRUE)
## ---- fig.height=5, fig.width=9, warning=TRUE, fig.cap="Sequences filtered with different filter lengths: a) Sequence with autocorrelation length equal to 10; b) Temporal autocorrelation of a); c) Sequence with autocorrelation length equal to 100: d) Temporal autocorrelation of c).", echo=FALSE----
p4 <- ggplot(data=data.frame(Time=time, Value=as.vector(moves.10)), aes(x=Time, y=Value)) +
geom_line(color="gray40", size=0.5) +
ggtitle("Time series") +
xlab("") +
theme(axis.line.x=element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank(), plot.margin=unit(c(0,0,0,0), "cm"))
p5 <- ggplot(data=acfToDf(moves.10, 200, 50), aes(x=lag, y=acf)) +
geom_hline(aes(yintercept = 0)) +
geom_hline(aes(yintercept = ci.max), color="red", linetype="dashed") +
geom_hline(aes(yintercept = ci.min), color="red", linetype="dashed") +
geom_segment(mapping = aes(xend = lag, yend = 0)) +
ggtitle("Temporal autocorrelation") +
xlab("") +
ylab("R-squared") +
theme(axis.line.x=element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank(), plot.margin=unit(c(0,0,0,0), "cm"))
p6 <- ggplot(data=data.frame(Time=time, Value=as.vector(moves.100)), aes(x=Time, y=Value)) +
geom_line(color="gray40", size=0.5) +
theme(plot.margin=unit(c(0,0,0,0), "cm"))
p7 <- ggplot(data=acfToDf(moves.100, 200, 50), aes(x=lag, y=acf)) +
geom_hline(aes(yintercept = 0)) +
geom_hline(aes(yintercept = ci.max), color="red", linetype="dashed") +
geom_hline(aes(yintercept = ci.min), color="red", linetype="dashed") +
geom_segment(mapping = aes(xend = lag, yend = 0)) +
xlab("Lag (years)") +
ylab("R-squared") +
theme(plot.margin=unit(c(0,0,0,0), "cm"))
plot_grid(p4, p5, p6, p7, labels = c("a", "b", "c", "d"), align = "v", nrow=2)
## ---- size="small"------------------------------------------------------------
moves.5000 <- filter(moves, filter=rep(1, 5000), circular=TRUE)
## ---- fig.height=2.5, fig.width=9, fig.cap="Sequence moves.5000 (a) and its temporal autocorrelation (b). In this case there is a large deviation between the required temporal autocorrelation (5000) and the outcome (2000).", echo=FALSE----
p8 <- ggplot(data=data.frame(Time=time, Value=as.vector(moves.5000)), aes(x=Time, y=Value)) +
geom_line(color="gray40", size=0.5) +
ggtitle("Time series") +
theme(plot.margin=unit(c(0,0,0,0), "cm"))
p9 <- ggplot(data=acfToDf(moves.5000, 5000, 50), aes(x=lag, y=acf)) +
geom_hline(aes(yintercept = 0)) +
geom_hline(aes(yintercept = ci.max), color="red", linetype="dashed") +
geom_hline(aes(yintercept = ci.min), color="red", linetype="dashed") +
geom_segment(mapping = aes(xend = lag, yend = 0)) +
ggtitle("Temporal autocorrelation") +
theme(plot.margin=unit(c(0,0,0,0), "cm")) +
ylab("R-square")
plot_grid(p8, p9, labels = c("a", "b"), align = "h")
## ---- fig.height=4, fig.width=9, fig.cap="Virtual driver and its temporal autocorrelation. Note that virtual driver B will not be used hereafter to simplify the explanation on the model dynamics."----
driver <- simulateDriverS(
random.seeds=c(60, 120),
time=1:10000,
autocorrelation.lengths = 600,
output.min=c(0, 0),
output.max=c(100, 100),
driver.names=c("A", "B"),
filename=NULL)
## ---- echo=FALSE--------------------------------------------------------------
parameters <- parametersDataframe(rows=1)
parameters[1,] <- c("Species 1", 50, 20, 2, 0.2, 0, 100, 10000, 1, 0, 50, 10, 0, 0, 600, 0)
parameters[, 2:ncol(parameters)] <- sapply(parameters[, 2:ncol(parameters)], as.numeric)
parameters.t <- data.frame(t(parameters))
parameters.t <- data.frame(parameters.t[-1,])
colnames(parameters.t) <- paste(parameters$label, sep="")
parameters.t$Parameter <- rownames(parameters.t)
rownames(parameters.t) <- NULL
parameters.t <- parameters.t[, c(ncol(parameters.t),1:(ncol(parameters.t)-1))]
#removing last two lines
parameters.t <- parameters.t[c(1:8, 10:11), ]
kable(parameters.t, caption="Parameters of a virtual species. Note that driver.B is ommited in order to simplify the explanation of the model.", booktabs = T) %>% kable_styling(latex_options = c("hold_position", "striped"))
## ---- size="small"------------------------------------------------------------
niche.A <- dnorm(x=0:100, mean=50, sd=10)
## ---- fig.height=2.5, fig.width=6, fig.cap="Ecological niche of the virtual species (blue) against the density (relative availability of values over time) of the driver (gray). Both densities have been scaled in the range [0, 1].", echo=FALSE----
#scaling suitability between 0 and 1
niche.A <- niche.A / max(niche.A)
#getting the density of the driver
driver.A = driver[driver$driver=="A", "value"]
density.driver.A = density(driver.A, from=min(driver.A), to=max(driver.A), n=101, bw=max(driver.A)/100)
density.driver.A.y = (density.driver.A$y - min(density.driver.A$y)) / (max(density.driver.A$y) - min(density.driver.A$y))
driver.A.range = seq(min(driver.A), max(driver.A), length.out = 101)
#dataframe for plot
plot.df = data.frame(Species=rep(paste(parameters[1, "label"], sep=""), 101),
Driver=c(rep("Driver A", 101)),
Driver.density.x=c(density.driver.A$x),
Driver.density.y=c(density.driver.A.y),
Value=driver.A.range,
Suitability=niche.A)
ggplot(data=plot.df, aes(x=Value, y=Suitability, group=Species)) +
geom_ribbon(data=plot.df, aes(ymin=0, ymax=Driver.density.y), color="gray80", fill="gray80", alpha=0.5) +
geom_ribbon(data=plot.df, aes(ymin=0, ymax=Suitability), alpha=0.8, colour=NA, fill="#4572A9") +
geom_line(data=plot.df, aes(x=Value, y=Driver.density.y), color="gray80", alpha=0.5) +
xlab("Driver values") +
ylab("Environmental suitability") +
theme(text = element_text(size=12), strip.background = element_rect(fill=NA), panel.spacing = unit(1, "lines"))
## ---- fig.height=5, fig.width=9, fig.cap="Driver and environmental suitability of the virtual taxa. Burn-in period is highlighted by a gray box in the Environmental suitability panel.", echo=FALSE----
#computing density
density.A <- dnorm(driver[driver$driver == "A", "value"], mean=50, sd=10)
#scaling to [0, 1]
suitability <- density.A / max(density.A)
burnin.suitability <- jitter(c(rep(1, parameters$maximum.age*5), seq(1, suitability[1], length.out = parameters$maximum.age*5)), amount=max(suitability)/100)
burnin.suitability[burnin.suitability < 0]<-0
burnin.suitability[burnin.suitability > 1]<-1
length.burnin.suitability<-length(burnin.suitability)
burnin.suitability <- c(burnin.suitability, suitability)
plot.df4 <- data.frame(Time=c(-length.burnin.suitability:-1, 1:(length(suitability))), Suitability=burnin.suitability, Period=c(rep("Burn-in", length.burnin.suitability), rep("Simulation", length(suitability))))
p1 <- ggplot(data=driver[driver$driver == "A", ], aes(x=time, y=value)) +
geom_line(size=0.5, color="gray40") +
ggtitle("Driver and environmental suitability") +
xlab("") +
ylab("Driver") +
coord_cartesian(xlim=c(-500, 10000)) +
theme(plot.margin = unit(c(0,0,-0.5,0), "cm"), axis.line.x=element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank())
p2 <- ggplot(data=plot.df4, aes(x=Time, y=Suitability)) +
geom_rect(aes(xmin=min(plot.df4$Time), xmax=0, ymin=0, ymax=Inf), inherit.aes=FALSE, fill="gray90") +
geom_line(size=0.5, color="#4572A9") +
xlab("Time (years)") +
ylab("Suitability") +
coord_cartesian(xlim=c(-500, 10000))+
theme(plot.margin = unit(c(0,0,0,0), "cm"))
plot_grid(p1, p2, ncol = 1, align = "v", rel_heights = c(1.1, 1))
## ---- fig.height=3, fig.width=6, fig.cap="Biomass vs. age curves resulting from different growth rates for a 400 years life-span.", echo=FALSE----
#objects to store results
growth.rate.vector <- biomass.vector <- age.vector <- vector()
#params
age <- 0:400
maximum.biomass <- 100
growth.rates <- c(0.025, 0.05, 0.1, 0.2, 0.4, 0.8)
#iterating through growth rates
for(growth.rate in growth.rates){
biomass.vector <- c(biomass.vector,maximum.biomass / (1 + maximum.biomass * exp(- growth.rate * age)))
growth.rate.vector <- c(growth.rate.vector, rep(growth.rate, length(age)))
age.vector <- c(age.vector, age)
}
plot.df3 <- data.frame(growth.rate=growth.rate.vector, age=age.vector, biomass=biomass.vector)
plot.df3$growth.rate <- as.factor(plot.df3$growth.rate)
p3 <- ggplot(data=plot.df3, aes(x=age, y=biomass, group=growth.rate, color=growth.rate)) +
geom_line(size=1) +
scale_color_viridis(option="D", discrete=TRUE, direction=-1) +
ggtitle("Biomass gain under different growth rates") +
ylab("Biomass (relative)") + xlab("Age (years)")
p3
rm(age, age.vector, biomass.vector, growth.rate, growth.rate.vector, growth.rates, p3, plot.df3)
## ---- fig.height=3, fig.width=6, fig.cap="Risk curve defining the probability of removal of a given individual as a function of its fractional age when maximum carrying capacity is reached.", echo=FALSE----
temp.df <- data.frame(Pm=NA, Age=0:100)
temp.df$Pm <- 1 - sqrt(temp.df$Age/100)
temp.df$Age <- temp.df$Age/100
ggplot(data=temp.df, aes(x=Age, y=Pm, color=Pm)) +
geom_line(size=1) +
scale_color_viridis(option="D", direction=-1) +
ggtitle("Probability of mortality when carrying capacity is reached") +
ylab("Removal probability") +
xlab("Age of the individual (as proportion to the maximum age)") +
theme(legend.position = "none")
rm(temp.df)
## ---- eval=FALSE, size="small"------------------------------------------------
#
# #parameters (1st line in dataframe "parameters")
# maximum.age <- parameters[1, "maximum.age"]
# reproductive.age <- parameters[1, "reproductive.age"]
# growth.rate <- parameters[1, "growth.rate"]
# carrying.capacity <- parameters[1, "carrying.capacity"]
# fecundity <- parameters[1, "fecundity"]
#
# #reproductive age to proportion
# reproductive.age <- reproductive.age / maximum.age
#
# #years scaled taking maximum.age as reference
# scaled.year <- 1/maximum.age
#
# #vector to store pollen counts
# pollen.count <- vector()
#
# #starting population
# population <- sample(seq(0, 1, by=scaled.year),
# 100,
# replace=TRUE)
#
# #iterating through suitability (once per year)
# #------------------------------------
# for(suitability.i in suitability){
#
# #AGING -----------------------------------------------
# population <- population + scaled.year
#
# #SENESCENCE ------------------------------------------
# #1 is the maximum age of ages expressed as proportions
# population <- population[population < 1]
#
# #LOCAL EXTINCTION AND RECOLONIZATION -----------------
# if (length(population) == 0){
#
# #local extinction, replaces population with a seedbank
# population <- rep(0, floor(100 * suitability.i))
#
# #adds 0 to pollen.count
# pollen.count <- c(pollen.count, 0)
#
# #jumps to next iteration
# next
# }
#
# #PLANT GROWTH ---------------------------------------
# #biomass of every individual
# biomass <- maximum.biomass /
# (1 + maximum.biomass *
# exp(- (growth.rate * suitability.i) *
# (population * maximum.age)
# )
# )
#
# #SELF-THINNING --------------------------------------
# #carrying capacity reached
# while(sum(biomass) > carrying.capacity){
#
# #removes a random individual based on risk curve
# individual.to.remove <- sample(
# x = length(population),
# size = 1,
# replace = TRUE,
# prob = 1 - sqrt(population) #risk curve
# )
#
# #removing individuals from population and biomass
# population <- population[-individual.to.remove]
# biomass <- biomass[-individual.to.remove]
#
# }#end of while
#
# #REPRODUCTION --------------------------------------
# #identifyies adult individuals
# adults <- population > reproductive.age
#
# #seeds (vector of 0s)
# #fractional biomass of adults * fecundity * suitability
# seeds <- rep(0, floor(sum((biomass[adults]/maximum.biomass) *
# fecundity) * suitability.i))
#
# #adding seeds to the population
# population <- c(population, seeds)
#
# #POLLEN OUTPUT -------------------------------------
# #biomass of adults multiplied by suitability
# pollen.count <- c(pollen.count,
# sum(biomass[adults]) * suitability.i)
#
# } #end of loop through suitability values
#
## ---- size="small", fig.height=6.5, fig.width=9, message=TRUE, warning=TRUE, error=TRUE, results="hide", fig.cap="Simulation outcome. Green shades represent different age groups (seedlings, saplings, and adults).", warning=FALSE, message=FALSE----
#simulate population based on parameters
simulation <- simulatePopulation(parameters=parameters[1, ],
drivers=driver)
#plotting simulation output
plotSimulation(simulation.output=simulation,
burnin=FALSE,
panels=c("Driver A",
"Suitability",
"Population",
"Pollen"),
plot.title="",
text.size=12,
line.size=0.4)
## ---- message=TRUE, warning=TRUE, echo=FALSE----------------------------------
parameters[2,] <- c("Species 2", 50, 20, 4, 0.3, 0, 100, 10000, 1, 0, 50, 15, 0, 0, 600, 600)
parameters[3,] <- c("Species 3", 50, 20, 6, 0.4, 0.5, 100, 10000, 1, 0, 50, 20, 0, 0, 600, 600)
parameters[, 2:ncol(parameters)] <- sapply(parameters[, 2:ncol(parameters)], as.numeric)
parameters.t <- data.frame(t(parameters))
parameters.t <- parameters.t[c(2:9, 11:12),]
names(parameters.t) <- paste(parameters$label, sep="")
parameters.t$Parameter <- rownames(parameters.t)
rownames(parameters.t) <- NULL
parameters.t <- parameters.t[, c(ncol(parameters.t),1:(ncol(parameters.t)-1))]
kable(parameters.t, caption="Parameters of the three virtual species.", booktabs = T) %>% kable_styling(latex_options = c("hold_position", "striped"))
## ---- fig.height=5, fig.width=9, fig.cap="Comparison of the pollen abundance and environmental suitability (same in all cases) for the three virtual species shown in Table 2 within the period 5600-6400. Species 2 has a higher fecundity than Species 1 (1 vs 10)", message=FALSE----
#simulating species 2 and 3 of the dataframe parameters
simulation.2 <- simulatePopulation(parameters=parameters,
species=c(2,3),
drivers=driver)
#adding the results to the previous simulation
simulation <- c(simulation, simulation.2)
rm(simulation.2)
#plotting the comparison for the time interval between 4000 and 5000.-
compareSimulations(simulation.output=simulation,
species = "all",
columns = c("Suitability", "Pollen"),
time.zoom = c(5600, 6400))
## ---- echo=FALSE--------------------------------------------------------------
parameters.test <- parametersDataframe(rows=1)
parameters.test[1,] <- c("Test 1", 4, 1, 0.55, 2, 0, 1, 30, 0.5, 0.5, 50, 10, 50, 10, 100, 100)
parameters.test[2,] <- c("Test 2", 3, 1, 0.5, 2, 0, 1, 30, 0.5, 0.5, 50, 10, 50, 10, 100, 100)
parameters.test[, 2:ncol(parameters.test)] <- sapply(parameters.test[, 2:ncol(parameters.test)], as.numeric)
parameters.t <- data.frame(t(parameters.test))
parameters.t <- data.frame(parameters.t[-1,])
names(parameters.t) <- c("Test 1", "Test 2")
parameters.t$Parameter <- rownames(parameters.t)
rownames(parameters.t) <- NULL
parameters.t <- parameters.t[, c(ncol(parameters.t),1:(ncol(parameters.t)-1))]
parameters.t <- parameters.t[1:(nrow(parameters.t)-2),]
kable(parameters.t, caption="Parameters of virtual taxa used to test model limits.")
rm(parameters.t)
## ---- message=FALSE-----------------------------------------------------------
simulation.test.1 <- simulatePopulation(
parameters=parameters.test,
driver.A=jitter(rep(50, 500), amount=4)
)
## ---- fig.height=3, fig.width=9, fig.cap="Pollen output of virtual taxa Test 1 and Test 2 for a 200 years time-window."----
compareSimulations(simulation.output=simulation.test.1,
columns="Pollen",
time.zoom = c(0, 200))
## ---- echo=FALSE--------------------------------------------------------------
parameters.test[3,] <- c("Test 3", 1000, 100, 0.5, 0.05, 0, 100, 10000, 0.5, 0.5, 50, 10, 50, 10, 100, 100)
parameters.test[4,] <- c("Test 4", 1000, 500, 0.5, 0.05, 0, 100, 10000, 0.5, 0.5, 50, 10, 50, 10, 100, 100)
parameters.test[5,] <- c("Test 5", 1000, 900, 0.5, 0.05,0, 100, 10000, 0.5, 0.5, 50, 10, 50, 10, 100, 100)
parameters.test[, 2:ncol(parameters.test)] <- sapply(parameters.test[, 2:ncol(parameters.test)], as.numeric)
parameters.t <- data.frame(t(parameters.test[3:5,]))
parameters.t <- data.frame(parameters.t[-1,])
names(parameters.t) <- c("Test 3", "Test 4", "Test 5")
parameters.t$Parameter <- rownames(parameters.t)
rownames(parameters.t) <- NULL
parameters.t <- parameters.t[, c(ncol(parameters.t),1:(ncol(parameters.t)-1))]
parameters.t <- parameters.t[1:(nrow(parameters.t)-2),]
kable(parameters.t, caption="Parameters of virtual taxa used to test model limits.")
rm(parameters.t)
## ---- message=FALSE-----------------------------------------------------------
simulation.test.2 <- simulatePopulation(
parameters=parameters.test,
species=c(3:5),
driver.A=jitter(rep(50, 2000), amount=4)
)
## ---- fig.height=3, fig.width=9, fig.cap="Pollen output of virtual taxa Test 1 and Test 2 for a 200 years time-window."----
compareSimulations(simulation.output=simulation.test.2,
columns="Pollen",
time.zoom = c(0, 800))
## ---- fig.height=3, fig.width=9, fig.cap="Virtual sediment accumulation rate."----
accumulation.rate <- simulateAccumulationRate(seed=140,
time=1:10000,
output.min=1,
output.max=50)
## ---- echo=FALSE--------------------------------------------------------------
kable(accumulation.rate[1:20, ], caption="Dataframe resulting from the function to generate virtual accumulation rates. Each group in the grouping column has as many elements as accumulation.rate the given group has.", booktabs = T) %>% kable_styling(latex_options = c("hold_position", "striped"))
## ---- size="small"------------------------------------------------------------
simulation.aggregated <- aggregateSimulation(
simulation.output=simulation,
accumulation.rate=accumulation.rate,
sampling.intervals=c(2, 6, 10)
)
## ---- fig.height=5, fig.width=9, fig.cap="Effect of applying accumulation rate and different sampling depth intervals to a section of the the annual data produced by the simulation (represented in the Legend by the label Annual). Note that the 10 cm resampling completely misses the whole high-suitability event in the Pollen panel, and barely registers it in the Suitability panel. Inter-decadal variability shown by the Annual data is completely lost even at 1 cm, the higher sampling resolution.", echo=FALSE----
#checking results
temp.list <- simulation.aggregated[1, 1:5]
names(temp.list) <- c("Annual", "1 cm", "2 cm", "6 cm", "10 cm")
compareSimulations(simulation.output=temp.list, columns = c("Suitability","Pollen"), time.zoom=c(5800, 6200))
## ---- fig.height=3, fig.width=6, fig.cap="Histogram of the time differences (in years) between consecutive samples for the outcome of aggregateSimulation when resampled at intervals of 6 centimeters on Species 1. It clearly shows how the data are not organized in regular time intervals, and therefore are unsuitable for analyses requiring regular time lags.", echo=FALSE----
#getting example data sampled at 2cm intervals
simulated.data = simulation.aggregated[[1, 3]]
#checking distribution of differences in age between consecutive samples
hist(simulated.data[2:nrow(simulated.data),"Time"] - simulated.data[1:(nrow(simulated.data)-1),"Time"], main="Age differences between consecutive samples", xlab="Age difference between consecutive samples", col=viridis(12, begin = 0, end=1))
## ---- size="small", eval=FALSE------------------------------------------------
#
# #getting example data sampled at 2cm intervals
# simulated.data = simulation.aggregated[[1, 3]]
#
#
# #span values to be explored
# span.values=seq(20/nrow(simulated.data),
# 5/nrow(simulated.data),
# by=-0.0005)
#
#
# #plotting the optimization process in real time
# x11(height=12, width=24)
#
# #iteration through candidate span values
# for(span in span.values){
#
# #function to interpolate the data
# interpolation.function = loess(
# Pollen ~ Time,
# data=simulated.data,
# span=span,
# control=loess.control(surface="direct"))
#
# #plot model fit
# plot(simulated.data$Pollen, type="l", lwd=2)
# lines(interpolation.function$fitted, col="red", lwd=2)
#
# #if correlation equals 0.9985 loop stops
# if(cor(interpolation.function$fitted,
# simulated.data$Pollen) >= 0.9985){break}
#
# }
#
# #gives time to look at result before closing the plot window
# Sys.sleep(5)
## ---- warning=FALSE, size="small"---------------------------------------------
simulation.interpolated <- toRegularTime(
x=simulation.aggregated,
time.column="Time",
interpolation.interval=10,
columns.to.interpolate=c("Pollen",
"Suitability",
"Driver.A")
)
## ---- echo=FALSE , fig.height=5, fig.width=9, message=TRUE, error=TRUE, warning=TRUE, fig.cap="Data aggregated using virtual accumulation rate and reinterpolated into a regular time grid of 10 years resolution."----
#getting the results for Species 1
temp.list <- simulation.interpolated[1, 1:5]
names(temp.list) <- c("Annual", "1 cm", "2 cm", "6 cm", "10 cm")
#plotting comparison
compareSimulations(simulation.output=temp.list, columns = c("Suitability","Pollen"), time.zoom=c(5800, 6200))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.