Как воспроизвести случайные эффекты в lme4 от SAS?

Я хочу запустить линейную смешанную модель для зависимой переменной DV, которая собирается под двумя разными Condition в трех разных Timepoint. Данные структурированы следующим образом:

   ## dput(head(RawData,5))
    structure(list(Participant = structure(c(2L, 2L, 2L, 2L, 4L), 
  .Label = c("Jessie", "James", "Gus", "Hudson", "Flossy", 
 "Bobby", "Thomas", "Alfie", "Charles", "Will", "Mat", "Paul", "Tim", 
  "John", "Toby", "Blair"), class = "factor"), 
 xVarCondition = c(1, 1, 0, 0, 1), 
 Measure = structure(c(1L, 2L, 3L, 4L, 1L), 
.Label = c("1", "2", "3", "4", "5", "6", "7", "8", 
"9", "10", "11", "12"), class = "factor"), 
Sample = structure(c(1L, 2L, 1L, 2L, 1L), 
.Label = c("1", "2"), class = "factor"), 
 Condition = structure(c(2L, 2L, 1L, 1L, 2L),
 .Label = c("AM", "PM"), class = "factor"),
 Timepoint = structure(c(2L, 2L, 2L, 2L, 1L), 
.Label = c("Baseline", "Mid", "Post"), class = "factor"),
 DV = c(83.6381348645853, 86.9813802115179, 69.2691666620429, 
 71.3949807856125, 87.8931998204771)), 
.Names = c("Participant", "xVarCondition", "Measure", 
   "Sample", "Condition", "Timepoint", "DV"), 
 row.names = c(NA, 5L), class = "data.frame")

Каждый Participant выполняет два испытания на Condition в течение трех Timepoint, как показано Measure; однако отсутствуют данные, поэтому не обязательно 12 уровней на участника. Столбец xVarCondition — это просто фиктивная переменная, которая включает 1 для каждой записи AM в Condition. Столбец Sample относится к 2 испытаниям для каждого Condition в каждом Timepoint.

Я пользователь R, но статистик — пользователь SAS, который считает, что код для модели должен быть таким:

proc mixed data=RawData covtest cl alpha=α
class Participant Condition Timepoint Measure Sample;
model &dep=Condition Timepoint/s ddfm=sat outp=pred residual noint;
random int xVarCondition xVarCondition*TimePoint*Sample 
          TimePoint/subject=Participant s;

Приведенный выше код SAS дает разумные ответы и работает отлично. Мы считаем, что результирующий синтаксис lme4 для приведенной выше модели будет следующим:

TestModel = lmer(DV ~ Condition + Timepoint + 
              (1 | Participant/Timepoint) +
              (0 + xVarCondition | Participant) +
              (1 | Participant:xVarCondition:Measure), data = RawData)

Однако при запуске этой модели я получаю следующую ошибку:

Error: number of levels of each grouping factor must be < number of observations

Правильно ли указаны случайные эффекты?


person user2716568    schedule 16.06.2016    source источник
comment
Я заметил, что оба ваших фиксированных эффекта, Condition и Timepoint, являются факторами. Вы уверены, что смешанная линейная модель — лучший подход в этом случае? Кроме того, я не понимаю разницы между xVarCondition и Condition.   -  person Adam Quek    schedule 16.06.2016
comment
Я считаю, что линейная смешанная модель уместна, поскольку нас интересует вариация внутри и между участниками. xVarCondition — это просто фиктивная переменная со значением 1 каждый раз, когда участник выполняет условие AM.   -  person user2716568    schedule 16.06.2016


Ответы (1)


Я не могу точно сказать из вашего описания, но, скорее всего, ваш термин Participant:xVarCondition:Measure создает группирующую переменную, которая имеет не более одного наблюдения на каждом уровне классификации, что сделает термин (1|Participant:xVarCondition:Measure) избыточным с термином остаточной ошибки, который всегда включается в модели lmer. Вы можете переопределить ошибку, если действительно хотите, включив

control=lmerControl(check.nobs.vs.nlev = "ignore")

в вашем вызове функции, но (если я правильно диагностировал проблему) это приведет к тому, что остаточная дисперсия и дисперсия Participant:xVarCondition:Measure будут совместно неидентифицируемыми. Такая неидентифицируемость обычно не вызывает проблем с остальной частью модели, но мне удобнее идентифицируемая модель (всегда есть вероятность, что такая неидентифицируемость приведет к числовым проблемам).

Аналогичный пример есть здесь.

Мое предположение можно проверить следующим образом:

ifac <- with(RawData,
   interaction(Participant,xVarCondition,Measure,drop=TRUE))
length(levels(ifac)) == nrow(RawData)
person Ben Bolker    schedule 16.06.2016
comment
Однако спасибо за аналогичный пример, даже когда я пытаюсь переопределить ошибку, я все равно получаю сообщение об ошибке. Могу ли я указать случайные эффекты, чтобы они соответствовали модели в SAS? Он работает гладко в SAS с разумными ответами. - person user2716568; 16.06.2016
comment
вы можете написать по адресу [email protected] (сегодня у меня может не быть на это времени). Также будет полезен воспроизводимый пример. - person Ben Bolker; 16.06.2016