Survival-Time Analysis of White Spruce during Spruce Budworm Defoliation

Mortality and defoliation (DF%) in 987 white spruce (Picea glauca (Moench) Voss) trees were followed from 1992 to 2003 during an outbreak of the spruce budworm Choristoneura fumiferana (Clem.) in 15 white-spruce-dominated uneven-aged stands in the Fort Nelson Forest District near Prince George, British Columbia. Four stands were aerially sprayed with Bacillus thuringiensis (Bt). Defoliation and mortality levels were elevated in non-sprayed stands. The relationship between defoliation and survival-times was captured in a Cox proportional hazard model with a defoliation stress index (DSI), diameter (DBH), crown class (CCL), a random stand effect, Bt-treatment, and number of years of exposure to stand-level defoliation (DYEAR) as predictors. The DSI, optimized for discrimination between survivors and non-survivors, is the discounted sum of five lagged DF% values. Survival probabilities were predicted with a maximum error of 0.02. Hazard rates increased by 0.06 for every one point increase in DSI. CCL and random stand effects were highly significant. Bt-treatment effects were fully captured by DSI, CCL, and DYEAR.


Introduction
The spruce budworm, Choristoneura fumiferana (Clem.), is the most widely distributed native defoliator of balsam fir, Abies balsamea (L.) Mill., Abies spp, and white spruce, Picea glauca (Moench) Voss, forests in North America.This budworm inhabits the boreal forest, extending from the Atlantic coast westward to Alaska and north to the Mackenzie River delta in the Northwest Territories.
Prolonged budworm outbreaks cause a significant increase in tree mortality and growth reduction in Canadian forests (Royama 1984, Solomon et al. 2003, Steinman and MacLean 1994).Past cutting practices and fire suppression may have exacerbated the impact of this defoliator (Blais 1983, Hummel andAgee 2003).
Most budworms share a similar life cycle (Sanders 1992).Eggs are laid in masses on host foliage in mid-summer.Newly hatched larvae move to sheltered locations on the branches where they establish hibernacula, enter obligatory diapause, and pass the winter.Larvae become active the following spring, feed first on pollen in male strobili or mine needles, but soon migrate to the new expanding shoots and form feeding shelters.For most species, feeding is completed in about six weeks, usually by early July.Mature caterpillars feeding on the current year's growth cause most damage.Outbreaks can last as long as 15 years and may be synchronized over a large geographic area (Royama 1984).Feeding by budworm larvae results in severed and partly consumed needles which, when dry, become red and give the trees a scorched look that is easily recognized in aerial photographs (Fig. 1).
The susceptibility of forest stands to budworm has been quantified as a function of species composition, site factors, stand structure, density, crown closure and climate (Alfaro et al. 2001, MacKinnon and MacLean 2003, MacLean and MacKinnon 1997).Impacts of defoliation on growth and yield have been researched in detail (Krause andMorin 1999, MacLean andMacKinnon 1996) and found to be substantial.Extensive records of tree defoliation status and stand-level mortality rates have allowed reliable forecasting of the expected (5-or 10-year) mortality (by tree size) in stands with a known defoliation history (Erdle and MacLean 1999, Piene et al. 2003, Steinman and MacLean 1994).
These impact studies provide the information required to intervene to prevent mortality and salvage dead timber.However, the trend towards a more holistic approach to forest management (Bergeron et al. 1999, Galindo-Leal and Bunnell 1995, Kohm and Franklin 1997) emphasizes the implications of the budworm within an ecosystem context.This approach requires a more detailed understanding of the spatial and temporal aspects of forest disturbances (Almquist et al. 2002, Schweiger andSterba 1997).For example, tracking the temporal trend in tree mortality during a spruce budworm outbreak cycle may be required for assessment of changes in stand structure and succession that could affect specific ecosystem functions, such as the supply of habitat for cavity-nesting birds.
While the general relationship between defolia- tion and tree mortality is well understood, there is a paucity of models that explicitly link the survival (death) times of a single tree during a defoliation cycle as a function of current and past defoliation levels.The distribution of survival times provides insight into the dynamics of the interaction between time, defoliation and mortality in 'real time', which is critical for the development of ecological process-oriented impact models.Volney's study of jack pine, Pinus banksiana Lamb., mortality during the first ten years of a budworm outbreak (Volney 1998) and a study on host-insect interactions by Nealis and Régenière (2004) are exceptions.Logistic (Monserud and Sterba 1999) or classification and regression models of tree mortality (Dobbertin and Brang 2001) optimize the prediction that a tree will survive during a fixed (constant) period of t years given a set of fixed predictor values but they do not predict the distribution of survival times during an outbreak cycle.To do that they would have to have timevarying coefficients and the estimation would be exceedingly complex due to the longitudinal character of the data (Hosmer and Lemeshow 1999).
In this study, we develop a Cox proportional hazard regression model for the survival times of white spruce trees exposed to defoliation by spruce budworm in the Prince George Forest Region of British Columbia, Canada.The effect of current and past tree defoliation values is captured by an optimized discriminatory defoliation stress index.The intended use of the model is for dynamic forecasting of mortality during an outbreak of budworm defoliation.

