0
$\begingroup$

So, I am trying to understand some odd results in one of my mixed-effects models.

I am fitting data from 50 individual units over 20 timepoints each. There is also a time varying covariate $C$ which represents one of three conditions (treatment coded). I'm using the nlme package in R.

The level one model is: $$y_{tj} = \beta_{0j} + \beta_{1j} Time_{tj} + \beta_2 C_{tj} + \beta_3 (Time \times C)_{tj} + \epsilon_{tj}$$

with a level two model of
$$\beta_{0j} = \pi_{00} + u_{0j}$$ $$\beta_{1j} = \pi_{10} + u_{1j}$$

Now, when I investigate the residuals using ACF() function I see a spike at lag one. Reading through Pinheiro and Bates (2000) Chapter 5, I've tried an AR(1) structure as well as an ARMA(1,1) structure for the residuals. Comparing the fit among these three the ARMA(1,1) clearly is the better fitting structure.

> anova(mod0, mod0.ar, mod0.arma)
          Model df      AIC      BIC    logLik   Test   L.Ratio p-value
mod0          1 12 5913.776 5973.023 -2944.888                         
mod0.ar       2 13 5702.923 5767.109 -2838.462 1 vs 2 212.85227  <.0001
mod0.arma     3 14 5689.686 5758.808 -2830.843 2 vs 3  15.23749   1e-04

However, the ACF() shows that the correlations increased with each successive step. Here are the ACF() results for the three models of interest: mod0 = 'no structure', mod0.ar = 'AR(1) structure', mod0.arma = 'ARMA(1,1)'

> cbind(ACF(mod0), ACF(mod0.ar), ACF(mod0.arma))
    lag       ACF   lag       ACF   lag       ACF
1    0  1.00000000   0  1.00000000   0  1.00000000
2    1  0.33602521   1  0.42870350   1  0.60794719
3    2  0.08305023   2  0.19884848   2  0.45676945
4    3 -0.09015774   3  0.03441100   3  0.34080482
5    4 -0.14064216   4 -0.02889772   4  0.28317209
6    5 -0.19556329   5 -0.09967618   5  0.20736343
7    6 -0.23900230   6 -0.17554920   6  0.12053133
8    7 -0.19520183   7 -0.16130212   7  0.09409834
9    8 -0.17774278   8 -0.16866118   8  0.05719358
10   9 -0.15304269   9 -0.18176084   9  0.01063171
11  10 -0.18854378  10 -0.21749268  10 -0.04218369
12  11 -0.17917976  11 -0.24312063  11 -0.11463608
13  12 -0.16505966  12 -0.25805307  12 -0.18200234
14  13 -0.03499787  13 -0.17128993  13 -0.17436168

Now, my understanding is that if I get the correct structure for the residuals then the ACF() should reflect white noise. This is clearly not what is happening here, have I unintentionally introduced structure among the residuals?

Any comments, links, suggestions are greatly appreciated.

UPDATE

Below are the results from running the ACF() using normalized residuals.

> cbind(ACF(mod0, resType='n'),ACF(mod0.ar, resType='n'),ACF(mod0.arma, resType='n'))
   lag         ACF lag         ACF lag          ACF
1    0  1.00000000   0  1.00000000   0  1.000000000
2    1  0.33602521   1  0.44946820   1  0.606817548
3    2  0.08305023   2  0.22057896   2  0.454756584
4    3 -0.09015774   3  0.03580304   3  0.331220675
5    4 -0.14064216   4 -0.06923082   4  0.231779900
6    5 -0.19556329   5 -0.12016149   5  0.167267715
7    6 -0.23900230   6 -0.17831475   6  0.093203247
8    7 -0.19520183   7 -0.17073701   7  0.067495864
9    8 -0.17774278   8 -0.20500130   8  0.008955958
10   9 -0.15304269   9 -0.24103752   9 -0.055438933
11  10 -0.18854378  10 -0.24494946  10 -0.086814640
12  11 -0.17917976  11 -0.23472983  11 -0.137555183
13  12 -0.16505966  12 -0.24058335  12 -0.186371822
14  13 -0.03499787  13 -0.15469957  13 -0.160794285
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

Your ACF() calls are not including the error correlation. Looking at ?ACF.lme we note argument resType. This is documented as follows

resType: an optional character string specifying the type of residuals
          to be used. If ‘"response"’, the "raw" residuals (observed -
          fitted) are used; else, if ‘"pearson"’, the standardized
          residuals (raw residuals divided by the corresponding
          standard errors) are used; else, if ‘"normalized"’, the
          normalized residuals (standardized residuals pre-multiplied
          by the inverse square-root factor of the estimated error
          correlation matrix) are used. Partial matching of arguments
          is used, so only the first character needs to be provided.
          Defaults to ‘"pearson"’.

The key bit is the last line. You want resType = "normalized" in each of your ACF() calls.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for the reply, though this bit does not really address the issue of the increase in the auto-correlations due to the different residual structures, but rather how they are scaled. I did however try it and found the same trend as in my prior post. $\endgroup$
    – JSteele
    Commented Aug 22, 2015 at 5:11
  • $\begingroup$ You were computing the wrong residuals. Anyway, if you still see strong autocorrelation after computing them correctly can you update the question with the real outputs? I find it hard to reconcile what your say in this comment and the anova() output. $\endgroup$
    – Gavin Simpson
    Commented Aug 22, 2015 at 5:17
  • $\begingroup$ You could be trading off modelling the fixed effects for modelling residual correlation, which may point to an identifiability issue with the autocorrelation and time-varying terms. But even that doesn't fit with the anova() output you show. $\endgroup$
    – Gavin Simpson
    Commented Aug 22, 2015 at 5:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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