Как подобрать модель случайных эффектов с субъектом как случайным в R?

Приведенные данные следующего вида

myDat = structure(list(Score = c(1.84, 2.24, 3.8, 2.3, 3.8, 4.55, 1.13, 
2.49, 3.74, 2.84, 3.3, 4.82, 1.74, 2.89, 3.39, 2.08, 3.99, 4.07, 
1.93, 2.39, 3.63, 2.55, 3.09, 4.76), Subject = c(1L, 1L, 1L, 
2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 
7L, 7L, 8L, 8L, 8L), Condition = c(0L, 0L, 0L, 1L, 1L, 1L, 0L, 
0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 
1L), Time = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L)), .Names = c("Score", 
"Subject", "Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))

Я хотел бы моделировать Score как функцию Subject, Condition и Time. Каждый (человеческий) предметный счет был измерен три раза, обозначенный переменным временем, поэтому я повторил меры.

Как я могу построить в R случайную модель эффектов с эффектами субъекта, установленными как случайные?

ДОБАВЛЕНИЕ: спросили, как я создал эти данные. Вы догадались, данные поддельные, поскольку день длинный. Оценка - это время плюс случайный шум, и в состоянии 1 добавляется точка к счету. Это поучительно, как типичная психологическая установка. У вас есть задача, в которой люди набирают очки с практикой (время) и наркотиками (условие == 1), что увеличивает баллы.

Вот некоторые более реалистичные данные для целей этого обсуждения. Теперь моделируемые участники имеют случайный уровень "умения", который добавляется к их баллам. Кроме того, теперь факторы являются строками.

myDat = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41, 
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01, 
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E", 
"F", "G", "H"), class = "factor"), Condition = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"), 
    Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM", 
    "2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject", 
"Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))

Посмотрите:

library(ggplot2)
qplot(Time, Score, data = myDat, geom = "line", group = Subject, colour = factor(Condition))

Ответы

Ответ 1

используя библиотеку nlme...

Отвечая на указанный вами вопрос, вы можете создать случайную модель смешанного эффекта intecept, используя следующий код:

> library(nlme)
> m1 <- lme(Score ~ Condition + Time + Condition*Time,
+ data = myDat, random = ~ 1 | Subject)
> summary(m1)
Linear mixed-effects model fit by REML
 Data: myDat 
       AIC      BIC    logLik
  31.69207 37.66646 -9.846036

Random effects:
 Formula: ~1 | Subject
         (Intercept)  Residual
StdDev: 5.214638e-06 0.3151035

Fixed effects: Score ~ Condition + Time + Condition * Time 
                   Value Std.Error DF  t-value p-value
(Intercept)    0.6208333 0.2406643 14 2.579666  0.0218
Condition      0.7841667 0.3403507  6 2.303996  0.0608
Time           0.9900000 0.1114059 14 8.886423  0.0000
Condition:Time 0.0637500 0.1575517 14 0.404629  0.6919
 Correlation: 
               (Intr) Condtn Time  
Condition      -0.707              
Time           -0.926  0.655       
Condition:Time  0.655 -0.926 -0.707

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.5748794 -0.6704147  0.2069426  0.7467785  1.5153752 

Number of Observations: 24
Number of Groups: 8

Перехват перехвата в основном равен 0, что указывает на отсутствие эффекта субъекта, поэтому эта модель не фиксирует между отношениями времени. Случайная модель перехвата редко является типом модели, которую вы хотите для повторного проектирования мер. Случайная модель перехвата предполагает, что корреляции между всеми точками времени равны. то есть корреляция между временем 1 и временем 2 является такой же, как между временем 1 и временем 3. При нормальных обстоятельствах (возможно, не в том, что генерирует ваши поддельные данные), мы ожидаем, что позднее будет меньше первого. Авторегрессивная структура, как правило, лучше всего подходит.

> m2<-gls(Score ~ Condition + Time + Condition*Time,
+ data = myDat, correlation = corAR1(form = ~ Time | Subject))
> summary(m2)
Generalized least squares fit by REML
  Model: Score ~ Condition + Time + Condition * Time 
  Data: myDat 
       AIC      BIC    logLik
  25.45446 31.42886 -6.727232

Correlation Structure: AR(1)
 Formula: ~Time | Subject 
 Parameter estimate(s):
       Phi 
-0.5957973 

Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)    0.6045402 0.1762743  3.429543  0.0027
Condition      0.8058448 0.2492895  3.232566  0.0042
Time           0.9900000 0.0845312 11.711652  0.0000
Condition:Time 0.0637500 0.1195452  0.533271  0.5997

 Correlation: 
               (Intr) Condtn Time  
Condition      -0.707              
Time           -0.959  0.678       
Condition:Time  0.678 -0.959 -0.707

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.6850557 -0.6730898  0.2373639  0.8269703  1.5858942 

Residual standard error: 0.2976964 
Degrees of freedom: 24 total; 20 residual

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

Приложение:

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

> library(nlme)
> m1 <- lme(Score ~ Condition + as.numeric(Time) + Condition*as.numeric(Time),
+ data = myDat, random = ~ 1 | Subject)
> summary(m1)
Linear mixed-effects model fit by REML
 Data: myDat 
       AIC      BIC    logLik
  38.15055 44.12494 -13.07527

Random effects:
 Formula: ~1 | Subject
        (Intercept)  Residual
StdDev:   0.2457355 0.3173421

Fixed effects: Score ~ Condition + as.numeric(Time) + Condition * as.numeric(Time) 
                                  Value Std.Error DF   t-value p-value
(Intercept)                    1.142500 0.2717382 14  4.204415  0.0009
ConditionYes                   1.748333 0.3842958  6  4.549447  0.0039
as.numeric(Time)               0.575000 0.1121974 14  5.124898  0.0002
ConditionYes:as.numeric(Time) -0.197500 0.1586710 14 -1.244714  0.2337
 Correlation: 
                              (Intr) CndtnY as.(T)
ConditionYes                  -0.707              
as.numeric(Time)              -0.826  0.584       
ConditionYes:as.numeric(Time)  0.584 -0.826 -0.707

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.44560367 -0.65018585  0.01864079  0.52930925  1.40824838 

Number of Observations: 24
Number of Groups: 8

Мы видим значительный эффект условия, указывающий, что условие "да" имеет более высокие оценки (примерно на 1,7) и значительный временный эффект, указывающий на то, что обе группы со временем увеличиваются. Поддерживая сюжет, мы не находим дифференциального влияния времени между двумя группами (взаимодействие). т.е. склоны одинаковы.

Ответ 2

Это не ответ на ваш вопрос, но вы можете найти эту визуализацию своих информационных данных.

library(ggplot2)
qplot(Time, Score, data = myDat, geom = "line", 
  group = Subject, colour = factor(Condition))

Data visulation

Ответ 3

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

m <- lmer( Score ~ Condition + Time + (1|Subject), data=myDat )

Чтобы увидеть случайные эффекты, вы можете просто использовать

ranef(m)

Как отметил Иан Феллоуз, ваши данные могут также иметь случайные компоненты Condition и Time. Вы можете проверить это с помощью другой модели. В приведенном ниже состоянии, времени и перехвате разрешено изменять случайным образом субъект. Он также оценивает их корреляции.

mi <- lmer( Score ~ Condition + Time + (Condition + Time|Subject), data=myDat )

и попробуйте

summary(mi)
ranef(mi)

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

Обратите внимание, что Subject - ваш фактор группировки, так что вы подходите под другие эффекты как случайные. Это не то, что вы тогда явно вписываете как предиктор.