Study Design
The annual extent of defoliation and mortality in 987 white spruce trees was monitored between 1992 and 2003 in 15 white spruce-dominated mixed stands of spruce, fir, and aspen in the Prince George Forest Region of British Columbia.The study area comprised 1.5 million ha near the town of Fort Nelson in the boreal white and black spruce biogeoclimatic zone (Meidinger and Lewis 1983).Geographically the area is between 122°00'00"W and 124°00'00"W longitude and latitudes 58°24'00"N and 59°36'00"N.Defoliation history for the area prior to plot establishment (1986 to 1991) was obtained from the annual aerial Forest Insect and Disease Survey (FIDS) maps of defoliation provided by the Canadian Forest Service.Stand-level estimates of budworm defoliation -light (1 to 25%), moderate (26-65%), and severe (66-100%) -obtained by observations by trained personnel from fixed-wing aircraft were available for the entire duration of the current outbreak cycle.Further details of the study area can be found in Burleigh et al. (2002).
Monitored trees were sampled along five systematically located 50-m-long survey lines established in each of the 15 stands in 1991.The survey lines were located 50 m from any stand boundary and with 50 m (or more) between any two lines.The design provides for a representative sample of the stand condition.At the time of plot establishment (spring of 1992) the diameter at breast height (DBH) and crown class (CCL) of dominant, co-dominant, intermediate, and suppressed trees was recorded for each tree located within 1 m from the survey line.There was an average of 11 trees per survey line (range 9-30), and 66 trees per stand.

Stand Treatment
Four of the stands had a history of Bt (Bacillus thuringiensis Berliner) spraying prior to 1992.Two stands were sprayed twice (1991 and 1992) and two were sprayed four times (1989)(1990)(1991)(1992) during the first years of budworm defoliation.The primary reason for the spray operation was to secure the production of seed for local reforestation use.In the current context of survival-time analysis it was deemed appropriate to collapse the Bt-treatments into a binary "BT2+" for the four Bt treated stands and "BT0" for the 11 non-Bt treated stands.The difference in mortality between stands treated twice and four times were too small to consider them as separate treatments.Total mortality between 1992 and 2003 in the two stands treated four times with Bt was 21 (1.1%) while the mortality in the twice treated stands were 14 (1.7%).A t-test under the null hypothesis of no difference failed to reject the null hypothesis of no difference (P = 0.22).Year specific mortalities in the two Bt treatments varied from 0 to 5 with an annual average difference of 1.4.Treated stands had, on average, a larger DBH (31 cm versus 23 cm, P < 0.000, two-sample t-test), and a higher proportion of dominant trees (44% versus 34%, P < 0.000, Chi-squared test).Other characteristics like ground vegetation and soils were deemed similar.The relative crown class distribution and DBH by treatment group is in Table 1.

