Jak dopasować model efektów losowych z Tematem jako losowym w R?

Question

z danych w następującej formie

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)) 

chciałbym modelować wynik w zależności od tematu, stan i czas. Ocenę każdego (ludzkiego) podmiotu mierzono trzykrotnie, co wskazuje zmienna Czas, więc powtarzam mi środki.

Jak mogę zbudować w R model losowy z efektami Temat dopasowany jako losowy?

ADDENDUM: Zapytano mnie, w jaki sposób wygenerowałem te dane. Zgadłeś, dane są fałszywe, ponieważ dzień jest długi. Wynik to czas plus losowy szum, a będąc w Warunku 1 dodaje punkt do wyniku. To pouczające, jak typowa konfiguracja Psych. Masz zadanie, w którym wynik jest lepszy dzięki praktyce (czas) i lekowi (warunek == 1), który poprawia wynik.

Oto kilka bardziej realistycznych danych do celów tej dyskusji. Teraz symulowani uczestnicy mają losowy poziom "umiejętności", który jest dodawany do ich wyników. Również czynniki są teraz łańcuchami.

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)) 

Zobacz:

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

Answer 1

wykorzystaniem biblioteki nlme ...

odpowiadając na wyrażoną pytanie, można stworzyć efekt losowy intecept model mieszany za pomocą następującego kodu:

> 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 

Wariancja punktu przecięcia jest w zasadzie 0, co wskazuje na brak efektu podmiotu, więc ten model nie jest odpowiednie dopasowanie czasu pomiędzy relacjami. Losowy model przechwytywania rzadko jest typem modelu, który ma być stosowany w projektach powtarzalnych pomiarów. Losowy model przechwytywania zakłada, że ​​korelacje między wszystkimi punktami czasowymi są równe. tj. korelacja między czasem 1 i czasem 2 jest taka sama jak między czasem 1 i czasem 3. W normalnych okolicznościach (być może nie tych generujących fałszywe dane) spodziewalibyśmy się, że później będzie mniej niż ten pierwszy. Automatyczna regresywna struktura jest zwykle lepszym rozwiązaniem.

> 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 

Twoje dane pokazują korelację -.596 między korelacją w punkcie czasowym, która wydaje się dziwna. zazwyczaj powinna istnieć co najmniej dodatnia korelacja między punktami czasowymi. W jaki sposób wygenerowano te dane?

Uzupełnienie:

Z nowych danych wiemy, że proces generujący dane jest odpowiednikiem modelu losowego przechwytującym (choć nie jest to najbardziej realistyczne dla podłużnego badania Wizualizacja pokazuje, że wpływ czasu wydaje. być dość liniowa, więc powinniśmy czuć się komfortowo, traktując ją jako zmienną numeryczną.

> 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 

Widzimy znaczący wpływ Stan, wskazując, że „tak” stan zwykle mają wyższe wyniki (około 1,7) i znaczący efekt czasowy, wskazujący, że obie grupy rosną w miarę upływu czasu, wspierając fabułę, my f nie ma różnicowego wpływu czasu między dwiema grupami (interakcja). tj. stoki są takie same.

Answer 2

To nie jest odpowiedź na twoje pytanie, ale ta wizualizacja danych może okazać się przydatna.

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

Answer 3

(używając lme4 biblioteki) ten pasuje do Twojego efekt losowy przedmiot, jak i również zmienna, że ​​losowe efekty są pogrupowane pod. W tym modelu efektem losowym jest przechwycenie zmieniające się w zależności od podmiotu.

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

aby zobaczyć, losowe efekty można po prostu użyć

ranef(m) 

Ian Fellows wspomniano, dane mogą również mieć losowych stanie i czasowe elementy. Możesz to przetestować z innym modelem. W tej warstwie poniżej Warunek, Czas i punkt przecięcia mogą różnić się losowo w zależności od tematu. Ocenia także ich korelacje.

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

i spróbuj

summary(mi) 
ranef(mi) 

Można również przetestować ten bez korelacji z osią, z interakcji pomiędzy stan i czas, i wielu innych modeli, aby zobaczyć, które najlepiej pasuje do Twoich danych oraz/lub teoria. Twoje pytanie jest trochę niejasne, ale te kilka poleceń powinno zacząć.

Pamiętaj, że tematem jest twój czynnik grupujący, dlatego dopasowujesz inne efekty do losowych. To nie jest coś, co następnie można wyraźnie określić jako predyktor.