adRes-package: Assessing whether phenotypic changes in response to climate...

In radchukv/adRes: Assessing Whether Phenotypic Changes in Response to Climate are Adaptive

Description

This R package aims at providing the data and documenting the analysis behind the results from the paper entitled 'Adaptive responses of animals to climate change are most likely insufficient' by Radchuk et al. Nature Communications (2019).

Details

This package has not been conceived for general use!

All main functions of this package contain a short documentation and examples which could be useful for those who try to understand our code. Type ls("package:adRes") for a list of all exported functions.

We recommend you to follow the examples below to reproduce the results of our analysis and understand the structure of our workflow.

You may also directly explore the files contained in the package after uncompressing the content of the *.tar.gz file (link available on the GitHub page https://github.com/radchukv/adRes). You can use the R function untar to extract the content of the tarball.

The package also contains several fitted models directly, to spare you the time required to fit them.

In the examples below, we provide the workflow leading the results presented in the paper.

References

Radchuk et al. (2019) Adaptive responses of animals to climate change are most likely insufficient". Nature Communications.

Examples

 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####              Settings for the workflow              ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 nb_cores <- 2L ## increase the number for using more cores
 digit <- 3L
 par_ini <- par(no.readonly = TRUE) ## default parameters for display



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####              Temperature: Condition 1              ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ## for PRCS
 mod_T_prcs <- fit_all(data = dat_Clim, temperature = TRUE,
                       precipitation = FALSE, phenology = TRUE,
                       morphology = TRUE, condition = '1',
                       nb_cores = nb_cores, rand_trait = FALSE,
                       fixed = NULL, digit = digit)

 ## for PRC
 ## Not run: 
 mod_T_prc <- fit_all(data = dat_Clim_prc, temperature = TRUE,
                      precipitation = FALSE, phenology = TRUE,
                      morphology = TRUE, condition = '1',
                      nb_cores = nb_cores, rand_trait = FALSE,
                      fixed = NULL, digit = digit)

 ## forest plot as in Fig. 2
 par(par_ini)
 plot_T_cond1 <- plot_forest(meta_obj1 = mod_T_prcs$meta_res,
                             meta_obj2 = NULL,
                             list_extra_meta_obj =
                             list(mod_T_prcs$meta_res, mod_T_prc$meta_res),
                             sort = c("Coord"),
                             increasing = FALSE,
                             labels = c(traits = FALSE,
                                        fitness = FALSE,
                                        country = TRUE,
                                        authors = FALSE))
 mtext('Across', side = 2, line = 3,
       at = -1, las = 2, cex = 0.9, col = 'black')
 mtext('studies', side = 2, line = 3,
       at = -2, las = 2, cex = 0.9, col = 'black')
 
## End(Not run)

 ## plots of raw data and extracted slopes per study, as in Suppl. Fig. S11
 par(par_ini); par(oma = c(3, 5, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Clim, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE, morphology = TRUE,
          condition = '1', id_to_do = c(1:20))

 ## in the full dataset there are 41 ids (21_41). In the publicly shared 36
 par(par_ini); par(oma = c(3, 5, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Clim, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE,
          morphology = TRUE, condition = '1',
          id_to_do = c(21:36))



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####        Temperature: Condition 2, phenology         ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ## for PRCS
 mod_phen_T_prcs <- fit_all(data = dat_Trait, temperature = TRUE,
                            precipitation = FALSE, phenology = TRUE,
                            morphology = FALSE, condition = '2',
                            nb_cores = nb_cores, rand_trait = FALSE,
                            fixed = NULL, digit = digit)

 ## fixed effects
 mod_phen_T_prcs_Type <- fit_all(data = dat_Trait,
                                 temperature = TRUE, precipitation = FALSE,
                                 phenology = TRUE, morphology = FALSE,
                                 condition = '2', nb_cores = nb_cores,
                                 rand_trait = FALSE, fixed = 'Trait_Cat',
                                 digit = digit)

 ## for PRC
 ## Not run: 
 mod_phen_T_prc <- fit_all(data = dat_Trait_prc,
                           temperature = TRUE, precipitation = FALSE,
                           phenology = TRUE, morphology = FALSE,
                           condition = '2', nb_cores = nb_cores,
                           rand_trait = FALSE, fixed = NULL,
                           digit = digit)

 ## fixed effects
 mod_phen_T_prc_Type <- fit_all(data = dat_Trait_prc,
                                temperature = TRUE, precipitation = FALSE,
                                phenology = TRUE, morphology = FALSE,
                                condition = '2', nb_cores = nb_cores,
                                rand_trait = FALSE, fixed = 'Trait_Cat',
                                digit = digit)

 mod_phen_T_prc_Taxon <- fit_all(data = dat_Trait_prc,
                                 temperature = TRUE, precipitation = FALSE,
                                 phenology = TRUE, morphology = FALSE,
                                 condition = '2', nb_cores = nb_cores,
                                 rand_trait = FALSE, fixed = 'Taxon',
                                 digit = digit)
 
## End(Not run)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####        Temperature: Condition 2, morphology        ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ## for PRCS
 mod_morph_T_prcs <- fit_all(data = dat_Trait, temperature = TRUE,
                             precipitation = FALSE, phenology = FALSE,
                             morphology = TRUE, condition = '2',
                             nb_cores = nb_cores, rand_trait = FALSE,
                             fixed = NULL, digit = digit)

 ## fixed effects
 mod_morph_T_prcs_Type <- fit_all(data = dat_Trait,
                                  temperature = TRUE, precipitation = FALSE,
                                  phenology = FALSE, morphology = TRUE,
                                  condition = '2', nb_cores = nb_cores,
                                  rand_trait = FALSE, fixed = 'Morph_type',
                                  digit = digit)

 ## for PRC
 ## Not run: 
   mod_morph_T_prc <- fit_all(data = dat_Trait_prc,
                              temperature = TRUE, precipitation = FALSE,
                              phenology = FALSE, morphology = TRUE,
                              condition = '2', nb_cores = nb_cores,
                              rand_trait = FALSE, fixed = NULL,
                              digit = digit)

   ## fixed effects
   mod_morph_T_prc_Type <- fit_all(data = dat_Trait_prc,
                                   temperature = TRUE, precipitation = FALSE,
                                   phenology = FALSE, morphology = TRUE,
                                   condition = '2', nb_cores = nb_cores,
                                   rand_trait = FALSE, fixed = 'Morph_type',
                                   digit = digit)

   mod_morph_T_prc_Taxon <- fit_all(data = dat_Trait_prc,
                                    temperature = TRUE, precipitation = FALSE,
                                    phenology = FALSE, morphology = TRUE,
                                    condition = '2', nb_cores = nb_cores,
                                    rand_trait = FALSE, fixed = 'Taxon',
                                    digit = digit)

   mod_morph_T_prc_Therm <- fit_all(data = dat_Trait_prc,
                                    temperature = TRUE, precipitation = FALSE,
                                    phenology = FALSE, morphology = TRUE,
                                    condition = '2', nb_cores = nb_cores,
                                    rand_trait = FALSE, fixed = 'Blood',
                                    digit = digit)


 ## forest plot as in Fig. 3
 par(par_ini)
 plot_T_cond2 <- plot_forest(meta_obj1 = mod_phen_T_prcs$meta_res,
                             meta_obj2 = mod_morph_T_prcs$meta_res,
                             list_extra_meta_obj =
                             list(mod_phen_T_prcs$meta_res,
                                  mod_morph_T_prcs$meta_res,
                                  mod_phen_T_prc$meta_res,
                                  mod_phen_T_prc_Taxon$meta_res,
                                  mod_morph_T_prc$meta_res,
                                  mod_morph_T_prc_Taxon$meta_res),
                             sort = c('Species', 'Trait_Categ_det',
                                      'PaperID'),
                             increasing = TRUE,
                             labels = c(traits = TRUE,
                                        fitness = FALSE,
                                        country = FALSE,
                                        authors = FALSE))
 mtext('Across studies, ', at = -0.6, side = 2,
       line = 4, las = 2, cex = 0.9, col = 'black')
 mtext('PRCS dataset', at = -2.2, side = 2,
        line = 4.3, las = 2, cex = 0.9, col = 'black')
 mtext('Across studies, ', at = -7.4, side = 2,
       line = 4, las = 2, cex = 0.9, col = 'black')
 mtext('PRC dataset ', at = -9, side = 2,
       line = 4.3, las = 2, cex = 0.9, col = 'black')
 
## End(Not run)

 ## plots of raw data and extracted slopes per study, as in Suppl. Fig. S12-S13
 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Trait, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE,
          morphology = FALSE, condition = '2',
          id_to_do = c(1:20))

 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Trait, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE,
          morphology = FALSE, condition = '2',
          id_to_do = c(21:38))

 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Trait, temperature = TRUE,
          precipitation = FALSE, phenology = FALSE,
          morphology = TRUE, condition = '2',
          id_to_do = c(1:4))



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####               Temperature: Condition 3             ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ## phenology
 mod_Sel_T_phen <- fit_all(data = dat_Sel, temperature = TRUE,
                           precipitation = FALSE, phenology = TRUE,
                           morphology = FALSE, condition = '3',
                           nb_cores = nb_cores, rand_trait = FALSE,
                           fixed = NULL, digit = digit)

 ## fixed effects
 mod_Sel_T_phen_Fitn <- fit_all(data = dat_Sel, temperature = TRUE,
                                precipitation = FALSE, phenology = TRUE,
                                morphology = FALSE, condition = '3',
                                nb_cores = nb_cores, rand_trait = FALSE,
                                fixed = 'Fitness_Categ', digit = digit)

 mod_Sel_T_phen_Gener <- fit_all(data = dat_Sel, temperature = TRUE,
                                 precipitation = FALSE, phenology = TRUE,
                                 morphology = FALSE, condition = '3',
                                 nb_cores = nb_cores, rand_trait = FALSE,
                                 fixed = 'GenerationLength_yr',
                                 digit = digit)
 ## morphology
 mod_Sel_T_morph <- fit_all(data = dat_Sel, temperature = TRUE,
                            precipitation = FALSE, phenology = FALSE,
                            morphology = TRUE, condition = '3',
                            nb_cores = nb_cores, rand_trait = FALSE,
                            fixed = NULL, digit = digit)
 ## fixed effects
 mod_Sel_T_morph_Fitn <- fit_all(data = dat_Sel, temperature = TRUE,
                                 precipitation = FALSE, phenology = FALSE,
                                 morphology = TRUE, condition = '3',
                                 nb_cores = nb_cores, rand_trait = FALSE,
                                 fixed = 'Fitness_Categ', digit = digit)

 mod_Sel_T_morph_Gener <- fit_all(data = dat_Sel, temperature = TRUE,
                                  precipitation = FALSE, phenology = FALSE,
                                  morphology = TRUE, condition = '3',
                                  nb_cores = nb_cores, rand_trait = FALSE,
                                  fixed = 'GenerationLength_yr',
                                  digit = digit)
 ## forest plot as in Fig. 4
 par(par_ini)
 plot_Sel_T_cond3 <- plot_forest(meta_obj1 = mod_Sel_T_phen$meta_res,
                         meta_obj2 = mod_Sel_T_morph$meta_res,
                         list_extra_meta_obj =
                         list(mod_Sel_T_phen$meta_res,
                              mod_Sel_T_morph$meta_res,
                              mod_Sel_T_phen_Fitn$meta_res),
                         sort = c('Species', 'Fitness_Categ', 'PaperID'),
                         increasing = TRUE,
                         labels = c(traits = TRUE,
                                    fitness = TRUE,
                                    country = FALSE,
                                    authors = FALSE))
 mtext('Across', side = 2, line = 4.5, at = -2.2,
       las = 2, cex = 0.85, col = 'black')
 mtext('studies', side = 2, line = 4.5, at = -3.4,
       las = 2, cex = 0.85, col = 'black')

 ## plot of raw data for supplementary, as in S14-S15
 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Sel, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE, morphology = FALSE,
          condition = '3', id_to_do = c(1:20))

 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Sel, temperature = TRUE,
          precipitation = FALSE, phenology = TRUE, morphology = FALSE,
          condition = '3', id_to_do = c(21:36))

 par(par_ini); par(oma = c(2, 2, 0, 0), mar = c(1, 2, 3, 1))
 plot_raw(data = dat_Sel, temperature = TRUE,
          precipitation = FALSE, phenology = FALSE, morphology = TRUE,
          condition = '3', id_to_do = c(1:9))

 ## change of selection over years, phenology
 mod_Sel_T_phen_year <- fit_all(data = dat_Sel, temperature = TRUE,
                                precipitation = FALSE, phenology = TRUE,
                                morphology = FALSE, condition = '3b',
                                nb_cores = nb_cores, rand_trait = FALSE,
                                fixed = NULL, digit = digit)
 ## change of selection over years, morphology
 mod_Sel_T_morph_year <- fit_all(data = dat_Sel, temperature = TRUE,
                                precipitation = FALSE, phenology = FALSE,
                                morphology = TRUE, condition = '3b',
                                nb_cores = nb_cores, rand_trait = FALSE,
                                fixed = NULL, digit = digit)

 par(par_ini)
 plot_Sel_T_overYrs <- plot_forest(meta_obj1 = mod_Sel_T_phen_year$meta_res,
                       meta_obj2 = mod_Sel_T_morph_year$meta_res,
                       list_extra_meta_obj =
                       list(mod_Sel_T_phen_year$meta_res,
                            mod_Sel_T_morph_year$meta_res),
                       sort = c('Species', 'Study_Authors',
                                 'Fitness_Categ'),
                       increasing = TRUE,
                       labels = c(traits = TRUE,
                                  fitness = TRUE,
                                  country = FALSE,
                                  authors = FALSE))
 mtext('Across', side = 2, line = 6, at = -1,
       las = 2, cex = 0.8, col = 'black')
 mtext('studies', side = 2, line = 6, at = -2.2,
       las = 2, cex = 0.8, col = 'black')



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
####  Assess whether responses are adaptive: products of slopes             #
####   from conditions 1 & 2 vs. the slopes from condition 3                #
####                       binomial test                                    #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 # assess condition 3: whether response is adaptive?
 ## extracting slopes
 all_slopes_phen <- merge_3slopes(meta_obj = mod_T_prcs,
                                meta_obj1 = mod_phen_T_prcs,
                                meta_obj2 = mod_Sel_T_phen)

