Come adattare un modello di effetti casuali con Soggetto come casuale in R?

Question

dati forniti del seguente modulo

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

vorrei modellare punteggio in funzione del soggetto, Stato e Ora. Ogni punteggio del soggetto (umano) è stato misurato tre volte, indicato dalla variabile Tempo, quindi ho ripetuto le misure.

Come posso creare in R un modello di effetti casuali con gli effetti del soggetto montati come casuali?

ADDENDUM: è stato chiesto come ho generato questi dati. Hai indovinato, i dati sono falsi come il giorno è lungo. Il punteggio è il tempo più il rumore casuale e l'essere in Condizione 1 aggiunge un punto a Punteggio. È istruttivo come un tipico setup psicologico. Hai un compito in cui il punteggio delle persone migliora con la pratica (tempo) e un farmaco (condizione == 1) che migliora il punteggio.

Ecco alcuni dati più realistici per gli scopi di questa discussione. Ora i partecipanti simulati hanno un livello di "abilità" casuale che viene aggiunto ai loro punteggi. Inoltre, i fattori sono ora stringhe.

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

Guardalo:

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

Answer 1

utilizzando la libreria NLME ...

rispondere alla tua domanda dichiarato, è possibile creare un modello di intecept effetto casuale misto utilizzando il seguente codice:

> 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 

La varianza di intercettazione è fondamentalmente 0, che indica no all'interno dell'effetto soggetto, quindi questo modello non è c aprendo bene la relazione tra il tempo. Un modello di intercettazione casuale è raramente il tipo di modello che si desidera per un disegno di misure ripetute. Un modello di intercettazione casuale presuppone che le correlazioni tra tutti i punti di tempo siano uguali. cioè la correlazione tra il tempo 1 e il tempo 2 è la stessa del tempo 1 e del tempo 3. In circostanze normali (forse non quelli che generano i tuoi dati falsi) ci aspetteremmo che il futuro sia inferiore al primo. Una struttura auto regressiva è di solito un modo migliore per andare.

> 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 

I dati mostrano un -.596 tra la correlazione del punto temporale, che sembra strano. normalmente dovrebbe esserci, come minimo, una correlazione positiva tra i punti temporali. Come sono stati generati questi dati?

addendum:

Con i nuovi dati sappiamo che il processo di generazione dei dati è equivalente ad un modello di intercettazione casuale (anche se questo non è il più realistico per uno studio longitudinale La visualizzazione mostra che l'effetto del tempo sembra. ad essere abbastanza lineare, quindi dovremmo stare tranquillo trattandolo come una variabile numerica.

> 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 

noi vediamo un significativo effetto Condizione, che indica che la condizione 'sì' tende ad avere punteggi più alti (da circa 1,7), e un significativo effetto temporale, che indica che entrambi i gruppi salgono nel tempo. Supportando la trama, abbiamo f ind senza effetto differenziale di tempo tra i due gruppi (l'interazione). cioè le pendenze sono le stesse.

Answer 2

Non è una risposta alla tua domanda, ma potresti trovare questa visualizzazione dei tuoi dati informativa.

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

Answer 3

(usando lme4 biblioteca) Questo si inserisce l'effetto soggetto come casuale e anche la variabile che i vostri effetti casuali sono raggruppati sotto. In questo modello l'effetto casuale è l'intercettazione che varia per soggetto.

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

Per vedere gli effetti casuali si può semplicemente utilizzare

ranef(m) 

come Ian Fellows accennato, i tuoi dati potranno anche avere componenti Condizione e tempo casuale. Puoi testarlo con un altro modello. In quella sotto Condizione, Tempo e l'intercettazione possono variare a caso per soggetto. Valuta anche le loro correlazioni.

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

e cercare

summary(mi) 
ranef(mi) 

Si potrebbe anche verificare questo senza correlazioni con l'intercetta, con interazioni tra Stato e ora, e numerosi altri modelli per vedere che meglio si adatta i vostri dati e/o il vostro teoria. La tua domanda è un po 'vaga, ma questi pochi comandi dovrebbero farti iniziare.

Nota che il soggetto è il tuo fattore di raggruppamento, quindi è quello che si adatta ad altri effetti come casuale sotto. Non è qualcosa che poi esplicitamente si adatta anche a un predittore.