Почему линейная смешанная модель работает в SAS и nlme, но не в lme4?

Question

Мои данные состоят из 20 предметов в контрольной группе и 20 в экспериментальной группе. Представляющий интерес DV представляет собой оценку изменения пиковой мощности, измеренную для каждого участника. Существует также фиктивная переменная xVarExp это включает 1 для субъектов только в экспериментальной группе. Я заинтересован в отдельных ответах, и дисперсия этих чисел - статистика, суммирующая это. Я также заинтересован в средствах каждой группы; Эксталь и Контроль.

Мои данные структурированы следующим образом:

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")

Статистик является пользователем SAS и использовал приведенный ниже код для получения разумных ответов:

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;

Я начинаю использовать R и lme4, поэтому я считаю, что код:

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

Тем не менее, я получаю следующее: Error: number of levels of each grouping factor must be < number of observations

Мне удалось заставить работать моделирование, используя приведенный ниже синтаксис в nlme, который соответствует выводу SAS:

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

Мои вопросы: 1) Почему модель работает в nlme, а не в lme4? и 2) Как мне соответствовать синтаксису SAS, чтобы модель работала в lme4?

Спасибо!

Answer 1

Пакет lme4 имеет некоторые встроенные проверки моделей, которые приводят к ошибкам. Если вам нужно приспособить необычную линейную смешанную модель с lmer, вы можете изменить игнорировать модель проверяет эту ошибку по умолчанию с помощью аргументов в lmerControl,

Чтобы учесть случайный эффект, который имеет то же количество уровней, что и термин остаточной ошибки, как в модели, которую вы подходите, вам нужно изменить check.nobs.vs.nlev а также check.nobs.vs.nRE в "ignore" по умолчанию "stop", Таким образом, модель, в которой требуется разная остаточная дисперсия для каждой группы, может выглядеть примерно так:

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

Однако, если вам нужна модель, которая допускает различную остаточную дисперсию для группы, то вы можете рассмотреть возможность использования gls от NLME. В gls Вы можете легко расширить линейную модель, чтобы учесть непостоянную дисперсию. Эта модель будет выглядеть

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

Эти две модели дают одинаковые оценки и стандартные ошибки для фиксированных эффектов:

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