all_slopes_morph <- merge_3slopes(meta_obj = mod_T_prcs,
                                 meta_obj1 = mod_morph_T_prcs,
                                 meta_obj2 = mod_Sel_T_morph)
 ## testing adaptation
 res_phen <- table(all_slopes_phen$slope.x * all_slopes_phen$slope.y > 0)

 (bc_phenT <- binom.confint(res_phen["TRUE"], sum(res_phen)))
 (bc_phenF <- binom.confint(res_phen["FALSE"], sum(res_phen)))
 binom.test(res_phen["TRUE"], sum(res_phen))

 res_morph <- table(all_slopes_morph$slope.x * all_slopes_morph$slope.y > 0)
 (bc_morphT <- binom.confint(res_morph["TRUE"], sum(res_morph)))
 (bc_morphF <- binom.confint(res_morph["FALSE"], sum(res_morph)))
 binom.test(res_morph["TRUE"], sum(res_morph))


 par(par_ini)
 mat <- matrix(c(1, 2, 3, 3), 2,2, byrow = TRUE)
 layout(mat, widths =  rep.int(1, ncol(mat)),
        heights = rep.int(1, nrow(mat)), respect = TRUE)
 plot_slopes_products(data = merge_3slopes(meta_obj = mod_T_prcs,
                      meta_obj1 = mod_phen_T_prcs,
                      meta_obj2 = mod_Sel_T_phen),
                      phenological = TRUE)

