sasmixed-models

Trying to convert SAS code to R (mixed effects model)


I have the following code in SAS:

PROC MIXED DATA = mydata;
  CLASS time;
  MODEL y = trt time / SOLUTION;
  REPEATED time / SUBJECT = patient_id TYPE = UN;
RUN;

I would like to run this in R, but the output does not match. I tried this

library(lme4)
model <- lmer(y ~ trt + time + (1 | patient_id), data = mydata)

I was expecting a p-value for trt of less than 0.05 but in R it is more than 0.1. Also, none of the other output matches. I assume it is something related to TYPE = UN, but I don't know how to do that with lme4. Maybe lme4 can't do that, in which case how do I achieve similar results ?


Solution

  • Your model is not a single random-intercepts one, but considers the outcome as multivariate normal. The R (and thereby V) part of your variance-covariance has as many distinct rows and columns as you have timepoints, a single random intercept would apply to all timepoints equally.

    This is identical in interpretation (but not computation or variance parametrization!) to having a random intercept per timepoint. Alternatives that get you the same result to within convergence precision should be:

    1. changing the SAS repeated statement into random, making the model actually random-effects using G-side parametrization with a single residual variance profiled out. This is just to show how different parametrizations match in their interpretation, it's a more complicated/slower computation to get to the same result but will do "exactly" what the R alternatives do.
    2. Using (time | patient_id) as the random effect in R (assuming time is a factor). Edit: after trying this out lme4 might only want to do this if you add control = lmerControl(check.nobs.vs.nRE = 'ignore') and will still produce convergence warnings while returning the same estimates as SAS, so #3 (or nlme) might be the better option for MMRMs in R.
    3. Using R's mmrm package which was purpose-built to reproduce these kind of models, the random effect becomes us(time | patient_id) (again assuming both are factors). This will also allow you to match several degrees of freedom approximations from SAS (if you're interpreting P-values you are using an appropriate degrees of freedom approximation, right?).

    What the SAS random and repeated statements do exactly is explained in detail in the documentation, a key point to remember when translating to R is that the implementations there usually only do actual random (G-side) effects. I've never done so myself but recall reading that nlme may be more flexible in this regard than lme4, though for models like this it doesn't matter since you can get them to match either way.

    Here's a blog (not mine) that looks at implementing MMRMs across different software packages as well.