Why does a linear mixed model work in SAS and nlme but not lme4?

Question

My data consists of 20 subjects in a control group and 20 in an experimental group. The DV of interest is a change score of peak power measured on each participant. There is also a dummy variable xVarExp that includes a 1 for subjects in the experimental group only. I am interested in individual responses and the variance of these numbers is the statistic summarising this. I am also interested in the means of each group; Exptal and Control.

My data is structured as follows:

structure(list(Subject = structure(1:40, .Label = c("Alex", "Ariel", 
"Ashley", "Bernie", "Casey", "Chris", "Corey", "Courtney", "Devon", 
"Drew", "Dylan", "Frances", "Gene", "Jaimie", "Jean", "Jesse", 
"Jo", "Jody ", "Jordan", "Kelly", "Kerry", "Kim", "Kylie", "Lauren", 
"Lee", "Leslie", "Lindsay", "Morgan", "Pat", "Reilly", "Robin", 
"Sage", "Sam", "Sidney", "Terry", "Tristan", "Vic", "Wil", "Wynn", 
"Zane"), class = "factor"), Group = structure(c(1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L), .Label = c("Control", "Exptal"), class = "factor"), 
    xVarExp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1), DV = c(3.3, -0.8, 2.7, 2.8, 0.6, 5.2, 
    1, 3.4, 1.3, -2.4, 8.5, 3.5, -1.9, 4.3, 1.2, -1.9, -0.6, 
    1.3, -2.6, -1, -3.7, 1.9, 4.6, 2.9, 7.2, -1.7, 4.2, 3.9, 
    -3.2, 9.9, 2.7, -1.7, 7.9, 8.1, 3.8, 2.8, 4.6, 0.8, 2.5, 
    4.1)), .Names = c("Subject", "Group", "xVarExp", "DV"), row.names = c(NA, 
-40L), class = "data.frame")

The statistician is a SAS user and has used the code below to obtain sensible answers:

title "Analyzing change scores";
proc mixed data=import plots(only)=StudentPanel(conditional) alpha=0.1 nobound;
class Subject Group;
model DV=Group/residual outp=pred ;
random xVarExp/subject=Subject;
lsmeans Group/diff=control("Control") cl alpha=0.1;
run;

I am beginning to use R and lme4, whereby I believe the code is:

Model1 <- lmer(DV ~ Group + (1|Subject/xVarExp), 
             data = RawData)

However, I receive the following: Error: number of levels of each grouping factor must be < number of observations

I managed to get the modelling working, using the syntax below, in nlme which matches the output of SAS:

Model2 <- lme(DV ~ Group, 
            random = ~ 1|xVarExp/Subject, data = RawData)

My questions are: 1) Why does the model work in nlme but not lme4? and 2) How can I match the SAS syntax to get the model going in lme4?

Thank you!

Yes, I simply did dput(RawData) within R. Is there something missing? — user2716568, Jun 30 '16 at 02:41
So you don't have repeated measures on your subjects and xVarExp is literally just a numeric version of Group? — Dason, Jun 30 '16 at 02:42
Each subject had a pre and post measurement, the difference between these is the change score/ DV of interest. I can upload the pre and post data if required. There was only one measurement for pre and one measurement for post. — user2716568, Jun 30 '16 at 02:44
But by differencing the pre and post measures you only have a single measurement for each subject. I don't really see why you think you need a mixed model for your data. If you kept the pre/post measures separate then you could use a mixed model but by analyzing the change scores directly you already eliminated any correlation within subjects that might cause an issue. — Dason, Jun 30 '16 at 02:46
Think of it this way.... At the moment how is your data any different than randomizing 20 people to a treatment and 20 to a control and taking a single observation from each subject? Just from looking at the data - it isn't any different. If you kept it in pre/post format then you it could make sense to talk about a mixed model so you could model the variability between subjects (since you would have more than one measurement from each) but after differencing you lose ALL of the ability to distinguish any sort of subject effect and it gets sucked into the overall error term. — Dason, Jun 30 '16 at 02:55
You may want to run `Model2 <- lm(DV ~ Group, data = RawData)` and compare with your `lme` and SAS inputs. I like to think of random effects models as a generalization of fixed effects models. A random effect model with only one observation per group will reduce to the fixed effect version of the model. It seems `lme` and SAS perform this reduction implicitly. `lme4` requires the reduction to be explicit. — Benjamin, Jun 30 '16 at 03:01
Thank you for the feedback. The statistician believes this is the correct SAS model/ linear mixed model, to analyse the change score. The xVarExp is included to estimate the extra error variance arising from individual responses. — user2716568, Jun 30 '16 at 03:42
How could they believe that the xVarExp variable is adding something that the group variable isn't already covering? Without having more than one measure on each subject (which you don't have in your data anymore) it's impossible to estimate the variance arising from the individual responses versus the pure error variance. I sincerely question your statistician's abilities... — Dason, Jun 30 '16 at 03:54
The statistician works on estimating an SD to represent individual differences and responses for athletic populations, as per: http://jap.physiology.org/content/jap/early/2015/02/18/japplphysiol.00098.2015.full.pdf — user2716568, Jun 30 '16 at 05:39
You can bypass this behavior in `lmer` via some of the `lmerControl` arguments. For example, you can change the checks for numbers of observations vs number of levels to `"ignore"` instead of stop. Such as `control = lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE="ignore")`. If you are trying to allow for different residual variances per group it may be that you will then want your random effects to then be something like `(xVarExp-1|Subject)`. — aosmith, Jun 30 '16 at 20:52
@aosmith Thank you so very much, that worked perfectly! The statistician will be happy. If you post your comment as an answer, I can accept it and up-vote. Cheers! — user2716568, Jul 01 '16 at 04:30

score 4 · Accepted Answer · answered Jul 01 '16 at 14:46

The lme4 package has some built-in model checks that are leading to errors. If you need to fit an unusual linear mixed model with lmer, you can change ignore the model checks that error by default via arguments in lmerControl.

To allow for a random effect that has the same number of levels as the residual error term like in the model you are fitting, you would need to change check.nobs.vs.nlev and check.nobs.vs.nRE to "ignore" from the default "stop". So a model where you want a different residual variance per group might look something like

Model1 <- lmer(DV ~ Group + (xVarExp-1|Subject), 
            data = RawData, control = lmerControl(check.nobs.vs.nlev = "ignore",
                                        check.nobs.vs.nRE="ignore"))

However, if the model you want is one that allows for a different residual variance per group then you might consider using gls from nlme. In gls you can easily extend the linear model to allow for nonconstant variance. That model would look like

Model2 <- gls(DV ~ Group, data = RawData, weights = varIdent(form = ~1|Group))

These two models give the same estimates and standard errors for the fixed effects:

summary(Model1)$coefficients
            Estimate Std. Error  t value
(Intercept)    1.395  0.6419737 2.172986
GroupExptal    1.685  1.0449396 1.612533

summary(Model2)$tTable
            Value Std.Error  t-value    p-value
(Intercept) 1.395 0.6419737 2.172986 0.03608378
GroupExptal 1.685 1.0449396 1.612533 0.11512145

Why does a linear mixed model work in SAS and nlme but not lme4?

1 Answers1