plot_slopes_products(data = merge_3slopes(meta_obj = mod_T_prcs,
                     meta_obj1 = mod_morph_T_prcs,
                     meta_obj2 = mod_Sel_T_morph),
                     phenological = FALSE, xlim = c(-0.1, 0.1),
                     ylim = c(-0.1, 0.1))

 par(par_ini); par(mar = c(4, 8, 3, 6))
 theplot_Adapt <- barplot(cbind(c(res_phen["TRUE"]/sum(res_phen),
                                  res_phen["FALSE"]/sum(res_phen)),
                                    c(res_morph["TRUE"]/sum(res_morph),
                                      res_morph["FALSE"]/sum(res_morph))),
                          names.arg = c("Phenology", "Morphology"),
                          las = 1, ylim = c(0, 1.1),
                          ylab = "Proportion of studies",
                          col = c("grey", "black"),
                          cex.lab = 1.3, cex.axis = 1.3, cex.names = 1.3)
 legend(x = 0.8, y = 1, horiz = TRUE, fill = c("grey", "black"),
        legend = c("adaptative", "maladaptive"), bg = 'white')
 arrows(x0 = theplot_Adapt[1], x1 = theplot_Adapt[1],
        y0 = bc_phenT[11, "lower"], y1 = bc_phenT[11, "upper"],
        code = 3, angle = 90, col = "grey30", lwd = 4, length = 0.1)
 arrows(x0 = theplot_Adapt[2], x1 = theplot_Adapt[2],
        y0 = bc_morphT[11, "lower"], y1 = bc_morphT[11, "upper"],
        code = 3, angle = 90, col = "grey30", lwd = 4, length = 0.1)
 mtext('C)', side = 3, line = 0, adj = 0, cex = 2)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####         Assess whether responses are adaptive         ####
 ####                     meta-analysis                     ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 dat_AllTPhen <- prepare_data(data = dat_All, temperature = TRUE,
                               phenology = TRUE, morphology = FALSE)
 dat_Adpt <- compute_adaptation(data = dat_AllTPhen)

 dat_Adpt$slope <- dat_Adpt$WMS * sign(dat_Adpt$prod12)
 dat_Adpt$SE_slope <- dat_Adpt$se_WMS

 mod_Adapt_T_Phen <- fit_meta(meta_data = dat_Adpt, fixed = NULL)



 ## morphology
 dat_AllTMorph <- prepare_data(data = dat_All, temperature = TRUE,
                               phenology = FALSE, morphology = TRUE)
 dat_AllTMorph <- dat_AllTMorph[! is.na(dat_AllTMorph$ScaledSE), ]
 dat_Adpt_Morph <- compute_adaptation(data = dat_AllTMorph)

 dat_Adpt_Morph$slope <- dat_Adpt_Morph$WMS * sign(dat_Adpt_Morph$prod12)
 dat_Adpt_Morph$SE_slope <- dat_Adpt_Morph$se_WMS

 mod_Adapt_T_Morph <- fit_meta(meta_data = dat_Adpt_Morph, fixed = NULL)

 # using selection dataset only: to identify whether the sign
 # is in the same direction as phenotypic change
 dat_SelT <- prepare_data(data = dat_Sel, temperature = TRUE,
                               phenology = TRUE, morphology = FALSE)
 dat_AllTPhen_Sel <- dat_All[dat_All$id %in% unique(dat_SelT$id), ]
 dat_Adpt_Sel <- compute_adaptation(data = dat_AllTPhen_Sel)

 dat_Adpt_Sel$slope <- dat_Adpt_Sel$WMS * sign(dat_Adpt_Sel$prod12)
 dat_Adpt_Sel$SE_slope <- dat_Adpt_Sel$se_WMS

 mod_Adapt_T_Sel <- fit_meta(meta_data = dat_Adpt_Sel, fixed = NULL)

 mod_Adapt_T_Sel_Sens <- fit_meta(meta_data =
                                  dat_Adpt_Sel[dat_Adpt_Sel$slope < 0.6, ],
                                  fixed = NULL)


 ## morphology
 dat_SelT_Morph <- prepare_data(data = dat_Sel, temperature = TRUE,
                               phenology = FALSE, morphology = TRUE)
 dat_AllMorph_Sel <- dat_All[dat_All$id %in% unique(dat_SelT_Morph$id), ]
 dat_AllMorph_Sel <- dat_AllMorph_Sel[! is.na(dat_AllMorph_Sel$ScaledSE), ]
 dat_Adpt_Morph_Sel <- compute_adaptation(data = dat_AllMorph_Sel)

 dat_Adpt_Morph_Sel$slope <- dat_Adpt_Morph_Sel$WMS *
                             sign(dat_Adpt_Morph_Sel$prod12)
 dat_Adpt_Morph_Sel$SE_slope <- dat_Adpt_Morph_Sel$se_WMS

 mod_Adapt_T_Morph_Sel <- fit_meta(meta_data = dat_Adpt_Morph_Sel, fixed = NULL)

 par(par_ini)
 plot_adapt_T_Sel <- plot_forest(meta_obj1 = mod_Adapt_T_Sel,
                                 meta_obj2 = mod_Adapt_T_Morph_Sel,
                                 list_extra_meta_obj =
                                   list(mod_Adapt_T_Sel, mod_Adapt_T_Morph_Sel),
                                 sort = c('Species', 'Study_Authors',
                                          'Fitness_Categ'),
                                 increasing = TRUE,
                                 labels = c(traits = TRUE,
                                            fitness = TRUE,
                                            country = FALSE,
                                            authors = FALSE),
                                 mar = c(7, 10, 2, 2))

 mtext('the direction of climate-driven', side = 1, line = 4.2, cex = 1.3)
 mtext('trait change', side = 1, line = 5.3, cex = 1.3)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####                   Precipitation: Condition 1              ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 ## for PRCS
 mod_Prec_prcs <- fit_all(data = dat_Clim, temperature = FALSE,
                          precipitation = TRUE, phenology = TRUE,
                          morphology = TRUE, condition = '1',
                          nb_cores = nb_cores, rand_trait = FALSE,
                          fixed = NULL, digit = digit)
 ## for PRC
 mod_Prec_prc <- fit_all(data = dat_Clim_prc, temperature = FALSE,
                         precipitation = TRUE, phenology = TRUE,
                         morphology = TRUE, condition = '1',
                         nb_cores = nb_cores, rand_trait = FALSE,
                         fixed = NULL, digit = digit)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####                Precipitation: Condition 2                 ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 ## for PRCS
 mod_phen_Prec_prcs <- fit_all(data = dat_Trait, temperature = FALSE,
                               precipitation = TRUE, phenology = TRUE,
                               morphology = FALSE, condition = '2',
                               nb_cores = nb_cores, rand_trait = FALSE,
                               fixed = NULL, digit = digit)

 ## for PRC
 mod_phen_Prec_prc <- fit_all(data = dat_Trait_prc, temperature = FALSE,
                              precipitation = TRUE, phenology = TRUE,
                              morphology = FALSE, condition = '2',
                              nb_cores = nb_cores, rand_trait = FALSE,
                              fixed = NULL, digit = digit)

 ## fixed effects
 mod_phen_Prec_prc_Taxon <- fit_all(data = dat_Trait_prc,
                                    temperature = FALSE,
                                    precipitation = TRUE,
                                    phenology = TRUE, morphology = FALSE,
                                    condition = '2', nb_cores = nb_cores,
                                    rand_trait = FALSE, fixed = 'Taxon',
                                    digit = digit)




 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####                  Precipitation: Condition 3               ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ## for phenology, because there is no data on morphological traits
 ## in the dataset on precipitation
 mod_Sel_Prec_phen <- fit_all(data = dat_Sel,
                              temperature = FALSE, precipitation = TRUE,
                              phenology = TRUE, morphology = FALSE,
                              condition = '3', nb_cores = nb_cores,
                              rand_trait = FALSE, fixed = NULL,
                              digit = digit)

 ## fixed effects
 mod_Sel_Prec_phen_Fitn <- fit_all(data = dat_Sel, temperature = FALSE,
                                   precipitation = TRUE, phenology = TRUE,
                                   morphology = FALSE, condition = '3',
                                   nb_cores = nb_cores, rand_trait = FALSE,
                                   fixed = 'Fitness_Categ', digit = digit)

 mod_Sel_Prec_phen_Gener <- fit_all(data = dat_Sel, temperature = FALSE,
                                    precipitation = TRUE, phenology = TRUE,
                                    morphology = FALSE, condition = '3',
                                    nb_cores = nb_cores, rand_trait = FALSE,
                                    fixed = 'GenerationLength_yr',
                                    digit = digit)

 ## supplementary plot of all effects for each condition (as Fig. S3)
 par(par_ini); par(oma = c(2,1,1,0), mfrow = c(1,3))

 ## panel a)
 plot_Prec_Cond1 <- plot_forest(meta_obj1 = mod_Prec_prcs$meta_res,
                                meta_obj2 = NULL,
                                list_extra_meta_obj =
                                list(mod_Prec_prcs$meta_res,
                                     mod_Prec_prc$meta_res),
                                sort = c("Coord"),
                                increasing = FALSE,
                                labels = c(traits = FALSE,
                                           fitness = FALSE,
                                           country = TRUE,
                                           authors = FALSE),
                                mar = c(4, 9, 2, 2))
 mtext('Across', side = 2, line = 5,
       at = -1.2, las = 2, cex = 0.9, col = 'black')
 mtext('studies', side = 2, line = 5,
       at = -1.8, las = 2, cex = 0.9, col = 'black')
 mtext('a)', side = 2, at = 11.1, cex = 1.5, las = 2, line = 1)

 ## panel b)
 plot_Prec_Trait_Cond2 <- plot_forest(meta_obj1 = mod_phen_Prec_prcs$meta_res,
                         meta_obj2 = NULL,
                         list_extra_meta_obj =
                         list(mod_phen_Prec_prcs$meta_res,
                              mod_phen_Prec_prc$meta_res),
                         sort = c('Species', 'Trait_Categ_det'),
                         increasing = FALSE,
                         labels = c(traits = TRUE,
                                    fitness = FALSE,
                                    country = FALSE,
                                    authors = FALSE),
                         mar = c(4, 9, 2, 2))
 mtext('Across', side = 2, line = 6,
       at = -1.2, las = 2, cex = 0.8, col = 'black')
 mtext('studies', side = 2, line = 6,
       at = -1.8, las = 2, cex = 0.8, col = 'black')
 mtext('b)', side = 2, at = 12.3, cex = 1.5, las = 2, line = 1.2)

 ## panel c)
 plot_Sel_Prec_cond3 <- plot_forest(meta_obj1 = mod_Sel_Prec_phen$meta_res,
                         meta_obj2 = NULL,
                         list_extra_meta_obj =
                         list(mod_Sel_Prec_phen$meta_res,
                              mod_Sel_Prec_phen_Fitn$meta_res),
                         sort = c('Species', 'Study_Authors',
                                  'Fitness_Categ'),
                         increasing = FALSE,
                         labels = c(traits = TRUE,
                                    fitness = TRUE,
                                    country = FALSE,
                                    authors = FALSE),
                         mar = c(4, 13, 2, 2))
 mtext('Across', side = 2, line = 8,
       at = -2, las = 2, cex = 0.8, col = 'black')
 mtext('studies', side = 2, line = 8,
       at = -3, las = 2, cex = 0.8, col = 'black')
 mtext('c)', side = 2, at = 13.3, cex = 1.5, las = 2, line = 1.2)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####                      Sensitivity analysis                           ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 ## excluding the study by Goodenough et al. (2011)
 mod_phen_T_prcs_noOutlier <- fit_all(data = dat_Trait[
                              dat_Trait$Study_Authors != 'Goodenough', ],
                              temperature = TRUE, precipitation = FALSE,
                              phenology = TRUE, morphology = FALSE,
                              condition = '2', nb_cores = nb_cores,
                              rand_trait = FALSE, fixed = NULL,
                              digit = digit)

 ## selection
 mod_Sel_T_phen_noOutlier <- fit_all(data = dat_Sel[
                             dat_Sel$Study_Authors != 'Goodenough', ],
                             temperature = TRUE, precipitation = FALSE,
                             phenology = TRUE, morphology = FALSE,
                             condition = '3', nb_cores = nb_cores,
                             rand_trait = FALSE, fixed = NULL,
                             digit = digit)

 ## fixed effects
 mod_Sel_T_phen_Fitn_noOutlier <- fit_all(data = dat_Sel[
                                  dat_Sel$Study_Authors != 'Goodenough', ],
                                  temperature = TRUE, precipitation = FALSE,
                                  phenology = TRUE, morphology = FALSE,
                                  condition = '3', nb_cores = nb_cores,
                                  rand_trait = FALSE, fixed = 'Fitness_Categ',
                                  digit = digit)

 mod_Sel_T_phen_Gener_noOutlier <- fit_all(data = dat_Sel[
                                   dat_Sel$Study_Authors != 'Goodenough', ],
                                   temperature = TRUE, precipitation = FALSE,
                                   phenology = TRUE, morphology = FALSE,
                                   condition = '3', nb_cores = nb_cores,
                                   rand_trait = FALSE,
                                   fixed = 'GenerationLength_yr',
                                   digit = digit)

 ## excluding the study on mammal
 mod_phen_T_prcs_noMammal <- fit_all(data = dat_Trait[
                             dat_Trait$Study_Authors != 'Plard_et_al', ],
                             temperature = TRUE, precipitation = FALSE,
                             phenology = TRUE, morphology = FALSE,
                             condition = '2', nb_cores = nb_cores,
                             rand_trait = FALSE, fixed = NULL,
                             digit = digit)

 ## phenology
 mod_Sel_T_phen_noMammal <- fit_all(data = dat_Sel[
                            dat_Sel$Study_Authors != 'Plard_et_al', ],
                            temperature = TRUE, precipitation = FALSE,
                            phenology = TRUE, morphology = FALSE,
                            condition = '3', nb_cores = nb_cores,
                            rand_trait = FALSE, fixed = NULL,
                            digit = digit)

 ## fixed effects
 mod_Sel_T_phen_Fitn_noMammal <- fit_all(data = dat_Sel[
                                 dat_Sel$Study_Authors != 'Plard_et_al', ],
                                 temperature = TRUE, precipitation = FALSE,
                                 phenology = TRUE, morphology = FALSE,
                                 condition = '3', nb_cores = nb_cores,
                                 rand_trait = FALSE, fixed = 'Fitness_Categ',
                                 digit = digit)

 mod_Sel_T_phen_Gener_noMammal <- fit_all(data = dat_Sel[
                                  dat_Sel$Study_Authors != 'Plard_et_al', ],
                                  temperature = TRUE, precipitation = FALSE,
                                  phenology = TRUE, morphology = FALSE,
                                  condition = '3', nb_cores = nb_cores,
                                  rand_trait = FALSE,
                                  fixed = 'GenerationLength_yr',
                                  digit = digit)

 ## excluding both the study on mammal and the study by Goodenough et al. (2011)
 mod_phen_T_prcs_noBoth <- fit_all(data = dat_Trait[!
                             dat_Trait$Study_Authors %in%
                             c('Goodenough', 'Plard_et_al'), ],
                             temperature = TRUE, precipitation = FALSE,
                             phenology = TRUE, morphology = FALSE,
                             condition = '2', nb_cores = nb_cores,
                             rand_trait = FALSE, fixed = NULL, digit = digit)

 mod_Sel_T_phen_noBoth <- fit_all(data = dat_Sel[!
                                  dat_Sel$Study_Authors %in%
                                  c('Goodenough', 'Plard_et_al'), ],
                                  temperature = TRUE, precipitation = FALSE,
                                  phenology = TRUE, morphology = FALSE,
                                  condition = '3', nb_cores = nb_cores,
                                  rand_trait = FALSE, fixed = NULL,
                                  digit = digit)

 ## fixed effects
 mod_Sel_T_phen_Fitn_noBoth <- fit_all(data = dat_Sel[!
                                 dat_Sel$Study_Authors %in%
                                 c('Goodenough', 'Plard_et_al'), ],
                                 temperature = TRUE, precipitation = FALSE,
                                 phenology = TRUE, morphology = FALSE,
                                 condition = '3', nb_cores = nb_cores,
                                 rand_trait = FALSE, fixed = 'Fitness_Categ',
                                 digit = digit)

