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

У меня есть фрейм данных с переменными subject, wd и group и переменная ответа value. Каждый испытуемый разделен на одну группу, и каждый будний день выполняет 7 измерений. Поскольку каждый субъект полностью вложен в группу, я хочу использовать модель вложенных случайных эффектов с subject и group, а также добавить третий случайный эффект для wd. В настоящее время я использую это для этого:

model = lmer(value ~ 1+ (1|wd) + (1|group) + (1|subject), 
             data = dframe, REML = 0)

Я нашел код, на котором основано это, на странице 40 этого руководства . Я использовал и REML = TRUE, и REML = 0. Однако когда я использую VarCorr(model)$variances, я получаю

Groups   Name        Std.Dev.  
subject  (Intercept) 94.9534363
wd       (Intercept) 42.5931401
group    (Intercept)  0.0015608
Residual              0.9589836

Эта групповая дисперсия конфликтует с кодом, который я использовал для генерации данных, который имеет групповые средние 36,9, 28,78 и -15,269. Когда я смотрю на «остатки» для предсказанных случайных эффектов (с использованием ranef) по сравнению с истинными случайными эффектами, я получаю остатки с очень высокой корреляцией с группой, в которой они находятся (если бы я моделировал residuals ~ group, значение R-квадрата было бы более 0,9 ).

Как правильно подогнать вложенную модель случайных эффектов в R? Я предпочитаю использовать lme4, но подойдет любой пакет.

Вот код, который я использовал для генерации данных:

library(dplyr)
generate_data <- function(n = 10, g = 3, seed = 1, mean.overall = 300,
                          sigma.g = 50, sigma.wd = 50, 
                          sigma.subject = 100, sigma. = 30) {
    set.seed(seed)
    means.wd = rnorm(7) * sigma.wd
    means.g = rnorm(g) * sigma.g
    means.subject = rnorm(n*g) * sigma.subject
    dframe = data.frame(subject = rep(1:(g*n), each = 7),
                        wd = rep(1:7, g*n), 
                        group = rep(1:g, each = (7*n)))
    dframe = mutate(dframe,
       value = mean.overall + means.wd[wd] +    
           means.subject[subject] + means.g[group] + rnorm(7*g*n),
       subject = factor(subject, levels = 1:(n*g)),
       wd = factor(wd), 
       group = factor(group, levels = 1:g))
    dframe$value = round(pmax(5,dframe$value))
    truefx = list(wd = means.wd, group = means.g, 
                  subject = means.subject)
    list(data = dframe, effects = truefx)
}

dframe = generate_data()$data

r lme4 modeling random-effects
person Max Candocia    schedule 15.05.2015    source источник
comment
У меня проблемы с чтением ваших данных, но я не уверен, в чем проблема. Вы бы проверили, работает ли это для вас?   -  person Aaron left Stack Overflow    schedule 15.05.2015
comment
Почему-то раньше он не работал в блоках кавычек, но теперь работает построчно с форматированием кода.   -  person Max Candocia    schedule 15.05.2015
comment
Можете ли вы показать код того, как вы сгенерировали эти данные? Знание истинной модели, которую вы смоделировали, скорее всего, прояснит, какую модель вы пытаетесь / должны соответствовать.   -  person aosmith    schedule 15.05.2015
comment
Я добавил функцию выше. Данные, которые я использую, используют параметры по умолчанию (включая начальное число) для данных.   -  person Max Candocia    schedule 15.05.2015
comment
Я взял на себя смелость заменить дамп dput кодом для его генерации, а также отредактировал вашу функцию так, чтобы ее можно было читать без прокрутки.   -  person Aaron left Stack Overflow    schedule 15.05.2015
comment
Я думаю, вы подходите к той модели, которой хотите. Ваши оценки будут в среднем приближаться к истине для многих смоделированных наборов данных (поэтому вам нужно удалить начальное значение из генерации данных). Этот пример хорошо показывает, насколько плохими могут быть оценки дисперсии, когда у нас очень мало уровней.   -  person aosmith    schedule 15.05.2015
comment
Да, @aosmith, хотя я не уверен, что с этими несколькими уровнями оценки дисперсии действительно будут в среднем приближаться к истине. Кроме того, если вас действительно интересуют оценки дисперсии, обычно считается, что REML обеспечивает более точные оценки дисперсии.   -  person Aaron left Stack Overflow    schedule 15.05.2015
comment
Меня в первую очередь интересуют случайные эффекты для каждого предмета. Я использовал дисперсию случайного эффекта, чтобы показать, что есть что-то проблематичное, когда оно очень мало.   -  person Max Candocia    schedule 15.05.2015
comment
Я согласен, @Aaron, очень вероятно, что в среднем он будет иметь низкое смещение и будет иметь особенно странное распределение. Вот что такого крутого в симуляциях: вы действительно можете увидеть, как возникают эти проблемы.   -  person aosmith    schedule 15.05.2015

Ответы (1)


arrow_upward
0
arrow_downward

Я предполагаю, что вам нужна группа как фиксированный эффект, так как есть только пара уровней. Также кажется ненужным объединение буднего дня внутри предмета, поскольку для каждой комбинации предмет / будний день есть только одно наблюдение. Если это так, все, что вам нужно, это

lmer(value ~ group + (1|subject), data = dframe)

Однако не ясно, действительно ли день недели вложен в тему; если бы для всех испытуемых были одни и те же будние дни, возможно, более подходящим был бы другой вариант. Однако это лучший вопрос для stats.stackexchange.com.

Если вы действительно хотите вложить их, что-то вроде этого может сработать.

lmer(value ~ (1|group/subject/wd), data = dframe)
person Aaron left Stack Overflow    schedule 15.05.2015
comment
В настоящее время я моделирую данные, и мне еще предстоит определить окончательные размеры / количество групп, поэтому я спрашиваю. - person Max Candocia; 15.05.2015
comment
Кроме того, я не собираюсь вкладывать переменную wd. Я, наверное, включу это как фиксированный эффект. Я пробовал запустить model = lmer(value ~ 1+(1|wd) + (1|group/subject),data = dframe), но все равно получаю те же проблемы. - person Max Candocia; 15.05.2015