In ggcostoya/tpcurves2: T-P Curves
knitr::opts_chunk$set(echo = TRUE)
devtools::load_all()
In this page I'll run two basic simulations representing the evolution of TPTs in a population to present the capabilities of the tpcurves2 package so far. These two simulations are aimed at replicating two scenarios described in the Logan & Cox 2020 Frontiers in Genetics paper, namely: an Specialist-generalist trade off (SGT) and the Thermodynamic effect (HIB) (a.k.a "Hotter is better").
Before moving forward with the simulations, I'll generate a basic population TPD which will be the starting point of all simulations:
# starting population TPD
pop <- gen_pop_tpd(N = 100, Topt = 25, CTmax = 28, CTmin = 23, Pmax = 10, Pmin = 0,
Error = 0.5, Samples = 10, Degree = 3)
I am considering an $N = 100$, which is a compromise since it might be a little too large to represent an actual isolated population but it just big enough to avoid extreme stochasticity.
Specialist-generalist trade off (SGT)
The first simulation focuses on the SGT scenario. In this case, I'll simulate how a population's TPC and TPTs respond to an increase in temperature variation. This scenario aims to represent the effects of selection due to increasingly extreme weather events, which is one of the expected threats from climate change.
First, I'll define a sequence of temperature ($Tm$) values. I'll do so using the gen_tmseq function included in this package. Using this function, I can specify from which to which value I want mean temperatures to shift in the sequence, how many samples should be taken and the amount of error. For the SGT example, the sequence will barely move from the $T_{opt}$ of the starting population but it will have a high Error value to insert a lot of temperature variability.
# SGT temperature sequence
sgt_tm <- gen_tmseq(Tmin = 25, Tmax = 25 + 10^-3, Nt = 50, Error = 0.3)
# modifiy data for plotting
sgt_tm_plot <- tibble(Temperature = sgt_tm, t = seq(1, length(sgt_tm), by = 1))
# plot
ggplot(data = sgt_tm_plot, aes( x = t, y = Temperature)) +
geom_line( col = "red", alpha = 0.5, lwd = 2) +
geom_point( col = "red", alpha = 0.55, size = 3) +
theme_classic() +
ggtitle("SGT Temperature Sequence") +
geom_hline(aes(yintercept = 25), lty = 2, lwd = 0.01)
Now, I'll run the function sim_mtime, an extension of the sim_time function discussed in detail in "Simulation Basics", that runs a simulation for multiple time periods. Before I run it, some considerations:
- I'll use the previously defined
sgt_tm vector of temperatures as the environmental temperature ($T_{e}$) for each interval.
- I consider the relationship between performance and survival to be $S_{p} = 0.95$. It seems reasonable that all performance is a good indicator for survival.
- I consider the relationship between performance and reproduction to be $R_p = 0.10$. Although a better performing individual will be more likely to reproduce, I am making it so that reproduction, like in nature, is not likely to happen every time.
- I keep mutation rate to a realistically low $\mu = 10^{-5}$ and I am limiting the mean amount of change in a trait to 1% it's current value with an sd of 0.5% the current trait value.
The code to run the function looks like:
sgt_sim <- sim_mtime(TPD = pop, Te = sgt_tm, Sp = 0.95, Rp = 0.1,
Mu = 10^-5, Mum = 0.01, Musd = 0.005,
Pmin = 0, Error = 0.1, Samples = 10, Degree = 3)
And results in the following transition in TPC shapes:
# preparing the overall plotting data
sgt_curves <- data.frame(Tm = NA, P = NA, L = NA, U = NA, t = NA) # prepare holding dataset
sgt_tpc <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent)) # nest columns on existing data
for(i in 1:nrow(sgt_tpc)){ # start loop to get each t's curve
which_time <- as.numeric(sgt_tpc$t[i]) # select time value
time <- sgt_tpc[i,] %>% unnest(cols = c(data)) # select a column and unnest that t's TPD
time_tpc <- get_tpc(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) # get the tpc for that t's TPD
time_tpc$t <- rep(which_time, nrow(time_tpc)) # add the time variable
sgt_curves <- rbind(sgt_curves, time_tpc) # bind data to holding dataset
}
sgt_curves <- sgt_curves[-c(1),]
# preparing the first and last data
target <- c(0,length(sgt_tm))
sgt_firstlast <- sgt_curves %>% filter(t %in% target)
# plotting the curves
ggplot(data = sgt_curves, aes(x = Tm, y = P)) +
geom_line(colour = "lightgrey", lwd = 0.5, alpha = 0.75) +
theme_classic() +
theme(legend.position = "none") +
geom_line(data = sgt_firstlast, aes(x = Tm, y = P, colour = factor(t)), lwd = 2) +
labs( x = "Temperature", y = "Performance", colour = "t")
As shown above, the shape of the population's TPC has changed significantly (In lightgrey all the other TPC shapes the simulation has gone through). On a first glance, $P_{max}$ seems to have moved upwards while $CT_{max}$ and $CT_{min}$ have not moved that much. My first observation would be to say that the population has adapted to a more extreme $T_m$ range through survival of those with better performances not those with wider TPCs. Let's take a look at the population's TPTs.
# prepare TPTs data
sgt_tpts <- data.frame(trait = NA, value = NA, t = NA)
sgt_tpts_data <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent))
for(i in 1:nrow(sgt_tpts_data)){
which_time <- as.numeric(sgt_tpts_data$t[i])
time <- sgt_tpts_data[i,] %>% unnest(cols = c(data))
time_tpts <- get_tpts(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0)
time_tpts$t <- rep(which_time, nrow(time_tpts))
sgt_tpts <- rbind(sgt_tpts, time_tpts)
}
sgt_tpts <- sgt_tpts[-c(1),]
# plot TPTs data
# T traits plot
ttraits <- c("Topt", "CTmax", "CTmin")
tsplot <- sgt_tpts %>% filter(trait %in% ttraits) %>%
ggplot(aes(x = t, y = value)) +
geom_line(aes(colour = factor(trait)), lwd = 1) +
geom_point(aes(colour = factor(trait))) +
theme_classic() +
theme(legend.position = "top") +
geom_hline(yintercept = c(sgt_tpts$value[1], sgt_tpts$value[2], sgt_tpts$value[3]), lty = 2) +
labs(x = "t", y = "TPT value", colour = "TPT")
# Pmax plot
pmaxplot <- sgt_tpts %>% filter(trait == "Pmax") %>%
ggplot(aes(x = t, y = value)) +
geom_line(lwd = 1) +
geom_point() +
theme_classic() +
geom_hline(yintercept = c(sgt_tpts$value[5]), lty = 2) +
labs(x = "t", y = "Pmax")
# arrange in a grid
grid.arrange(tsplot, pmaxplot, ncol = 2)
Although the y-axis is the same scale it is clear that the only trait that has shifted upwards is $P_{max}$. Another important bit of information is the fluctuation in population size:
# prepare dataset
sgt_n <- data.frame(N = NA, t = NA)
sgt_n_data <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent))
for(i in 1:nrow(sgt_n_data)){
which_time <- as.numeric(sgt_n_data$t[i])
popsize <- sgt_n_data[i,] %>% unnest(cols = c(data)) %>% nrow() %>% as.numeric()
time_popsize <- data.frame(N = popsize/10, t = which_time)
sgt_n <- rbind(sgt_n, time_popsize)
}
sgt_n <- sgt_n[-c(1),]
# plot data
ggplot(data = sgt_n, aes(x = t, y = N)) +
geom_line(lwd = 1.5, alpha = 0.5) +
geom_point(size = 2, alpha = 0.5) +
theme_classic() +
labs(x = "t", y = "Population size")
Every time the simulation is ran the population size assumes a different dynamic over time. One thing I notice every time is that, in all cases, extreme temperature events result in extreme population size drops, whereas more stable periods result in stable growth.
Thermodynamic effect / Hotter is better (HIB)
The second simulation focuses on the HIB scenario. This time I'll simulate a population's response to a steady, stable increase of 1 C in mean temperature during 50 time periods. Once again, I'll use gen_tmseq to generate a temperature sequence:
# SGT temperature sequence
hib_tm <- gen_tmseq(Tmin = 25, Tmax = 26, Nt = 50, Error = 0.05)
# modifiy data for plotting
hib_tm_plot <- tibble(Temperature = hib_tm, t = seq(1, length(hib_tm), by = 1))
# plot
ggplot(data = hib_tm_plot, aes( x = t, y = Temperature)) +
geom_line( col = "red", alpha = 0.5, lwd = 2) +
geom_point( col = "red", alpha = 0.55, size = 3) +
theme_classic() +
ggtitle("HIB Temperature Sequence") +
geom_hline(aes(yintercept = 25), lty = 2, lwd = 0.01)
To simulate the survival and reproduction of the initial population I'll use the exact same parameter values as in the SGT example:
hib_sim <- sim_mtime(TPD = pop, Te = hib_tm, Sp = 0.95, Rp = 0.1,
Mu = 10^-5, Mum = 0.01, Musd = 0.005,
Pmin = 0, Error = 0.1, Samples = 10, Degree = 3)
Let's see how the TPC, TPTs and Population size changed.
# TPC plot
# preparing the overall plotting data
hib_curves <- data.frame(Tm = NA, P = NA, L = NA, U = NA, t = NA) # prepare holding dataset
hib_tpc <- hib_sim %>% nest(data = c(ID, Tm, P, Parent)) # nest columns on existing data
for(i in 1:nrow(hib_tpc)){ # start loop to get each t's curve
which_time <- as.numeric(hib_tpc$t[i]) # select time value
time <- hib_tpc[i,] %>% unnest(cols = c(data)) # select a column and unnest that t's TPD
time_tpc <- get_tpc(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) # get the tpc for that t's TPD
time_tpc$t <- rep(which_time, nrow(time_tpc)) # add the time variable
hib_curves <- rbind(hib_curves, time_tpc) # bind data to holding dataset
}
hib_curves <- hib_curves[-c(1),]
# preparing the first and last data
target <- c(0,length(hib_tm))
hib_firstlast <- hib_curves %>% filter(t %in% target)
# plotting the curves
hib_tpcplot <- ggplot(data = hib_curves, aes(x = Tm, y = P)) +
geom_line(colour = "lightgrey", lwd = 0.5, alpha = 0.75) +
theme_classic() +
theme(legend.position = "none") +
geom_line(data = hib_firstlast, aes(x = Tm, y = P, colour = factor(t)), lwd = 2) +
labs( x = "Temperature", y = "Performance", colour = "t")
# TPTs plot
# prepare TPTs data
hib_tpts <- data.frame(trait = NA, value = NA, t = NA)
hib_tpts_data <- hib_sim %>% nest(data = c(ID, Tm, P, Parent))
for(i in 1:nrow(hib_tpts_data)){
which_time <- as.numeric(hib_tpts_data$t[i])
time <- hib_tpts_data[i,] %>% unnest(cols = c(data))
time_tpts <- get_tpts(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0)
time_tpts$t <- rep(which_time, nrow(time_tpts))
hib_tpts <- rbind(hib_tpts, time_tpts)
}
hib_tpts <- hib_tpts[-c(1),]
# plot TPTs data
# T traits plot
ttraits <- c("Topt", "CTmax", "CTmin")
hib_tsplot <- hib_tpts %>% filter(trait %in% ttraits) %>%
ggplot(aes(x = t, y = value)) +
geom_line(aes(colour = factor(trait)), lwd = 1) +
geom_point(aes(colour = factor(trait))) +
theme_classic() +
theme(legend.position = "top") +
geom_hline(yintercept = c(hib_tpts$value[1], hib_tpts$value[2], hib_tpts$value[3]), lty = 2) +
labs(x = "t", y = "TPT value", colour = "TPT")
# Pmax plot
hib_pmaxplot <- hib_tpts %>% filter(trait == "Pmax") %>%
ggplot(aes(x = t, y = value)) +
geom_line(lwd = 1) +
geom_point() +
theme_classic() +
geom_hline(yintercept = c(hib_tpts$value[5]), lty = 2) +
labs(x = "t", y = "Pmax")
# N plot
# prepare dataset
hib_n <- data.frame(N = NA, t = NA)
hib_n_data <- hib_sim %>% nest(data = c(ID, Tm, P, Parent))
for(i in 1:nrow(hib_n_data)){
which_time <- as.numeric(hib_n_data$t[i])
popsize <- hib_n_data[i,] %>% unnest(cols = c(data)) %>% nrow() %>% as.numeric()
time_popsize <- data.frame(N = popsize/10, t = which_time)
hib_n <- rbind(hib_n, time_popsize)
}
hib_n <- hib_n[-c(1),]
# plot data
hib_nplot <- ggplot(data = hib_n, aes(x = t, y = N)) +
geom_line(lwd = 1.5, alpha = 0.5) +
geom_point(size = 2, alpha = 0.5) +
theme_classic() +
labs(x = "t", y = "Population size")
grid.arrange(hib_tpcplot, hib_nplot, hib_tsplot, hib_pmaxplot, ncol = 2, nrow = 2)
For now, in most simulation runs, the shape of the TPC shifts as the HIB theory predicts, upwards and to the right. Nonetheless, the population undergoes extreme and unrealistic population size shifts that I can't explain. I hope to explain it and make it more realistic in future iterations.
Potential Improvements
-
Improve the simulation speed significantly by avoiding loops and implementing vectorization.
-
Improve visualization of the curve shift. Ideas are welcome.
-
Improve overall realism of the simulations.
ggcostoya/tpcurves2 documentation built on Jan. 1, 2021, 2:19 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.
knitr::opts_chunk$set(echo = TRUE) devtools::load_all()
In this page I'll run two basic simulations representing the evolution of TPTs in a population to present the capabilities of the tpcurves2 package so far. These two simulations are aimed at replicating two scenarios described in the Logan & Cox 2020 Frontiers in Genetics paper, namely: an Specialist-generalist trade off (SGT) and the Thermodynamic effect (HIB) (a.k.a "Hotter is better").
Before moving forward with the simulations, I'll generate a basic population TPD which will be the starting point of all simulations:
# starting population TPD pop <- gen_pop_tpd(N = 100, Topt = 25, CTmax = 28, CTmin = 23, Pmax = 10, Pmin = 0, Error = 0.5, Samples = 10, Degree = 3)
I am considering an $N = 100$, which is a compromise since it might be a little too large to represent an actual isolated population but it just big enough to avoid extreme stochasticity.
Specialist-generalist trade off (SGT)
The first simulation focuses on the SGT scenario. In this case, I'll simulate how a population's TPC and TPTs respond to an increase in temperature variation. This scenario aims to represent the effects of selection due to increasingly extreme weather events, which is one of the expected threats from climate change.
First, I'll define a sequence of temperature ($Tm$) values. I'll do so using the gen_tmseq function included in this package. Using this function, I can specify from which to which value I want mean temperatures to shift in the sequence, how many samples should be taken and the amount of error. For the SGT example, the sequence will barely move from the $T_{opt}$ of the starting population but it will have a high Error value to insert a lot of temperature variability.
# SGT temperature sequence sgt_tm <- gen_tmseq(Tmin = 25, Tmax = 25 + 10^-3, Nt = 50, Error = 0.3)
# modifiy data for plotting sgt_tm_plot <- tibble(Temperature = sgt_tm, t = seq(1, length(sgt_tm), by = 1)) # plot ggplot(data = sgt_tm_plot, aes( x = t, y = Temperature)) + geom_line( col = "red", alpha = 0.5, lwd = 2) + geom_point( col = "red", alpha = 0.55, size = 3) + theme_classic() + ggtitle("SGT Temperature Sequence") + geom_hline(aes(yintercept = 25), lty = 2, lwd = 0.01)
Now, I'll run the function sim_mtime, an extension of the sim_time function discussed in detail in "Simulation Basics", that runs a simulation for multiple time periods. Before I run it, some considerations:
- I'll use the previously defined
sgt_tmvector of temperatures as the environmental temperature ($T_{e}$) for each interval. - I consider the relationship between performance and survival to be $S_{p} = 0.95$. It seems reasonable that all performance is a good indicator for survival.
- I consider the relationship between performance and reproduction to be $R_p = 0.10$. Although a better performing individual will be more likely to reproduce, I am making it so that reproduction, like in nature, is not likely to happen every time.
- I keep mutation rate to a realistically low $\mu = 10^{-5}$ and I am limiting the mean amount of change in a trait to 1% it's current value with an sd of 0.5% the current trait value.
The code to run the function looks like:
sgt_sim <- sim_mtime(TPD = pop, Te = sgt_tm, Sp = 0.95, Rp = 0.1, Mu = 10^-5, Mum = 0.01, Musd = 0.005, Pmin = 0, Error = 0.1, Samples = 10, Degree = 3)
And results in the following transition in TPC shapes:
# preparing the overall plotting data sgt_curves <- data.frame(Tm = NA, P = NA, L = NA, U = NA, t = NA) # prepare holding dataset sgt_tpc <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent)) # nest columns on existing data for(i in 1:nrow(sgt_tpc)){ # start loop to get each t's curve which_time <- as.numeric(sgt_tpc$t[i]) # select time value time <- sgt_tpc[i,] %>% unnest(cols = c(data)) # select a column and unnest that t's TPD time_tpc <- get_tpc(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) # get the tpc for that t's TPD time_tpc$t <- rep(which_time, nrow(time_tpc)) # add the time variable sgt_curves <- rbind(sgt_curves, time_tpc) # bind data to holding dataset } sgt_curves <- sgt_curves[-c(1),] # preparing the first and last data target <- c(0,length(sgt_tm)) sgt_firstlast <- sgt_curves %>% filter(t %in% target) # plotting the curves ggplot(data = sgt_curves, aes(x = Tm, y = P)) + geom_line(colour = "lightgrey", lwd = 0.5, alpha = 0.75) + theme_classic() + theme(legend.position = "none") + geom_line(data = sgt_firstlast, aes(x = Tm, y = P, colour = factor(t)), lwd = 2) + labs( x = "Temperature", y = "Performance", colour = "t")
As shown above, the shape of the population's TPC has changed significantly (In lightgrey all the other TPC shapes the simulation has gone through). On a first glance, $P_{max}$ seems to have moved upwards while $CT_{max}$ and $CT_{min}$ have not moved that much. My first observation would be to say that the population has adapted to a more extreme $T_m$ range through survival of those with better performances not those with wider TPCs. Let's take a look at the population's TPTs.
# prepare TPTs data sgt_tpts <- data.frame(trait = NA, value = NA, t = NA) sgt_tpts_data <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent)) for(i in 1:nrow(sgt_tpts_data)){ which_time <- as.numeric(sgt_tpts_data$t[i]) time <- sgt_tpts_data[i,] %>% unnest(cols = c(data)) time_tpts <- get_tpts(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) time_tpts$t <- rep(which_time, nrow(time_tpts)) sgt_tpts <- rbind(sgt_tpts, time_tpts) } sgt_tpts <- sgt_tpts[-c(1),] # plot TPTs data # T traits plot ttraits <- c("Topt", "CTmax", "CTmin") tsplot <- sgt_tpts %>% filter(trait %in% ttraits) %>% ggplot(aes(x = t, y = value)) + geom_line(aes(colour = factor(trait)), lwd = 1) + geom_point(aes(colour = factor(trait))) + theme_classic() + theme(legend.position = "top") + geom_hline(yintercept = c(sgt_tpts$value[1], sgt_tpts$value[2], sgt_tpts$value[3]), lty = 2) + labs(x = "t", y = "TPT value", colour = "TPT") # Pmax plot pmaxplot <- sgt_tpts %>% filter(trait == "Pmax") %>% ggplot(aes(x = t, y = value)) + geom_line(lwd = 1) + geom_point() + theme_classic() + geom_hline(yintercept = c(sgt_tpts$value[5]), lty = 2) + labs(x = "t", y = "Pmax") # arrange in a grid grid.arrange(tsplot, pmaxplot, ncol = 2)
Although the y-axis is the same scale it is clear that the only trait that has shifted upwards is $P_{max}$. Another important bit of information is the fluctuation in population size:
# prepare dataset sgt_n <- data.frame(N = NA, t = NA) sgt_n_data <- sgt_sim %>% nest(data = c(ID, Tm, P, Parent)) for(i in 1:nrow(sgt_n_data)){ which_time <- as.numeric(sgt_n_data$t[i]) popsize <- sgt_n_data[i,] %>% unnest(cols = c(data)) %>% nrow() %>% as.numeric() time_popsize <- data.frame(N = popsize/10, t = which_time) sgt_n <- rbind(sgt_n, time_popsize) } sgt_n <- sgt_n[-c(1),] # plot data ggplot(data = sgt_n, aes(x = t, y = N)) + geom_line(lwd = 1.5, alpha = 0.5) + geom_point(size = 2, alpha = 0.5) + theme_classic() + labs(x = "t", y = "Population size")
Every time the simulation is ran the population size assumes a different dynamic over time. One thing I notice every time is that, in all cases, extreme temperature events result in extreme population size drops, whereas more stable periods result in stable growth.
Thermodynamic effect / Hotter is better (HIB)
The second simulation focuses on the HIB scenario. This time I'll simulate a population's response to a steady, stable increase of 1 C in mean temperature during 50 time periods. Once again, I'll use gen_tmseq to generate a temperature sequence:
# SGT temperature sequence hib_tm <- gen_tmseq(Tmin = 25, Tmax = 26, Nt = 50, Error = 0.05)
# modifiy data for plotting hib_tm_plot <- tibble(Temperature = hib_tm, t = seq(1, length(hib_tm), by = 1)) # plot ggplot(data = hib_tm_plot, aes( x = t, y = Temperature)) + geom_line( col = "red", alpha = 0.5, lwd = 2) + geom_point( col = "red", alpha = 0.55, size = 3) + theme_classic() + ggtitle("HIB Temperature Sequence") + geom_hline(aes(yintercept = 25), lty = 2, lwd = 0.01)
To simulate the survival and reproduction of the initial population I'll use the exact same parameter values as in the SGT example:
hib_sim <- sim_mtime(TPD = pop, Te = hib_tm, Sp = 0.95, Rp = 0.1, Mu = 10^-5, Mum = 0.01, Musd = 0.005, Pmin = 0, Error = 0.1, Samples = 10, Degree = 3)
Let's see how the TPC, TPTs and Population size changed.
# TPC plot # preparing the overall plotting data hib_curves <- data.frame(Tm = NA, P = NA, L = NA, U = NA, t = NA) # prepare holding dataset hib_tpc <- hib_sim %>% nest(data = c(ID, Tm, P, Parent)) # nest columns on existing data for(i in 1:nrow(hib_tpc)){ # start loop to get each t's curve which_time <- as.numeric(hib_tpc$t[i]) # select time value time <- hib_tpc[i,] %>% unnest(cols = c(data)) # select a column and unnest that t's TPD time_tpc <- get_tpc(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) # get the tpc for that t's TPD time_tpc$t <- rep(which_time, nrow(time_tpc)) # add the time variable hib_curves <- rbind(hib_curves, time_tpc) # bind data to holding dataset } hib_curves <- hib_curves[-c(1),] # preparing the first and last data target <- c(0,length(hib_tm)) hib_firstlast <- hib_curves %>% filter(t %in% target) # plotting the curves hib_tpcplot <- ggplot(data = hib_curves, aes(x = Tm, y = P)) + geom_line(colour = "lightgrey", lwd = 0.5, alpha = 0.75) + theme_classic() + theme(legend.position = "none") + geom_line(data = hib_firstlast, aes(x = Tm, y = P, colour = factor(t)), lwd = 2) + labs( x = "Temperature", y = "Performance", colour = "t") # TPTs plot # prepare TPTs data hib_tpts <- data.frame(trait = NA, value = NA, t = NA) hib_tpts_data <- hib_sim %>% nest(data = c(ID, Tm, P, Parent)) for(i in 1:nrow(hib_tpts_data)){ which_time <- as.numeric(hib_tpts_data$t[i]) time <- hib_tpts_data[i,] %>% unnest(cols = c(data)) time_tpts <- get_tpts(Tm = time$Tm, P = time$P, Degree = 3, Pmin = 0) time_tpts$t <- rep(which_time, nrow(time_tpts)) hib_tpts <- rbind(hib_tpts, time_tpts) } hib_tpts <- hib_tpts[-c(1),] # plot TPTs data # T traits plot ttraits <- c("Topt", "CTmax", "CTmin") hib_tsplot <- hib_tpts %>% filter(trait %in% ttraits) %>% ggplot(aes(x = t, y = value)) + geom_line(aes(colour = factor(trait)), lwd = 1) + geom_point(aes(colour = factor(trait))) + theme_classic() + theme(legend.position = "top") + geom_hline(yintercept = c(hib_tpts$value[1], hib_tpts$value[2], hib_tpts$value[3]), lty = 2) + labs(x = "t", y = "TPT value", colour = "TPT") # Pmax plot hib_pmaxplot <- hib_tpts %>% filter(trait == "Pmax") %>% ggplot(aes(x = t, y = value)) + geom_line(lwd = 1) + geom_point() + theme_classic() + geom_hline(yintercept = c(hib_tpts$value[5]), lty = 2) + labs(x = "t", y = "Pmax") # N plot # prepare dataset hib_n <- data.frame(N = NA, t = NA) hib_n_data <- hib_sim %>% nest(data = c(ID, Tm, P, Parent)) for(i in 1:nrow(hib_n_data)){ which_time <- as.numeric(hib_n_data$t[i]) popsize <- hib_n_data[i,] %>% unnest(cols = c(data)) %>% nrow() %>% as.numeric() time_popsize <- data.frame(N = popsize/10, t = which_time) hib_n <- rbind(hib_n, time_popsize) } hib_n <- hib_n[-c(1),] # plot data hib_nplot <- ggplot(data = hib_n, aes(x = t, y = N)) + geom_line(lwd = 1.5, alpha = 0.5) + geom_point(size = 2, alpha = 0.5) + theme_classic() + labs(x = "t", y = "Population size") grid.arrange(hib_tpcplot, hib_nplot, hib_tsplot, hib_pmaxplot, ncol = 2, nrow = 2)
For now, in most simulation runs, the shape of the TPC shifts as the HIB theory predicts, upwards and to the right. Nonetheless, the population undergoes extreme and unrealistic population size shifts that I can't explain. I hope to explain it and make it more realistic in future iterations.
Potential Improvements
-
Improve the simulation speed significantly by avoiding loops and implementing vectorization.
-
Improve visualization of the curve shift. Ideas are welcome.
-
Improve overall realism of the simulations.
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.