dfoliatR : An R package for detection and analysis of insect defoliation signals in tree rings

.


Introduction
Variation in the width and morphology of annual radial growth rings in trees permits dating and quantification of past forest insect defoliator outbreaks. Defoliation can be distinguished from climate-and other disturbance-related influences by comparing ring-width or other annually-resolved features in the wood of host species to that of non-host species or annually-resolved climate records. The effect of defoliation on radial growth of trees has been recognized since the 1860s, and used to reconstruct outbreak regimes since the 1950s (Blais, 1954;Alfaro et al., 1982;Lynch, 2012). It was not until the 1980s, however, that precise dendrochronological techniques were applied for inferring defoliation events and reconstructing defoliator outbreak regimes (Swetnam et al., 1985;Speer, 2010;Lynch, 2012). The first studies (Swetnam et al., 1985;Lynch, 1989, 1993) focused on developing historical outbreak reconstructions of western spruce budworm (WSBW; Choristoneura freemani (Lepidoptera: Tortricidae); previously known as C. occidentalis). The methodology has since been successfully applied to a wide range of defoliator species, most of which are conifer herbivores, and has evolved in sophistication for a wide range of ecosystem situations (Lynch, 2012).
The main dendrochronological tool for inferring, dating, and characterizing defoliator outbreaks from tree-ring records has been the software routine OUTBREAK (Swetnam et al., 1985;Holmes and Swetnam, 1986;Swetnam and Lynch, 1989). OUTBREAK computes indices (described later in detail) of suppressed growth by subtracting a detrended and standardized climate series (a "control" chronology) from individual host-tree detrended and standardized radial growth series after the host and non-host series have been brought to a common variance. The non-host chronology usually consists of a site chronology developed from non-host tree species growing on a climate-sensitive site, but a gridded climate data point series, like the North American Drought Atlas (Cook and Krusic, 2004) also suffices. If the host and non-host species respond similarly to climate (which can and should be tested), the derived series retains variability that the host and non-host series do not have in common, generally the insect signal and some unexplained variability (noise). The user defines a rule base specifying the magnitude and duration that a period of indexed growth suppression must meet or surpass for a period of suppressed growth to be inferred as a defoliation event at the tree level. Rule bases are derived from the user's knowledge of insect and host ecologies, and from consideration of the likelihood and relative importance of Type I and II errors.
Though powerful, OUTBREAK is outdated and increasingly difficult to use in modern computing environments. It was written in FORTRAN V with inherently severe restrictions, as RAM and disk space were limited at that time (256 kb and 10 MB, respectively) and FORTRAN conventions imposed very strict formatting, file naming, and output conventions. The program lacks a graphical interface or capabilities, forcing users to import generated text files into spreadsheets or other software to assess results and perform analyses. Furthermore, OUTBREAK can only handle one test at a time, creating barriers to batch operation and a large burden for researchers with datasets including multiple sites. We developed dfoliatR (Guiterman et al., 2020) as an R-and dplR-based library to overcome these issues.
dfoliatR adds to a growing suite of dendrochronology packages in the R computing environment (R Core Team, 2019). Stemming from the dplR library (Bunn, 2008) that enables R users to read and write an array of tree-ring data formats, standardize ring width series, build and evaluate chronologies, and perform quality control (to name a few), one can now also measure ring widths from scanned images of prepared samples (Lara et al., 2015;Shi et al., 2019), conduct and check crossdating (Bunn, 2010), analyze sub-annual anatomical features (Rathgeber et al., 2011;Campelo et al., 2016), and perform many analytical tests (Zang and Biondi, 2015;Jevšenak and Levanič, 2018). Tools for assessing stand dynamics and disturbance analyses are under rapid development, with new packages for assessing release events (TRADER: Altman et al., 2014), metrics of growth resilience (pointRes: van der Maaten-Theunissen et al., 2015), and fire history (burnr: Malevich et al., 2018). The key objective of dfoliatR is to provide tools to identify and analyze insect defoliation and outbreak events by building on the methods employed by OUTBREAK. It capitalizes on the robust software already available in R by using dplR data formats for incoming tree-ring series and providing output data formats embodied by the tidyverse (Wickham et al., 2019) that include efficient data manipulation (dplyr: Wickham et al., 2020) andgraphics (ggplot2: Wickham, 2016).
In this paper, we describe the statistical methods employed by dfoliatR, compare results to those produced by OUTBREAK, and present an example analysis including test data sets and script. Users need not have much experience in R to replicate the analyses and graphics as presented. The R code below is executable in an R session once the required libraries are installed and loaded. Support documentation in addition to this paper is provided within the package via standard help menus and on the package website (https://chguiterman. github.io/dfoliatR/), which includes up-to-date vignettes that describe various routines. Code to generate a preprint of this manuscript, including the R scripts and tabular and graphical output is available from https://github.com/chguiterman/dfoliatR_paper.

Overview of the software
The dfoliatR library requires two sets of tree-ring data to infer defoliation and outbreak events: • Standardized ring-width series for individual trees of the host species. • A standardized tree-ring chronology from a local non-host species, or a climate reconstruction.
Users can develop these data sets in the software of their choosing, such as dplR or ARSTAN (Cook and Holmes, 1996). It is important that the host-tree data include only one tree-ring series per tree. dplR (via the dplR::treeMean() function) and dpl versions of ARSTAN have options for averaging multiple sample series into a tree-level series. At the heart of dfoliatR lies two functions: defoliate_trees() and outbreak(). These identify defoliation events on individual trees ( Fig. 1) and then composite across multiple trees to infer stand or site level outbreak events (Fig. 2).

Identifying defoliation of trees
The defoliate_trees() function is the point of entry to the dfoliatR library. It performs two processes, removing climate-related growth signals from the host-tree series and identifying tree-level defoliation events. The climatic or non-defoliation signals in each host-tree series are characterized by a non-host chronology or climate reconstruction. dfoliatR removes the non-defoliation signal by subtracting the non-host series from each host-tree series, which generates a residual index. In OUTBREAK, this residual index was termed the "corrected index." We call it the "growth suppression index" (GSI). The GSI is calculated the same as in OUTBREAK for each host tree as where H and NH are the host-tree series and the non-host chronology, in year i, respectively. Only the common period between the host-tree series and the non-host chronology are used in Eq. (1). The host and nonhost chronologies are brought to common variance by scaling the nonhost chronology by its mean (NH) and multiplying by the ratio of host and non-host standard deviations ( σH σNH ), which approximates the variance of the host tree series.
Negative departures in the normalized GSI (NGSI, or GSI converted to z-scores) that surpass user-specified thresholds in duration and magnitude are defined as defoliation events. As in OUTBREAK, the lowest NGSI value in the particular sequence being assessed must reach the magnitude threshold. The default setting is − 1.28 (NGSI is in units of standard deviation), which was previously determined to be representative of WSBW effects (Swetnam and Lynch, 1989) and is commonly used for other species (see Lynch, 2012). The year with the lowest value is termed the "year of maximum departure" and becomes a central point in time for assessing other thresholds before being included as a defoliation event. If the year of maximum departure is higher than the threshold (i.e., NGSI lowest > − 1.28), the sequence being assessed is omitted from the event results. Event duration is assessed by examining sequences of negative NGSI (for which one or more values exceeded the magnitude threshold) before and after the year of maximum departure. Each defoliation event is allowed one single-year positive excursion on each side of the year of maximum departure. Duration is computed across the entire sequence that may include these two positive excursions. As in OUTBREAK, the user specifies a duration threshold (minimum number of years) for a departure sequence to be inferred as a defoliation event. The default threshold is eight years, as is commonly used in WSBW studies (Swetnam and Lynch, 1989). If the sequence is shorter than the duration threshold, the sequence is omitted from the event results (i.e., both thresholds must be met). Researchers can, and should, adjust the duration and magnitude parameters accordingly and critically evaluate the results, as insect species vary in the length of their outbreaks and the degree to which they can suppress tree growth. OUTBREAK provides two sets of default values, those for WSBW, which typically has lengthy outbreaks, and ones for Douglas-fir tussock moth (Orgyia pseudotsugata (Lepidoptera: Tortricidae)) of three years duration with − 1.28 departure threshold that may be suitable for more eruptive species.
Like OUTBREAK, users are provided an option to suspend the duration threshold at the recent end of the series in cases where an outbreak event is known to be ongoing. This should be used if the user has direct knowledge of defoliation at the site during the sampling campaign. The advantage of allowing potentially short, series-end events is that it allows a current event to be included in return-interval estimates, and can aid in identifying the start-year for the current defoliation event or outbreak.
Diverging from OUTBREAK, dfoliatR includes an option allowing users to extend defoliation events on individual trees by bridging between sequential events (Fig. 3). In cases where two defoliation events are separated by a single year, bridging will link them into a single event. This option was added to dfoliatR during the testing phase of development, when we realized that OUTBREAK deliberately omits sequential, or back-to-back events, even when both events surpass the magnitude and duration thresholds. Instead, OUTBREAK will select the one sequential event with the lowest negative departure year. In every case we assessed (described below) we felt that the OUTBREAK-omitted defoliation events should have been maintained and recorded. Due in large part to reconstructions using OUTBREAK (see papers cited by Lynch (2012)), we now know considerably more about forest defoliator outbreak regimes than we did in the 1980s when OUTBREAK was under development. We think that two or more prolonged events separated by a single year should in some situations be considered a single event. This is particularly relevant to WSBW and spruce budworm (C. fumiferana), for which multiple outbreak regime reconstructions, as well as other research and forest health observations, show that outbreaks can be very long (Schmitt et al., 1984;Sanders et al., 1985;Brookes et al., 1987, and many later publications). Often the greatest growth suppression occurs late in the outbreak due to lag effects between defoliation and radial growth, and to cumulative effects accrued on a tree's resources (Brubaker, 1978;Alfaro et al., 1982;Wickman, 1986;Lynch, 1989, 1993;Mason et al., 1997;Axelson et al., 2014). We urge caution in using the bridging option, however, because it may not be appropriate for all studied insects, such as in situations where impacted stands barely recover from one outbreak before another begins, as with pine processionary caterpillars (Thaumetopoea pityocampa (Lepidoptera: Thaumetopoeidae)) (Carus, 2004(Carus, , 2009 or where outbreaks are known to be very short, such as larch budmoth (Zeiraphera diniana Gn.) in the European Alps (Esper et al., 2007).

Inferring outbreak events
Defoliation of one or a few trees does not constitute an outbreak. To determine when defoliation becomes an outbreak event, dfoliatR composites the individual tree defoliation series into a site-level chronology with the outbreak() function. Users have options to define the number and/or the proportion of trees required for an event to be considered an outbreak. Three parameters control whether a defoliation event constitutes an outbreak: the minimum number of trees available, the minimum number of trees recording defoliation, and the percent of trees recording defoliation. The first allows the researcher to make a judgment call as to the confidence ascribed to reduced sample depth toward the ends of their chronologies, thus compensating for the "fading record problem" (Swetnam et al., 1999). The second two parameters adjust the scale of defoliation considered to be an outbreak. Absolute numbers of trees and percentages can be applied separately or in conjunction, following filtering conventions in tree-ring fire history studies (Malevich et al., 2018). We urge users to carefully consider the choice of absolute numbers in situations where the number of trees represented in the series varies with time, or the choice of percentages when sample size is small.

Approach
We tested dfoliatR against OUTBREAK by comparing NGSI to OUTBREAK's normalized corrected indices for individual trees and years, defoliation status for individual trees and years, and percentage of trees recording outbreaks at the site level. Our tests used standardized ring-width data from eight host-tree sites spanning the range of WSBW.  (Ryerson et al., 2003) are listed on the Y-axis, dotted lines represent the series length for each tree, and colored segments show periods of defoliation. The colors of defoliation segments represent its severity, for which users can define cut-off values to determine severe-moderate-minor defoliation intensities. The default break points for severity classes are the mean and first quartile for event NGSI values. Colors and other features of the graphic can be adjusted using ggplot2 parameters, as shown below. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.
Default dfoliatR graphic for site-level outbreak events, produced by the plot_outbreak() function. These plots summarize the data in Fig. 1 for the DMJ site. Top panel shows the sample depth, the middle panel shows the mean GSI with inferred outbreak events filled, and the bottom panel shows the percent of trees defoliated, used to identify time periods of inferred outbreaks.
The sites were sampled in British Columbia (Axelson et al., 2015), Wyoming (Axelson et al., 2018), Colorado (Ryerson et al., 2003), and New Mexico (Swetnam and Lynch, 1993). These host data were compared to non-host chronologies from the original studies, but we made no effort here to replicate the reconstructions or analyses of those studies.
We detrended host data for both dfoliatR and OUTBREAK using ARSTAN (version 6.1) with cubic smoothing splines (50% frequency response on 100-150 year wavelengths depending on the site). In both dfoliatR and OUTBREAK we used event thresholds of − 1.28 normalized indices, 8 years duration, and allowed for events at the end of series in seven of eight sites that had known outbreaks at the time of sampling. We found it necessary to be consistent in how we detrended and what software we employed (e.g., ARSTAN vs. dplR) because subtle differences in standardized ring-width indices generated between the programs transferred into differences between dfoliatR and OUTBREAK. In the end, we chose to only use the standardization output files from ARSTAN, which are easily read into R (and then dfoliatR) using the dplR package.
The R code to replicate our comparisons is available from https://gi thub.com/chguiterman/dfoliatR_paper.

Findings
Across the 43,280 ring-width indices from 222 trees included in our evaluation, we found that dfoliatR and OUTBREAK compute identical growth suppression indices at 0.00 precision. We expected this outcome because both programs apply Eq. (1) to calculate disturbance indices. At the tree-level, the programs identified 11,530 total index years with defoliation. The programs agreed on 97.9% of the years, leaving 927 "difference" years in which only one program identified defoliation on an individual tree. The differences included 102 events on 85 trees. We carefully inspected each of these events in the full context of each tree's ring-series, and categorized the differences as follows: • Series-end events (40% of the total) in which OUTBREAK included "truncated outbreaks" (for seven sites) at the end of each series. In dfoliatR, this option is controlled by the "series_end_events" parameter to defoliate_trees(). In OUTBREAK, the option appears while changing the duration parameter (option 3). When selected, OUTBREAK will include any sequences of negative indices at the beginning and the end of each tree series as a defoliation event, without consideration of either duration or magnitude thresholds. In dfoliatR, the duration threshold is omitted and the magnitude threshold is retained in series-end-events. Each of the 13 events included in these differences did not meet the "max_reduction" parameter (− 1.28 NGSI) in dfoliatR and were excluded. In two cases, OUTBREAK included events at the beginning of the series where dfoliatR does not allow truncated events. In four cases, OUTBREAK omitted only the last year of the series because the index was positive, but dfoliatR allowed this single positive excursion. Finally, there were two cases in which dfoliatR omitted possible events because it had already included a positive NGSI excursion after the "max_reduction" year, and since it will only allow one excursion on either side of the max year, the events were omitted due to short duration. • Sequential events (36%) in which OUTBREAK omitted back-to-back events that occur one year prior to, or one year following an identified event. When this occurs, OUTBREAK selects the one event sequence with the lowest negative index year (e.g., Fig. 3). On two trees, OUTBREAK omitted two of three sequential events. While inspecting these differences, we added an option to defoliate_trees() that would "bridge" between sequential events (that each surpass the magnitude and duration thresholds) into single, long events. We felt that this was ecologically justified, especially for studies of WSBW, because outbreaks are known to be of long duration and tree-ring reconstructions have shown that outbreaks may persist for as long as 30-50 years at the site level. • Undetermined differences (22%) occurred in cases where OUTBREAK omitted events without clear cause that dfoliatR correctly identified as defoliations. • Rounding differences (2%) in the indices either omitted or cut short events on two trees. In both cases the indices were very close to zero, and the difference was less than the precision of the raw data measurement.
At the site level, OUTBREAK and dfoliatR produce similar time series of percent trees defoliated (Fig. 4), which forms the basis for inferring outbreak occurrence, intensity, and duration. In nearly all sitelevel comparisons, dfoliatR included either more events or it inferred a longer duration outbreak. These differences arise from the inclusion of tree-level events by dfoliatR that were omitted by OUTBREAK (see note on sequential events above). Thus, in dfoliatR, there were a greater number of trees experiencing defoliation during outbreak periods, or outbreaks were represented by a single tree when there was low sample depth.
This comparison revealed what we believe are shortcomings in how OUTBREAK identifies defoliation events on individual trees. In every one of the 102 cases we inspected, we felt that dfoliatR provided a more biologically and statistically appropriate assessment of defoliation, translating to more robust inferences of outbreak events and associated statistics at the site level.

Availability and installation
The dfoliatR library is provided free and open source from the Comprehensive R Archive Network (CRAN; https://cran.r-project.org/). To install dfoliatR from CRAN use Fig. 3. Examples of inferred defoliation events on individual trees. For each tree (DMJ26 and EFK22), OUTBREAK and dfoliatR identify most of the same events, but there is one added event (in blue) that was omitted by OUTBREAK. These were omitted because they were separated by a single year of positive normalized growth suppression index (NGSI) and OUTBREAK selected the one event with the lowest maximum departure value. dfoliatR provides an option to bridge these sequential events into single long events that may better represent the duration of defoliation given the insect and sites under consideration. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) In each R session, dfoliatR can be loaded via Development versions of dfoliatR are available on GitHub and installed using the devtools library, Issues, bug reports, and ideas for improving dfoliatR can be posted to https://github.com/chguiterman/dfoliatR/issues. As an open source library, we welcome and encourage community involvement in future development. The best ways to contribute to dfoliatR are through standard GitHub procedures or by contacting the corresponding author.

Example usage
Once dfoliatR is loaded into an R session (via library(dfo-liatR)) users can access two sets of tree-ring data to aid in exploring the functions, graphics, and outputs. Each data set consists of individual host-tree series and a local non-host chronology. The host-tree series were standardized using 128-year splines with a 50% frequency response, while the non-host ring-width data were standardized using 150-year splines with a 50% frequency response and then averaged via Tukey's biweight robust mean procedure. Host trees from Demijohn Peak site (DMJ; 2902 m asl) in the San Juan Mountains of southern Colorado include Douglas-fir (Pseudotsuga menziesii) compared against a local non-host ponderosa pine (Pinus ponderosa) chronology (Ryerson et al., 2003). The East Fork site (EFK; 2580 m asl) in the Jemez Mountains of north-central New Mexico includes Douglas-fir and white fir (Abies concolor) host trees and a ponderosa pine non-host chronology (Swetnam and Lynch, 1993).
With dfoliatR loaded, the datasets are accessible using the data () function. The data object names are prefixed by their site codes. For instance, the dmj_* objects come from the DMJ site and include the host-tree series (dmj_h), the non-host chronology (dmj_nh), the defoliation series (dmj_defol), and the outbreak series (dmj_obr). The same suite of data are available for EFK using the efk_* prefix.
In our example scripts below, ## and # denote user comments, per standard R coding, which are colored in brown. Text in blue denotes functions; black are loaded objects, and green are quoted variables and links. Values or other information provided after equal signs are filenames and parameters provided for this example, and in actual use would be replaced with user-specified information. In this example "dmj_h" and "dmj_nh" are the individual-tree host series and non-host site chronology files for the Demijohn site, thresholds are set at 8 years and − 1.28 standard deviations, bridging is used, series-end events are included in the interval computations, and comprehensive results information is not included in the output.

Tree-level defoliation events
The function defoliate_trees() performs the GSI indexing procedure on each host-tree series and then identifies defoliation events. Fig. 4. Comparison of reconstructed western spruce budworm outbreaks computed by dfoliatR and OUTBREAK. Input parameters were identical between programs. Differences arise because dfoliatR will identify and record more defoliation events on individual trees.

C.H. Guiterman et al.
The result is a long-format (stacked) data frame with five variables: "year", "series", "gsi", "ngsi", and "defol_status." The "defol_status" column indicates whether that year has defoliation or not, with a set of factors that include "nd" for non-defoliation year, "defol" for a defoliation year, "max_defol" for the year of maximum suppression (that acts as the basis for individual events), "bridge_defol" to identify years that link subsequent events (only one is present at DMJ), and "series_end_defol" to identify defoliation at the present-end of the series.
Selecting list_output = TRUE in defoliate_trees() provides a list-object of data frames, each with an rwl object that combines the host tree and non-host series and the other columns created by defo-liate_trees(). This option is not used by subsequent functions in dfoliatR, but researchers can examine it to check the results of the GSI calculation (Eq. (1)), such as the non-host series after scaling to a common variance with a particular host-tree series.
The results of running defoliate_trees() can be assessed through graphical and table outputs. The function get_defol_events () will provide a list of every defoliation event for every tree, with the corresponding mean "ngsi" value. A summary table of the results for each tree is produced by defol_stats() (Table 1).
The plot_defol() function produces a "ggplot" graphics object with line segments showing the measured sequence of each series and a filled segment for each identified defoliation event (Fig. 1). The defoliation segments are colored by their relative severity based on their average NGSI value. By default, plot_defol() will calculate the average NGSI for all identified events, and assign severity based on the mean and first quartile of the averages. "Severe" events have a mean NGSI above the overall average event-period NGSI. "Moderate" events fall between the mean and first quartile. "Minor" events fall below the first quartile. Users can re-define the breaks to suit their needs via the "breaks" parameter in plot_defol().
These output functions aid in assessing the sensitivity of input parameters to defoliate_trees(), including the duration and magnitude thresholds for identifying defoliation events. Using plot_defol () also provides a direct assessment of the between-tree variability in defoliation.

Site-level events
To infer outbreak events at the site level, the function outbreak() composites tree-level defoliation series into a single chronology, with input parameters that control thresholds in the number and proportions of trees recording a defoliation event.
Input parameters to outbreak() include "filter_min_series" to control the chronology cut-off points with regard to sample depth, "fil-ter_min_defol" and "filter_perc" to control the minimum number and percent of trees recording a defoliation event in a given year. outbreak () produces a new data frame with eight variables: "year", "num_defol", "percent_defol", "num_max_defol", "mean_gsi", "mean_ngsi", and "out-break_status." All of these variables are populated regardless of an inferred outbreak event, providing a continuous outbreak reconstruction. The "num_max_defol" variable counts the number of trees recording their maximum defoliation in a given year. The "mean_gsi" and "mean_ngsi" variables provide averages of these indices across all available trees. Finally, the "outbreak_status" column shows if an Table 1 Tree-level tabular output provided by the defol_stats() function for the DMJ example site. Note that these calculations exclude the ongoing "series-end" events as selected in defoliate_trees(). outbreak event is inferred ("outbreak") or not ("not_obr"), and whether it represents an ongoing series-end event ("se_outbreak"). The default plotting function to visualize results from outbreak() is plot_outbreak(). It creates a three-panel graph showing the sample depth, mean site-level chronology, and percent of trees recording a defoliation over time (Fig. 2).
Inferred outbreak events are shown in the middle panel of Fig. 2 as the filled-in spaces. Users can change the time series in this panel with the "disp_index" parameter, choosing between the mean NGSI (the default) or GSI.
A summary table of the inferred outbreak events is generated by the outbreak_stats() function ( Table 2). The table provides a range of summary statistics, including the start and end years of each outbreak event, along with the corresponding duration, the number and percent of trees in defoliation at the start of the event ("n_df_start" and "perc_df_start", respectively), the maximum number of trees recording the outbreak event during a single year ("max_df_obr"), the year corresponding to that peak ("yr_max_df"), the year with the lowest mean NGSI during the event ("yr_min_ngsi"), and the minima of mean GSI and mean NGSI indices during the event.
Saving the results of outbreak_stats() (the dmj_obr_stats object above) provides an array of options for assessing metrics of the insect outbreak regime. For example, taking the first year of each outbreak event, we can calculate the duration of years between outbreaks, via the diff() function in R. The average of those differences, calculated via mean() is the mean return interval of reconstructed outbreak events at the DMJ site.

Conclusions
The dfoliatR package provides dendroecologists with tools to infer, quantify, analyze, and visualize tree-ring growth suppression events and to reconstruct forest insect defoliator outbreak regimes. It is built on the long-accepted host to non-host comparison methodology used in the 1980s FORTRAN program OUTBREAK (Swetnam et al., 1985;Lynch, 1989, 1993). Our evaluation of the two programs revealed that dfoliatR excels in identifying defoliation events on single trees, providing researchers with more consistent and biologically-justifiable results. dfoliatR provides easier control of the rule base for suppression thresholds, additional output tables, and high-quality and customizable graphics. These features allow users to compare insect outbreak regimes of different tree species or geographic regions, evaluate sample-size considerations, examine a multitude of relevant insect disturbance questions, and more readily evaluate the potential for Type I and II errors in their results. Finally, dfoliatR operates in the open source R environment that is stable across computing platforms and is under active development and maintenance by a large and growing community.
Using dfoliatR requires standardized ring-width measurements from insect host trees and either an indexed tree-ring chronology from local non-host trees or suitable climate chronology. It performs an indexing procedure to remove the climatic signal represented in the nonhost chronology from the host-tree series. It then infers defoliation events in individual trees. Site-level analyses identify outbreak events that synchronously affect a user-defined number or proportion of the host trees. Functions are provided for summary statistics and graphics of tree-and site-level series. The package produces publication-quality plots, and tabulates growth suppression indices and tree-and site-level outbreak event statistics for user-defined post-processing needs, including those suitable for charting and tabulating landscape-and regional-level results.
dfoliatR adds a new option for dendroentomology to combine, or "bridge," sequential tree-level defoliation events into single events. In practice, we suggest that researchers carefully evaluate if bridging is ecologically applicable to their study system and insect ecology, and to carefully explore the data before deciding whether or not to use this option. It is probably not appropriate for insects with high-frequency, high-severity outbreaks, such as processionary caterpillars. Outbreak reconstructions of insects for which the interval is notably longer than typical outbreak duration, such as Douglas-fir tussock moth, are unlikely to be significantly affected. Species for which outbreak duration or individual tree resilience to defoliation varies considerably, or which may chronically infest trees or sites, or alternate between chronic, outbreak, and minimal activity states such as conifer-feeding Choristoneura, present more complicated challenges, and the researcher should use discretion with the bridge option.
dfoliatR adds to the on-going open source software development for dendrochronological methods (e.g., Bunn, 2008;Brewer, 2014;Brewer and Guiterman, 2016). The R environment enables automation of analyses, allowing input/output processes to become routine, enables efficient sensitivity analyses, and empowers batch processing of large multi-site projects. It also facilitates additional statistical analyses, such as spectral analyses and superposed epoch analyses (e.g., Malevich et al., 2018), with easy transfer from dfoliatR and dplR to other libraries in R. Source code for dfoliatR is available in the Comprehensive R Archive Network (CRAN) and GitHub https://github.com/chguiterman /dfoliatR with updated descriptions and helpful vignettes on the package website https://chguiterman.github.io/dfoliatR/. Researchers wishing to contribute to the further development of dfoliatR are encouraged to do so via the GitHub repository.