[R-sig-ME] Including an autocorrelation term dramatically reduces random variance (lme)

[David Villegas Ríos]

Dear list,

My dataset includes 285 individuals for which I have measured one trait
over three months. The traits has been measured at two different time
scales: monthly (maximum 3 replicates per fish) and daily (maximum 90
replicates per fish). For one third of the fish, the trait was measured in
2012, for another third in 2013 and for the last third in 2014. But I
measured always in the same three months: june, july and august.

My objective is to use a mixed modelling approach to estimate repeatability
(Vrandom/(Vrandom+Vresidual)) of this trait. I want to account for three
fixed factors: area (two sites), total length (tl) and time. For
comparative purposes, I want to do it for the two temporal scales, i.e.,
one model using month and therefore three replicates per ID, and a
different model using (julian) day and 90 replicates per ID.

My attempts so far for the model with month:

MODEL 1: Mixed model without autocorrelation term

lme1=lme(loghr~area+tl+factor(month2),random=~1|fish,data=hrm3,method="REML")

> summary(lme1)
Linear mixed-effects model fit by REML
 Data: hrm3
       AIC      BIC    logLik
  1498.771 1531.537 -742.3854

Random effects:
 Formula: ~1 | fish
                (Intercept) Residual
StdDev:   0.4649101 0.4790193

Fixed effects: loghr ~ area + tl + factor(month2)
                            Value         Std.Error       DF    t-value
   p-value
(Intercept)           -1.1688029 0.15189536 514  -7.694790    0.0000
areatvedestrand -0.9943160 0.06490892 283  -15.318635  0.0000
tl                         -0.0034115 0.00308169 283  -1.107031    0.2692
factor(month2)7 -0.2242556 0.04074810 514  -5.503462     0.0000
factor(month2)8 -0.1269165 0.04245662 514  -2.989322     0.0029

Correlation:
                            (Intr)   artvds   tl         fc(2)7
areatvedestrand -0.224
tl                         -0.943  0.019
factor(month2)7 -0.121  0.003 -0.011
factor(month2)8 -0.113 -0.006 -0.011  0.475

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-2.7792032 -0.5127146 -0.1031226  0.3935509  4.3119462

Number of Observations: 802
Number of Groups: 286

The estimation of Repeatability from this model would be
0.46/(0.46+0.47)=0.485
However if I extract the residuals from this model
res=resid(lme,type="normalized") and do acf(residuals) the plot shows a
high autocorrelation in the first 15 lags. So I tried the following:

MODEL 2: same model including an autocorrelation term

lme3=update(lme1,correlation=corAR1(form=~month2))

> summary(lme)
Linear mixed-effects model fit by REML
 Data: hrm3
       AIC     BIC    logLik
  1461.013 1498.46 -722.5066

Random effects:
 Formula: ~1 | fish
              (Intercept)       Residual
StdDev: 0.0005272246 0.6819065

Correlation Structure: ARMA(1,0)
 Formula: ~month2 | fish
 Parameter estimate(s):
     Phi1
0.6097222
Fixed effects: loghr ~ area + tl + factor(month2)
                            Value         Std.Error      DF    t-value
  p-value
(Intercept)           -1.1584548 0.15740075 514  -7.359906  0.0000
areatvedestrand -1.0019165 0.06730508 283  -14.886194 0.0000
tl                         -0.0035497 0.00319449 283  -1.111192   0.2674
factor(month2)7 -0.2212453 0.03628563 514  -6.097327   0.0000
factor(month2)8 -0.1275377 0.04735955 514  -2.692968   0.0073
 Correlation:
                            (Intr) artvds tl     fc(2)7
areatvedestrand -0.227
tl                         -0.944  0.022
factor(month2)7 -0.104  0.003 -0.009
factor(month2)8 -0.124 -0.005 -0.013  0.611

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.62204450 -0.68663130 -0.08503995  0.58177543  3.79290975

Number of Observations: 802
Number of Groups: 286

If I plot acf(residual) now the plot looks much nicer (almost all the
autocorrelation has been corrected)

As you can see, the fact of including the autocorrelation term dramatically
reduces the random variance from 0.46 to 0.0005, which in turn reduces the
Repeatability from 0.485 to 5.98e-07.

NOW: if repeat this exercise working at a daily scale (90 replicates per
ID), the reduction in random variance and repeatability is much smaller
(from 0.46 to 0.16), but very noticeable still.

My questions are:
- How should I interpret this? Is the autocorrelation term in model 2
taking most of the random variance previously explained by the individual
ID in model 1?
- Model 2 yields a very low value of repeatability, however model selection
is telling me than the random effect (individual ID) has to be included in
the model. I interpret this as the repeatability being significantly
different from zero, even it is really low. Correct?
- With a different trait, I observe the same reduction after including a
"weights" argument. And with another trait, I observe a similar reduction
when I pass from an autocorrelation structure corAR1 to a correlation
structure corARMA (2,2). Any thought?
- Could I say that including explanatory variables (fixed part of the
model) will reduce the residual variance, and including terms in the random
part of the model (random factors, autocorrelations terms, weights, ...)
will reduce the random variance?

Many thanks in advance,

David

	[[alternative HTML version deleted]]