CODEWITHSUNDEEP
rsaslme4nlme
user2716568
user2716568

Reputation: 1946

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

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!

Upvotes: 0

Views: 753

Answers (1)

aosmith
aosmith

Reputation: 36076

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

Upvotes: 5

Related Questions