I want to apply a GAMM with R to this time series but I am not sure how to handle the station P18, as shown in the figure below.

enter image description here

If I shrink the dataset to the point where P18 ends (i.e. left side of the dashed line) I would loose many reversing trends in P07,P14,P22,P2O and so on. The question is thus: can I use the whole dataset including P18? I am not sure what the GAMM is doing with such a gap. what would be if there were more gaps?

any idea is welcome!

  • $\begingroup$ How's that? ;-) $\endgroup$ – Gavin Simpson Jul 6 '11 at 15:48
  • $\begingroup$ To me, that looks like you have very different trends, so how are you using GAMM? What model do you have in mind? I would probably fit an AM with a different smooth trend for each site. Is that what you want? $\endgroup$ – Gavin Simpson Jul 6 '11 at 15:50
  • $\begingroup$ I have also a set of nine potentially influencing predictor variables like temperature, salinity, distance to land, etc that might influence the respones (activity of harbor porpoises per day in %). I thought GAMM would be most useful in order to find the influencing factors. so it is not too much about the trend in this case, but more about finding the effects that explain it. many thanks for putting in the figure! do you have another suggestion instead a GAMM? ah...the Generalized approach only because there is no nonlinear lmm otherwise (?). it looked partly similiar to Zuur et al. 2009 $\endgroup$ – Jens Jul 6 '11 at 15:57
  • $\begingroup$ There are a lot of examples in Zuur et al 2009. Care to be specific - I have a copy on my desk... I'm not clear what the Q is now. Is Julian day one of the covariates you want to use. Take a look at my Answer and see if that helps and perhaps follow-up there with comments? $\endgroup$ – Gavin Simpson Jul 6 '11 at 16:14

If the covariates (not including Julian day) are sufficient to model the response in and of themselves, then the problem of missing data in the one site is irrelevant. For example, if the model were

mod <- gam(porpoise ~ s(temperature) + s(salinity) + s(dist.to.land) +
                      s(site, bs = "re"), data = foo))

Which is an AM with a random effect term for site employed using a random effect s(foo, bs = "re") spline, then we are saying that, irrespective of site, there is a functional, additive relationship between the porpoise acitivity and the 3 covariates. The random effect for site just says that the mean activity per site is allowed to differ, i.e. some sites can have more activity than others in general. You can get the same with gamm() but that does involve a bit more heavy lifting - though if your covariates don't fully model the temporal dependence in the data then you'll need gamm() and a correlation structure for the residuals.

In that case, the lack of data after a certain date in one or more sites is irrelevant because you are not using time itself to model the response, it is the combination of the covariates that is being modelled.

By the way, the generalised in GAM and GAMM relates to the generalisation of the model to non-Gaussian errors. It is the additive bit, the A in the acronyms, that gives the non-linear bit. You are making things difficult for yourself by expressing activity as a % variable because those usually need to be modelled using things like beta regression which are not GLMs and thus are not GAMs. It might be better to work with actual events, counts in other words, and not a %.

  • $\begingroup$ I did not get the re(site) running, is the function gam from mgcv, gam, gamlss? I did not find it, but it looks interesting. However, there is quite some temporal autocorrelation present. but yes, I also have count data which are very overdispersed and I was not yet sure how to get the overdispersion corrected in the GAMM (family=quasipoisson?). Chapter 18 of Zuur comes closest. Temporal autocorrelation + gaps in measurement. so I still feel gamm() with a quasipoisson and arma(p,q) could help. but what about this gap? $\endgroup$ – Jens Jul 6 '11 at 16:33
  • $\begingroup$ re() is in newer versions of mgcv it allows you to add simple random effects to an additive model without having to use the heavy-lifting of a full mixed-model setup (lme() and/or glmmPQL() with gamm() or lmer() via the gamm4 package). The basic reasoning is that you can vie the models fitted by GAM in a linear mixed model way, the splines are split into fixed and random effects parts, so the re() smoother is just a way of doing it easily, for simple random effects within gam(). It won't be as efficient as lme/lmer etc because they are designed to fit MMs, however. $\endgroup$ – Gavin Simpson Jul 6 '11 at 16:48
  • $\begingroup$ aaaah, it worked with the following: s(site, bs="re") thanks! looks indeed nice, but still the serial correlation is there...so back to gamm. $\endgroup$ – Jens Jul 6 '11 at 17:14
  • $\begingroup$ Ooops. Sorry. My bad, yes, it should be bs = "re" $\endgroup$ – Gavin Simpson Jul 6 '11 at 17:48
  • 1
    $\begingroup$ OK, with gamm(), if my description of the problem is what you want to do - you can ignore the gap at that site. Make sure you nest the AR or ARMA correlation within site. And be prepared to wait around for a while for the model to fit... ;-) $\endgroup$ – Gavin Simpson Jul 6 '11 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.