Defoliation and Tree Status
Defoliation during the period of maximum defoliation visibility (July-August) was assessed by visual examination of trees using binoculars (Shore et al. 1988).In the first year of assessment each tree crown was divided into thirds (top, middle, bottom) and each third was assigned an estimated cumulative percent of foliage missing from all foliage age cohorts due to current and past budworm feeding (DF%).Thereafter, once a year, individual tree crown sections were reassessed to determine the percent of budworm defoliation of current-year foliage.A weighted average of the current-year's and last-year's defoliation value was then computed to obtain an estimate of the current defoliation value by crown third.Current-year foliage in the top, middle, and bottom section was given a weight of 0.3, 0.25, and 0.25, respectively.The weights reflect the relative contribution of new foliage to the overall foliage pool (Alfaro et al. 2001).The same observer recorded defoliation for 10 out of the 11 years of observation.
At the time of defoliation assessment the tree status (alive/dead) was also recorded.A tree with no live foliage was declared dead and the year of dying as the last year the tree displayed live foliage.Due to logistical constraints no defoliation assessments were made in 1993 and 1998, only tree status (dead/alive) was recorded during those years.

Imputation of Defoliation Values
A complete defoliation record for each tree would greatly facilitate the computational aspects and interpretation of survival-time analysis (Hosmer and Lemeshow 1999).Missing records can be treated as periods when the observational unit 'left' the study or, alternatively, be replaced by imputed values (Rubin 1987).We opted for the latter for sake of convenience.Missing DF% values for 1993 and 1998 were imputed as the DF% value for 1992 viz. 1997 plus a random proportion (p) of the apparent change ΔDF% in DF% between 1992 and 1994 viz.1997 and 1999.The distribution of p was triangular with a mean of 0.5.Based on sequences of 3, 4, and 5 years of observed DF%, we infer that imputation errors should be less than 10% in most (90%) cases.

Defoliation Time-Lines
Aerial defoliation records for the area indicate that the current outbreak of spruce budworm defoliation began in 1986.In eight of the 15 monitored stands the first stand-level defoliation was noted in 1986.In the remaining seven stands defoliation began in 1988 (1 stand), 1989 (4 stands), and 1990 (2 stands).
The year of onset of defoliation in a stand was defined as the first year when more than 10% of the trees in a stand have a defoliation value of 10% or higher.For stands with an onset of defoliation prior to plot establishment the year of onset was obtained from the aerial survey data.
The stand-level differences in year of onset means that the length of exposure to budworm defoliation also differs among stands.To remove this heterogeneity and facilitate the survival-time analysis we created a stand-level time-variable called "Defoliation Year" (or DYEAR) which measures the number of years of exposure to budworm defoliation.All trees in a given stand will have the same DYEAR value.DYEAR is zero in any year prior to the onset of defoliation.In this study DYEAR covers a range from 0 to 17 years of exposure.

Defoliation Stress Index (DSI)
Defoliation values above a certain threshold is known to strain a tree's growth, maintenance, and defense (Piene et al. 2003).Continued defoliation aggravates the situation and may lead to a premature death of the tree (Solomon et al. 2003).
A combination of current and past defoliation values is therefore an indicator of the cumulative stress on a tree.For the purpose of predicting mortality a generic defoliation stress index (DSI) is developed.The DSI is designed for optimum discrimination between trees that died and trees that survived in a given DYEAR.DSI for tree i in DYEAR = t is a weighted sum of current and past DF% by crown section.The effect of past DF% on current stress is assumed to depreciate at a compound rate with a constant annual discount factor λ (0 < λ < 1).A high λ means that the effect of past DF% is discounted faster than with a lower λ value.The number of immediate past years of DF% values to include in the index is k, with k to be determined but constrained to be less than or equal to five years.The limit imposed on k is largely determined by the length of the available defoliation time series.Past DF% values for years prior to the first observations in 1992 were needed to compute some DSI values.They were obtained by linear interpolation between the DF% value recorded for 1992 and the presumed value of 0% in the last year before a stand level onset of defoliation.
The relative contribution of each crown third to the overall stress index of a tree is determined by three weight factors (0 ≤ w c ≤ 1, c = {top, middle, bottom}) summing to one.With these notational conventions the stress index of tree i in year t is Note, a discount factor of λ = 1 implies that only the current defoliation value enters the index.The estimates of the index parameters were those that maximized, on a logistic scale, the average (over all DYEAR) DSI-difference between trees that died and those that survived in a given DYEAR.The solution was

