У меня есть несбалансированный набор данных повторных измерений для анализа, и я прочитал, что большинство статистических пакетов обрабатывают это с помощью ANOVA (т.е. сумма квадратов типа III) неверно. Поэтому я хотел бы использовать модель смешанных эффектов для анализа этих данных. Я много читал о смешанных моделях R, но я все еще очень Rплохо знаком с моделями со смешанным эффектом и не очень уверен, что все делаю правильно. Обратите внимание, что я пока не могу полностью отделиться от «традиционных» методов, и мне все еще нужны и специальные тесты.$p$

Я хотел бы знать, имеет ли смысл следующий подход или я делаю что-то ужасно неправильное. Вот мой код:

# load packages
library(lme4)
library(languageR)
library(LMERConvenienceFunctions)
library(coda)
library(pbkrtest)

# import data
my.data <- read.csv("data.csv")

# create separate data frames for each DV & remove NAs
region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1))

# output summary of data
data.summary <- summary(region.data)

# fit model
# "time" is a factor with three levels ("t1", "t2", "t3")
region.lmer <- lmer(dv ~ time + (1|subject), data=region.data)

# check model assumptions
mcp.fnc(region.lmer)

# remove outliers (over 2.5 standard deviations)
rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5)
region.data <- rm.outliers$data
region.lmer <- update(region.lmer)

# re-check model assumptions
mcp.fnc(region.lmer)

# compare model to null model
region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data)
region.krtest <- KRmodcomp(region.lmer, region.lmer.null)

# output lmer summary
region.lmer.summary <- summary(region.lmer)

# run post hoc tests
t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data)
t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data)
t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

# Get mcmc mean and 50/95% HPD confidence intervals for graphs
# repeated three times and stored in a matrix (not shown here for brevity)
as.numeric(t1.pvals$fixed$MCMCmean)
as.numeric(t1.pvals$fixed$HPD95lower)
as.numeric(t1.pvals$fixed$HPD95upper)
HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
    HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)

У меня есть несколько конкретных вопросов:

1. Это правильный способ анализа моделей со смешанными эффектами? Если нет, то что я должен делать вместо этого.
2. Достаточно ли хороши выходные графики критики mcp.fnc для проверки предположений модели или я должен предпринять дополнительные шаги?
3. Я понимаю, что для того, чтобы смешанные модели были достоверными, данные должны учитывать предположения о нормальности и гомоскедастичности. Как определить, что является «приблизительно нормальным», а что нет, взглянув на графики критики, сгенерированные mcp.fnc? Мне просто нужно почувствовать это, или это предписанный способ делать вещи? Насколько надежны смешанные модели в отношении этих предположений?
4. Мне нужно оценить различия между тремя временными точками для ~ 20 характеристик (биомаркеров) субъектов в моем образце. Подходит и тестирует ли отдельные модели для каждой приемлемой, пока я сообщаю обо всех проведенных тестах (значимых или нет), или мне нужна какая-либо форма коррекции для множественных сравнений.

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

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

1. повторно использовать фактор времени для получения сравнений t1-t2, t2-t3 и t1-t3 с pvals.fnc (из пакета languageR)
2. сравните мою смешанную модель с нулевой моделью, используя приблизительный F-тест, основанный на подходе Кенварда-Роджера (с использованием пакета pbkrtest), а не тест отношения правдоподобия (потому что я читал, что Кенварда-Роджера лучше рассматривать прямо сейчас)
3. Используйте пакет LMERConvenienceFunctions, чтобы проверить предположения и удалить выбросы (потому что я читал, что смешанные модели очень чувствительны к выбросам)

Ваш вопрос (ы) немного "большой", поэтому я начну с некоторых общих комментариев и советов.

Некоторое фоновое чтение и полезные пакеты

Вам, вероятно, следует взглянуть на некоторые вводные руководства по использованию смешанных моделей, я бы порекомендовал начать с Baayen et al (2008) - автор Baayen languageR. Барр и др. (2013) обсуждают некоторые проблемы со структурой случайных эффектов, и Бен Болкер является одним из текущих разработчиков lme4. Baayen и др. Особенно хороши для ваших вопросов, потому что они проводят много времени, обсуждая использование тестов начальной загрузки / перестановки (материал позади mcp.fnc).

Литература

• Baayen, RH, D. Davidson и D. Bates (2008). Моделирование смешанных эффектов со скрещенными случайными эффектами для предметов и предметов . Журнал памяти и языка, 59 (390–412).
• Barr, Dale J, R. Levy, C. Scheepers und HJ Tily (2013). Структура случайных эффектов для проверки подтверждающей гипотезы: Сохраняйте ее максимальной . Журнал памяти и языка, 68: 255–278.
• Bolker, Benjamin M, ME Brooks, CJ Clark, SW Geange, JR Poulsen, MHH Stevens und J.-SS White (2009). Обобщенные линейные смешанные модели: практическое руководство по экологии и эволюции . Trends Ecol Evol, 24 (3): 127–35.

У Флориана Йегера также есть куча материала о практической стороне смешанных моделей в блоге его лаборатории .

Существует также ряд полезных пакетов R для визуализации и тестирования смешанных моделей, таких как lmerTestи effects. effectsПакет особенно приятно , потому что она позволяет построить линейные регрессии и доверительные интервалы происходит за кулисами.

Хорошая посадка и сравнение моделей

$p$lmerTest

anova()mer$\chi^2$$\chi^2$$p$-значения для прямого сравнения двух моделей. Недостатком этого является то, что не сразу понятно, насколько хорошо вы подходите.

$t$summary()$|t|>2$fixef()

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

Многократные сравнения

Я отвечу на ваш вопрос №4 чуть более прямо. Ответ в основном такой же, как и в случае с традиционной ANOVA, но, к сожалению, это место, где существует большая неопределенность для большинства исследователей. (Это то же самое, что и традиционная ANOVA, потому что и линейные модели смешанных эффектов, и ANOVA основаны на общей линейной модели, просто у смешанных моделей есть дополнительный термин для случайных эффектов.) Если вы предполагаете, что временные окна создают Разница и хотите сравнить влияние времени, вы должны включить время в вашей модели. Между прочим, это также даст вам удобный способ оценить, изменило ли время разницу: есть ли основной (фиксированный) эффект для времени? Недостатком этого маршрута является то, что ваша модель станет намного более сложной и единственной «супер» модели со временем как параметром, вероятно, потребуется больше времени для вычисления, чем тремя меньшими моделями без времени как параматера. Действительно, классический учебный пример для смешанных моделейsleepstudy который использует время в качестве параметра.

$t$foreachlme4$\chi^2$

Если ваши характеристики являются зависимой переменной, то вам все равно придется вычислять разные модели, а затем вы можете использовать AIC и BIC для сравнения результатов.