# How to model time-series data with nested random and fixed effects?

I'm analysing PAM fluorescence data from an experimental set-up that I duplicated from an earlier experiment with a missing control. That's why I haven't given the statistics of the experiment much (if any) thought in advance.

The set-up consisted of 8 containers with peat moss (Sphagnum magellanicum), divided over 4 treatments, so that each treatment was performed in duplicate. At regular (weekly) intervals, over the course of 3 months, I performed life PAM fluorescence measurements on a number of capitula (growth tops) in each container to determine a kinetic response curve for each of these capitula.

To minimize intraleaf (in my case, intracapitula) variance, ideally, PAM fluorescence measurements would have been repeated for the same leaf every week in the 3-month time series, but for practical reasons, my AOIs (areas of interest) for the fluorescence meter where located on different capitula every week. This is also my first subquestion: can I consider measurements at different time points in the same container as repeated measures, or would this only be valid if I had been measuring the same AOIs every week? And does this depend on whether I aggregate the measured values of the different AOIs per container before further analysis?

After nightfall, once every week, for 5–7 AOIs in each container, I determined a kinetic curve, for which the PAM software performs 20 measurements. The first measurement represents the dark-adapted fluorescence values, after which an actinic light source (at a wavelength that can facilitate photosynthesis) is started for the 19 remaining measurements. From the start of the kinetic curve (the dark adapted $\phi_{PSII}$ values), I determine $F_v/F_m$ and from the end of the curve (the flat part), I determine $\text{mean}(\phi_{PSII})$. $\phi_{PSII}$ and $F_v/F_m$ measure the quantum yield of photosystem II and the max. efficiency of photosystem II, respectively; $F_v/F_m = \phi_{PSII}$ in a dark-adapted state.

I'm interested in building two models, one in which the response (dependent) variable is $\phi_{PSII}$ and one in which it is $F_v/F_m$. The (independent) predictor variables are:

• AOI (factor): a number between 1–6;
• Container (factor): a number between 1–8;
• Treatment: (factor): a number between 1–4; and
• DaysTreated (integer): the number of days since the treatments began.

My guess is that I should treat AOI and Container as random effects variables, with AOI nested in Container and Container nested in the fixed effect variable Treatment. DaysTreated, then, would be my continuous predictor (covariate). For $\phi_{PSII}$, I would model this in R like this:

library(nlme)
YII_m1 <- lme(mean_YII ~ DaysTreated * Treatment,
random = ~1 | Container / AOI,
method = "ML",
data = fluor_aoi)
# fluor_aoi is a data-frame in which each AOI kinetic curve is
# aggregated into one row, where mean_YII = mean( YII[15:19] )
# and FvFm = YII[1]


I'm not sure if this is the most parsimious model. To find out, I want to try different models with different fixed effects but all with the same random effects. anova.lme() warned me that comparing between these models is a no-go when using the default method (method = "REML"), which is why I use method = "ML".

anova(YII_m1, # ~ DaysTreated * Treatment
YII_m2, # ~ DaysTreated:Treatment + Treatment
YII_m3, # ~ DaysTreated:Treatment + DaysTreated
YII_m4, # ~ DaysTreated:Treatment
YII_m5, # ~ DaysTreated + Treatment
YII_m6, # ~ DaysTreated
YII_m7  # ~ Treatment
)

Model df       AIC       BIC   logLik   Test  L.Ratio p-value
YII_m1     1 11 -2390.337 -2340.578 1206.168
YII_m2     2 11 -2390.337 -2340.578 1206.168
YII_m3     3  8 -2390.347 -2354.158 1203.173 2 vs 3  5.99019  0.1121
YII_m4     4  8 -2390.347 -2354.158 1203.173
YII_m5     5  8 -2366.481 -2330.293 1191.241
YII_m6     6  5 -2363.842 -2341.224 1186.921 5 vs 6  8.63915  0.0345
YII_m7     7  7 -2264.868 -2233.203 1139.434 6 vs 7 94.97389  <.0001


I would have liked it if the best fit was model 2 with the fixed effects formula ~ DaysTreated:Treatment + Treatment, because my expectation at the onset of my experiment was to see a decline in Sphagnum vitality, but only for some of the treatments and hopefully not in the controls. (The acclimatization period was very long, hoping that any effects on the mosses of the new (greenhouse) environment would have flattened out by the onset of the treatments.)

Edit 2015-May-1: First I compared only 6 models; model 4 was missing from my initial question. Also, I forgot to factorize treatment, so that instead of model 2, now, different models give the ‘best fit’.

Anyway, so far (unless you tell me otherwise), I feel I can continue to use model 2, which also best fits the visual observation that 4 of the 8 containers where doing very badly at the end of the experiment while the other 4 seemed to do ok.

anova(YII_m2)
numDF denDF  F-value p-value
(Intercept)           1   620 526.9698  <.0001
Treatment             3     4   5.0769  0.0753
DaysTreated:Treatment 4   620  36.4539  <.0001


An ANCOVA test on model 2 reveals that only the interaction between DaysTreated and Treatment is significant, which makes sense to me, given that the containers started out in roughly the same condition after acclimatization. There was visible difference between containers in the same treatments, but that should have been taken care of by correcting for the random error effect.

Mean $\phi_{PSII}$ plotted per container shows that there have been differences between containers within the same treatment from the onset of the treatments:

Now that I've made an attempt at constructing and testing a somewhat decent model (which I'd love to receive criticism on), I'd like to perform a multiple pairwise comparison to find out which treatments diverge significantly from each other over time, but I have no idea what is the proper way to approach this.

Also, I want to try a linear correlation, but again, I'm clueless as to how. Is there an appropriate way to integrate this in my model or should I try to model a regression per treatment?

Please forgive the ignorance in my approach and my questions. I'm a BSc student whose statistical background mainly consists of a brief entry-level course, followed by a recipe-level R course. RTFM comments are definitely welcome, as long as they include a link to TFM.