2
$\begingroup$

I have $Y$ measurements per several Subjects and I'm studying impact of factor on $Y$ measurements. I've fit a lognormal mixed model with a random interaction, but I'm finding autoregressive dependence on residuals. However, some exploration suggests different temporal dependence per subject.

Can lme() incorporate different autoregressive dependence by subject? That is, can lme() implement a different AR order per subject?

I could do this fairly easily by dividing the residuals by subject and using the function ar(). But I'm hoping there's a more cohesive alternative in R.

Improve this question
$\endgroup$
10
  • $\begingroup$ Question about programming are better suited for stack overflow. $\endgroup$
    – user25658
    Sep 16, 2013 at 22:06
  • 3
    $\begingroup$ @BabakP This isn't a question about programming; just because one uses R, which also happens to be a programming language, doesn't make such question programming ones. This is clearly about the use of statistical software. $\endgroup$
    – Gavin Simpson
    Sep 16, 2013 at 22:34
  • $\begingroup$ @Gavin perhaps programming was not the right word, but i still think the question is off topic since it is basically asking can this software do what I want. I would like to hear from a moderator on this. $\endgroup$
    – user25658
    Sep 16, 2013 at 22:42
  • 2
    $\begingroup$ @Gavin We may be miscommunicating, but this is not the place to find out in detail where my posts could be clarified. Your meta post is a good place to continue this conversation. Here I would just like to mention that I am not defending the initial comment and I have not voted to close this question. In the future, though, I hope you will be a little more restrained and think harder before publicly accusing people of lack of knowledge, of writing "confused messes," etc.: that unconstructive approach just makes you sound petty and mean. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 17:52
  • 1
    $\begingroup$ @whuber Not to prolong this comment thread, but my point was that there are differing opinions on what is or is not OT, yours was the only Mod opinion I found when I wrote the comment, my only point here was to say that I don't agree with your position entirely in the thread I linked to, and finally I didn't mean to put words in your mouth or imply you would have closed this Q&A, apologies for that. $\endgroup$
    – Gavin Simpson
    Sep 17, 2013 at 18:07

1 Answer 1

Reset to default
4
$\begingroup$

The specific answer is no, this is not how lme() works. corAR1() will estimate a single additional parameter $\phi$ which is the AR(1) coefficient applied within any specified nesting via argument form. In other words, the same AR(1) is then assumed for, in your case, each subject.

The same applies to corARMA(); a single estimate of of the AR and MA parameters is made, across all levels of any grouping factor. Hence if you specify an ARMA(2,1), the same ARMA(2,1) is assumed to operate within the levels of the grouping factor.

You could investigate this for yourself by fitting a model with an AR correlation structure for residuals nested within subject corAR1(form = time | subject) and note the change in model degrees of freedom compared with a model fitted without the corAR1().

Specifically, when you specify say corAR1(form = time | subject), you are asking lme() to estimate a single parameter $\phi$ which is then applied within the groups defined by subject and that there is zero correlation between groups (subjects).

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for your response. I was aware of the standard implementation of lme fitting on ARMA type model for within group temporal dependance by potentially assuming independent subjects (an adequate assumption in my case). There's several ways around it, but wanted to make sure there wasn't a cohesive way to construct this model $\endgroup$
    – Rob
    Sep 17, 2013 at 0:33
  • 1
    $\begingroup$ @Rob Nope, only option I know is to go Bayesian and cook the correlations yourself. Or you could look at adapting the current corStruct methods/functions to do what you want, but that would be quite involved especially if you aren't familiar with nlme; these are extensible by the user and can be quite flexible, e.g. the mgcv package. $\endgroup$
    – Gavin Simpson
    Sep 17, 2013 at 2:52
  • $\begingroup$ A Bayesian implementation would not be necessary. What I did is that I included a lagged version of the response as a covariate. Granted, the estimate of this coefficient will be biased, but as long as the other assumptions of the model hold, the estimator will be fairly consistent and efficient. Besided, the intention of the model was for interpretation of a main covariate and not prediction. The other model assumptions were fine, and the temporal dependence had disappeared in the ACF plot (in the previous model it displayed exponential decay). $\endgroup$
    – Rob
    Sep 18, 2013 at 14:57
  • $\begingroup$ Strangely, acf plots of the previous model (with no lag response) by subject suggested different AR orders per Subject. When fitting an lme using corAR1() the temporal association was not being modeled properly. But when I used a lagged version of the response, the ACF plot did not suggest temporal dependance. A different slope for the lagged response didn't dramatically improve on the original lagged response model. By accounting for the temporal dependance this way, I have effectively being able to focus on the main covariate of interest. $\endgroup$
    – Rob
    Sep 18, 2013 at 15:06
  • $\begingroup$ @Rob Good point about explicitly modelling using the lag1 covariate - I'd overlooked that. I was thinking more of using the Bayesian tools and something convenient to allow computation of different within group AR(1)s. Perhaps you should add your modelling lag1 as a Answer to your Question? $\endgroup$
    – Gavin Simpson
    Sep 18, 2013 at 15:43

Your Answer

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

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