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.

```
rm(list=ls())
library(rtrim)
#> Welcome to rtrim version 2.3.0; Type ?`rtrim-package` to get started.
#>
#> Attaching package: 'rtrim'
#> The following object is masked from 'package:stats':
#>
#> heatmap
```

```
data(skylark2) # use extended version of the Skylark dataset
summary(skylark2)
#> site year count habitat deposition
#> site_01: 8 Min. :1984 Min. : 1.00 dunes:232 Min. :1.000
#> site_02: 8 1st Qu.:1986 1st Qu.: 2.00 heath:208 1st Qu.:2.000
#> site_03: 8 Median :1988 Median : 5.00 Median :3.000
#> site_04: 8 Mean :1988 Mean : 12.55 Mean :2.818
#> site_05: 8 3rd Qu.:1989 3rd Qu.: 11.00 3rd Qu.:3.000
#> site_06: 8 Max. :1991 Max. :131.00 Max. :4.000
#> (Other):392 NA's :238
#> weight
#> Min. : 1.000
#> 1st Qu.: 1.000
#> Median :10.000
#> Mean : 5.745
#> 3rd Qu.:10.000
#> Max. :10.000
#>
```

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.

```
idx <- which(names(skylark2)=="year") # rename year->season
names(skylark2)[idx] <- "season"
count_summary(skylark2, year_col="season") # show that it works
#> Total number of sites 55
#> Sites without positive counts (0):
#> Number of observed zero counts 0
#> Number of observed positive counts 202
#> Total number of observed counts 202
#> Number of missing counts 238
#> Total number of counts 440
```

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).

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())`

.

```
z1 <- trim(count ~ site + year, data=skylark2, model=3, serialcor=TRUE, overdisp=TRUE)
summary(z1)
#> Call:
#> tools::buildVignettes(dir = ".", tangle = TRUE)
#>
#> Model : 3
#> Method : GEE (Convergence reached after 11 iterations)
#>
#> Coefficients:
#> time add se_add mul se_mul
#> 1 1 0.00000000 0.0000000 1.0000000 0.00000000
#> 2 2 -0.32017834 0.1054679 0.7260195 0.07657174
#> 3 3 -0.16867813 0.1054130 0.8447808 0.08905090
#> 4 4 -0.18965668 0.1083411 0.8272431 0.08962440
#> 5 5 -0.08241227 0.1070220 0.9208922 0.09855569
#> 6 6 0.02080919 0.1058650 1.0210272 0.10809105
#> 7 7 0.09969167 0.1082296 1.1048302 0.11957536
#> 8 8 0.15580476 0.1107841 1.1685980 0.12946210
#>
#> Overdispersion : 1.3672
#> Serial Correlation : 0.3024
#>
#> Goodness of fit:
#> Chi-square = 191.40, df=140, p=0.0026
#> Likelihood Ratio = 194.80, df=140, p=0.0015
#> AIC (up to a constant) = -85.20
```

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.

A possible improvement of the model for a better fit might be the
inclusion of the *Habitat* covariate.

```
z2 <- trim(count ~ site + year + habitat, data=skylark2, model=3, serialcor=TRUE, overdisp=TRUE)
summary(z2)
#> Call:
#> tools::buildVignettes(dir = ".", tangle = TRUE)
#>
#> Model : 3
#> Method : GEE (Convergence reached after 12 iterations)
#>
#> Coefficients:
#> covar cat time add se_add mul se_mul
#> 1 baseline 0 1 0.0000000 0.0000000 1.0000000 0.0000000
#> 2 baseline 0 2 -0.2165232 0.1991433 0.8053138 0.1603728
#> 3 baseline 0 3 -0.3781379 0.2303445 0.6851360 0.1578173
#> 4 baseline 0 4 -0.4981991 0.2437072 0.6076240 0.1480823
#> 5 baseline 0 5 -0.7391541 0.2565409 0.4775177 0.1225028
#> 6 baseline 0 6 -0.5212781 0.2408642 0.5937612 0.1430158
#> 7 baseline 0 7 -0.6366106 0.2633641 0.5290827 0.1393414
#> 8 baseline 0 8 -0.7215137 0.2651985 0.4860160 0.1288907
#> 9 habitat 2 1 0.0000000 0.0000000 1.0000000 0.0000000
#> 10 habitat 2 2 -0.1444846 0.2323601 0.8654682 0.2011003
#> 11 habitat 2 3 0.2771464 0.2553396 1.3193596 0.3368847
#> 12 habitat 2 4 0.3865051 0.2681110 1.4718279 0.3946133
#> 13 habitat 2 5 0.7747286 0.2789551 2.1700031 0.6053334
#> 14 habitat 2 6 0.6449194 0.2642595 1.9058334 0.5036346
#> 15 habitat 2 7 0.8588036 0.2856443 2.3603350 0.6742164
#> 16 habitat 2 8 1.0307734 0.2882726 2.8032330 0.8080952
#>
#> 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
```

Now, the \(p\)-value of the likelihood ratio is (just slightly) above the classical threshold value of 0.05, and we decide to accept this model.

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:

```
wald(z3)
#> Wald test for significance of covariates
#> Covariate W df p
#> habitat 21.54981 7 0.003036166
#>
#> Wald test for significance of changes in slope
#> Changepoint Wald_test df p
#> 1984 10.27482731 2 0.005872859
#> 1985 9.17718172 2 0.010167175
#> 1986 3.08043508 2 0.214334471
#> 1987 1.53703040 2 0.463701061
#> 1988 1.63613306 2 0.441284040
#> 1989 0.88658542 2 0.641919282
#> 1990 0.01468516 2 0.992684312
```

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**.

```
gof(z3)
#> 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
```

```
LR3 <- gof(z3)$LR$LR # Get raw LR info for run 4
df3 <- gof(z3)$LR$df
gof(z4)
#> 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
```

```
LR4 <- gof(z4)$LR$LR # idem for run 3
df4 <- gof(z4)$LR$df
# Test the differece by using the fact that the difference of two LR measures is
# asymptotically Chi^2 distributed
LR <- abs(LR4 - LR3)
df <- abs(df4 - df3)
p <- 1 - pchisq(LR, df=df) # Use Chi-squared distribution
p
#> [1] 0.9997119
```

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%.

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.

```
index(z4, which="both")
#> time fitted se_fit imputed se_imp
#> 1 1984 1.0000000 0.00000000 1.0000000 0.00000000
#> 2 1985 0.7530579 0.06553658 0.7372610 0.06663569
#> 3 1986 0.7915817 0.06477129 0.8303956 0.07292113
#> 4 1987 0.8369110 0.06696343 0.8179245 0.07314308
#> 5 1988 0.8895273 0.07203043 0.8859143 0.07976879
#> 6 1989 0.9499771 0.07998826 0.9628009 0.08812600
#> 7 1990 1.0188773 0.09098601 1.0268504 0.09643874
#> 8 1991 1.0969220 0.10530053 1.1098103 0.10876964
```

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:

```
index(z4, which="both", covars=TRUE)
#> covariate category time fitted se_fit imputed se_imp
#> 1 Overall (none) 1984 1.0000000 0.00000000 1.0000000 0.00000000
#> 2 Overall (none) 1985 0.7530579 0.06553658 0.7372610 0.06663569
#> 3 Overall (none) 1986 0.7915817 0.06477129 0.8303956 0.07292113
#> 4 Overall (none) 1987 0.8369110 0.06696343 0.8179245 0.07314308
#> 5 Overall (none) 1988 0.8895273 0.07203043 0.8859143 0.07976879
#> 6 Overall (none) 1989 0.9499771 0.07998826 0.9628009 0.08812600
#> 7 Overall (none) 1990 1.0188773 0.09098601 1.0268504 0.09643874
#> 8 Overall (none) 1991 1.0969220 0.10530053 1.1098103 0.10876964
#> 9 habitat dunes 1984 1.0000000 0.00000000 1.0000000 0.00000000
#> 10 habitat dunes 1985 0.7641004 0.13933319 0.7791479 0.17017428
#> 11 habitat dunes 1986 0.7070608 0.12049989 0.7022059 0.17788473
#> 12 habitat dunes 1987 0.6542792 0.10988682 0.6448394 0.17652462
#> 13 habitat dunes 1988 0.6054377 0.10615237 0.5533376 0.17018093
#> 14 habitat dunes 1989 0.5602421 0.10724676 0.5814696 0.18256289
#> 15 habitat dunes 1990 0.5184204 0.11109150 0.5270963 0.17883059
#> 16 habitat dunes 1991 0.4797206 0.11609263 0.4820668 0.17468761
#> 17 habitat heath 1984 1.0000000 0.00000000 1.0000000 0.00000000
#> 18 habitat heath 1985 0.7486689 0.07298069 0.7203724 0.07707145
#> 19 habitat heath 1986 0.8251756 0.07643505 0.8820811 0.09258138
#> 20 habitat heath 1987 0.9095004 0.08191091 0.8877115 0.09474977
#> 21 habitat heath 1988 1.0024425 0.09019392 1.0200076 0.10685876
#> 22 habitat heath 1989 1.1048823 0.10206243 1.1165518 0.11809890
#> 23 habitat heath 1990 1.2177904 0.11821918 1.2283488 0.12999561
#> 24 habitat heath 1991 1.3422366 0.13928361 1.3629135 0.14869062
```

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.

One may try to extend the model further by also incorporating the second covariate ‘deposition’ in the model. This covariate is a measure for the amount of acidic aerial deposition. This time, the time-effects model (model 3) cannot be estimated due to lack of data in particular years:

```
check_observations(skylark2, model=3, covars=c("habitat","deposition"))
#> $errors
#> $errors$deposition
#> year deposition
#> 1 1990 1
#> 2 1991 1
#>
#>
#> $sufficient
#> [1] FALSE
```

The linear trend model with covariates can still be estimated.

```
z5 <- trim(count ~ site + year + habitat+deposition, data=skylark2, model=2,
serialcor=TRUE, overdisp=TRUE)
#> Warning: Serial correlation is very low (rho=0.162); consider disabling it.
```

```
summary(z5)
#> 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 1991 -0.0547486466 0.30356576 0.9467231 0.28739271
#> 2 habitat 2 1984 1991 0.1631245102 0.04752131 1.1771833 0.05594129
#> 3 deposition 2 1984 1991 -0.0399174384 0.30595598 0.9608688 0.29398355
#> 4 deposition 3 1984 1991 -0.0734776978 0.30615962 0.9291569 0.28447032
#> 5 deposition 4 1984 1991 -0.0002787989 0.30637534 0.9997212 0.30628993
#>
#> Overdispersion : 1.1557
#> Serial Correlation : 0.1620
#>
#> Goodness of fit:
#> Chi-square = 164.11, df=142, p=0.0987
#> Likelihood Ratio = 147.99, df=142, p=0.3482
#> AIC (up to a constant) = -136.01
```

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.

```
z6 <- trim(count ~ site + year + habitat, data=skylark2, model=2, changepoints="auto",
serialcor=TRUE, overdisp=TRUE, weights="weight")
idx = index(z6, "fitted", covars=TRUE)
plot(idx)
```

The separate indices for Dunes and Heathland remain similar, of course, but due to the weighting the overall index decreases from 1985 onwards.