mod_Sel_T_phen_Gener_noBoth <- fit_all(data = dat_Sel[!
                                   dat_Sel$Study_Authors %in%
                                   c('Goodenough', 'Plard_et_al'), ],
                                   temperature = TRUE, precipitation = FALSE,
                                   phenology = TRUE, morphology = FALSE,
                                   condition = '3', nb_cores = nb_cores,
                                   rand_trait = FALSE,
                                   fixed = 'GenerationLength_yr',
                                   digit = digit)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####              Sensitivity to the inclusion of abundance             ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 dat_PhenT_ab <- prepare_data(data = dat_Trait, temperature = TRUE,
                              precipitation = FALSE, phenology = TRUE,
                              morphology = FALSE)
 dat_PhenT_ab <- dat_PhenT_ab[! is.na(dat_PhenT_ab$ScaledPop), ]

 # to exclude studies with < 11 years
 num_years <- dat_PhenT_ab %>%
              dplyr::group_by(., id) %>%
              dplyr::summarise(num = dplyr::n())

 num_years <- num_years[num_years$num > 10, ]
 dat_PhenT_ab <- dat_PhenT_ab[dat_PhenT_ab$id %in% num_years$id, ]

 EfSizes_phen_ab <- extract_effects_all_ids(data = dat_PhenT_ab, condition = '2b',
                                           nb_cores = nb_cores)
 EfSizes_phen_ab$slope <- EfSizes_phen_ab$slope_clim
 EfSizes_phen_ab$SE_slope <- EfSizes_phen_ab$se_slope_clim
 meta_phen_abund <- fit_meta(EfSizes_phen_ab, fixed = NULL,
                            rand_trait = FALSE, digit = 3)


 ## and a model without abundance on this subset of data
 mod_phen_T_prcs_subs <- fit_all(data = dat_PhenT_ab, temperature = TRUE,
                                 precipitation = FALSE, phenology = TRUE,
                                 morphology = FALSE, condition = '2',
                                 nb_cores = nb_cores, rand_trait = FALSE,
                                 fixed = NULL, digit = digit)

 ## plot as Supplementary Fig. S14 (for all slopes)
 par(par_ini)
 plot_abund_effects(meta_data = EfSizes_phen_ab)



 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####               Supplementary: Funnel plots                          ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 ## funnel plots for PRCS, as Supplementary Fig. S16
 par(par_ini); par(mfrow = c(2, 3))
 plot_funnel(meta_obj = mod_T_prcs)
 mtext('a)', line = 1, adj = 0, cex = 1.7)

 plot_funnel(meta_obj = mod_phen_T_prcs)
 mtext('b)', line = 1, adj = 0, cex = 1.7)

 plot_funnel(meta_obj = mod_morph_T_prcs)
 mtext('c)', line = 1, adj = 0, cex = 1.7)

 plot_funnel(meta_obj = mod_Sel_T_phen, model = 'lm')
 mtext('d)', line = 1, adj = 0, cex = 1.7)

 plot_funnel(meta_obj = mod_Sel_T_morph)
 mtext('e)', line = 1, adj = 0, cex = 1.7)


 ## funnel plots for PRC, as Supplementary Fig. S17
 ## run only if the models were fit
 ## Not run: 
   par(par_ini); par(mfrow = c(1,3))

   plot_funnel(meta_obj = mod_T_prc)
   mtext('a)', line = 1, adj = 0, cex = 1.7)

   plot_funnel(meta_obj = mod_phen_T_prc)
   mtext('b)', line = 1, adj = 0, cex = 1.7)

   plot_funnel(meta_obj = mod_morph_T_prc)
   mtext('c)', line = 1, adj = 0, cex = 1.7)