Probability of Survival and Hazard Rates
The probability of a tree surviving the first t DYEAR of an exposure to a stand-level defoliation was computed as the Kaplan-Meier estimator (Cleves et al. 2004, p. 93) ˆ( ) where n j is the number of live trees at risk (exposed) at DYEAR = t j , and d j is the number of trees that died (failed) at DYEAR = t j .The product is over all observed DYEARs less than or equal to t. Survival probabilities in Eq. 2 were computed separately for the BT0 and the BT2+ stands.Confidence intervals of ˆ( ) The intensity of the mortality process (hazard rate) in year DYEAR = t (haz(t)) is the difference between the cumulative hazard rate of that year (H(t)) and that of the preceding year (H(t -1)) where H(t) is the Nelson-Aalen estimator of the cumulative hazard (Cleves et al. 2004 with an estimated variance of (5) Apparent random fluctuations in the discrete approximation to the continuous hazard rate were smoothed by a Gaussian kernel with bandwidth of 1.25 (Silverman 1986).An estimator of the standard error of a hazard rate in DYEAR = t was obtained as d t × n t -2 .

Cox Proportional Hazard Regression
The effects of DSI, crown class, DBH (in 1991), Bt treatment, and a random stand effect on survival and hazard of an individual tree (i) in DYEAR = t was estimated with a Cox tree-level proportional hazard regression model (Cox 1972) where haz 0 (t) is the base-line hazard function, DBH' i the deviation in DBH of the ith tree from the mean DBH of trees with a co-dominant crown class in the BT0 stands (average = 25.9 cm), δ(Bt) is a Bt treatment indicator taking the value of 1 if the Bt treatment of the ith tree is BT2+ and 0 otherwise, and ν s(i) is a random effect of the stand s (s = 1,....,15) in which the ith tree is located.It is assumed that the distribution of the stand effects is Log-Gamma with a mean of 0 and a variance Log [φ] to be estimated from the data.The random effect introduces a positive correlation of survival times of trees in a stand since all trees in the stand 'share' the same random effect.The cause behind a stand effect cannot be determined in this retrospective study; essentially, a random stand effect is the part of the non-explained hazard residual that is common to all trees from a single stand.In Eq. 6 CCL was coded as -1 for dominant, 0 for co-dominant, 2 for intermediate and 3 for suppressed.Accordingly, the base-line hazard function is for a non Bt-treated co-dominant tree with a DBH of 26 cm and a DSI value of 0. Regression coefficients in Eq. 6 were estimated by a penalized maximum partial likelihood method treating estimates of ν s(i) and φ as known).Variance and covariance of the coefficients were obtained from the inverse of the negative information matrix (Kalbfleisch and Prentice 2002).
We chose Efron's method for handling tied failure times (Efron 1977).Alternatives, such as, for example, Breslow's, partial, exact-partial, or conditional logistic methods gave virtually identical results.Schoenfeld residuals and Martingale residuals (Hosmer and Lemeshow 1999) were used to assess goodness of fit and model assumptions.

Results
The average DF% value of a live tree in the BT0 stands fluctuated between 48 and 64 (mean 56) in the 2 to 17 years since exposure to stand-level defoliation (Fig. 2).Note, the first DF% entry in Fig. 2 is for DYEAR = 2 since no observations were made for the first 2 DYEARs.A similar offset applies throughout if not otherwise indicated.As the regional budworm population appears to wane in DYEAR = 15, tree-level values of DF% begin to decrease.Defoliation values in BT2+ stands were, as expected, consistently lower (14-36%, mean 22%), especially in the first 6-8 DYEARs.
With an exception of DYEAR 11 and 12 the differences in mean defoliation values between BT0 and BT2+ stands were statistically significant (P < 0.01, two-sample t-tests with a Bonferoni-type adjustment of the significance aimed at controlling the family-wise Type I error rate to 0.05 or less (Miller 1980)).During DYEAR 1-7 the defoliation values in the mid and lower crown sections of live trees were 5-10% higher than in the top crown section.However, these differences faded after 7 years of exposure.Stand effects within each Bt treatment group were significant (P < 0.01, F-test of square-root transformed DF% values).The pooled among-stand standard deviation of DF% was about 5.
During the first 7 DYEARs over 90% of the trees in the BT0 stands had one or more crown thirds with DF% ≥ 10; in later years only 70-80% met this criterion.BT2+ stands followed a different trend in DF%.During the first 4 DYEARs, a time of active Bt control, only about one third of the trees had any appreciable defoliation; in DYEARs 5-12 defoliation was visible in about 80% of the previously treated trees, and subsequently DF% declined to the levels of the first 4 years.The proportion of trees with DF% over 10 in the BT0 stands was (in all DYEARs except for the 11 th and the 12 th ) significantly higher than in the BT2+ stands (two-sample t-test, with Bonferoni-type adjustment of significance levels).
The estimated discount rate of 0.60 for the contribution of past defoliation to the defoliation stress index implies a fairly rapid write-down of past defoliation values as indicators of current stress.The discount rate was significantly less than 1.00 (P < 0.000, t-test).In contrast, the number of past years to include in the index reached the maximum of 5. Thus, the optimum value of k is greater.However, given the discount rate the truncation had only a negligible effect.DF% values from the middle and bottom crown thirds do not appear to carry any additional information relevant to mortality beyond what is captured by the topmost third.The high colinearity between section-specific DF% values (correlations were in excess of 0.95) confirmed the redundancy of the middle and bottom crown defoliation values.The maximum possible value of DSI i,t is 166.Individual tree-level values of DSI were computed for the years 1992 to 2003.
The mean defoliation stress index (DSI) of live trees in BT0 stands remained high (75-85) over time (Fig. 3) with a local maximum of 92 attained in DYEAR = 15.Trends in DSI show, as expected, less variation than trends in DF% due to the averaging effect of the index.BT2+ stands enjoyed significantly lower DSI values throughout.During DYEARs 1-9 their DSI was about 50 points below those in BT0 stands.During the next 3 years this treatment-related gap narrowed to about 30 points as the initial benefit of Bt spraying faded.The among-stand standard deviation in DSI averaged 7 points in both groups.
During the 11 years of observation a total of 281 (36%) of the 770 trees monitored in the BT0 stands died, predominantly due to severe defoliation.Comparable numbers in the BT2+ stands were 35 (16%) dead out of a total of 217 observed trees.The difference of 20% in overall mortality rates was highly significant (P < 0.01, two-sample t-test).Annual mortality rates in BT0 stands peaked at 7-9% annually in DYEARs 9-12.Annual losses in BT2+ stands fluctuated between 0 and 4%.Treatment-induced differences in annual culling rate were statistically significant in DYEARs 9, 12, and 13 only (P < 0.04, twosample t-test with Bonferoni-type adjustment of significance levels).For trees that died during the study period the DSI value, for the last year in which they were recorded as alive, was about double the average of survivors (135 versus 69).All trees that died had a stress index in excess of  100 in the last year they were recorded as alive.
The probability that a tree in a BT0 viz. a BT2+ stand survives t years of exposure to budworm is displayed in Fig. 4 along with the 95% confidence bands for the Kaplan-Meier estimators.During the first four DYEARs survival is high (> 98%).Then, between DYEARs 5 and 7 the survival probability drops to 92% for the BT0 trees and to 96% in the BT2+ trees.In DYEARs 8 to 12 the survival probability of BT0 trees drops to 66% while that of BT2+ trees drops to just 85% (P < 0.01, two sample t-test).During the last 4 DYEARs we see a drop of approximately 1% per year.A Log-rank test of equality of survivor functions (Hosmer and Lemeshow 1999) rejected the null hypothesis of no difference (χ 2 = 26.8,P < 0.01).Alternative tests produced similar results.
Temporal trends in the intensity of the mortality process (hazard rate) are in Fig. 5.As expected, the BT0 and BT2+ curves begin and end at about the same level.Whereas the hazard for BT2+ trees remains nearly constant throughout the time of exposure, that of BT0 trees first increases, then remains more or less constant between DYEARs 10 and 12, and finally drops to values estimated for the initial years of exposure.
DSI, CCL, and stand effects, but not DBH nor Bt per se, exerted a significant influence on the hazard rate of failure (death).Statistics of the coefficients in Cox's proportional hazard regression model are summarized in Table 2.A global test of the proportional hazards assumption based on Schoenfeld residuals did not detect any major violation (P = 0.14).Nor did an inspection of Martingale residuals reveal any non-quantified time trend in the covariates.None of the tabled results are sensitive to the way tied failure times have been handled.The regression coefficient for DSI suggests that a one-point increase in DSI increases the hazard by an average of about 6%, albeit with a downward trend as DYEAR increases.
An effect of crown class on hazard was apparent.Suppressed trees had a mortality rate that was about 50% higher than otherwise expected under the null hypothesis of no crown class effect.Conversely, dominant trees died at about one-third the rate expected under the same null hypothesis.A change of one crown class shifted hazard by about 25%.
No isolated Bt treatment effect on hazard (P = 0.18, Chi-squared likelihood ratio test) was  identified after the effects of DSI and CCL were accounted for.Nor did we find any significant interactions between CCL and DSI (P > 0.55, Chi-squared likelihood ratio test).Random stand effects explained a significant part of the tree-level variation in hazard rates.Predictions of P surv for an individual tree in DYEARs 2 to17 indicated a satisfactory fit to the empirical Kaplan-Meier estimators (Fig. 6).The fit to BT0 data was especially favorable with a median absolute deviation (MAD) of just 0.005 compared to a MAD of 0.019 for BT2+.The nearly 4:1dominance of BT0 observations explains, by and large, this differential in goodness-of-fit.Mortality (1 -P surv ) predicted for base-line trees with no defoliation stress (DSI = 0), a co-dominant crown class, and a diameter of 25.9 cm was consistently low (< 2%).
DSI does not capture all the stress of a tree; time since exposure is also needed to describe time trends in hazard rates.Hazard rates predicted for a co-dominant tree with an average DBH of 25.9 cm is illustrated in Fig. 7 for constant DSI values of 60, 100, and 140.The hazard for a given DSI value first increases and then decreases with DYEAR.A hazard rate above 0.50 is only predicted for DSI ≥ 140.

Discussion
Variable time-line survival data can be difficult to analyse outside the framework of survival-time analysis (Hosmer and Lemeshow 1999).We chose the Cox proportional hazard regression model as it requires no assumption about the baseline hazard rate.An accelerated failure-time model, in which the impact of covariates is measured through their effect on time would have produced very similar results if we assumed that conditioned on the covariates survival-times follow a known distribution (Kalbfleisch and Prentice 2002).A distinct shortcoming of survival-time analysis is the current limitation to a single random effect, when a multi-tiered model with random effects of stands, survey line, and trees would appear desirable.Methods for overcoming this limitation are emerging (Guo and Carlin 2004).
The tendency of the western spruce budworm to repeatedly defoliate the same trees for many years in a row (median 13) combined with a strong yearto-year correlation of defoliation levels (0.60) and moderate year-to-year fluctuation in defoliation values suggests that a current defoliation value would be a reasonable predictor of future defoliation values during an outbreak cycle.We exploited this opportunity in our imputations of missing defoliation values.
Defoliation, years of exposure to defoliation and crown-class were confirmed as highly significant predictors of survival times (Volney 1998) and, thus, mortality (Alfaro et al. 1999, Erdle and MacLean 1999, MacLean and MacKinnon 1996, Steinman and MacLean 1994).A tendency of budworms to drop from overstorey to understorey trees creates a disproportionate defoliation 'stress' on small-crowned trees.A probable cause of the reported relationship between crown class and hazard.The gradual rise in hazard rates to a peak 11-13 years after the onset of stand-level defoliation matches the general time profile of mortality rates during a spruce budworm defoliation outbreak (Alfaro et al. 1999, Solomon et al. 2003, Steinman and MacLean 1994).Defoliation stress index and crown class do not explain, fully, the time trends in observed mortality.The effect of current and past defoliation values is modified by the time of exposure.
The shape of the predicted DYEAR-profile of hazard has also been reported for shorter duration (2-6 years) defoliation cycles (Cedervind and Langstrom 2003, Langstrom et al. 2001, Piene et al. 2003, Volney 1998).Volney (1998) found the mortality force (hazard) to be 19 times higher for a heavily defoliated jack pine (DF% > 75) than for a lightly defoliated tree (DF% < 25%); our results suggest that a very similar ratio applies to white spruce.
Sprays with Bacillus thuringiensis appear to have had lasting benefit.An extended "protection" of Bt-sprayed trees may at first seem unlikely.Persistent DF% values in Bt-treated stands indicate, however, that budworm populations were present and not completely wiped out by Bt.We hypothesize that the apparent lasting effect may simply be due to the (repeatedly) sprayed trees having larger, healthier crowns, which, everything else being equal, may sustain less damage for a given budworm population density (MacKinnon and MacLean 2003).
The role of past defoliation values on survival times, although understood in general terms as an added stress factor augmenting mortality rates (Cedervind and Langstrom 2003, Pedersen 1998, Steinman and MacLean 1994), is difficult to quantify accurately.Our defoliation stress index was optimized for discrimination of a binary mortality event.Thus the carry-over of discounted past defoliation values confirmed that while the current defoliation value is the most important predictor of stress, information on past defoliation values and the number of years of exposure is also needed to improve prediction of hazard and survival times.The significant effect of exposure time on hazard (and mortality) for a fixed set of predictor values illustrated the importance of exposure time.We can only speculate on the nature of this interaction.Secondary stressors (bark beetles, fungi) may play an important role for accelerating mortality (Alfaro et al. 2001, Pedersen 1998).
Significant stand-level effects on mortality are to be expected.Stem density, competition, species mixture, age, crown closure, and site quality have been reported to influence mortality (Alfaro et al. 2001, MacKinnon and MacLean 2003, MacLean and MacKinnon 1997).Hence, an among-stand variation in any of these factors will manifest itself as a stand effect in our study.Accurate stand-level predictions of survival times is thus only possible if one can predict the stand effect.The non-replicated and retrospective nature of our study precludes such a prediction.
Survival-time and hazard-rate models obtained from this study lend themselves to integration with stand-level process-oriented models of budworm defoliation outbreaks and their effect on individual trees (Crookston 1991, Greenbank et al. 1980, Hassell et al. 1999, Krause and Morin 1999, Piene et al. 2003, Royama 1984, Volney and Cerezke 1992).Single time-step predictions of defoliation status of individual trees are then used, along with other tree attributes, to predict the DSI and then the chance that a tree will survive the current year.Model trees can then be culled in a manner consistent with these predictions.This process is repeated for another year and so on.
Survival-time predictions are important for stand-level management.Thresholds for intervention and modification of thinning operations can be worked out once the expected outcome of a defoliation cycle on mortality is quantified.While our results are specific to the spruce budworm we expect that a survival-time analysis of shorter-duration defoliation cycles caused by different insects will be equally useful in generating a process-oriented quantitative assessment of hazard rates, survival, and mortality during an insect outbreak.

Fig. 2 .
Fig. 2. Average tree defoliation values (DF%) in Bacillus thuringiensis treated (BT2+) and non-treated (BT0) stands of white spruce plotted against years of exposure to stand-level budworm defoliation (DYEAR).Vertical bars delineate the interval of a mean ±1 × the among-stand standard deviation.

Fig. 3 .
Fig. 3. Average defoliation stress index (DSI) in Bacillus thuringiensis treated (BT2+) and non-treated (BT0) white spruce stands plotted against years of exposure to stand-level defoliation (DYEAR) Vertical bars delineates the interval of a mean ±1 × the among-stand standard deviation.

Fig. 4 .
Fig. 4. The average tree-level probability of surviving (P surv ) DYEARs of exposure to budworm defoliation in Bacillus thuringiensis treated (bottom) and non-treated (top) stands.Dashed lines indicate the 95% confidence interval of the Kaplan-Meier estimators.

Fig. 5 .
Fig. 5. Mortality hazard rates (haz.)plotted as a function of DYEAR (years of exposure to stand-level defoliation) in Bacillus thuringiensis treated (BT2+) and non-treated (BT0) stands.Shown hazard rates have been smoothed with a Gaussian kernel.68% confidence bands are indicated by thin dashed lines.

Fig. 6 .
Fig. 6.Predicted and empirical Kaplan-Meier (triangles and stars) probabilities of a tree surviving (P surv ) DYEARs of exposure to defoliation in Bacillus thuringiensis treated (BT2+) and non-treated (BT0) stands of white spruce.Topmost nearly horizontal and overlapping lines depict treatment specific estimates of base-line survival for co-dominant trees with DSI = 0 and a class average DBH of 25.9 cm.

Fig. 7 .
Fig. 7.An illustration of the impact of the significant interactions between defoliation stress index values (DSI) and DYEAR on hazard rates (haz.) of death at three constant values of DSI.

Table 1 .
Relative crown class distribution (%) and DBH by Bt treatment group.