GAMM with fragmented spatial data - presence/absence response

Question

Quite stuck so hope someone could help me in the right direction.

Some background for the data: For the last 15 years data have been collected on a specific species during the period from early summer to late autumn in three different geographical areas. This have cumulated to quite a large dataset. An overview of the number of individuals investigated per week in two of the areas:

table(data$year, data$week, data$area)

Area 1:

         26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44
  2002    0    0    0    0    0    0    0    0  931 1147 1142 1275 1238 1164 1295 1388 1270  707  138
  2003    0    0    0    0    0    0  595  378  499  782 1283 1319  977  820  874 1167  793  762  271
  2004    0    0    0    0    0    0    0    0    0    0  597  646  660  609  552  113   62    0    0
  2007    0    0   58   57   95  261  170  202  249  266  155  252  173  178   94   76   47   92   55
  2008    0    0    0    0    0    0  145  178  169  200  209  171  154  104  144   35   45    0    0
  2009    0    0    0    0    0    0   43   53   44   69   34   56   44   30   62   59    0    0    0
  2010    0    0   72   55  171  136  132  122  103   78  124   93   23   86    0    0    0    0    0
  2012   67    0   78   41   93   87   51  216  278  328  177  235  274  300  306  282  318  125   99
  2013    0   27   89  109   73   76  160  184  188  208  243  188  204  176  170  208  199  156  164
  2014    0    0    0    0    0    0    0    0   52   78  108  109  132  142  125  110  129   84   65
  2015    0    0    0    0   59   78  125  130  134  144  130  127  128  192   44   89   53   33   41

Area 2:

      26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44
  2002    0    0    0    0    0  216  417  480  549  668  758  811  625  827  603  104    0    0    0
  2003   88  101  681  320  362  110  512   80  504    0  331   54  289  237    0    0    0    0    0
  2004    0    0    0    0    0    0    0    0    0    0  320  352  435  420  392    0    0    0    0
  2007    0  137  167  189  277  227  222  241  262  264  157  205   45   81   44   75   70    0    0
  2008    0    0   45   50   81   97   88   82  113   75   79    0    0    0    0    0    0    0    0
  2009    0    0  131  144  221  168  189  190  184  222  153  204  150   56    0    0    0    0    0
  2010    0    0   71  139  123  197  155  149  133  129  158  190   49   61   53    0    0    0    0
  2012    0    0    0    0    0   59  109  202  140  122  157  107  197  192   45    0    0    0    0
  2013    0    0    0    0   82   27  151   91   77   85  122   99  140  108    0    0    0    0    0
  2014    0    0    0  104   90   95   91  113   53  155  154  125  176   68   65   60   56   56   82
  2015    0    0    0    0  171  187  173   57   74   77  170  174  173  173   67   79  200    0    0

As you can see my data have quite a lot of missing values both in the temporal and spatial scale. And also varying number of individuals have been investigated different years and in different areas (due to annual variation in staff, funding, weather etc.). In the table I have dropped some years where sampling were only conducted in one or two of the areas.

One of the parameters investigated for each (and every) individual is a presence/absence variable (parasite). I am hoping to fit a GAM with week as a smoothing term and area as a factor with the binary parasite as a response. I.e. I want to see if the timing of when there is a peak in parasites and if this varies between areas. (Experience tells us it is). A simple version of my model (without all the explanatory variables) is:

model<-gamm4(parasite~s(week, by=area)+area+s(size), family=binomial, data=data)

(I have included the parameter size just to emphasize that I also have explanatory variables which relate to the individual)

So my main question

How could I include year in my model? I am not really interested in the annual variation but mainly the overall trend of the response in one area and the difference between areas. This is why I have "pooled" the years. But I know should include year as some sort of random variable (factor), but not sure how..

I am also prepared to hear that my data is too fragmented in the temporal and/or spatial scale to use a GAM, but I am also very grateful for any input on the right statistical approach for my question.

Answer 1

You could add s(Year) + s(Year, by = area, m = 1) which will generate a common trend (s(Year)) and area-specific smooth departures from this common trend. The m = 1 in the area-specific term use a penalty on the first derivative of the spline and hence penalises departures from a flat, horizontal line, which when coupled with the common trend models how the Year effect in each area departs from the overall/polled Year effect.

If you wanted to treat Year as smooth but also truly random, you could do one of:

1. s(Year) + s(Year, area, bs = "fs", m = 1), or
2. s(Year, area, bs = "fs")

The first form is just the "true-random-effect" of the by version I showed above (common/population effect plus area-specific departures). The second form fits a random intercept and smooth for each area directly.

If you included Year as a standard random effect (random intercept for Year), this would be very close to fitting a spline in Year given that you're representing the splines as random effects in this model anyway.

(I'm presuming you want gamm4() because you have some complex random effects that need lme4? Other wise look at random effect splines in gam() via the bs = "re" basis and with method = "REML".)