Извлечь средства обработки из объекта lmer и рассчитать погрешности

Question

[Я подробно расскажу об эксперименте, который у меня есть для фона - Я понимаю метод lmer s, просто неясно, как извлечь некоторые значения, которые мне нужны/от руки, поэтому я разместил это на SO, а не на CV. Я надеюсь, что это было правильное место для публикации!]

data are here.

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

Есть 6 блоков. В каждом блоке есть 2 графика и два подзаголовка на каждом участке. Лечение 1 имеет два уровня (A и B) и применяется на уровне участка: в каждом блоке имеется один участок, получающий Лечение 1 уровня A, и один приемный уровень 1 лечения B.

Лечение 2 применяется на уровне подзаголовка а также имеет два уровня (C и D): на каждом участке есть один подзабор, получающий обработку. Два уровня A и один подзабор, получающий обработку 2 уровня B.

Эксперимент проводился в течение 3 лет. Мне интересно, как каждая комбинация двух методов лечения влияет на мою зависимую переменную (DV).

Как таковой у меня есть 4 комбинации лечения:

TMT1A:TMT2C 

TMT1B:TMT2C 

TMT1A:TMT2D 

TMT1b:TMT2D 

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

The lmer с для каждого года выглядит следующим образом:

m2011<- lmer (DV2011~ TMT1*TMT2 + (1|Block/TMT1)) 
m2012<- lmer (DV2012~ TMT1*TMT2 + (1|Block/TMT1)) 
m2013<- lmer (DV2013~ TMT1*TMT2 + (1|Block/TMT1)) 

Для графического представления изменения этих средств лечения в течение долгого времени, я хочу, чтобы извлечь средства лечения для каждого уровня каждого обращения (см четыре уровня выше) на каждый год и нанести их на каждый год эксперимента, сродни the example in this post

Мне интересно, можно ли извлечь средства для обработки для четырех различных комбинаций лечения (таких как перечисленные выше) от объекта lmer? Или они должны быть рассчитаны вручную?

Один из способов, который я решил сделать, - это создать другой фактор, который представляет собой 4 комбинации обработки (см. Столбец «TMT1x2» в вставленных данных). Затем я мог бы запускать следующую модель за каждый год:

m2011<- lmer (DV2011~ TMT1x2 + (1|Block/TMT1)) 

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

Кроме того, если мне нужно вычислить средства лечения вручную, кто-нибудь знает, как это можно сделать, учитывая уровни гнездования в моем эксперименте?

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

Если кто имеет представление о том, все это было бы весьма признателен!

Answer 1

Альтернатива, которая использует функции в пакете languageR. Я звоню в ваш набор данных df.

library(lme4) 
library(languageR) 
library(ggplot2) 

# fit model 
# n.b. I don't claim that this is a sensible model 
# It is just used to demonstrate the plot 
mod <- lmer(DV ~ TMT1 * TMT2 + (1|Block), data = df) 

# create MCMC matrix 
mcmc <- pvals.fnc(mod, nsim = 1000, withMCMC = TRUE) 
# pval.fnc also calculates MCMC-based p-values and HPD confidence intervals, 
# and plot the posterior distributions of the parameters 

# plot using plotLMER.fnc 
# in addition, set withList = TRUE to create a list of data frames with plot data 
# which can be used for a (possibly prettier) plot in ggplot 
ll <- plotLMER.fnc(mod, withList = TRUE, pred = "TMT1", 
       intr = list(
       "TMT2", 
       c("C", "D"), 
       "end", 
       list(c("red", "blue"), rep(1, 2))), 
       addlines = TRUE, 
       mcmcMat = mcmc$mcmc) 

# here follows additional steps to plot using ggplot 

# convert list to data frame 
df <- do.call(rbind, ll$TMT1) 

# rename 
names(df)[names(df) == "Levels"] <- "TMT1" 

# add TMT2 
df$TMT2 <- rep(c("C", "D"), each = 2) 

# plot using ggplot 
dodge <- position_dodge(width = 0.1) 
ggplot(data = df, aes(x = TMT1, y = Y, col = TMT2, group = TMT2)) + 
    geom_point(position = dodge, size = 3) + 
    geom_errorbar(aes(ymax = upper, ymin = lower, width = 0.1), position = dodge) + 
    geom_line(position = dodge) + 
    ylab("DV") + 
    theme_classic() 

