# Introduction

I am analysing temporal population data on the amphipod Orchestia gammarellus. At several moments each year, all animals were collected from a small spot, and several life history traits were measured. For these traits, I want to know whether there are seasonality and/or trends. I have data since mid 2014.

In this example, I focus on the proportion of reproductively active females in each sample. This peaks every year around July, so I would expect to find seasonality.

The data have been collected irregularly, so I thought time series analysis to be inappropriate. Rather, I have read good things about Generalised Additive Models (GAMs), and am trying to use them in R. No luck so far (see below).

# Methods

## Study system

Animals are collected from the grazed salt marsh on the island of Schiermonnikoog, a barrier island in the north of the Netherlands. This is the western part of the salt marsh, which is grazed by cattle. The cattle eat primarily Elytrigia atherica, but shun the tougher Juncus gerardii. This creates patches of "high vegegation" J. gerardii habitat, in which the soil is loose, amidst a sea of "low vegetation" E. atherica habitat, where the soil is compacted by trampling.

In comparison, the eastern part of the salt marsh is ungrazed, and there are no mosaics of vegetation types. Here the soil is much looser, and densities of amphipods much higher.

## Animals

The amphipod Orchestia gammarellus is a keystone species in European salt marshes. Through burrowing small tunnels, it changes soil nutrient levels, stimulates soil development, promotes soil oxygen conditions and impacts vegetation succession. However, on the compacted, grazed soil, this behaviour is limited to patches of "high vegetation".

Like isopods, amphipods have a brood pouch in which the juveniles, and are live-bearing.

## Study design

On a plot of about 50 m2 on the grazed salt marsh, animals are collected ideally once every four (in summer) or 6 (in winter). A cylinder (diameter 15 cm) is placed on a piece of cow dung in the low vegetation, which is a hot spot in low vegetation patches, and all animals are collected in 70% ethanol.

All animals are inspected under a microscope. We record gender (or juvenile, for which this is not possible), the number of antennal segments (an indicator of age, as one is added every moult), and length of the first body segment (pereon 1, as an indicator of body size). For females, we record the stage of the brood pouch (we can identify four levels, of which only the fourth is reproductively active), and, for fully reproductive females, brood size (if present) and the stage of the brood pouch (we can identify three levels).

# Statistical models

One of the variables I want to look into is the proportion of fully reproductive females in the samples, for which I want to know whether there are seasonality and trend. Here are the data:

brood_pouch <- data.frame(
date = c("2015-07-08","2015-08-20","2015-10-07","2015-12-10","2016-02-26",
"2016-04-29","2016-06-26","2016-07-12","2016-08-13","2016-09-29",
"2016-11-29","2017-03-03","2017-04-25","2017-06-17","2017-07-11",
"2017-08-13","2017-10-03","2017-11-29","2018-03-14","2018-04-29",
"2018-06-17","2018-12-07"),
n = c(101, 57, 75, 95, 118, 203, 197, 134, 77, 175, 314, 236, 64, 171, 103,
288, 49, 61, 133, 51, 75, 154),
proportion = c(0.900990099009901, 0.210526315789474, 0, 0, 0,
0.108374384236453, 0.903553299492386, 0.985074626865672,
0.246753246753247, 0.0171428571428571, 0, 0, 0.03125,
0.906432748538012, 0.300970873786408, 0.197916666666667,
0.0612244897959184, 0, 0, 0.176470588235294, 0.986666666666667,
0.00649350649350649)
)


As can be seen in the figure below, reproductively active females are abundant in summer but absent in winter.

Were it a generalised linear model, I would create separate variables of year and month, and create the following model:

glm(proportion ~ year * month, data = brood_pouch,
weight = n_total, family = binomial)


where weight is the number of individuals that make up each proportion.

With a generalised additive model, I would expect to create a similar model, where in this case there would be clear seasonality but no clear trend.

I have been looking at the gam package and the prophet package, but they seem to operate on single data points per date.

