Парный t-критерий как частный случай линейного смешанного моделирования

20

Мы знаем, что парное t- тестирование - это всего лишь частный случай одностороннего повторного измерения (или внутри субъекта) ANOVA, а также линейной модели смешанного эффекта, которую можно продемонстрировать с помощью функции lme () пакета nlme в R как показано ниже.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Когда я запускаю следующий парный t-тест:

t.test(x1, x2, paired = TRUE)

Я получил этот результат (вы получите другой результат из-за генератора случайных чисел):

t = -2.3056, df = 9, p-value = 0.04657

С помощью подхода ANOVA мы можем получить тот же результат:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Теперь я могу получить тот же результат в IME со следующей моделью, предполагая положительно определенную симметричную корреляционную матрицу для двух условий:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

Или другая модель, предполагающая составную симметрию для корреляционной матрицы двух условий:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

С парным t-тестом и односторонним повторным измерением ANOVA я могу записать традиционную модель среднего значения ячейки как

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

где i индексирует условие, j индексирует субъект, Y ij - переменная отклика, µ - постоянная для фиксированного эффекта для общего среднего значения, α i - фиксированный эффект для условия, β j - случайный эффект для субъекта, следующего за N (0, σ p 2 ) (σ p 2 - дисперсия совокупности), а ε ij - остаточное число после N (0, σ 2 ) (σ 2 - дисперсия внутри объекта).

Я думал, что приведенная выше модель среднего значения ячейки не подходит для моделей Ime, но проблема в том, что я не могу придумать разумную модель для двух подходов Ime () с предположением о структуре корреляции. Причина в том, что модель Ime, кажется, имеет больше параметров для случайных компонентов, чем предложенная выше модель среднего значения ячейки. По крайней мере, модель lme обеспечивает точно такое же значение F, степени свободы и значение p, что не может сделать gls. Более конкретно, gls дает неправильные DFs из-за того, что он не учитывает тот факт, что у каждого субъекта есть два наблюдения, что приводит к сильно раздутым DFs. Скорее всего, модель lme чрезмерно параметризована при определении случайных эффектов, но я не знаю, что это за модель и каковы параметры. Так что проблема до сих пор не решена для меня.

bluepole
источник
2
Не совсем уверен, что вы спрашиваете. Модель, которую вы записали, является именно моделью для модели случайных эффектов; корреляционная структура индуцируется случайным эффектом.
Аарон - Восстановить Монику
@ Аарон: случайный эффект βj в модели среднего значения ячейки должен следовать за N (0, σp2). Я путаюсь, как этот термин (только с одним параметром σp2) связан со структурой корреляции, определяемой либо сложной симметрией, либо простой симметричной матрицей в модели Ime?
bluepole
Когда вы вычисляете корреляцию между двумя наблюдениями по одному и тому же предмету, корреляция будет sigma_p ^ 2 / (sigma_p ^ 2 + sigma ^ 2), потому что они имеют один и тот же бета_j. См. Пинейро / Бейтс, стр. 8. Кроме того, модель случайного эффекта, как вы написали, эквивалентна сложной симметрии; другие корреляционные структуры являются более сложными.
Аарон - Восстановить Монику
@ Аарон: Спасибо! Я уже читал об этом книгу Пиньейру / Бейтса и до сих пор не могу выяснить особенности случайных эффектов. Более релевантные страницы кажутся примером на стр.160-161. Кроме того, вывод случайных эффектов из lme () с предположением о сложной симметрии, похоже, не согласуется с корреляцией σp2 / (σp2 + σ2) в модели среднего значения ячейки. Все еще сбит с толку по поводу структуры модели.
bluepole
Ну, почти эквивалентно сложной симметрии; в CS корреляция может быть отрицательной, но не со случайными эффектами. Возможно, в этом ваша разница. См. Stats.stackexchange.com/a/14185/3601 для подробной информации.
Аарон - Восстановить Монику

Ответы:

16

Эквивалентность моделей можно наблюдать, рассчитав корреляцию между двумя наблюдениями от одного и того же человека следующим образом:

YяJзнак равноμ+αя+βJ+εяJβJ~N(0,σп2)εяJ~N(0,σ2)Соv(YяК,YJК)знак равноСоv(μ+αя+βК+εяК,μ+αJ+βК+εJК)знак равноСоv(βК,βК)знак равноσп2σ 2 p / ( σ 2 p + σ 2 )Вaр(YяК)знак равноВaр(YJК)знак равноσп2+σ2 , поэтому соотношение равно .σп2/(σп2+σ2)

