How to fit different ARMA models by group in R? 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. 
 A: 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).
