Introduction
As an example of the use of rtrim, counts of the Skylark (Alauda arvensis) will be analysed (the data are obtained from the Breeding Bird Monitoring Scheme in the Netherlands of Sovon and Statistics Netherlands). A first view of the overall structure of the data can be obtained from base R functions.
A more specific overview of the data can be obained by running the rtrim command count_summary. This function expects the presence of columns names count, year and site. If one or more of the actual data columns have different names, these can be specified.
In this case, we find that the Skylark dataset contains counts for 55 sites in 8 years (1984–1991). Of these 440 Site by Year combinations 202 were observed and the other 238 were missing. Two covariates are included: Habitat, which distinguishes between Dunes and Heathland sites and Deposition, which indicates the amount of acidic aerial deposition (This dataset was collected in the 1990’s when acidification was a prominent theme in ecological research).
Initial model estimation
To analyse these data with rtrim, we start with a model with time effects (model 3), ignoring the Habitat covariate. Model 3 is chosen because it makes no assumption about how population changes over time. Year effects are strictly independent of each other. A quick overview of the model results can be obtained by running summary() and plot(overall()).
Output from summary() includes:
- The call to
trim used to estimate the model.
- Details on the model, estimation method, and the numerical solution.
- The estimated model parameters (from both the additive and the multiplicative perspective), and the associated standard errors.
- Estimations of overdispersion and serial correlation parameters, if applicable.
- Several model goodness-of-fit measures.
The goodness-of-fit test (Likelihood Ratio) for this model amounts 194.8, with 140 degrees of freedom and p<0.05, which implies that the model has to be rejected.
Model simplification
The advantage of Model 3 is, as argued above, the absence of any assumptions regarding the temporal trend. This, however, comes at a price: Postive counts are required for all individual years to allow estimation of the model parameters. So, this model cannot be used for cases where one or more years are missing. Furthermore, the model is far from being parsimonious. Even if the Skylark population follows a perfectly theoretical trend with constant population increase or decrease, each year is assigned it’s own growth parameters, even if these are identical to last year’s.
For both reasons it may be preferable to replace model 3 by model 2 (piecewise linear), especially because in one extreme case these models are equivalent. This is the case when all years are treated as change points, and each year the trend changes into a different one. Let’s first check this.
z3 <- trim(count ~ site + year + habitat, data=skylark2, model=2, changepoints="all",
serialcor=TRUE, overdisp=TRUE)
summary(z3)
#> Call:
#> tools::buildVignettes(dir = ".", tangle = TRUE)
#>
#> Model : 2
#> Method : GEE (Convergence reached after 11 iterations)
#>
#> Coefficients:
#> covar cat from upto add se_add mul se_mul
#> 1 baseline 0 1984 1985 -0.21652346 0.1991430 0.8053136 0.1603726
#> 2 baseline 0 1985 1986 -0.16161496 0.2207074 0.8507687 0.1877709
#> 3 baseline 0 1986 1987 -0.12006156 0.2194797 0.8868658 0.1946491
#> 4 baseline 0 1987 1988 -0.24095494 0.2260048 0.7858770 0.1776120
#> 5 baseline 0 1988 1989 0.21787599 0.2249234 1.2434329 0.2796772
#> 6 baseline 0 1989 1990 -0.11533261 0.2180151 0.8910697 0.1942667
#> 7 baseline 0 1990 1991 -0.08490346 0.2329769 0.9186010 0.2140128
#> 8 habitat 2 1984 1985 -0.14448463 0.2323599 0.8654682 0.2011001
#> 9 habitat 2 1985 1986 0.42163114 0.2480259 1.5244461 0.3781021
#> 10 habitat 2 1986 1987 0.10935908 0.2335696 1.1155628 0.2605616
#> 11 habitat 2 1987 1988 0.38822335 0.2389049 1.4743590 0.3522316
#> 12 habitat 2 1988 1989 -0.12980911 0.2366361 0.8782631 0.2078287
#> 13 habitat 2 1989 1990 0.21388409 0.2301313 1.2384791 0.2850127
#> 14 habitat 2 1990 1991 0.17197025 0.2458819 1.1876425 0.2920198
#>
#> Overdispersion : 1.1616
#> Serial Correlation : 0.2265
#>
#> Goodness of fit:
#> Chi-square = 154.50, df=133, p=0.0979
#> Likelihood Ratio = 159.64, df=133, p=0.0575
#> AIC (up to a constant) = -106.36
and indeed this results in a similar model fit (although parameter values are different, of course).
The graphical display of the time-totals suggests that after an initial decline in counts, Skylark population recovers with approximately the same rate. One could either just argue if this recovery starts in 1985 or in 1987, and to what extent the recovery rate is ‘constant’, or one can look at the model statistics tfor a more objective analysis. In this case, we look at the Wald statistics associated with the Habitat covariate, and the individual changepoints:
The first test shows that there is a significant (at the 5% level) effect of the Habitat covariate on slopes (or year indices), showing that the slopes (year indices) for Dunes are different from those for Heathland. The tests for the significance of changes in slopes show that the only significant changes are for the years 1984, which means that the slope between 1984 and 1985 is different from zero, and 1985, which means that the slope between 1985 and 1986 is different from the slope between 1984 and 1985. This suggests that it should be possible to describe these data with a model with less than the full set of seven changepoints. To investigate this possibility, the stepwise procedure for selection of changepoints can be used by including stepwise=TRUE in the call to trim():
z4 <- trim(count ~ site + year + habitat, data=skylark2, model=2, changepoints="all",
stepwise=TRUE, serialcor=TRUE, overdisp=TRUE)
summary(z4)
#> Call:
#> tools::buildVignettes(dir = ".", tangle = TRUE)
#>
#> Model : 2
#> Method : GEE (Convergence reached after 11 iterations)
#>
#> Coefficients:
#> covar cat from upto add se_add mul se_mul
#> 1 baseline 0 1984 1985 -0.26905602 0.18234931 0.7641004 0.13933319
#> 2 baseline 0 1985 1991 -0.07758256 0.04105089 0.9253506 0.03798646
#> 3 habitat 2 1984 1985 -0.02040236 0.20676977 0.9798044 0.20259392
#> 4 habitat 2 1985 1991 0.17488185 0.04372458 1.1911055 0.05208058
#>
#> Overdispersion : 1.1265
#> Serial Correlation : 0.2283
#>
#> Goodness of fit:
#> Chi-square = 161.09, df=143, p=0.1431
#> Likelihood Ratio = 160.76, df=143, p=0.1471
#> AIC (up to a constant) = -125.24
wald(z4)
#> Wald test for significance of covariates
#> Covariate W df p
#> habitat 18.50555 2 9.584531e-05
#>
#> Wald test for significance of changes in slope
#> Changepoint Wald_test df p
#> 1984 10.99440 2 0.0040982293
#> 1985 14.64858 2 0.0006593285
Not surprisingly, this results in a model with only two changepoints left, at 1984 and 1985.
The difference between the models of this run (z4) and the previous (z3) can be tested by comparing their Likelihood Ratio’s, see also Section 2.5.
Since \(p \gg 0.05\) the \(H_0\) hypothesis that model z4 is a submodel of z3 in the sense that z4 can be obtained from z3 by setting some of z4 parameters to 0, cannot be rejected at the \(alpha=0.05\) level. In other words, both models are practically equivalent. The model z4, however, is the most sparse model, as shown by Akaike’s Information Criterion.
Concerning model z4, the Wald-test for the significance of the effects of the covariate on the slope parameters shows that this effect is very significant (p=0.0001) and the Wald-tests for the significance of changes in slope shows that both changes (at 1984 and 1985) are, as expected, also very significant (p is 0.004 and 0.0007, respectively).
The slope (in the additive parameterization) for a site is the sum of the constant term and the effects for the covariate values for that site. The effect for the first category of a covariate is zero and omitted from the output of summary() and coefs(). Thus, sites with covariate value 1 (Dunes) have slope -0.269 between 1984 and 1985 and -0.078 from 1985 onwards. The corresponding multiplicative parameters show that for Dunes there is a sharp decrease (the multipicative coefficient of 0.76 corresponds to \((1-0.76)\times100=24\%\) decrease) between 1984 and 1985 and a much smaller annual decrease (0.93, equivalent to 7% decrease) from 1985 to 1991. For sites with covariate value 2 (Heathland) the slope between 1984 and 1985 is \(-0.269 - 0.020 = -0.289\) corresponding with a multiplicative effect of \(0.764 \times 0.980 = 0.75\) which is only slightly different from the effect for Dunes for this time period. Apparently, the significant effect of the covariate is mainly determined by the trend from 1985 onwards. The parameters show indeed that Skylark populations increase in Heathland, while they decrease in Dunes. The slope is 0.097 (additive) and 1.10 multiplicative, corresponding to an annual increase of 10%.
Indices and time-totals
Model-based and imputed overall indices (based on the time-totals for all sites) can be obtained from TRIM output by calling index() or totals(). By default, only imputed indices or time-totals are returned. Model-based indices or time-totals can be added by using the which="both" option.
Indices can also be plotted:
In this plot the solid red line connects the indices for the individual years. In this case, the first year, 1984, is chosen as base year. Standard errors for the indices are shown using a transparent band and white ‘error-bars’.
In the last trim run, habitat was used as a covariate. Indices for covariates can also be computed, by setting the covars flag:
Indices are collected in a single dataframe, but can be easily separated by using e.g. subset().
Again, indices for the covariate categories can be plotted without much effort:
The model based indices reflect the strong decrease from 1984 to 1985 and the smaller decrease from 1985 onwards for Dunes (habitat category 1) and the similar decrease from 1984 to 1985 and the increase from that year onwards for Heathland (category 2). The overall model based indices are between the indices for Dunes and Heathland and show much less change over time than when Dunes and Heathland are treated separately. The imputed indices are very similar to the model based indices with the exception that the imputed index for 1986 is larger than the model based index for that year.
Weighting
So far, the overall indices are the indices that correspond with the time totals summed over all sites. The next run shows the results if the sites in Dunes are weighted 10 times.
The separate indices for Dunes and Heathland remain similar, of course, but due to the weighting the overall index decreases from 1985 onwards.