## End(Not run)

 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####           Supplementary: temperature over years VS                ####
 ####                    study duration VS first year                   ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 ## make sure the model is fitted before running this model
 ## Not run: 
   plot_eff_dur_firstY(data = dat_Clim_prc, meta_obj = mod_T_prc)
 
## End(Not run)

 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
 ####              Prepare Supplementary Tables                          ####
 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

 # a list of models
 # can only be run if all the models were fitted,
 # some of which are time-consuming to fit
 ## Not run: 
   mod_list <- list(mod_T_prcs, mod_phen_T_prcs, mod_phen_T_prcs_Type,
                  mod_morph_T_prcs, mod_morph_T_prcs_Type, mod_Sel_T_phen,
                  mod_Sel_T_phen_Fitn, mod_Sel_T_phen_Gener, mod_Sel_T_morph,
                  mod_Sel_T_morph_Fitn, mod_Sel_T_morph_Gener, mod_Prec_prcs,
                  mod_phen_Prec_prcs, mod_Sel_Prec_phen, mod_Sel_Prec_phen_Fitn,
                  mod_Sel_Prec_phen_Gener, mod_T_prc, mod_phen_T_prc,
                  mod_phen_T_prc_Taxon, mod_phen_T_prc_Type, mod_morph_T_prc,
                  mod_morph_T_prc_Taxon, mod_morph_T_prc_Type,
                  mod_morph_T_prc_Therm, mod_Prec_prc, mod_phen_Prec_prc,
                  mod_phen_Prec_prc_Taxon)

   ## A table with the effect sizes and their SE
   tab_efSizes_ST1(model_list = mod_list)

   ## A table with LRT statistics
   tab_LRT_ST2(model_list = mod_list)
 
