This Rmarkdown accompanies the Perspective article for Conservation Physiology authored by Arnold, Walker, Wishart & Guja.
It runs simulations microclimate at a specified location and time period. It then conducts Thermal Load Sensitivity (TLS) analyses for key life stage transitions in a generic plant.
Disclaimer: This is intended as a demonstration of the potential for the TLS approach to simulate thermal vulnerability at different developmental stages. The values used in simulation are based on unpublished data for a range of species and are intended as plausible values rather than empirical values of a particular species.
Contact Pieter Arnold (pieter.arnold@anu.edu.au) for any further information.
1) Packages
Install and load the required packages.
Code
# Load required libraries after installation of missing packages. Some developer versions available only via github.# library(devtools) # remotes::install_github("ropensci/rnoaa")# remotes::install_github("ilyamaclean/microclima")# remotes::install_github("mrke/NicheMapR")pacman::p_load(tidyverse, # Tidy tools ggpubr, # For combined plotting cowplot, # For combined plotting NicheMapR, # For microclimate and ectotherm biophysical leaf models parallel) # For parallel computation# Colourblind friendly colour palettecbp2 <-c("#000000", "#E69F00", "#56B4E9", "#009E73","#F0E442", "#0072B2", "#D55E00", "#CC79A7")
2) Load functions
Load the required functions for TLS analysis. Note that several damage_accum functions have been developed in collaboration with Daniel Noble, and these may include some modification for specific purpose here. Further discussion and related TLS code can be accessed here: https://github.com/daniel1noble/thermal_tol_interactions
Code
# Function source: Noble et al. 2025 https://doi.org/10.32942/X2BP8N # Calculates cumulative heat damage through time with temperature-dependent repair (Sharpe–Schoolfield form), tracking the damage state over time and computing instantaneous hazard and mortality probability at each time step. The threshold is defined such that when cumulative damage reaches 100 * threshold, mortality is 50% (e.g., threshold = 0.5 for LT50).damage_repair_accum_state <-function(time, temp, intercept, slope,threshold =0.5, TA, TAL, TAH, TL, TH, TREF, kdot_ref) {if (!is.numeric(temp)) stop("'temp' must be numeric (°C).")if (length(time) !=length(temp)) {stop("'time' and 'temp' must have the same length.") }if (length(time) <2L) {stop("Need at least two time points to compute damage increments.") }# Kelvin temps TK <- temp +273.15# Interval lengths and step temperatures (max over interval) dt <-diff(time) # minutes temp_st <-pmax(temp[-1], temp[-length(temp)]) TK_st <-pmax(TK[-1], TK[-length(TK)]) n <-length(dt)# Damage rate per interval (without repair) damage_rate <-100/ (10^(slope * temp_st + intercept))# Sharpe–Schoolfield repair rate per interval (per minute) repair_rate <-vapply(seq_len(n), function(i) { T_K <- TK_st[i] num <- kdot_ref *exp(TA * (1/ TREF -1/ T_K)) den <-1+exp(TAL * (1/ T_K -1/ TL)) +exp(TAH * (1/ TH -1/ T_K)) num / den }, numeric(1))# Raw increments for bookkeeping damage_i <- damage_rate * dt repair_inc <- repair_rate * dt# State variables net_damage_inc <-numeric(n) # can be negative cum_damage_state <-numeric(n) # damage STATE with repair hazard_inst <-numeric(n) # instantaneous hazard from state mortality_prob <-numeric(n) D50 <-100* threshold D_prev <-0for (i inseq_len(n)) {# net increment in damage state this interval net_inc <- (damage_rate[i] - repair_rate[i]) * dt[i] # can be < 0 D_new <-max(D_prev + net_inc, 0) # repair can undo damage# instantaneous hazard given CURRENT damage state h_i <-log(2) * (D_new / D50)# interpret this as "mortality if tested now" mort_i <-1-exp(-h_i) net_damage_inc[i] <- net_inc cum_damage_state[i] <- D_new hazard_inst[i] <- h_i mortality_prob[i] <- mort_i D_prev <- D_new } out <-data.frame(damage_i = damage_i,repair_inc = repair_inc,net_damage_inc = net_damage_inc,cum_damage_repair = cum_damage_state, # current damage statehazard = hazard_inst, # from state, not cumulativemortality_prob = mortality_prob, # can go up or downtemp = temp_st,time_min = time[-1] )return(out)}# Function source: Noble et al. 2025 https://doi.org/10.32942/X2BP8N # Calculates a cumulative percentage increase through time, where each interval’s increment depends on the temperature.damage_accum <-function(time, temp, intercept, slope, Tc, threshold =0.5) { output =data.frame(damage_i =NA, cum_damage =NA, mortality_prob =NA)# Calculate the damage in a given incrementfor (i in1:(length(time) -1)) {if(temp[i] > Tc) {# Calculate the accumulated damage increment for each time step acc = (100* (time[i +1] - time[i])) / (10^(slope *max(temp[i +1], temp[i]) + intercept)) output[i, "damage_i"] = acc }else { acc =0 output[i, "damage_i"] = acc } }# Sum the increments cumulatively to get accumulated heat damage through time output$cum_damage[1] <- output$damage_i[1]for (j in2:nrow(output)) { output$cum_damage[j] <- output$cum_damage[j -1] + output$damage_i[j] }# Convert HI to survival probability at each time step given a threshold value used (50%, or 80% usually) output$mortality_prob = (output$cum_damage * threshold) /100# Now check that mortality cum damage and probabilities do not exceed 1 (100%) output$mortality_prob[output$mortality_prob >1] <-1 output$cum_damage[output$cum_damage >100] <-100 output$temp <- temp[1:(length(time) -1)] output$time_min <- time[1:(length(time) -1)]return(output)}# Function source: Noble et al. 2025 https://doi.org/10.32942/X2BP8N # Computes a temperature-dependent rate using an Arrhenius rise with both low- and high-temperature inactivation terms (Schoolfield form).damage_repair_accum <-function(time, temp, intercept, slope, threshold =0.5, Tc, TA, TAL, TAH, TL, TH, TREF, kdot_ref) {if (!is.numeric(temp)) stop("'temp' must be numeric (°C).")for (nm inc("TA", "TAL", "TAH", "TL", "TH", "TREF", "kdot_ref")) {if (!is.numeric(get(nm)) ||length(get(nm)) !=1L) {stop(sprintf("'%s' must be a numeric scalar.", nm)) } }## Convert input temperatures to Kelvin TK <-273.15+ temp # evaluation temperatures (K)# Scalars below are already in Kelvin per argument definitions:# TL, TH, TREF# Step 1: Get the damage without repair for each time step damage <-damage_accum(time = time, temp = temp, intercept = intercept, slope = slope, threshold = threshold, Tc = Tc)# Step 2: Calculate how much repair occurs at each time step based on temperature damage_repair_i <-c()for(i in1:(length(TK)-1)){ repair_rate <- (exp(TA*(1/TREF-1/max(TK[i +1], TK[i]))/(1+exp(TAL*(1/max(TK[i +1], TK[i])-1/TL))+exp(TAH*(1/TH-1/max(TK[i +1], TK[i])))))*kdot_ref) # per minute repair rate damage$damage_repair_i[i] <- (repair_rate * (time[i+1] - time[i])) # time in minutes on original scale }# Step 3: Calculate net damage after repair at each time step. If repair produces a negative value, reset to 0 because there is no damage to repair if it's negative.for (j in1:nrow(damage)) { net_damage_repair_i <- damage$damage_i[j] - damage$damage_repair_i[j] damage$net_damage_repair_i[j] <- net_damage_repair_i }# Step 4: Sum the increments cumulatively to get accumulated heat damage through time. Always reset to 0 of the repair produces a negative value because there is no damage to repair if it's negative. This prevents a buildup of negative values over time. damage$cum_damage_repair[1] <- damage$net_damage_repair_i[1]for (j in2:nrow(damage)) { cum_damage_repair <- damage$cum_damage_repair[j -1] + damage$net_damage_repair_i[j] cum_damage_repair <-ifelse(cum_damage_repair <0, 0, cum_damage_repair) damage$cum_damage_repair[j] <- cum_damage_repair }# Convert heat injury to survival probability at each time step given a threshold value used damage$mortality_prob = (damage$cum_damage_repair * threshold) /100# Now check that mortality cum damage and probabilities do not exceed 1 (100%) damage$mortality_prob[damage$mortality_prob >1] <-1 damage$cum_damage_repair[damage$cum_damage_repair >100] <-100# Set all values after the first instance of 100% mortality to 100% mortality first_max_index <-which(damage$cum_damage_repair ==100)[1]if (!is.na(first_max_index)) { damage$cum_damage_repair[first_max_index:nrow(damage)] <-100 damage$mortality_prob[first_max_index:nrow(damage)] <-1 }return(damage)}# Function Source: Faber, Ørsted, Ehlers (2024) https://doi.org/10.1093/jxb/erae096 # Calculate the accumulated damage based on time in minutesAcc_function <-function(x, y, intercept, slope) { output <-data.frame(z =NA, cumsum =NA) last_row <-nrow(data.frame(x = x, y = y)) x <-c(x, x[last_row]) y <-c(y, y[last_row])for (i in1:length(x) -1) {#print(paste0("round nr. ", i)) acc <-max(0, min(100, (100* (x[i +1] - x[i])) / (10^(slope *max(y[i +1], y[i]) + intercept)))) output[i, 1] <- acc } output$cumsum <-cumsum(output$z) z <- output$cumsum z <-append(z, 0) z[z >100] <-100# Set values exceeding 100% to 100% z <- z[-length(z)]return(z)}# Function source: Arnold et al. 2025 Global Change Biology https://doi.org/10.1111/gcb.70315# Arrhenius repair function ArrFunc5 developed by Michael KearneyArrFunc5 <-function(x,TA,TAL,TAH,TL,TH,TREF,kdot_ref) {exp(TA * (1/TREF-1/(273.15+x)))/(1+exp(TAL * (1/(273.15+x)-1/TL)) +exp(TAH * (1/TH-1/(273.15+x)))) * kdot_ref}# Helper function to run damage_repair_accum_state simulation for an individual.run_individual_sim <-function(i, alpha_vec, beta_vec) { result <-damage_repair_accum_state(time = times, temp = temps.min, intercept = alpha_vec[i],slope = beta_vec[i], threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref) result$individual <- i result$alpha <- alpha_vec[i] result$beta <- beta_vec[i]return(result)}# Helper function to generate individual parameters for a life stage using alpha/beta variationgenerate_individuals <-function(stage_params, seed =123) {set.seed(seed) n <- stage_params$n_individuals# Calculate alpha and beta from ctmax and z alpha_mean <--stage_params$ctmax / stage_params$z beta_mean <-1/ stage_params$z# Generate alpha and beta with variation alpha_values <-rnorm(n, mean = alpha_mean, sd = stage_params$alpha_sd) beta_values <-rnorm(n, mean = beta_mean, sd = stage_params$beta_sd)# Back-transform to ctmax and z z_bt <-1/ beta_values ctmax_bt <--alpha_values * z_btdata.frame(individual =1:n,life_stage = stage_params$name,alpha = alpha_values,beta = beta_values,z = z_bt,ctmax = ctmax_bt)}# Helper function to calculate TDT (Thermal Death Time) curve for one individualcalculate_tdt <-function(ctmax, z, time_range =c(1, 1000)) { time2 <-log10(seq(time_range[1], time_range[2], 1)) ht <- ctmax + (z * time2)data.frame(log_time = time2,time_min =10^time2,temp = ht)}# Helper function to calculate damage and repair rates across life stages with specified and shared parameterscalculate_damage_repair <-function(stage_params, shared_params) {# Temperature sequence Tbs <-seq(5, 60, 0.05501)# Generate Arrhenius functions for repair rates TempCorr0 <-ArrFunc5(Tbs, shared_params$TA, shared_params$TAL, shared_params$TAH, shared_params$TL, shared_params$TH, shared_params$TREF, shared_params$kdot0) TempCorr1 <-ArrFunc5(Tbs, shared_params$TA, shared_params$TAL, shared_params$TAH, shared_params$TL, shared_params$TH, shared_params$TREF, shared_params$kdot1)# Calculate damage rate using mean alpha and beta alpha_mean <--stage_params$ctmax / stage_params$z beta_mean <-1/ stage_params$z time_seq <-seq(1, 1000, 1)# Calculate accumulation of damage damageacc <-Acc_function(x = time_seq, y = Tbs, intercept = alpha_mean, slope = beta_mean) damage <-c(diff(damageacc), NA) damage[which.max(damage):length(damage)] <-NA# Cap damage beyond maxdata.frame(Tbs = Tbs,TempCorr0 = TempCorr0,TempCorr1 = TempCorr1,damage = damage,life_stage = stage_params$name)}
3) Extract microclimate data
Extract microclimate data for point of interest with specified longitude and latitude. Run once and save object as it can take a while to run. Requires email address to access micro_silo() database.
Code
# Extract microclimate data for points of interest. Run once for now and save objectpoints <-c(150.737902, -33.618746)# 50 cm height aboveground for pollen stage # locs50cm <- NicheMapR::micro_silo(loc = c(points[1], points[2]), email = "INSERT_EMAIL_ADDRESS",# Usrhyt = 0.5, # 50cm for pollen# minshade = 0,# dstart = "01/01/2017", dfinish = "31/05/2020")# saveRDS(locs50cm, "./output/micro_locs50cm.RDS")# 0 cm height (soil surface, for imbibed seeds) is included in these microclimate simulations# 1 cm height aboveground for seedling stage # locs1cm <- NicheMapR::micro_silo(loc = c(points[1], points[2]), email = "INSERT_EMAIL_ADDRESS",# Usrhyt = 0.01, # 1cm for seedling# minshade = 0,# dstart = "01/01/2017", dfinish = "31/05/2020") # saveRDS(locs1cm, "./output/micro_locs1cm.RDS")# 1.2 m height aboveground for air.# locs1.2m <- NicheMapR::micro_silo(loc = c(points[1], points[2]), email = "INSERT_EMAIL_ADDRESS",# Usrhyt = 1.2, # 1.2m for air# minshade = 0,# dstart = "01/01/2017", dfinish = "31/05/2020") # saveRDS(locs1.2m, "./output/micro_locs1.2m.RDS")# Read in the RDS objects containing the microclimatemicro1cm <-readRDS("./output/micro_locs1cm.RDS")summary(micro1cm)
# Save meteorological, soil, and plant components of the microclimate simulation as separate data framesmetout1cm <-as.data.frame(micro1cm$metout)soilpot1cm <-as.data.frame(micro1cm$soilpot)plant1cm <-as.data.frame(micro1cm$plant)metout50cm <-as.data.frame(micro50cm$metout)soilpot50cm <-as.data.frame(micro50cm$soilpot)plant50cm <-as.data.frame(micro50cm$plant)metout1.2m <-as.data.frame(micro1.2m$metout)soilpot1.2m <-as.data.frame(micro1.2m$soilpot)plant1.2m <-as.data.frame(micro1.2m$plant)
4) Run NicheMapR ectotherm models
Run NicheMapR ectotherm models in leaf mode to simulate how the microclimate environment interacts biophysically with the organism at different life stages.
Code
# Set up NicheMapR ectotherm model in leaf mode for the different plant life stages # Tutorial for leaf mode uses of ectotherm model is available here: https://mrke.github.io/models/Plant-Models# And full publication Kearney & Leigh (2024) Methods in Ecology and Evolution# https://doi.org/10.1111/2041-210X.14373# Seed functional traitsseed_Ww_g <-0.3# wet weight, gseed_shape <-1# 0=plate, 1=cylinder, 2=ellipsoidseed_alpha_L <-0.5# solar absorptivity of the seed (-)seed_epsilon_L <-0.97# emissivity of the seed (-)seed_postur <-3# orientation to sun, 1=perpendicular, 2=parallel, 3=horizontal, 0=intermediate between perpendicular and parallelseed_fatosk <-0.5# radiation configuration factor to sky (-)seed_fatosb <-0.5# radiation configuration factor to substrateseed_g_vs_ab <-0.0# leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s seed_g_vs_ad <-0.0# leaf vapour conductance, adaxial (top of leaf), mol/m2/s# Seedling functional traitsseedling_Ww_g <-0.5# wet weight, gseedling_shape <-2# 0=plate, 1=cylinder, 2=ellipsoidseedling_shape_b <-0.0025# ratio of b axis:a axis for ellipsoidseedling_shape_c <-0.1176# ratio of c axis:a axis for ellipsoidseedling_alpha_L <-0.5# solar absorptivity of the leaf (-)seedling_epsilon_L <-0.97# emissivity of the leaf (-)seedling_postur <-0# orientation to sun, 1=perpendicular, 2=parallel, 3=horizontal, 0=intermediate between perpendicular and parallelseedling_fatosk <-0.5# radiation configuration factor to sky (-)seedling_fatosb <-0.5# radiation configuration factor to substrateseedling_g_vs_ab <-0.2# leaf vapour conductance, abaxial (bottom of leaf), mol/m2/sseedling_g_vs_ad <-0.0# leaf vapour conductance, adaxial (top of leaf), mol/m2/s# Pollen (i.e. flower) functional traitspollen_Ww_g <-1# wet weight, gpollen_shape <-0# 0=plate, 1=cylinder, 2=ellipsoidpollen_alpha_L <-0.7# solar absorptivity of the flower (-)pollen_epsilon_L <-0.97# emissivity of the flower (-)pollen_postur <-0# orientation to sun, 1=perpendicular, 2=parallel, 3=horizontal, 0=intermediate between perpendicular and parallelpollen_fatosk <-0.5# radiation configuration factor to sky (-)pollen_fatosb <-0.5# radiation configuration factor to substratepollen_g_vs_ab <-0.0# leaf vapour conductance, abaxial (bottom of leaf), mol/m2/s NApollen_g_vs_ad <-0.0# leaf vapour conductance, adaxial (top of leaf), mol/m2/s NA# Check shade time-series# par(mfrow = c(2,1))# plot(micro1cm$shadplant[,3], type = "l")# plot(micro50cm$shadplant[,3], type = "l")micro <- micro1cm# Simulate leaf (seed at 1cm height) with ectotherm modelseed_ecto1cm <-ectotherm(leaf =1, # turn on leaf modelive =0, # no behavioural thermoregulationpct_cond =0, # no conductionWw_g = seed_Ww_g,shape = seed_shape,alpha_max = seed_alpha_L,alpha_min = seed_alpha_L,postur = seed_postur,fatosk = seed_fatosk,fatosb = seed_fatosb,epsilon = seed_epsilon_L,g_vs_ab = seed_g_vs_ab,g_vs_ad = seed_g_vs_ad,metout = micro$metout,soil = micro$soil, soilpot = micro$soilpot,soilmoist = micro$soilmoist, humid = micro$humid,tcond = micro$tcond, shadmet = micro$shadmet,shadsoil = micro$shadsoil, shadpot = micro$shadpot,shadmoist = micro$shadmoist,shadhumid = micro$shadhumid,shadtcond = micro$shadtcond)# Simulate leaf (seedling at 1cm height) with ectotherm modelseedling_ecto1cm <-ectotherm(leaf =1, # turn on leaf modelive =0, # no behavioural thermoregulationpct_cond =0, # no conductionWw_g = seedling_Ww_g,shape = seedling_shape,shape_b = seedling_shape_b,shape_c = seedling_shape_c,alpha_max = seedling_alpha_L,alpha_min = seedling_alpha_L,postur = seedling_postur,fatosk = seedling_fatosk,fatosb = seedling_fatosb,epsilon = seedling_epsilon_L,g_vs_ab = seedling_g_vs_ab,g_vs_ad = seedling_g_vs_ad,metout = micro$metout,soil = micro$soil, soilpot = micro$soilpot,soilmoist = micro$soilmoist, humid = micro$humid,tcond = micro$tcond, shadmet = micro$shadmet,shadsoil = micro$shadsoil, shadpot = micro$shadpot,shadmoist = micro$shadmoist,shadhumid = micro$shadhumid,shadtcond = micro$shadtcond)# Simulate pollen with ectotherm modelmicro <- micro50cmpollen_ecto50cm <-ectotherm(leaf =1, # turn on leaf modelive =0, # no behavioural thermoregulationpct_cond =0, # no conductionWw_g = pollen_Ww_g,shape = pollen_shape,alpha_max = pollen_alpha_L,alpha_min = pollen_alpha_L,postur = pollen_postur,fatosk = pollen_fatosk,fatosb = pollen_fatosb,epsilon = pollen_epsilon_L,g_vs_ab = pollen_g_vs_ab,g_vs_ad = pollen_g_vs_ad,metout = micro$metout,soil = micro$soil, soilpot = micro$soilpot,soilmoist = micro$soilmoist, humid = micro$humid,tcond = micro$tcond, shadmet = micro$shadmet,shadsoil = micro$shadsoil, shadpot = micro$shadpot,shadmoist = micro$shadmoist,shadhumid = micro$shadhumid,shadtcond = micro$shadtcond,)
5) Plot microclimate data
Plot the microclimate conditions at the scale relevant to the different life stages for a set range of dates.
Code
# We can calculate VPD in addition to temperature if that is of interest.# VPD equation from: https://www.omnicalculator.com/biology/vapor-pressure-deficit # Define date rangedate_min <-"2018-09-20"date_max <-as.Date(date_min) +10# Partition microenvironments to life stages# Seeds (0cm height)environ_seed <-as.data.frame(seed_ecto1cm$soil)environ_seed$dates <- micro1cm$datesenviron_seed$vpd_1cm <- (0.61078*exp(17.2694* micro1cm$metout[,3] / (micro1cm$metout[,3] +237.3))) * (1- micro1cm$metout[,5]/100)environ_seed$soilmoist0cm <- micro1cm$soilmoist[,3] # 0cm soil depth (surface)# Seedlings (1cm height)environ_seedling <-as.data.frame(seedling_ecto1cm$environ)environ_seedling$dates <- micro1cm$datesenviron_seedling$vpd_1cm <- (0.61078*exp(17.2694* micro1cm$metout[,3] / (micro1cm$metout[,3] +237.3))) * (1- micro1cm$metout[,5]/100)environ_seedling$soilmoist2.5cm <- micro1cm$soilmoist[,4] # 2.5cm soil depthenviron_seedling$Tair <- micro1.2m$metout[,"TALOC"]# Pollen (50cm height)environ_pollen <-as.data.frame(pollen_ecto50cm$environ)environ_pollen$dates <- micro50cm$datesenviron_pollen$vpd_50cm <- (0.61078*exp(17.2694* micro50cm$metout[,3] / (micro50cm$metout[,3] +237.3))) * (1- micro50cm$metout[,5]/100)environ_pollen$soilmoist50cm <- micro50cm$soilmoist[,10] # 50cm soil depth, not used further in this simulation# Rainfall data for the date rangerain_seedling <-as.data.frame(seedling_ecto1cm$rainfall)colnames(rain_seedling) <-"rainfall"rain_seedling$dates <- micro1cm$dates2rain_seedling$dates <-as.POSIXct(paste0(micro1cm$dates2, " 00:00:00"), format ="%Y-%m-%d %H:%M:%S")rain_seedling2 <- rain_seedling %>%filter(dates > date_min & dates < date_max)# Filter environmental data to the relevant temperatures and rename themenviron_seed2 <- environ_seed %>%filter(dates > date_min & dates < date_max) %>% dplyr::select(dates, D0cm, D2.5cm) %>%pivot_longer(cols =c(D0cm, D2.5cm),names_to ="Temperatures",values_to ="Temperature") %>%mutate(Temperatures =ifelse(Temperatures =="D0cm", "Seed (0 cm)", "Seed (-2.5 cm)"))# Filter environmental data to the relevant temperatures and rename themenviron_seedling2 <- environ_seedling %>%filter(dates > date_min & dates < date_max) %>% dplyr::select(dates, TA, TC, Tair) %>%pivot_longer(cols =c(TA, TC, Tair),names_to ="Temperatures",values_to ="Temperature") %>%mutate(Temperatures =case_when( Temperatures =="TA"~"Air (1 cm)", Temperatures =="TC"~"Seedling (1 cm)", Temperatures =="Tair"~"Air (1.2 m)",TRUE~ Temperatures))# Filter environmental data to the relevant temperatures and rename themenviron_pollen2 <- environ_pollen %>%filter(dates > date_min & dates < date_max) %>% dplyr::select(dates, TA, TC) %>%pivot_longer(cols =c(TA,TC),names_to ="Temperatures",values_to ="Temperature") %>%mutate(Temperatures =ifelse(Temperatures =="TA", "Air (50 cm)", "Pollen (50 cm)"))# VPD at different heightsvpd_seedling2 <- environ_seedling %>%filter(dates > date_min & dates < date_max) %>% dplyr::select(dates, vpd_1cm) %>%pivot_longer(cols =c(vpd_1cm),names_to ="VPDs",values_to ="VPD") %>%mutate(VPDs =ifelse(VPDs =="vpd_1cm", "VPD seedling (1 cm)", "VPD pollen (50 cm)"))# Temp and rain plottemp_rain_plot <-ggplot() +geom_col(data = rain_seedling2, aes(y = rainfall*5, x = dates), fill ="skyblue2") +geom_line(data =subset(environ_seedling2, Temperatures=="Air (1.2 m)"), aes(y = Temperature, x = dates, col = Temperatures), show.legend =TRUE) +geom_line(data =subset(environ_seed2, Temperatures=="Seed (0 cm)"), aes(y = Temperature, x = dates, col = Temperatures), show.legend =TRUE) +geom_line(data =subset(environ_seedling2, Temperatures=="Seedling (1 cm)"), aes(y = Temperature, x = dates, col = Temperatures), show.legend =TRUE) +geom_line(data =subset(environ_pollen2, Temperatures=="Pollen (50 cm)"), aes(y = Temperature, x = dates, col = Temperatures), show.legend =TRUE) +labs(x ="") +scale_colour_manual(values = cbp2[c(1,8,2,4)], name ="", labels =c(bquote(italic(T)[air]~(1.2~m)), bquote(italic(T)[pollen]~(50~cm)),bquote(italic(T)[seed]~(0~cm)),bquote(italic(T)[seedling]~(1~cm)))) +scale_y_continuous(name =bquote(Temperature~(degree*C)), breaks =c(0,10,20,30,40,50,60),sec.axis =sec_axis(~./5, name ="Rainfall (mm)")) +scale_x_datetime(date_breaks = ("2 days"), date_labels ="%d %b") +theme_bw() +theme(legend.position ="inside",legend.position.inside =c(0.6,0.8),legend.direction ="vertical",legend.background =element_rect(fill ="transparent", colour =NA)) +guides(colour =guide_legend(nrow =3, byrow =TRUE))temp_rain_plot
6) Damage accumulation and repair
For demonstration purposes, run a damage accumulation and damage-repair function across temperatures.
Code
# Universal repair settings# Set Arrhenius model and repair rate parameters TA <-14065# Arrhenius temperature, KTL <-10.5+273.15# Arrhenius temperature lower threshold, KTH <-28.5+273.15# Arrhenius temperature upper threshold, KTAL <-50000# Arrhenius temperature lower, KTAH <-100000# Arrhenius temperature upper, KTREF <-20+273.15# reference temperature, Kkdot_ref <-0.05# repair rate per minute at ref temp 20 deg C# Temperature-time model sequencectmax <-42.5# critical temperature for 1 min durationz <--2.5# thermal sensitivity (decline in tolerated temperature for each 10-fold increase in duration)time2 <-log10(seq(1,1000,1))ht <- ctmax + (z*time2)sim <-data.frame(cbind(time2,ht))Tbs <-seq(5, 60, 0.05501)kdot0 <-0kdot1 <-0.05# Generate Arrhenius functions for repair ratesTempCorr0 <-ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot0)TempCorr1 <-ArrFunc5(Tbs, TA, TAL, TAH, TL, TH, TREF, kdot1)# Accumulated damage function, then convert to a damage ratedamage <-NULLdamageacc <-Acc_function(x =10^time2, y = Tbs, intercept =-(ctmax/z), slope = (1/z))for (i in1:length(damageacc)) { damage[i] <- damageacc[i+1] - damageacc[i] }damage[which.max(damage):1000] <-NA# cap damage rate beyond accumulated maximum# Format dataframe for plottingTempCorrDat <-data.frame(cbind(Tbs, TempCorr0, TempCorr1, damage))TempCorrDatL <- TempCorrDat %>%pivot_longer(cols =c(TempCorr0, TempCorr1), names_to =c("repair"), values_to ="repair_value")TempCorrDatL <-subset(TempCorrDatL, Tbs >6)TempCorrDatLrepair <-subset(TempCorrDatL, repair!="TempCorr0") time2 <-log10(seq(1,1000,1))htdat <-do.call(rbind, lapply(ht, data.frame))sim2 <-data.frame(cbind(time2,htdat))colnames(sim2) <-c("time", "temp")# Plot net damage modelnet_plot <-ggplot(data.frame(TempCorrDatL), aes(x = Tbs, y = damage-repair_value, colour = repair)) +geom_path(aes(group =rev(repair))) +geom_hline(yintercept =0, linetype =2) +scale_x_continuous(limits =c(5,50)) +scale_y_continuous(limits =c(-0.2,0.2)) +scale_colour_manual(name ="", values =c("#E69F00", "#009E73"), labels =c("No repair", "Repair")) +annotate("text", x =6, y =-0.05, label ="Repair \noutweighs \ndamage", hjust =0) +annotate("text", x =40, y =0.05, label ="Damage \noutweighs \nrepair", hjust =0) +annotate("text", x =18, y =-0.17, label ="Permissive", hjust =0) +annotate("text", x =36, y =-0.17, label ="Stressful", hjust =0) +geom_segment(aes(x =30, y =-0.15, xend =20, yend =-0.15),arrow =arrow(length =unit(0.2, "cm"), type ="closed"), colour ="#56B4E9") +geom_segment(aes(x =32, y =-0.15, xend =44, yend =-0.15),arrow =arrow(length =unit(0.2, "cm"), type ="closed"), colour ="#D55E00") +geom_segment(aes(x =36, y =0.15, xend =36.8, yend =0.2),arrow =arrow(length =unit(0.2, "cm"), type ="closed"), colour ="black") +geom_segment(x =32, xend =32, y =-0.25, yend =0, linetype =2, colour ="black") +labs(y =bquote(Net~damage~rate~("%"~min^-1)), x =bquote(Temperature~(degree*C))) +theme_bw() +theme(legend.position ="inside", legend.position.inside =c(0.2, 0.8),legend.background =element_rect(fill ="transparent", colour =NA))net_plot
7) TLS for pollen
Simulate thermal load sensitivity for pollen (single individual)
Code
# Pollen parameters for modellingpar(mfrow =c(1,1))T_pollen <- environ_pollen2 %>% dplyr::filter(Temperatures=="Pollen (50 cm)")ta_pollen <- T_pollen$Temperaturepollen_temps.min <-spline(ta_pollen, n =length(ta_pollen-1)*60)$y pollen_times <-1:length(pollen_temps.min)plot(pollen_temps.min, type ="l")
Code
pollen_ctmax <-42.5# critical temperature for 1 min durationpollen_z <--3# thermal sensitivity (decline in tolerated temperature for each 10-fold increase in duration)pollen_alpha <--pollen_ctmax/pollen_z # transformation to alpha when linear model has time on y-axis, for damage accumulationpollen_beta <-1/ pollen_z # transformation to beta when linear model has time on y-axis, for damage accumulationpollen_dam_sim <-damage_repair_accum_state(time = pollen_times, temp = pollen_temps.min, intercept = pollen_alpha,slope = pollen_beta, threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref)par(mfrow =c(3,1))plot(pollen_dam_sim$temp, type ="l")plot(pollen_dam_sim$net_damage_inc, type ="l")plot(pollen_dam_sim$mortality_prob, type ="l")
8) TLS pollen population
Apply the TLS simulation for pollen across a ‘population’
Code
# Setup parametersn_individuals <-100pollen_ctmax <-42.5pollen_z <--3pollen_alpha <--pollen_ctmax / pollen_zpollen_beta <-1/ pollen_zpollen_alpha_sd <-0.07# standard deviation for alpha simulationpollen_beta_sd <-0.05# standard deviation for beta simulationset.seed(123)pollen_alpha_values <-rnorm(n_individuals, mean = pollen_alpha, sd = pollen_alpha_sd)set.seed(123)pollen_beta_values <-rnorm(n_individuals, mean = pollen_beta, sd = pollen_beta_sd)pollen_bt_z <- (1/ pollen_beta_values) # check back-transformed values to see if sensiblepollen_bt_ct <--pollen_alpha_values * pollen_bt_z # check back-transformed values to see if sensibletemps.min <- pollen_temps.mintimes <- pollen_times# Run in paralleln_cores <-detectCores() -1cl <-makeCluster(n_cores)# Export necessary objects to clusterclusterExport(cl, c("damage_repair_accum_state", "times", "temps.min","pollen_alpha_values", "pollen_beta_values", "TA", "TAL", "TAH", "TL", "TH", "TREF", "kdot_ref"))# Runpollen_results <-parLapply(cl, 1:n_individuals, run_individual_sim,alpha_vec = pollen_alpha_values, beta_vec = pollen_beta_values)stopCluster(cl)# Combine resultspollen_data_all <-bind_rows(pollen_results)# Or for large n, you might want to store summaries rather than full time seriespollen_results_summary <-data.frame(individual =integer(),alpha =numeric(),beta =numeric(),max_mortality =numeric(),time_to_50pct_mortality =numeric(),final_cum_damage =numeric())for(i in1:n_individuals) { result <-damage_repair_accum_state(time = pollen_times, temp = pollen_temps.min, intercept = pollen_alpha_values[i],slope = pollen_beta_values[i], threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref )# Calculate summary statistics max_mort <-max(result$mortality_prob, na.rm =TRUE) time_50 <-if(any(result$mortality_prob >=0.5)) { result$time_min[which(result$mortality_prob >=0.5)[1]] } else {NA } final_damage <-tail(result$cum_damage_repair, 1) pollen_results_summary <-rbind(pollen_results_summary, data.frame(individual = i,alpha = pollen_alpha_values[i],beta = pollen_beta_values[i],pollen_z = (1/ pollen_beta_values[i]),pollen_ct =-pollen_alpha_values[i] * pollen_z,max_mortality = max_mort,time_to_50pct_mortality = time_50,final_cum_damage = final_damage ))if(i %%100==0) cat(sprintf("Completed %d/%d\n", i, n_individuals))}
Completed 100/100
Code
# pollen_results_summary
9) Plot TLS for pollen
Plot pollen TLS simulation
Code
# Plot damage/repair trajectories for all individualspollen_dam_plot <-ggplot(pollen_data_all, aes(x = time_min, y = net_damage_inc, group = individual)) +geom_line(alpha =0.2, color = cbp2[8]) +geom_hline(yintercept =0, colour ="black", linetype =2) +coord_cartesian(ylim =c(-0.2, 5)) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Net damage") +theme_bw()# Plot Mortality trajectories for all individualspollen_mort_plot <-ggplot(pollen_data_all, aes(x = time_min, y = mortality_prob, group = individual)) +geom_line(alpha =0.2, color = cbp2[8]) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Mortality probability") +theme_bw()pollen_plots <-plot_grid(pollen_dam_plot, pollen_mort_plot, ncol =1, nrow =2, rel_heights =c(1, 1), align ="v", labels =c("a", "b"))pollen_plots
10) TLS for seed
Simulate thermal load sensitivity for seed (single individual)
Code
# seed parameters for modellingpar(mfrow =c(1,1))T_seed <- environ_seed2 %>% dplyr::filter(Temperatures=="Seed (0 cm)")ta_seed <- T_seed$Temperatureseed_temps.min <-spline(ta_seed, n =length(ta_seed-1)*60)$y seed_times <-1:length(seed_temps.min)plot(seed_temps.min, type ="l")
Code
seed_ctmax <-48.5# critical temperature for 1 min durationseed_z <--2# thermal sensitivity (decline in tolerated temperature for each 10-fold increase in duration)seed_alpha <--seed_ctmax/seed_z # transformation to alpha when linear model has time on y-axis, for damage accumulationseed_beta <-1/ seed_z # transformation to beta when linear model has time on y-axis, for damage accumulationseed_dam_sim <-damage_repair_accum_state(time = seed_times, temp = seed_temps.min, intercept = seed_alpha,slope = seed_beta, threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref)par(mfrow =c(3,1))plot(seed_dam_sim$temp, type ="l")plot(seed_dam_sim$net_damage_inc, type ="l")plot(seed_dam_sim$mortality_prob, type ="l")
11) TLS for seed population
Apply the TLS simulation for seed across a ‘population’
Code
# Setup parametersn_individuals <-100seed_ctmax <-54seed_z <--5seed_alpha <--seed_ctmax / seed_zseed_beta <-1/ seed_zseed_alpha_sd <-0.07# standard deviation for alpha simulationseed_beta_sd <-0.05# standard deviation for beta simulationset.seed(123)seed_alpha_values <-rnorm(n_individuals, mean = seed_alpha, sd = seed_alpha_sd)set.seed(123)seed_beta_values <-rnorm(n_individuals, mean = seed_beta, sd = seed_beta_sd)seed_bt_z <- (1/ seed_beta_values) # check back-transformed values to see if sensibleseed_bt_ct <--seed_alpha_values * seed_bt_z # check back-transformed values to see if sensibletemps.min <- seed_temps.mintimes <- seed_times# Run in paralleln_cores <-detectCores() -1cl <-makeCluster(n_cores)# Export necessary objects to clusterclusterExport(cl, c("damage_repair_accum_state", "times", "temps.min","seed_alpha_values", "seed_beta_values", "TA", "TAL", "TAH", "TL", "TH", "TREF", "kdot_ref"))# Runseed_results <-parLapply(cl, 1:n_individuals, run_individual_sim,alpha_vec = seed_alpha_values, beta_vec = seed_beta_values)stopCluster(cl)# Combine resultsseed_data_all <-bind_rows(seed_results)# Or for large n, you might want to store summaries rather than full time seriesseed_results_summary <-data.frame(individual =integer(),alpha =numeric(),beta =numeric(),max_mortality =numeric(),time_to_50pct_mortality =numeric(),final_cum_damage =numeric())for(i in1:n_individuals) { result <-damage_repair_accum_state(time = seed_times, temp = seed_temps.min, intercept = seed_alpha_values[i],slope = seed_beta_values[i], threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref )# Calculate summary statistics max_mort <-max(result$mortality_prob, na.rm =TRUE) time_50 <-if(any(result$mortality_prob >=0.5)) { result$time_min[which(result$mortality_prob >=0.5)[1]] } else {NA } final_damage <-tail(result$cum_damage_repair, 1) seed_results_summary <-rbind(seed_results_summary, data.frame(individual = i,alpha = seed_alpha_values[i],beta = seed_beta_values[i],seed_z = (1/ seed_beta_values[i]),seed_ct =-seed_alpha_values[i] * seed_z,max_mortality = max_mort,time_to_50pct_mortality = time_50,final_cum_damage = final_damage ))if(i %%100==0) cat(sprintf("Completed %d/%d\n", i, n_individuals))}
Completed 100/100
Code
# seed_results_summary
12) Plot TLS for seed
Plot seed TLS simulation
Code
# Plot damage/repair trajectories for all individualsseed_dam_plot <-ggplot(seed_data_all, aes(x = time_min, y = net_damage_inc, group = individual)) +geom_line(alpha =0.2, color = cbp2[2]) +geom_hline(yintercept =0, colour ="black", linetype =2) +coord_cartesian(ylim =c(-0.2, 5)) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Net damage") +theme_bw()# Plot Mortality trajectories for all individualsseed_mort_plot <-ggplot(seed_data_all, aes(x = time_min, y = mortality_prob, group = individual)) +geom_line(alpha =0.2, color = cbp2[2]) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Mortality probability") +theme_bw()seed_plots <-plot_grid(seed_dam_plot, seed_mort_plot, ncol =1, nrow =2, rel_heights =c(1, 1), align ="v", labels =c("a", "b"))seed_plots
13) TLS for seedling
Simulate thermal load sensitivity for seedling (single individual)
Code
# Seedling parameters for modellingpar(mfrow =c(1,1))T_seedling <- environ_seedling2 %>% dplyr::filter(Temperatures=="Seedling (1 cm)")ta_seedling <- T_seedling$Temperatureseedling_temps.min <-spline(ta_seedling, n =length(ta_seedling-1)*60)$y seedling_times <-1:length(seedling_temps.min)plot(seedling_temps.min, type ="l")
Code
seedling_ctmax <-47.5# critical temperature for 1 min durationseedling_z <--4.3# thermal sensitivity (decline in tolerated temperature for each 10-fold increase in duration)seedling_alpha <--seedling_ctmax/seedling_z # transformation to alpha when linear model has time on y-axis, for damage accumulationseedling_beta <-1/ seedling_z # transformation to beta when linear model has time on y-axis, for damage accumulationseedling_dam_sim <-damage_repair_accum_state(time = seedling_times, temp = seedling_temps.min, intercept = seedling_alpha,slope = seedling_beta, threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref)par(mfrow =c(3,1))plot(seedling_dam_sim$temp, type ="l")plot(seedling_dam_sim$net_damage_inc, type ="l")plot(seedling_dam_sim$mortality_prob, type ="l")
14) TLS for seed population
pply the TLS simulation for seedling across a ‘population’
Code
# Setup parametersn_individuals <-100seedling_ctmax <-47.5seedling_z <--4.3seedling_alpha <--seedling_ctmax / seedling_zseedling_beta <-1/ seedling_zseedling_alpha_sd <-0.06# standard deviation for alpha simulationseedling_beta_sd <-0.03# standard deviation for beta simulationset.seed(123)seedling_alpha_values <-rnorm(n_individuals, mean = seedling_alpha, sd = seedling_alpha_sd)set.seed(123)seedling_beta_values <-rnorm(n_individuals, mean = seedling_beta, sd = seedling_beta_sd)seedling_bt_z <- (1/ seedling_beta_values) # check back-transformed values to see if sensibleseedling_bt_ct <--seedling_alpha_values * seedling_bt_z # check back-transformed values to see if sensibletemps.min <- seedling_temps.mintimes <- seedling_times# Run in paralleln_cores <-detectCores() -1cl <-makeCluster(n_cores)# Export necessary objects to clusterclusterExport(cl, c("damage_repair_accum_state", "times", "temps.min","seedling_alpha_values", "seedling_beta_values", "TA", "TAL", "TAH", "TL", "TH", "TREF", "kdot_ref"))# Runseedling_results <-parLapply(cl, 1:n_individuals, run_individual_sim,alpha_vec = seedling_alpha_values, beta_vec = seedling_beta_values)stopCluster(cl)# Combine resultsseedling_data_all <-bind_rows(seedling_results)# Or for large n, you might want to store summaries rather than full time seriesseedling_results_summary <-data.frame(individual =integer(),alpha =numeric(),beta =numeric(),max_mortality =numeric(),time_to_50pct_mortality =numeric(),final_cum_damage =numeric())for(i in1:n_individuals) { result <-damage_repair_accum_state(time = seedling_times, temp = seedling_temps.min, intercept = seedling_alpha_values[i],slope = seedling_beta_values[i], threshold =0.5, TA = TA, TAL = TAL, TAH = TAH, TL = TL, TH = TH, TREF = TREF, kdot_ref = kdot_ref )# Calculate summary statistics max_mort <-max(result$mortality_prob, na.rm =TRUE) time_50 <-if(any(result$mortality_prob >=0.5)) { result$time_min[which(result$mortality_prob >=0.5)[1]] } else {NA } final_damage <-tail(result$cum_damage_repair, 1) seedling_results_summary <-rbind(seedling_results_summary, data.frame(individual = i,alpha = seedling_alpha_values[i],beta = seedling_beta_values[i],seedling_z = (1/ seedling_beta_values[i]),seedling_ct =-seedling_alpha_values[i] * seedling_z,max_mortality = max_mort,time_to_50pct_mortality = time_50,final_cum_damage = final_damage ))if(i %%100==0) cat(sprintf("Completed %d/%d\n", i, n_individuals))}
Completed 100/100
Code
#seedling_results_summary
15) Plot TLS for seedling
Plot seedling TLS simulation
Code
# Plot damage/repair trajectories for all individualsseedling_dam_plot <-ggplot(seedling_data_all, aes(x = time_min, y = net_damage_inc, group = individual)) +geom_line(alpha =0.2, color = cbp2[4]) +geom_hline(yintercept =0, colour ="black", linetype =2) +coord_cartesian(ylim =c(-0.2, 5)) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Net damage") +theme_bw()# Plot Mortality trajectories for all individualsseedling_mort_plot <-ggplot(seedling_data_all, aes(x = time_min, y = mortality_prob, group = individual)) +geom_line(alpha =0.2, color = cbp2[4]) +stat_summary(aes(group =1), fun = mean, geom ="line", colour ="black", linewidth =1) +labs(x ="Time (min)",y ="Mortality probability") +theme_bw()seedling_plots <-plot_grid(temp_rain_plot, seedling_mort_plot, ncol =1, nrow =2, rel_heights =c(1, 1), align ="v", labels =c("a", "b"))seedling_plots
16) Plot heat failure probability across life stages
Plot heat failure probability (or mortality) across all life stages
Code
mort_plots <-ggplot() +geom_line(data = seedling_data_all, aes(x = time_min, y = mortality_prob, group = individual), alpha =0.2, color = cbp2[2]) +geom_line(data = seed_data_all, aes(x = time_min, y = mortality_prob, group = individual), alpha =0.2, color = cbp2[8]) +geom_line(data = pollen_data_all, aes(x = time_min, y = mortality_prob, group = individual), alpha =0.2, color = cbp2[4]) +stat_summary(data = seedling_data_all, aes(x = time_min, y = mortality_prob, group =1), fun = mean, geom ="line", colour = cbp2[2], linewidth =1) +stat_summary(data = seed_data_all, aes(x = time_min, y = mortality_prob, group =1), fun = mean, geom ="line", colour = cbp2[8], linewidth =1) +stat_summary(data = pollen_data_all, aes(x = time_min, y = mortality_prob, group =1), fun = mean, geom ="line", colour = cbp2[4], linewidth =1) +labs(x ="Time (min)",y ="Heat failure probability") +scale_y_continuous(limits =c(0,1)) +theme_bw()# save as 9 x 6.5 in plot_grid(temp_rain_plot, mort_plots, nrow =2, labels ="auto", align ="hv")
17) Plot TDT curves and damage-repair plots
Plot the TDT curves of each individual and population of each life stage on a single plot and combine with the damage-repair plot simulation across temperature.
Code
# Variance parameters alpha_sd <-0.02# default used = 0.02beta_sd <-0.01# default used = 0.01# Define parameters for each life stagelife_stages <-list(pollen =list(name ="Pollen",ctmax =mean(pollen_bt_ct, na.rm =TRUE),z =mean(pollen_bt_z, na.rm =TRUE),alpha_sd = alpha_sd,beta_sd = beta_sd,n_individuals =100,color = cbp2[8] ),seed =list(name ="Seed",ctmax =mean(seed_bt_ct, na.rm =TRUE),z =mean(seed_bt_z, na.rm =TRUE),alpha_sd = alpha_sd,beta_sd = beta_sd,n_individuals =100,color = cbp2[2] ),seedling =list(name ="Seedling",ctmax =mean(seedling_bt_ct, na.rm =TRUE),z =mean(seedling_bt_z, na.rm =TRUE),alpha_sd = alpha_sd,beta_sd = beta_sd,n_individuals =100,color = cbp2[4] ))# Shared parameters (same for all life stages)shared_params <-list(TA =298.15,TAL =283.15,TAH =318.15,TL =273.15,TH =323.15,TREF =298.15,kdot0 =0,kdot1 =0.05,kdot_ref =0.05)# Generate individual parameters for each life stageall_individuals <-bind_rows(lapply(life_stages, generate_individuals))cat("Generated parameters for:\n")
Generated parameters for:
Code
for(stage innames(life_stages)) { n <- life_stages[[stage]]$n_individuals stage_data <- all_individuals %>%filter(life_stage == life_stages[[stage]]$name)cat(sprintf(" %s: %d individuals\n", life_stages[[stage]]$name, n))cat(sprintf(" Original: CTmax = %.2f, z = %.2f\n", life_stages[[stage]]$ctmax, life_stages[[stage]]$z))cat(sprintf(" Mean simulated: CTmax = %.2f (SD=%.2f), z = %.2f (SD=%.2f)\n", mean(stage_data$ctmax), sd(stage_data$ctmax),mean(stage_data$z), sd(stage_data$z)))cat(sprintf(" Alpha: mean = %.3f (SD=%.3f), Beta: mean = %.3f (SD=%.3f)\n\n",mean(stage_data$alpha), sd(stage_data$alpha),mean(stage_data$beta), sd(stage_data$beta)))}
Pollen: 100 individuals
Original: CTmax = 44.01, z = -3.10
Mean simulated: CTmax = 43.91 (SD=1.34), z = -3.10 (SD=0.09)
Alpha: mean = 14.184 (SD=0.018), Beta: mean = -0.323 (SD=0.010)
Seed: 100 individuals
Original: CTmax = 58.98, z = -5.45
Mean simulated: CTmax = 58.82 (SD=3.21), z = -5.43 (SD=0.30)
Alpha: mean = 10.825 (SD=0.018), Beta: mean = -0.185 (SD=0.010)
Seedling: 100 individuals
Original: CTmax = 48.82, z = -4.41
Mean simulated: CTmax = 48.69 (SD=2.14), z = -4.40 (SD=0.19)
Alpha: mean = 11.060 (SD=0.018), Beta: mean = -0.228 (SD=0.010)
