specification of mixed effects model with two levels of repeated measures (in R)

Question

My colleagues and I conducted a study of the effects of an experimental translocation on the movement and activity patterns of common brushtail possums in New Zealand. This involved first capturing 12 individuals (6 males and 6 females), fitting them with GPS collars, and releasing each study animal within its home range. After seven days, we re-captured the animals, fitted them with new GPS collars, and then moved them to a common release site that was well outside of all of their home ranges. For both deployments the GPS collars recorded location data at 5-min intervals for an 11-h period during the night when possums are active. Our response variables were duration of nightly active periods, total distance moved per night, mean nightly speed and several other metrics that were descriptive of variation in movement behaviour. Our data are a bit messy (unbalanced) because we didn’t get observations from all animals each night of the two 7-d sampling periods (for various reasons), but basically look like this for each animal, with multiple response variables:

1. Up to seven nights of movement data prior to translocation;
2. Translocation event;
3. Up to seven nights of movement data after translocation.

Our major research questions were: 1. Do males and females differ in the responses to translocation? 2. How do responses to translocation (as measured by activity, movement, etc) change with respect to the time (day) since the translocation event? 3. The interaction between 1 and 2 above.

As far as I can tell, these data are best analysed with a mixed-effects model with REML, because of their repeated-measures nature (both before and after translocation), missing values, and combination of fixed (sex) and random factors. From what I’ve read (Zuur’s book), classical repeated-measures ANOVA is inappropriate for several reasons.

My question is thus: what is an appropriate mixed-effects model formulation for these data in R using the ‘lme4’ package? There is a before/after effect (with respect to translocation), a repeated-measures effect (seven sequential days of data for each GPS collar deployment), and again the fixed effect of sex. I am confused on how to nest the data properly to incorporate the two different levels of temporal autocorrelation (i.e., before/after translocation and then the time series of each GPS deployment).

Any help with what is the correct model code for analysis in R would be hugely appreciated!!

Answer 1

Our major research questions were: 1. Do males and females differ in the responses to translocation? 2. How do responses to translocation (as measured by activity, movement, etc) change with respect to the time (day) since the translocation event? 3. The interaction between 1 and 2 above.

I believe the basic structure here would be something like:

• fixed effects: ~ sex*day*I(day>0), where day is centred on the translocation day (i.e. the translocation occurs on day 0, pre-translocation days are negative, post-translocation days are positive). This allows for different slopes before and after as well as a discontinuity at the translocation time (different intercepts).
• random effects: ideally ~day*I(day>0)|individual to allow for varying slopes and intercepts in both periods among individuals. This may be slightly too ambitious (estimating a 3x3 covariance matrix from 12 individuals = 6 parameters), but it's not insane.
• you may want to use lme to more easily allow for autocorrelation via correlation=corCAR1(~day) (or something like that); this uses a continuous-time first-order autoregressive model. The continuous-time part in particular (in contrast to the more usual discrete-time corAR1) is useful when there are missing values in the data set or the values are otherwise unevenly spaced. You might be able to do this via the lme4ord package, but that's definitely at the 'alpha' (i.e. experimental) release stage.

You probably want to think about transforming the response. One useful trick is to find the residuals from your initial fit and then try MASS::boxcox(r~1). A log-transform is usually a good initial guess unless you have zeros in your data.

Some other more advanced stuff that might be interesting but is computationally harder and possibly beyond the scope of this data set:

