Сравнение смешанной модели (субъект как случайный эффект) с простой линейной моделью (субъект как фиксированный эффект)

Я заканчиваю анализ большого набора данных. Я хотел бы взять линейную модель, использованную в первой части работы, и переоснастить ее, используя линейную смешанную модель (LME). LME будет очень похожим, за исключением того, что одна из переменных, используемых в модели, будет использоваться в качестве случайного эффекта. Эти данные получены из многих наблюдений (> 1000) в небольшой группе субъектов (~ 10), и я знаю, что моделирование эффекта субъекта лучше выполнять как случайный эффект (это переменная, которую я хочу сместить). Код R будет выглядеть так:

my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Все работает отлично, и результаты очень похожи. Было бы неплохо, если бы я мог использовать что-то вроде RLRsim или AIC / BIC, чтобы сравнить эти две модели и решить, какая из них наиболее подходящая. Мои коллеги не хотят сообщать о LME, потому что нет легкодоступного способа выбора «лучше», хотя я думаю, что LME - более подходящая модель. Какие-либо предложения?

Это нужно добавить к ответу @ ocram, потому что он слишком длинный, чтобы оставлять комментарии. Я бы отнесся A ~ B + Cк вашей нулевой модели, чтобы вы могли оценить статистическую значимость Dслучайного перехвата на уровне вложенных моделей. Как указал Окрам, условия регулярности нарушаются, когда , и тестовая статистика отношения правдоподобия (LRT) не обязательно будет асимптотически распределена χ 2 . Решение, которому меня учили, состояло в том, чтобы загрузить LRT (чье распределение при начальной загрузке, скорее всего, не будет ) параметрически и вычислить p-значение начальной загрузки следующим образом:$H_0: \sigma^2 = 0$$\chi^2$$\chi^2$

library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

Доля загруженных LRT более экстремальна, чем наблюдаемая LRT, является p-значением.

$0$$H_0:variance = 0$$H_1: variance > 0$...

РЕДАКТИРОВАТЬ

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