В многоуровневой модели, каковы практические значения оценки параметров корреляции случайных эффектов, а не оценки?

В многоуровневой модели, каковы практические и связанные с интерпретацией последствия оценки, а не оценки параметров корреляции случайных эффектов? Практическая причина спрашивать это состоит в том, что в структуре Лмера в R нет реализованного метода для оценки p-значений с помощью методов MCMC, когда оценки делаются в модели корреляций между параметрами.

Например, взглянув на этот пример (приведенные ниже части), каковы практические последствия М2 по сравнению с М3. Очевидно, что в одном случае P5 не будет оцениваться, а в другом - будет.

Вопросов

1. По практическим соображениям (стремление получить значение p с помощью методов MCMC) может потребоваться подгонка модели без корреляции между случайными эффектами, даже если P5 по существу не равен нулю. Если это сделать, а затем оценить p-значения с помощью метода MCMC, интерпретируются ли результаты? (Я знаю, что @Ben Bolker ранее упоминал, что «комбинирование тестирования значимости с MCMC немного статистически непоследовательно, хотя я понимаю необходимость сделать это (получение доверительных интервалов более приемлемо)» , так что если это заставит вас спать лучше ночью притворяюсь, я сказал доверительные интервалы.)
2. Если не удается оценить P5, это то же самое, что утверждать, что оно равно 0?
3. Если P5 действительно ненулевой, то каким образом влияют оценочные значения P1-P4?
4. Если P5 действительно ненулевой, то каким образом влияют оценки ошибки для P1-P4?
5. Если P5 действительно не равен нулю, то каким образом интерпретации модели не могут содержать ошибки в P5?

Заимствуя из ответа @Mike Lawrence (те, кто более осведомлен, чем я, могут заменить его полной нотацией модели, я не совсем уверен, что смогу сделать это с разумной точностью):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(оценки P1 - P4)

М3: V1 ~ (1+V3|V2) + V3(оценки P1-P5)

Параметры, которые могут быть оценены:

P1 : глобальный перехват

P2 : Случайный эффект перехватывает для V2 (то есть для каждого уровня V2 отклонение перехвата этого уровня от глобального перехвата)

P3 : единая глобальная оценка эффекта (наклона) V3

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

P5 : корреляция между отклонениями перехвата и отклонениями V3 по уровням V2

Ответы, полученные из достаточно большого и широкого моделирования вместе с сопровождающим кодом на R с использованием lmer, будут приемлемыми.

Рассмотрим данные сна, включенные в lme4. Бейтс обсуждает это в своей онлайн- книге о lme4. В главе 3 он рассматривает две модели данных.

а также

В исследовании приняли участие 18 субъектов, изученных в течение 10 дней без сна. Время реакции рассчитывали в начале и в последующие дни. Существует четкий эффект между временем реакции и продолжительностью лишения сна. Есть также значительные различия между предметами. Модель А допускает возможность взаимодействия между случайным перехватом и эффектами наклона: представьте, скажем, что люди с плохим временем реакции более остро страдают от последствий лишения сна. Это предполагает положительную корреляцию в случайных эффектах.

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

Как видно из изображения, длительное время реакции связано с большей потерей производительности. Корреляция, использованная для моделирования, составила 0,58.

Я смоделировал 1000 образцов, используя метод имитации в Ime4, на основе установленных значений моих искусственных данных. Я подгоняю M0 и Ma к каждому и смотрю на результаты. Исходный набор данных содержал 180 наблюдений (по 10 для каждого из 18 субъектов), и моделируемые данные имеют такую ​​же структуру.

Суть в том, что разница очень мала.

1. Фиксированные параметры имеют одинаковые значения в обеих моделях.
2. Случайные эффекты немного отличаются. Для каждой моделируемой выборки имеется 18 случайных эффектов перехвата и 18 наклонов. Для каждой выборки эти эффекты вынуждены прибавлять к 0, что означает, что среднее различие между двумя моделями равно (искусственно) 0. Но различия и ковариации различаются. Медиана ковариации при МА составила 104, против 84 при М0 (фактическое значение 112). Дисперсии наклонов и пересечений были больше при M0, чем MA, по-видимому, чтобы получить дополнительное пространство для маневра, необходимое в отсутствие параметра свободной ковариации.
3. Метод ANOVA для lmer дает F-статистику для сравнения модели наклона с моделью только со случайным перехватом (без эффекта из-за лишения сна). Очевидно, что это значение было очень большим для обеих моделей, но обычно оно было (но не всегда) больше для МА (среднее значение 62 против среднего значения 55).
4. Ковариация и дисперсия фиксированных эффектов различны.
5. Примерно в половине случаев он знает, что МА правильно. Среднее значение p для сравнения M0 с MA составляет 0,0442. Несмотря на наличие значимой корреляции и 180 сбалансированных наблюдений, правильная модель будет выбрана только примерно в половине случаев.
6. Прогнозируемые значения отличаются для двух моделей, но очень незначительно. Средняя разница между прогнозами составляет 0, с SD 2,7. Sd самих предсказанных значений составляет 60,9

Так почему это происходит? @ Ганг обоснованно предположил, что отсутствие возможности корреляции приводит к некоррелированию случайных эффектов. Возможно, так и должно быть; но в этой реализации допускается корреляция случайных эффектов, а это означает, что данные могут выводить параметры в правильном направлении, независимо от модели. Неправильность неправильной модели проявляется в вероятности, поэтому вы можете (иногда) различать две модели на этом уровне. Модель смешанных эффектов в основном подгоняет линейные регрессии к каждому субъекту под влиянием того, что, по мнению модели, должно быть. Неправильная модель приводит к подгонке менее вероятных значений, чем при правильной модели. Но параметры, в конце концов, зависят от соответствия фактическим данным.

Вот мой несколько неуклюжий код. Идея заключалась в том, чтобы подогнать данные исследования сна, а затем создать имитированный набор данных с теми же параметрами, но с большей корреляцией для случайных эффектов. Этот набор данных был передан в simulate.lmer () для имитации 1000 образцов, каждый из которых подходил в обоих направлениях. После того, как я соединил подгоненные объекты, я мог извлечь различные особенности подгонки и сравнить их, используя t-тесты или что-то еще.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }

Плацидия уже предоставила подробный ответ, используя смоделированные данные, основанные на sleepstudyнаборе данных. Вот еще один (менее строгий) ответ, который также использует sleepstudyданные.

Мы видим, что можно повлиять на предполагаемую корреляцию между случайным пересечением и случайным наклоном, «сдвигая» случайную переменную предиктора. Посмотрите на результаты моделей fm1и fm2ниже:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

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

Де-соотнесение случайной перехватывать-наклонную ковариацию для LMMS обсуждаются в конспектах доктора Джека Вайса здесь . Вайс отмечает, что уменьшение дисперсионной корреляции таким образом иногда может помочь сходимости моделей, между прочим.

В приведенном выше примере изменяется случайная корреляция (параметр «P5»). Частично обращаясь к Q3 ОП, мы видим из вышеприведенного вывода, что:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected