Nothing
#' Generates drivers for \code{\link{simulatePopulation}}.
#'
#' @description Wrapper of \code{simulateDriver} to generate several drivers with different autocorrelation lengths, and return a long format dataframe to be used as input for \code{\link{simulatePopulation}}. It also produces a plot of the generated drivers. \strong{Important}: note that the variable \code{time} runs from left to right, with lower values representing older samples.
#'
#'@usage simulateDriverS(
#' random.seeds=c(60, 120),
#' time=1:10000,
#' autocorrelation.lengths=c(200, 600),
#' output.min=c(0,0),
#' output.max=c(100, 100),
#' driver.names=c("A", "B"),
#' filename=NULL)
#'
#' @param random.seeds vector of integers, seeds to be used by \code{set.seed}.
#' @param time integer, or numeric vector of integers with constant intervals. If a single integer is provided, a time sequence is generated from 0 to the given integer as \emph{seq(0, time, by = 1)}. Default is 1:10000.
#' @param autocorrelation.lengths vector of integers, represents the lengths of the convolution filters to be used to impose a particular temporal structure on each driver. Default is 100, equivalent to a filter composed by a hundred of ones.
#' @param output.min numeric vector, minimum values of the output time series. Used as input for \strong{rescaleVector}. Default is 0.
#' @param output.max numeric vector, maximum values of the output time series. Used as input for \strong{rescaleVector}. Default is 100.
#' @param driver.names character vector, with labels to be used to identify the drivers.
#' @param filename character string, name of output pdf file. If NULL or empty, no pdf is produced.
#'
#' @details It is recommended to use \code{time} vectors with a time step of 1 between consecutive values when the output is to be used as input for \code{\link{simulatePopulation}}, which considers annual time-steps while simulating virtual pollen curves. Initial random sequence is generated by \code{rnorm}. Desired temporal autocorrelation are approximate, but deviation becomes higher if \code{autocorrelation.length} is larger than half the length of \code{time}. Consequently, the function limits \code{autocorrelation.length} to \code{length(time)/2}.
#'
#' @author Blas M. Benito <blasbenito@gmail.com>
#'
#' @return A long format dataframe (see dataset \code{\link{drivers}}) with the following columns:
#'
#' \itemize{
#' \item \emph{time}: integer.
#' \item \emph{driver}: character, values are \code{A} and \code{B}
#' \item \emph{autocorrelation.length}: numeric, default values are 200, 600, and 1800.
#' \item \emph{value}: numeric, value of the driver for the given \emph{time}.
#' }
#'
#' @seealso \code{\link{drivers}}, \code{\link{simulateDriver}}, \code{\link{simulatePopulation}}, \code{\link{drivers}}
#'
#' @examples
#'
#'drivers <- simulateDriverS(
#' random.seeds=c(60, 120),
#' time=1:10000,
#' autocorrelation.lengths=c(200, 600),
#' output.min=c(0,0),
#' output.max=c(100, 100),
#' driver.names=c("A", "B"),
#' filename=NULL
#')
#'
#'str(drivers)
#'
#' @export
simulateDriverS <- function(random.seeds=c(60, 120),
time=1:10000,
autocorrelation.lengths=c(200, 600),
output.min=c(0,0),
output.max=c(100, 100),
driver.names=c("A", "B"),
filename=NULL){
#TESTING INPUT DATA
#number of driver names
n.names <- length(driver.names)
#driver.names
if(!is.character(driver.names)){stop("The argument driver.names should be a character vector.")}
#random.seeds
if(!is.numeric(random.seeds)){
random.seeds <- 1:n.names
} else {
random.seeds <- as.integer(random.seeds)
}
if(length(random.seeds) < n.names){
random.seeds <- random.seeds[1]:(random.seeds[1]+n.names)
}
#output.min and output.max
if(length(output.min) < n.names){
output.min <- rep(output.min[1], n.names)
}
if(length(output.min) > n.names){
output.min <- output.min[1:n.names]
}
if(length(output.max) < n.names){
output.max <- rep(output.min[1], n.names)
}
if(length(output.max) > n.names){
output.max <- output.max[1:n.names]
}
for(i in 1:n.names){
if(output.max[i] < output.min[i]){
temp.min <- output.max[i]
output.max[i] <- output.min[i]
output.min[i] <- temp.min
}
}
#time
if(!is.numeric(time)){stop("time should be a numeric vector, try 1:1000.")}
if(length(time)==1){time <- 1:floor(time)}
#autocorrelation.length
if(!is.numeric(autocorrelation.lengths)){stop("autocorrelation.length should be a numeric vector.")}
#data ranges and random seed for each drivers.2k
data.ranges <- data.frame(driver=driver.names, output.min, output.max, random.seeds)
#dataframes to store drivers
drivers <- drivers.temp <- data.frame(time=numeric(), driver=character(), autocorrelation.length=numeric(), value=numeric())
autocorrelation <- data.frame(lag=numeric(), acf=numeric(), ci.max=numeric(), ci.min=numeric(), driver=character(), autocorrelation.length=numeric())
#looping through drivers and memory lengths
for(driver in driver.names){
for(autocorrelation.length in autocorrelation.lengths){
#FILLING drivers.temp
#---------------------------------------
#fill drivers.temp with time and grouping
drivers.temp[max(time), ] <- NA
drivers.temp$time <- time
#fill with parameter values
drivers.temp$driver <- rep(driver, nrow(drivers.temp))
drivers.temp$autocorrelation.length <- rep(autocorrelation.length, nrow(drivers.temp))
#fill values of the driver
simulated.driver <- simulateDriver(
random.seed = data.ranges[data.ranges$driver==driver, "random.seeds"],
time = time,
autocorrelation.length = autocorrelation.length,
output.min = data.ranges[data.ranges$driver==driver, "output.min"],
output.max = data.ranges[data.ranges$driver==driver, "output.max"]
)
#getting the driver values
drivers.temp$value <- rescaleVector(x = simulated.driver, new.max = 100, new.min = 0)
#FILLING autocorrelation.temp
#---------------------------------------
#computing acf
autocorrelation.temp <- acfToDf(x = simulated.driver, lag.max = max(max(autocorrelation.lengths)), length.out = 100)
autocorrelation.temp$driver <- rep(driver, nrow(autocorrelation.temp))
autocorrelation.temp$autocorrelation.length <- rep(autocorrelation.length, nrow(autocorrelation.temp))
#merging with main dataframes
drivers <- rbind(drivers, drivers.temp)
autocorrelation <- rbind(autocorrelation, autocorrelation.temp)
}
}
#PLOTTING OUTPUT
p.drivers <- ggplot(data = drivers, aes(x = time, y = value, color = driver)) +
geom_line() +
scale_color_viridis(discrete = TRUE, begin = 0.2, end = 0.6) +
facet_wrap(driver ~ autocorrelation.length, ncol = 1, scales = "free_y") +
xlab("Time (years)") +
ylab("") +
ggtitle("Virtual drivers") +
theme(plot.margin = unit(c(0.5, -1, 0.5, 0.5), "cm"),
legend.position = "none",
panel.background = element_blank()) +
cowplot::background_grid(major = "none", minor = "none")
p.density <- ggplot(data = drivers, aes(value, fill = driver, colour = driver)) +
geom_density(aes(y = ..scaled..), alpha = 0.5) +
scale_color_viridis(discrete = TRUE, begin = 0.2, end = 0.6) +
scale_fill_viridis(discrete = TRUE, begin = 0.2, end = 0.6) +
facet_wrap(driver ~ autocorrelation.length, ncol = 1, scales = "free") +
xlab("Years") +
ylab("") +
xlab("") +
coord_flip() +
theme(plot.margin = unit(c(0.5, 0.5, 0.5, 0), "cm"),
axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
strip.background = element_blank(),
strip.text.x = element_blank(),
legend.position = "none",
panel.background = element_blank()) +
cowplot::background_grid(major = "none", minor = "none")
p.acfs <- ggplot(data = autocorrelation, aes(x = lag, y = acf, color = driver)) +
scale_color_viridis(discrete = TRUE, begin = 0.2, end = 0.6) +
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)) +
ylab("") +
xlab("Lag (years)") +
facet_wrap(driver ~ autocorrelation.length, ncol = 1, scales = "free_y") +
ggtitle("Temporal autocorrelation") +
theme(plot.margin = unit(c(0.5, 0.5, 0.5, -1), "cm"),
legend.position = "none",
panel.background = element_blank()) +
cowplot::background_grid(major = "none", minor = "none")
print(cowplot::plot_grid(p.drivers, NULL, p.density, NULL, p.acfs, align = "h", ncol = 5, rel_widths = c(1, 0, 0.3, 0, 1)))
if(!is.null(filename) & is.character(filename)){
ggsave(width = 12, height = (1.5*(length(driver.names) * length(autocorrelation.lengths))), filename = paste(filename, ".pdf", sep = ""))
}
return(drivers)
}
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.