Ошибка при попытке запустить lmer() в R

Итак, вот моя проблема. У меня есть набор данных в R, на котором мне нужно запустить модель смешанных эффектов. Вот код:

data <- read.csv("D:/blahblah.csv")
analysis.data <- lmer(intdiff ~ stress_limit * word_position * follows + (1|speaker), data)
summary(analysis.data)

Когда я пытаюсь запустить скрипт, он возвращает следующую ошибку:

 Error in mer_finalize(ans) : Downdated X'X is not positive definite, 15.

Я отследил ошибку до параметра «follows», потому что, когда я просто использую stress_limit и word_position, он работает нормально. Если это поможет, данные в "follows" - это только 3 строки: n или l, согласная, гласная. Я пытался заменить пробелы на _, но безуспешно. Есть ли что-то во внутренней работе функции lmer(), что препятствует использованию в этом случае «follows»? Любая помощь будет здорово!

Для получения дополнительной информации: intdiff содержит числовые значения, stress_limit — это строки (с ударением или без ударения), а позиция слова — это также строки (Word Medial или Word Initial).

РЕДАКТИРОВАТЬ: Вот пример данных, который воспроизводит ошибку:

structure(list(intdiff = c(11.45007951, 12.40144758, 13.47898367, 
6.279497762, 18.19461897, 16.15539707), word_position = structure(c(2L, 
2L, 2L, 1L, 1L, 1L), .Label = c("Word Initial", "Word Medial"
), class = "factor"), follows = structure(c(4L, 4L, 4L, 1L, 2L, 
4L), .Label = c("Consonant", "n or l", "Pause", "Vowel"), class = "factor"), 
stress_limit = structure(c(2L, 1L, 1L, 2L, 2L, 2L), .Label = c("Stressed", 
"Unstressed"), class = "factor"), speaker = structure(c(2L, 
2L, 2L, 2L, 2L, 2L), .Label = c("f11r", "f13r", "f15a", "f16a", 
"m09a", "m10a", "m12r", "m14r"), class = "factor")), .Names = c("intdiff", 
"word_position", "follows", "stress_limit", "speaker"), row.names = c(NA, 
6L), class = "data.frame")

Я также попробовал функцию lme(), но она вернула эту ошибку:

Error in MEEM(object, conLin, control$niterEM) : 
Singularity in backsolve at level 0, block 1

Код в моем исходном посте — это точный код, который я использую, за исключением вызова библиотеки (lme4), поэтому я не упускаю никакой информации, которую могу придумать.

Моя версия R 2.15.2.


person Shakesbeery    schedule 01.02.2013    source источник
comment
Сколько строк содержит ваш фактический фрейм данных data?   -  person Sven Hohenstein    schedule 02.02.2013
comment
Фрейм данных содержит около 1110 строк. Однако данные предсказуемы, учитывая выборку.   -  person Shakesbeery    schedule 02.02.2013
comment
Содержат ли ваши переменные-предикторы все возможные комбинации предела ударения, положения слова и следов, или некоторые из них отсутствуют (потому что они неосуществимы или вам не довелось их измерить)? Верно ли with(data,all(table(stress_limit,word_position,follows)>=1))? (Это превращается в вопрос статистики, а не вопрос программирования...)   -  person Ben Bolker    schedule 03.02.2013
comment
Результат ложный, поэтому я предполагаю, что действительно существуют не все возможные комбинации. Означает ли это, что модель смешанных эффектов вообще нельзя запустить? Или я могу учесть это при запуске функции?   -  person Shakesbeery    schedule 03.02.2013


Ответы (1)


Трудно сказать наверняка без воспроизводимого примера: Как сделать отличный R воспроизводимый пример?

Но, предположим: такого рода проблемы обычно возникают из-за коллинеарности в матрице дизайна. Центрирование вашего непрерывного предиктора (intdiff) может помочь. Вы также можете напрямую изучить матрицу дизайна

X <- model.matrix( ~ stress_limit * word_position * follows, data)

Коллинеарность между парами: cor(X). К сожалению, у меня нет первоначальных советов по обнаружению мульти-коллинеарности (т.е. не между парами, а между комбинациями >2 предикторов), хотя вы можете изучить инструменты для вычисления коэффициенты инфляции дисперсии (например, library("sos"); findFn("VIF")).

В качестве перекрестной проверки lme также должен иметь возможность обрабатывать вашу модель:

library(nlme)
lme(intdiff ~ stress_limit * word_position * follows, 
   random=~1|speaker, data=data)

Когда я запускаю ваши тестовые данные в разрабатываемой версии lme4 (доступно на github), я получаю Error in lmer(intdiff ~ stress_limit * word_position * follows + (1 | : rank of X = 5 < ncol(X) = 12. С другой стороны, с таким небольшим набором входных данных (6 наблюдений) невозможно подобрать 12 параметров. Немного сложнее сказать, где именно ваша проблема. Все ли 12 комбинаций ваших трех переменных действительно встречаются в ваших данных? Если какие-то отсутствуют, то вам нужно следовать советам, данным в справке разрабатываемой версии:

В отличие от некоторых более простых сред моделирования, таких как «lm» и «glm», которые автоматически обнаруживают идеально коллинеарные переменные-предикторы, «[gn]lmer» не может работать с матрицами проектирования менее полного ранга. Например, в случае моделей с взаимодействиями, которые имеют ненаблюдаемые комбинации уровней, пользователь должен определить новую переменную (например, создать «ab» в данных из результатов «droplevels(interaction(a,b)). )').

В частности, вы можете подогнать эту модель следующим образом:

data <- transform(data,
       allcomb=interaction(stress_limit,word_position,follow,drop=TRUE))
lme(intdiff ~ allcomb, random=~1|speaker, data=data)

Это даст вам односторонний ANOVA, рассматривающий уникальные комбинации уровней, которые фактически присутствуют в данных, как категории. Вам придется выяснить для себя, что они означают.

В качестве альтернативы можно уменьшить количество взаимодействий в модели до тех пор, пока вы не получите набор, в котором нет пропущенных комбинаций; если вам повезет, (stress_limit+word_position+follow)^2 (все двусторонние взаимодействия) будут работать, но вам, возможно, придется уменьшить модель еще больше (например, stress_limit + word_position*follow).

Другой способ проверить это — использовать lm() в предложенных вами моделях и убедиться, что в оценочных коэффициентах нет значений NA.

Главное, что вы потеряете в этих способах, это удобство/легкость интерпретации, потому что параметры недостающих комбинаций все равно не могли быть оценены по данным...

person Ben Bolker    schedule 02.02.2013
comment
Во-первых, большое спасибо за ответ! Я не очень понимаю аспект коллинеарности проблемы, возможно, вы могли бы объяснить? Я погуглил этот термин, но это не помогло многое прояснить (начинающий статистик). В противном случае я улучшил имеющуюся выше информацию, которая, возможно, поможет нам решить эту проблему? - person Shakesbeery; 02.02.2013
comment
@Shakesbeery, да, в этом сложность. Инструменты могут делать только то, что позволяют данные :-( . Напоминает мне о моем первом знакомстве с анализом наименьших квадратов. Я весело понял, что (сократив данные лаборатории физики бакалавриата), даже если базовая модель была y ~ x ^ 2, я может получить гораздо более крутую подгонку к полиному 8-го порядка. - person Carl Witthoft; 02.02.2013
comment
Фантастика, это был чрезвычайно полезный и информативный ответ! Думаю, я могу взять это отсюда. Еще раз примите мою самую искреннюю благодарность. - person Shakesbeery; 03.02.2013