Обратите внимание, что модели, однако, не совсем эквивалентны, так как модель случайного эффекта заставляет корреляцию быть положительной. Модель CS и модель t-test / anova - нет.

РЕДАКТИРОВАТЬ: Есть два других различия, а также. Во-первых, модели CS и случайного эффекта предполагают нормальность для случайного эффекта, а модель t-test / anova - нет. Во-вторых, модели CS и случайного эффекта подбираются с использованием максимального правдоподобия, в то время как анова подбирается с использованием средних квадратов; когда все сбалансировано, они согласятся, но не обязательно в более сложных ситуациях. Наконец, я бы с осторожностью использовал значения F / df / p из различных подгонок в качестве показателей того, насколько модели согласуются; см. знаменитую стяжку Дуга Бейтса на df's для более подробной информации. (КОНЕЦ РЕДАКТИРОВАНИЯ)

Проблема с вашим Rкодом заключается в том, что вы не указали структуру корреляции должным образом. Вам нужно использовать glsс corCompSymmкорреляционной структурой.

Создайте данные, чтобы получить эффект объекта:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

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

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

Стандартные ошибки из модели случайных эффектов:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

И корреляция и остаточное отклонение от модели CS:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

И они равны тому, что ожидается:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Другие корреляционные структуры обычно не соответствуют случайным эффектам, а просто указывают желаемую структуру; Одним из распространенных исключений является модель случайного эффекта AR (1) +, которая имеет случайный эффект и корреляцию AR (1) между наблюдениями одного и того же случайного эффекта.

РЕДАКТИРОВАТЬ 2: Когда я подхожу к трем вариантам, я получаю точно такие же результаты, за исключением того, что gls не пытается угадать df для рассматриваемого срока.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(Перехват здесь другой, потому что при кодировании по умолчанию это не среднее значение для всех предметов, а среднее значение для первого предмета.)

Также интересно отметить, что более новый lme4пакет дает те же результаты, но даже не пытается вычислить p-значение.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283
Аарон - Восстановить Монику
источник
Спасибо еще раз за помощь! Я знаю эту часть с точки зрения модели среднего значения. Однако со следующим результатом lme () с составной симметрией: Случайные эффекты: Формула: ~ x - 1 | subj Структура: составная симметрия StdDev xx1 1.1913363 xx2 1.1913363 Corr: -0.036 Остаточный 0.4466733. Я до сих пор не могу согласовать эти цифры со средней моделью ячейки. Может быть, вы можете помочь мне разобраться в этих числах?
bluepole
Кроме того, есть мысли о формулировке модели с другими корреляционными структурами, такими как простая симметричная матрица?
bluepole
Я вижу! Я должен был прочитать ваш ответ в другой ветке более внимательно. Я думал об использовании gls () раньше, но не смог выяснить спецификацию корреляции. Интересно, что lme () со сложной структурой симметрии для случайного эффекта по-прежнему отображает то же значение t, но, похоже, дисперсии для случайных эффектов не могут быть интерпретированы напрямую. Я очень ценю вашу помощь!
bluepole
После некоторой второй мысли я чувствую, что мое первоначальное замешательство все еще не разрешено. Да, gls можно использовать для демонстрации корреляционной структуры и среднеквадратичных ромов, но модель под ней не совсем совпадает с парным t-тестом (или ANOVA с односторонним повторным измерением в целом), и такая оценка далее поддерживаются неверные DFs и p-значение из gls. Напротив, моя команда lme с составной симметрией обеспечивает те же значения F, DF и p. Единственное, о чем я озадачен, это то, как параметризована модель lme, как указано в моем первоначальном посте. Любая помощь там?
bluepole
Не уверен, как вам помочь. Не могли бы вы написать, что вы думаете о двух разных моделях? Что-то не так в том, как вы думаете об одном из них.
Аарон - Восстановить Монику
3

Вы также можете рассмотреть возможность использования функции mixedв пакете afexдля возврата значений p с приближением Кенварда-Роджера df, которое возвращает идентичные значения p в виде парного t-теста:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

Или

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")
Том Венселерс
источник