Начальные значения смешанной модели для lme4

Я пытаюсь подогнать под смешанную модель с помощью функции lmer из пакета lme4. Однако я не понимаю, что следует вводить в параметр start. Моя цель - использовать простую линейную регрессию, чтобы использовать оцененные там коэффициенты в качестве начальных значений для смешанной модели.

Допустим, моя модель следующая:

linear_model = lm(y ~ x1 + x2 + x3, data = data)
coef = summary(linear_model)$coefficients[- 1, 1] #I remove the intercept
result = lmer(y ~ x1 + x2 + x3 | x1 + x2 + x3, data = data, start = coef)

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

Тогда я получаю ошибку следующего вида:

Error during wrapup: incorrect number of theta components (!=105) #105 is the value I get from the real regression I am trying to fit.

Я пробовал много разных решений, пытаясь предоставить список и назвать эти значения theta, как я видел на некоторых форумах.

Также код Github проверяет, подходит ли длина, но я не могу найти, к чему он относится:

# Assign the start value to theta
if (is.numeric(start)) {
        theta <- start
}

# Check the length of theta
length(theta)!=length(pred$theta)

Однако я не могу найти, где определено pred$theta, и поэтому не понимаю, откуда взялось это значение 105.

Любая помощь ?


person Yohan Obadia    schedule 08.09.2016    source источник
comment
Если вы этого не сделаете, он предоставит начальное значение. Просто оставь это.   -  person G. Grothendieck    schedule 08.09.2016
comment
Я хочу предоставить один, чтобы ускорить вычисления. Идея заключалась в том, чтобы запустить линейную регрессию и использовать полученные коэффициенты в качестве начальных значений.   -  person Yohan Obadia    schedule 08.09.2016
comment
Он уже выполняет внутренние вычисления, чтобы получить хорошие начальные значения, особенно если ваши единственные смешанные члены имеют форму 1 | x. Обратите внимание, что есть виньетка lme4 о производительности.   -  person G. Grothendieck    schedule 08.09.2016
comment
Это скорее форма y ~ x | Икс. Как он решает, какие начальные значения использовать?   -  person Yohan Obadia    schedule 08.09.2016
comment
Возможно, вам лучше попробовать поиграть с критериями остановки, а не с начальными значениями. Снова посмотрите виньетку.   -  person G. Grothendieck    schedule 13.09.2016


Ответы (1)


Несколько моментов:

  • lmer на самом деле явно не соответствует ни одному из коэффициентов фиксированного эффекта; они профилированы так, что они решаются неявно на каждом этапе процесса нелинейного оценивания. Оценка включает только нелинейный поиск параметров дисперсии-ковариации. Это подробно описано (довольно технически) в одной из виньеток lme4 (уравнения 30-31, стр. 15). Таким образом, предоставление начальных значений для коэффициентов фиксированного эффекта невозможно и бесполезно ...
  • glmer точно соответствует коэффициентам фиксированных эффектов как часть нелинейной оптимизации (как @Grothendieck обсуждает в комментариях), если nAGQ>0 ...
  • это, по общему признанию, довольно неясно, но начальные значения для theta параметров (единственные, которые явно оптимизированы в lmer fits) равны 0 для недиагональных элементов фактора Холецкого, 1 для диагональных элементов: это закодировано здесь
   ll$theta[] <- is.finite(ll$lower) # initial values of theta are 0 off-diagonal, 1 on

... где вам нужно знать, что в восходящем направлении значения вектора lower были закодированы так, что элементы вектора theta, соответствующие диагональным элементам, имеют нижнюю границу 0, недиагональные элементы имеют нижнюю границу -Inf ; это эквивалентно запуску с единичной матрицы для масштабированной ковариационной матрицы (т. е. матрицы дисперсии-ковариации параметров случайных эффектов, деленной на остаточную дисперсию) или дисперсии случайных эффектов -ковариационная матрица (сигма ^ 2 I).

Если у вас есть несколько случайных эффектов и большие матрицы дисперсии-ковариации для каждого из них, все может стать немного сложным. Если вы хотите восстановить начальные значения, которые lmer будут использовать по умолчанию, вы можете использовать lFormula() следующим образом:

library(lme4)
ff <- lFormula(Reaction~Days+(Days|Subject),sleepstudy)
(lwr <- ff$reTrms$lower)
## [1]    0 -Inf    0
ifelse(lwr==0,1,0)  ## starting values
## [1] 1 0 1

Для этой модели у нас есть единственная ковариационная матрица случайных эффектов 2x2. Параметры theta соответствуют фактору Холецкого нижнего треугольника этой матрицы в порядке столбцов, поэтому первый и третий элементы являются диагональными, а второй элемент - недиагональным.

  • Меня беспокоит тот факт, что у вас 105 theta параметров; подгонка такой большой модели со случайными эффектами будет чрезвычайно медленной и потребует огромного количества данных для надежной подгонки. (Если вы знаете, что ваша модель имеет смысл и у вас достаточно данных, вы можете изучить более быстрые варианты, например, использовать пакет MixedModels Дуга Бейтса для Джулии или возможно glmmTMB, который может масштабироваться лучше, чем lme4 для задач с большими theta векторами ...)
  • формула вашей модели, y ~ x1 + x2 + x3 | x1 + x2 + x3, кажется очень странной. Я не могу понять никакого контекста, в котором имело бы смысл иметь те же переменные, что и термины со случайным эффектом, и группирующие переменные в одной модели!
person Ben Bolker    schedule 09.09.2016
comment
Пусть y будет рейтингом, который человек i дает продукту, который в виде комбинации атрибутов x обозначается манекенами x1, x2, x3. Мы хотим зафиксировать коэффициенты, связанные с этими манекенами, принимая во внимание, что каждый человек i не имеет одинакового рейтингового шаблона, и у некоторых респондентов будет очень мало продуктов, которые они будут оценивать, что затрудняет использование только их ответа для определения своей собственной модели . - person Yohan Obadia; 09.09.2016
comment
оценка вариации между людьми в ответ на предикторы x1, x2 и x3 (а также ответы в среднем по населению) будет определяться формулой y~x1+x2+x3+ (x1+x2+x3|i) ... - person Ben Bolker; 09.09.2016
comment
Да, извините, я именно этим и занимаюсь. Спасибо за исправление. - person Yohan Obadia; 09.09.2016
comment
OK. По-прежнему будет очень сложно подогнать ковариационно-дисперсионную матрицу 14 * 14 ... сколько испытуемых, сколько всего наблюдений? - person Ben Bolker; 09.09.2016
comment
В примере, который я вам привел, у нас есть 735 наблюдений с 50 респондентами и 13 манекенами, которые используются для калибровки модели. Но считайте это просто примером, поскольку мы разрабатываем инструмент, который люди будут использовать, предоставляя собственные данные. Таким образом, у нас не будет никакого контроля над количеством респондентов, наблюдений или предикторов. - person Yohan Obadia; 10.09.2016
comment
Я могу почти гарантировать вам, что оценка ковариационной матрицы 14x14 на основе (приблизительно) этого количества респондентов и общего количества наблюдений не сработает. Вы могли бы рассмотреть возможность создания составно-симметричной или блочно-составной-симметричной модели в nlme ... - person Ben Bolker; 13.09.2016