прогноз для модели Лмера с вложенными случайными эффектами

Question

Я сейчас пытаюсь помочь коллеге и просто не могу найти решение. Поэтому я надеюсь, что кто-то еще сможет нам помочь.

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

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

m <- lmer(weight ~ design +(1|study/species), data=dataset)

Я пытался делать прогнозы для разных видов, но с помощью общего исследования c: я создал новый data.table new.dt, который содержит уникальные комбинации design-видов-комбинаций исходного набора данных, и добавил столбец для отчета.

new.dt <- unique(dataset[,.(design, species))
new.dt$study <- "xyz"

Затем я использовал функцию предсказания и разрешил новые уровни.

new.dt$p <- predict(m, newdata=new.dt, re.form= NULL, allow.new.levels=TRUE) 

Я не получаю ошибку, но я получаю одно и то же предсказание для каждого вида в дизайне.

Есть ли способ сохранить исходные уровни одной части вложенного случайного эффекта и сделать другую часть новым уровнем?

Заранее спасибо!

ОБНОВЛЕНИЕ - рабочий пример : Эта проблема не зависит от набора данных.

library(data.table)
library(lme4)

dt <- data.table(expand.grid(design=c("a", "b"), species=c("x", "y", "z"), report=c("1", "2", "3"), count=seq(1, 10, 1)))
dt$weight <- 0

dt[species=="x"]$weight <- rnorm(60, 70, 10)
dt[species=="y"]$weight <- rnorm(60, 80, 15)
dt[species=="z"]$weight <- rnorm(60, 90, 20)

dt[design=="a"]$weight <- dt[design=="a"]$weight- 0.1*dt[design=="a"]$weight

dt[report=="1"]$weight <- dt[report=="1"]$weight+0.15*dt[report=="1"]$weight
dt[report=="2"]$weight <- dt[report=="2"]$weight-0.15*dt[report=="1"]$weight

m <-lmer(weight~design+(1|report/species), data=dt)

dt.pred <- unique(dt[,c(1:2)])
dt.pred$report<- "xyz"
dt.pred$pred<-predict(m, newdata=dt.pred, re.form= NULL, allow.new.levels=TRUE) 

Answer 1

«Сходство» происходит от того факта, что вы устанавливаете re.form = NULL или эквивалентно re.form = ~ 0.

Линейная модель со смешанным эффектом моделирует Y|beta,b ~ intercept + X %*% beta + Z %*% b + e, и, задав re.form = NULL, вы задаете определение Z %*% b = 0 во время прогнозирования. Поскольку это случайная часть вашей модели (т.е. (1|report/species)), вы удаляете случайный эффект species и report.

В смешанных моделях вы бы назвали этот вид предсказания «безусловным предсказанием» (или предельным предсказанием) [в то время как на практике оно более псевдо-безусловное]. Часто он используется в моделях, в которых случайный эффект содержит individual. В этом случае, когда вы наблюдаете за новым человеком, у вас есть неизвестный случайный эффект, но в зависимости от вашего исследования вас может заинтересовать только «систематический c» или «фиксированный» эффект (т.е. ходил ли человек на работу до получения сбил машину? он на велосипеде?). Здесь имеет смысл взглянуть только на безусловный / маргинальный эффект.

Сказал иначе, установив re.form = NULL Вы говорите Z %*% b = 0. Поскольку вид является частью Z с весами b, вы не можете видеть специфичный для вида c эффект в своем прогнозе.

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

Ps. Пакет data.table имеет функцию, эквивалентную expand.grid, называемую CJ, которая для больших наборов будет немного быстрее и потребляет больше памяти.

Answer 2

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

Вот пример с вашими данными, больше примеров можно найти в этой виньетке .

library(data.table)
library(lme4)
#> Loading required package: Matrix

dt <- data.table(expand.grid(design=c("a", "b"), species=c("x", "y", "z"), report=c("1", "2", "3"), count=seq(1, 10, 1)))
dt$weight <- 0

dt[species=="x"]$weight <- rnorm(60, 70, 10)
dt[species=="y"]$weight <- rnorm(60, 80, 15)
dt[species=="z"]$weight <- rnorm(60, 90, 20)

dt[design=="a"]$weight <- dt[design=="a"]$weight- 0.1*dt[design=="a"]$weight

dt[report=="1"]$weight <- dt[report=="1"]$weight+0.15*dt[report=="1"]$weight
dt[report=="2"]$weight <- dt[report=="2"]$weight-0.15*dt[report=="1"]$weight

m <-lmer(weight~design+(1|report/species), data=dt)

library(ggeffects)
ggpredict(m, "design")
#> 
#> # Predicted values of weight
#> # x = design
#> 
#> x | Predicted |   SE |         95% CI
#> -------------------------------------
#> a |     72.64 | 6.57 | [59.77, 85.52]
#> b |     82.66 | 6.57 | [69.78, 95.54]
#> 
#> Adjusted for:
#> * species = 0 (population-level)
#> *  report = 0 (population-level)

ggpredict(m, c("design", "report"), type = "re")
#> 
#> # Predicted values of weight
#> # x = design
#> 
#> # report = 1
#> 
#> x | Predicted
#> -------------
#> a |     80.78
#> b |     90.80
#> 
#> # report = 2
#> 
#> x | Predicted
#> -------------
#> a |     64.91
#> b |     74.92
#> 
#> # report = 3
#> 
#> x | Predicted
#> -------------
#> a |     72.24
#> b |     82.26
#> 
#> Adjusted for:
#> * species = 0 (population-level)

plot(ggpredict(m, c("design", "report"), type = "re"))
#> Loading required namespace: ggplot2

^{Создано в 2020-02-07 пакетом Представления (v0.3.0)}