R: generalised additive model on proportional data

Question

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.

Answer 1

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.