• consider the possibility of differences between sexes in either within-individual variance (easy via weights=varIdent(~1|sex)) or among-individual variance (harder)
• 'stack' the data and fit the responses as single multivariate observations (e.g. see here rather than a variety of univariate observations
• consider GAMs for the temporal trends (e.g. the gamm4 package)

Answer 2

Following is a simple example of 6 individuals (3 males and 3 females) followed for 3 days after translocation:

> mydf
   id gender daynum value
1   1      M      1     2
2   2      M      1     3
3   3      M      1     5
4   4      F      1     6
5   5      F      1     4
6   6      F      1     2
7   1      M      2     3
8   2      M      2     5
9   3      M      2     6
10  4      F      2     8
11  5      F      2     9
12  6      F      2     5
13  1      M      3     3
14  2      M      3     2
15  3      M      3     1
16  4      F      3     3
17  5      F      3     2
18  6      F      3     4

It may be best to omit values before translocation since they form random effect which should be taken care of by including id as random effect:

> library(lme4)
> library(lmerTest)
> summary(lmer(value~gender*daynum+(1|id), data=mydf))


Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['merModLmerTest']
Formula: value ~ gender * daynum + (1 | id)
   Data: mydf

REML criterion at convergence: 69.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4940 -0.5824 -0.1393  0.4241  1.9245 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.000    0.000   
 Residual             4.813    2.194   
Number of obs: 18, groups:  id, 6

Fixed effects:
               Estimate Std. Error      df t value Pr(>|t|)   
(Intercept)      5.7778     1.9349 14.0000   2.986  0.00982 **
genderM         -1.1111     2.7364 14.0000  -0.406  0.69084   
daynum          -0.5000     0.8957 14.0000  -0.558  0.58550   
genderM:daynum  -0.1667     1.2667 14.0000  -0.132  0.89719   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gendrM daynum
genderM     -0.707              
daynum      -0.926  0.655       
gendrM:dynm  0.655 -0.926 -0.707
> 

Alternatively, all values before translocation can be averaged to form a baseline value (day number 0):

   id gender daynum value
1   1      M      0   5.0
2   2      M      0   6.0
3   3      M      0   4.0
4   4      F      0   5.3
5   5      F      0   4.0
6   6      F      0   5.0
7   1      M      1   2.0
8   2      M      1   3.0
9   3      M      1   5.0
10  4      F      1   6.0
11  5      F      1   4.0
12  6      F      1   2.0
13  1      M      2   3.0
14  2      M      2   5.0
15  3      M      2   6.0
16  4      F      2   8.0
17  5      F      2   9.0
18  6      F      2   5.0
19  1      M      3   3.0
20  2      M      3   2.0
21  3      M      3   1.0
22  4      F      3   3.0
23  5      F      3   2.0
24  6      F      3   4.0

And use day number after translocation as a factor to compare them with baseline average:

> summary(lmer(value~gender*factor(daynum)+(1|id), data=mydf))
Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['merModLmerTest']
Formula: value ~ gender * factor(daynum) + (1 | id)
   Data: mydf

REML criterion at convergence: 65.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6265 -0.6971  0.0000  0.6971  1.3942 

Random effects:
 Groups   Name        Variance             Std.Dev.     
 id       (Intercept) 0.000000000000002625 0.00000005123
 Residual             2.057916666666664618 1.43454406230
Number of obs: 24, groups:  id, 6

Fixed effects:
                        Estimate Std. Error      df t value  Pr(>|t|)    
(Intercept)               4.7667     0.8282 16.0000   5.755 0.0000295 ***
genderM                   0.2333     1.1713 16.0000   0.199    0.8446    
factor(daynum)1          -0.7667     1.1713 16.0000  -0.655    0.5221    
factor(daynum)2           2.5667     1.1713 16.0000   2.191    0.0436 *  
factor(daynum)3          -1.7667     1.1713 16.0000  -1.508    0.1510    
genderM:factor(daynum)1  -0.9000     1.6565 16.0000  -0.543    0.5944    
genderM:factor(daynum)2  -2.9000     1.6565 16.0000  -1.751    0.0991 .  
genderM:factor(daynum)3  -1.2333     1.6565 16.0000  -0.745    0.4673    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gendrM fct()1 fct()2 fct()3 gM:()1 gM:()2
genderM     -0.707                                          
fctr(dynm)1 -0.707  0.500                                   
fctr(dynm)2 -0.707  0.500  0.500                            
fctr(dynm)3 -0.707  0.500  0.500  0.500                     
gndrM:fc()1  0.500 -0.707 -0.707 -0.354 -0.354              
gndrM:fc()2  0.500 -0.707 -0.354 -0.707 -0.354  0.500       
gndrM:fc()3  0.500 -0.707 -0.354 -0.354 -0.707  0.500  0.500

It may not be correct to keep day number as numeric since the change from 0 to 1 and then to further days may not be linear. Change from day 1 after translocation to further days is more likely to be linear hence day number was kept as numeric if average baseline (day 0) is not included as in first method here.

Above is a practical answer since you asked "what is an appropriate mixed-effects model formulation for these data in R using the ‘lme4’ package?" In both methods, you get results regarding all your questions: role of gender, role of day number after translocation and interaction between these two.