Answer 2

Я думаю, что то, о чем вы просите, это некоторая форма predict(), для которой нет метода по умолчанию для класса mer в lme4 (по крайней мере, версия на CRAN). Однако вы можете использовать ez::ezPredict.

library(ez) 
library(ggplot2) 
to_predict <- expand.grid(TMT1=c("A","B"), TMT2=c("C","D")) 
t_means <- rbind(ezPredict(m2011, to_predict=to_predict, boot=F), ezPredict(m2012, to_predict=to_predict, boot=F), ezPredict(m2013, to_predict=to_predict, boot=F)) 
t_means$YEAR = rep(2011:2013, each = 4) 
ggplot(t_means, aes(x=YEAR, y=value, color=TMT1:TMT2)) + geom_point() + geom_line() 

Эта функция имеет некоторые дополнительные функции, которые могут оказаться полезными, например, обеспечить загрузочные значения.

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

mm = unique(model.matrix(m2011)) 
Y_bar <- c(mm%*%fixef(m2011), mm%*%fixef(m2012), mm%*%fixef(m2013)) 
ggplot(t_means, aes(x=YEAR, y=Y_bar, color=TMT1:TMT2)) + geom_point() + geom_line() 

Я не совсем конечно, что вы подразумеваете под «рассчитать средства обработки ... учет уровней гнездования в моем эксперименте». Случайные эффекты в смешанной модели структурированы, нормально распределенные отклонения от эффектов уровня населения (фиксированные эффекты). Может быть поучительно рассмотреть оценки случайного эффекта ranef(m2011) и соответствующую конструктивную матрицу m2011@Zt.

Итак, если все, что вы хотите построить, это средства лечения уровня населения, вы можете просто работать с оценками фиксированных эффектов fixef(m2011) и матрицей с фиксированным эффектом model.matrix(m2011), как указано выше. Если вы хотите включить некоторую меру неопределенности вокруг прогнозов на уровне населения или хотите предсказания для каждого блока/сюжета/подзаголовка, вам нужно будет работать как с случайным, так и с фиксированным эффектом. Предлагаю начать, посмотрев на http://glmm.wikidot.com/faq под заголовком «Прогнозы и/или доверительные интервалы (или предсказания) по прогнозам».

EDIT 8/26/2013:

Вы могли бы рассмотреть bootMer() в разрабатываемой версии lme4 для (параметрический бутстрапируемых) доверительных интервалов предсказания, которые должны включать неопределенность в дисперсиях случайных эффектов и будет работать с GLMM (см., Например, this thread).

Идея состоит в том, чтобы моделировать из интересующей модели, переоборудовать моделируемые значения и рассчитать статистику, представляющую интерес, из модели refit. Вы можете пройти самостоятельно шаги, с simulate() и refit():

t_sim <- apply(simulate(model, 999), 2, function(x) combn(unique(model.matrix(model))%*%fixef(refit(model, x)), 2, diff)) 

, который генерирует 999 бутстраповских повторений попарных различий между средствами лечения, которые вы могли бы использовать quantile() на (или любой вкус бутстраповского доверительного интервала вы хотите) :

apply(t_sim, 1, function(.) quantile(., c(0.975, 0.025))) 

Answer 3

С lme4 теперь я думаю bootMer(), вероятно, лучший способ пойти, потому что приходится неопределенности всех видов в модели. Однако для некоторых классов проблем bootMer() не может работать из-за того, сколько времени может потребоваться для каждой модели.Для этих больших проблем существует R-пакет, называемый merTools, который предоставляет метод predictInterval для использования arm::sim для учета неопределенности фиксированных и случайных эффектов, а также остаточной ошибки модели. Он довольно прост в использовании и намного быстрее для обеспечения прогнозов в тех случаях, когда модель занимает много времени. Он имеет хороший охват интервала предсказания, созданного bootMer(), для задач, где дисперсия в отношении между случайными эффектами достаточно хорошо определена.

Чтобы использовать его, вы бы просто пойти:

library(merTools) 
preds <- predictInterval(m2011, newdata = myData, level = 0.95, n.sims = 1000) 

Есть несколько настраиваемых параметров другого пользователя, но результат предсказания объекта, аналогичного тому, который получают при запросе интервалов прогнозирования от lm - данные в три столбца .frame с колонками fit, lwr и upr.