Search code examples
rsasmixed-modelsrandom-effectsautocorrelation

How do I suppress random effect but retain AR(1) correlation structure in R's lme mixed model


I am trying to re-create some SAS mixed model analyses using the lme function of R. The following SAS code:

proc mixed data = df noclprint covtest;
   class patid visno;
   model va = cst va0 cst0 / solution;
   random intercept / subject = patid;
   repeated visno / type = sp(pow)(month) subject = patid;
run;

postulates that every PATID gets a random intercept, and, in addition, for the repeated observations of each PATID, there is a set of errors with an AR(1)-type correlation structure (but with unequal time intervals) that gets added to them.

This analysis can apparently be re-created in R like this:

fitBoth <- lme(fixed = va ~ CST + cst0 + va0, data = muggeo,
           random = ~1 | PATID,
           correlation  = corAR1(form = ~ month | PATID))

I say it can be re-created because the numerical estimates of the parameters appear to be identical,

In SAS, I can drop the random intercept statement from the above SAS code, and SAS believes that the control statements are sensible, and estimates a different model (one without the random intercept).

But it looks like this doesn't work in R. If I drop the random statement from the R code, I get the following error message:

Error in lme.formula ... incompatible formulas for groups in 'random'
and 'correlation'

Is there a way to re-create the second SAS analysis (i.e. without the random statement) in lme?

Many thanks


Solution

  • As @user20650 suggests, you need to use gls ("generalized least squares") rather than lme ("linear mixed effects") if you want to fit a model with heteroscedasticity and/or correlation but no random effects. Something like

    fitBoth <- gls(va ~ CST + cst0 + va0, data = muggeo,
               correlation  = corAR1(form = ~ month | PATID))
    

    (note that the first argument of gls() is called model, not fixed as in lme) should work (but I can't test since you haven't provided a reproducible example) ... this does work:

    data("sleepstudy", package = "lme4")
    library(nlme)
    gls(Reaction ~ Days, data = sleepstudy,
       correlation = corAR1(form = ~Days | Subject))