• There is a python package for irregular time series traces.readthedocs.io/en/latest (not tested) – kjetil b halvorsen Mar 5 '19 at 10:29
• Thanks for the suggestion. I have no experience with python and all my data are in R already, so I would prefer to use R. I realised I can also use the list of all females with their reproductive state (0 or 1) with family = binomial, in which case I can skip the entire weight factor. That seems to work. I still need to figure out what to model then, because date is highly significant, but I do not know how to decompose. – Raoul Van Oosten Mar 6 '19 at 9:06

One way to do this with a GAM is to set up a model that decompose time into different components. In this case, you want a seasonal component, potentially a trend even if the data don't seem to have one, and then some remainder, which we will assume for now is independent noise such that given the model we have random draws from a binomial distribution with expectation given by the fitted values of the model.

Here is one way to achieve this using mgcv.

First you need to convert the date variable into something useful; initially a Date classed object and then we process that to get information about the day of year (week of year would also work), the year, which we also cast as a factor for one of the models we'll look at.

brood_pouch <- transform(brood_pouch,
date = as.Date(as.character(date)))
brood_pouch <- transform(brood_pouch,
doy   = as.numeric(format(date, '%j')),
year  = as.numeric(format(date, '%Y')),
fyear = as.factor(format(date, '%Y')))


The simplest reasonable model might be

m1 <- gam(proportion ~ fyear + s(doy), data = brood_pouch, method = 'REML',
weights = n, family = binomial)


which has a seasonal smooth plus a different mean effect for each year. Extending this model to include separate seasonal smooths for each year, we have:

m2 <- gam(proportion ~ fyear + s(doy, by = fyear, k = 5), data = brood_pouch,
method = 'REML', weights = n, family = binomial)


And the full blown model with the long term trend modelled smoothly and the seasonal smooth allowed to vary with the trend would be given by

m3 <- gam(proportion ~ te(year, doy, k = c(3,6)), data = brood_pouch,
method = 'REML', weights = n, family = binomial)


Of these three, there is a clear preference for m2:

> AIC(m1, m2 ,m3)
df      AIC
m1 11.23595 169.6490
m2 14.99039 114.4300
m3 14.43998 130.7556


but there are some issues with that model that a proper analysis might try to tackle.

You can also model the trend as a s(Date), but that won't decompose the data into seasonal effects plus change between years, and you'll need to increase k to get a flexible enough spline to model all the data in a single shot.

• Thank you so much! I would be inclined to run anova.gam(m2), which, according to the manual, is comparable to drop1(). This shows a clear year effect and clear doy effect for all years. Am I right to conclude that (1) there is clear seasonality in the proportion of reproductively active females in each year, and (2) there is a clear difference in the proportion of reproductively active females between years? Does the clear year effect mean there is a trend, like in time series analysis, or merely that there are differences between years? – Raoul Van Oosten Mar 7 '19 at 16:15
• You need to be careful with comparisons using anova() when models are fitted using 'REML' because these models have different "fixed" effects. You should refit the models with method = 'ML' I think if you want to do this comparison. I think you can claim 1), but I don't think you can be so bold on 2). All that the model with separate smooths per year is telling you in the summary is that it is unlikely to have gotten the smooth you did under the null hypothesis that each smooth was a flat line at 0. – Gavin Simpson Mar 7 '19 at 16:35
• I'd want to look at differences between smooths (which can be done using the model, eg: a blog post of mine), if I you wanted to make judgements about how clear the differences were; in the shape of the seasonality, when in particular etc. – Gavin Simpson Mar 7 '19 at 16:38
• As for "trend", the statistical definition is just that there are differences in the mean, or level, over time. This model suggests that there are, but we've only been able to model that trend as a different mean per year. So there is a statistical trend, but that might not mean what your readers/reviewers think it means. – Gavin Simpson Mar 7 '19 at 16:38
• Clear. So the year effect says there is a clear difference in the proportion between any of the years (classic anova). Would converting it to an ordinal factor (2015<2016<2017<2018) and finding a clear difference not suggest a trend? If not, I could show [i.postimg.cc/Kc7BFgMD/screenshot-1552031847.png](this figure) to clarify there is a clear year effect, but no trend. – Raoul Van Oosten Mar 8 '19 at 14:31