Model selection analysis examining population patterns of ladybeetles Harmonia population dynamics

Yearly average per-trap captures of Harmonia axyridis from 1994-2013 from the Michigan State University’s Kellogg Biological Station (KBS) were used to model population regulation processes in this species (Landis and Gage 2014). Data collection occurred as described in Bahlai et al. (2013) at KBS in southwestern Michigan (42°24’N, 85°24’W) at 288 m elevation. Coccinellid surveillance data were collected starting in 1989 at the KBS Long Term Ecological Research (LTER) site as part of the Main Cropping System Experiment (MCSE) and nearby forest site monitoring. See http://lter.kbs.msu.edu/research/long-term-experiments/main-cropping-system-experiment/ and http://lter.kbs.msu.edu/research/long-term-experiments/successional-and-forest-sites/ for a detailed experimental design. Because the sampling period varied from year to year, data were culled at day of the year 240 to minimize the effect of sampling period on the results.

To determine if regime shifts occurred during the study period and to identify break points in the time series, if present, an iterative model-selection approach was used. First, a density-dependent model was fitted to the entire time series of ladybeetle abundance values 1994-2013. Next, the time series was subdivided into two or three subsets, and the same model was fitted to data in the subsets. To identify break points in the dynamic regime, all possible breakpoints between subsets of the time series were tested, on the condition that any subset should at least have three data points to avoid overfitting. Finally, the data from the first and last subset were combined and fitted with one parameter set to determine whether the initial and final part of the time series were characterized by the same dynamic law. This scenario was tested because initial results suggested that Harmonia population dynamics had returned to previous patterns by the end of the study (see Results). Models were ranked according to AIC.

The model fitting was conducted with least squares using minpack.lm (Elzhov et al. 2013) in R 3.0.3 (R Development Core Team 2014). Two alternative models were tested for all break combinations: the Ricker and Logistic population models (Turchin 2003). Both Logistic and Ricker models are discrete time models: N(t+1), the population at time t+1, as a function of N(t), where N(t) is the population density of Harmonia observed at our site in year t. Both models contain the parameters K (the carrying capacity) and r (intrinsic rate of increase), observed for Harmonia during that time period.

AIC was used to determine the model type and break point combination with the greatest support from the data. AIC values associated with fitting each data subset was summed for each break point combination to give a total AIC by break point combination and model type. When AIC values differed by more than two units, the model with the lower AIC value was considered to have a better fit.

Sensitivity analyses- We tested the sensitivity of our model selection analysis procedure to variation in the sampling period from year to year by running the model selection algorithm on data which were left unculled and on data culled day of the year 220 and found model selection results were identical to when data were culled at day 240, as in the original analysis, although variability increased in some models tested. Thus, we used data culled at day 240 with yearly sampling periods ranging from 8 to 19 weeks, with a mean sampling period of 13 weeks for the final analysis.

To account for week-to-week variations in captures, we also computed the integral of the population curve for each growing season.