## End(Not run)
 # a list of models
 mod_list_sens <- list(mod_phen_T_prcs_noOutlier,
                       mod_Sel_T_phen_noOutlier,
                       mod_Sel_T_phen_Fitn_noOutlier,
                       mod_Sel_T_phen_Gener_noOutlier,
                       mod_phen_T_prcs_noMammal,
                       mod_Sel_T_phen_noMammal,
                       mod_Sel_T_phen_Fitn_noMammal,
                       mod_Sel_T_phen_Gener_noMammal,
                       mod_phen_T_prcs_noBoth,
                       mod_Sel_T_phen_noBoth,
                       mod_Sel_T_phen_Fitn_noBoth,
                       mod_Sel_T_phen_Gener_noBoth)

 tab_efSizes_ST4(model_list = mod_list_sens)

 ## a table for the LRT statistics and the variation due to random effects
 ## for the sensitivity analysis
 tab_LRT_ST5(model_list = mod_list_sens)

 ## a table for heterogeneity
 ## make sure all models are fitted before running this
 ## Not run: 
  mod_heterogen <- list(mod_T_prcs, mod_phen_T_prcs, mod_morph_T_prcs,
                        mod_Sel_T_phen, mod_Sel_T_morph, mod_Prec_prcs,
                        mod_phen_Prec_prcs, mod_Sel_Prec_phen, mod_T_prc,
                        mod_phen_T_prc, mod_morph_T_prc, mod_Prec_prc,
                        mod_phen_Prec_prc, mod_phen_T_prcs_noOutlier,
                        mod_Sel_T_phen_noOutlier, mod_phen_T_prcs_noMammal,
                        mod_Sel_T_phen_noMammal, mod_phen_T_prcs_noBoth,
                        mod_Sel_T_phen_noBoth)
  tab_heterog(model_list = mod_heterogen)
 
## End(Not run)

radchukv/adRes documentation built on June 1, 2019, 7:05 p.m.