termine di errore aov() in R: qual è la differenza bw Errore (id) e errore (id/timevar) specifica?

Question

Qual è la differenza tra le specifiche della formula aov(depvar~timevar+Error(id)) e aov(depvar~timevar+Error(id/timevar))? Queste due varianti producono risultati leggermente diversi.

La stessa domanda è stata posta una volta qui: https://stats.stackexchange.com/questions/60108/how-to-write-the-error-term-in-repeated-measures-anova-in-r Tuttavia, mi piacerebbe ripeterlo con un esempio più appropriato.

Ecco un esempio che ho creato:

var=rep(NA,180) 
id=rep(1:20,each=180/20) 
group=rep(rep(1:2,each=9),180/(9*2)) 
time1=rep(rep(1:3,each=3),180/(3*3)) 
time2=rep(c(8,15,20),180/3) 

var[group==1&time1==1&time2==8]=runif(10,105,115) 
var[group==2&time1==1&time2==8]=runif(10,105,115) 
var[group==1&time1==1&time2==15]=runif(10,95,105) 
var[group==2&time1==1&time2==15]=runif(10,95,105) 
var[group==1&time1==1&time2==20]=runif(10,85,95) 
var[group==2&time1==1&time2==20]=runif(10,85,95) 

var[group==1&time1==2&time2==8]=runif(10,95,105) 
var[group==2&time1==2&time2==8]=runif(10,95,105) 
var[group==1&time1==2&time2==15]=runif(10,85,95) 
var[group==2&time1==2&time2==15]=runif(10,75,85) 
var[group==1&time1==2&time2==20]=runif(10,75,85) 
var[group==2&time1==2&time2==20]=runif(10,65,75) 

var[group==1&time1==3&time2==8]=runif(10,95,105) 
var[group==2&time1==3&time2==8]=runif(10,95,105) 
var[group==1&time1==3&time2==15]=runif(10,85,95) 
var[group==2&time1==3&time2==15]=runif(10,75,85) 
var[group==1&time1==3&time2==20]=runif(10,75,85) 
var[group==2&time1==3&time2==20]=runif(10,65,75) 

df=data.frame(id,var,group,time1,time2) 
df$id=factor(df$id) 
df$group=factor(df$group) 
df$time1=factor(df$time1) 
df$time2=factor(df$time2) 

Performing AOV() su questo ottiene risultati leggermente diversi a seconda di errore() specifica termine:

termine solo per una volta:

> summary(aov(var~time1+Error(id),data=df)) 
Error: id 
     Df Sum Sq Mean Sq F value Pr(>F) 
     Residuals 19 958.4 50.44    
Error: Within 
     Df Sum Sq Mean Sq F value Pr(>F)  
     time1  2 7538 3769 30.41 6.72e-12 *** 
     Residuals 158 19584  124   

> summary(aov(var~time1+Error(id/time1),data=df)) 
Error: id 
     Df Sum Sq Mean Sq F value Pr(>F) 
     Residuals 19 958.4 50.44    
Error: id:time1 
     Df Sum Sq Mean Sq F value Pr(>F)  
     time1  2 7538 3769 211.5 <2e-16 *** 
     Residuals 38 677  18     
     --- 
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
Error: Within 
     Df Sum Sq Mean Sq F value Pr(>F) 
     Residuals 120 18907 157.6  

O per entrambi i termini di tempo (non digitare l'output qui per motivi di spazio, è possibile controllare da soli):

summary(aov(var~group*time1*time2+Error(id/(group*time1*time2)),data=df)) 
summary(aov(var~group*time1*time2+Error(id),data=df)) 

Perché succede? Quale variante è corretta?

Answer 1

Ecco un post sul blog che ti aiuterà a capire cosa significhi sotto la sezione "Effetti casuali in ANOVA classica".

Dal blog, ecco un riepilogo di cosa significa "dividere" nel termine Error.

aov(Y ~ Error(A), data=d)    # Lone random effect 
aov(Y ~ B + Error(A/B), data=d)   # A random, B fixed, B nested within A 
aov(Y ~ (B*X) + Error(A/(B*X)), data=d) # B and X interact within levels of A 

Quindi dalla tua domanda,

aov(depvar~timevar+Error(id/timevar)) 

significa avere un effetto casuale da id ma poi correggere timevar con timevar nidificato all'interno id livelli contro

aov(depvar~timevar+Error(id)) 

che è solo prendendo id come gli effetti casuali senza vincoli su altre variabili.

Fonte

: http://conjugateprior.org/2013/01/formulae-in-r-anova/

This potrebbe rivelarsi utile pure, che è un codice che va oltre l'analisi della varianza che ha alcune raccomandazioni su di apprendimento ANOVA.

Answer 2

La differenza tra aov(depvar~timevar+Error(id)) e aov(depvar~timevar+Error(id/timevar)) indica se si include o meno timevar come un effetto casuale.

Nota che c'è più di un modo per includere una variabile come effetto casuale. Puoi anche utilizzare aov(depvar~timevar+Error(id*timevar)) o aov(depvar~timevar+Error(id + timevar)). Ognuno di questi significa qualcosa di molto diverso, ma può essere fonte di confusione perché spesso ti danno risultati simili se applicati allo stesso set di dati, a causa dei vincoli dei dati stessi.

La barra / utilizzata in aov() indica nesting. Quando si utilizza /, R lo espande automaticamente all'effetto principale della variabile inferiore più l'interazione tra il fondo e la parte superiore. Ad esempio, A/B si espande automaticamente a A + A:B.Questo è simile al modo in cui A*B si espande automaticamente a A + B + A:B, ma con l'annidamento, la variabile nel nido non appare mai al di fuori del proprio nido (vale a dire non può esserci alcun effetto principale di B da solo).

Si può vedere questa espansione accadendo nella vostra uscita:

> summary(aov(var~time1+Error(id/time1))) 

Error: id 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 1 52.24 52.24    

Error: id:time1 
     Df Sum Sq Mean Sq 
time1 1 4291 4291 

Error: Within 
      Df Sum Sq Mean Sq F value Pr(>F) 
time1  1 1239 1238.7 10.19 0.00167 ** 
Residuals 176 21399 121.6     
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

I Error termini denotano effetti casuali. Si noti che si ottiene uno per l'effetto principale di id perché è la base del nido, e uno per l'interazione tra id e time1, perché time1 è nidificato all'interno id (si ottiene anche un effetto Error per Within che è la durata residua di base per il modello, cioè l'effetto casuale delle singole osservazioni stesse).

Quindi qual è l'approccio corretto per i dati?

Dipende da 1) in che modo i dati sono effettivamente strutturati e 2) quale modello si intende eseguire. Nota: non esiste un test definitivo da eseguire sui dati per determinare la struttura o il modello corretto; questo è un esercizio di pensiero piuttosto che uno computazionale.

Nei modelli di esempio che hai fornito, voi hanno un risultato var, e poi quelli che sembrano essere variabili di raggruppamento group e id, e poi due variabili temporali time1 e time2. Ogni id è solo in un gruppo, non in entrambi i gruppi, suggerendo che l'id è nidificata all'interno del gruppo.

> table(group, id) 
    id 
group 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
    1 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 
    2 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 0 9 

supporrò che id riferisce a un singolo partecipante, e le misure su 9 time1 e time2 sono entro soggetti test di ciascun partecipante (cioè ogni partecipante è stato misurato 9 volte su var, quindi questo è un repeated measures design).

Per rendere più concreto, diciamo var è un punteggio su qualche compito risolvere il problema, e time1 e time2 sono i minuti partecipanti sono autorizzati a studiare il problema e la quantità di tempo sono dato per completare il problema, rispettivamente, . Poiché i valori time1 e time2 vengono attraversati, ciascun partecipante completa l'attività 9 volte, in base a ciascuna combinazione di circostanze.

> table(time1, time2) 
    time2 
time1 8 15 20 
    1 20 20 20 
    2 20 20 20 
    3 20 20 20 
> table(time1, time2, id) 
, , id = 1 

    time2 
time1 8 15 20 
    1 1 1 1 
    2 1 1 1 
    3 1 1 1 

, , id = 2 

    time2 
time1 8 15 20 
    1 1 1 1 
    2 1 1 1 
    3 1 1 1 
(output truncated) 

I partecipanti sono testati in gruppi, con la metà dei partecipanti nel gruppo 1 e l'altra metà nel gruppo 2. Forse lo studio è stato eseguito in aule, e il gruppo 1 è una classe e gruppo 2 è la seconda classe . Probabilmente, l'identità di gruppo non è in realtà una variabile di interesse, ma non dovremmo lasciarla fuori dal modello perché potrebbe esserci qualche variazione di disturbo dovuta alle differenze tra i gruppi. Ad esempio, forse la prima aula ha una luce migliore, dando a tutti i membri del gruppo 1 e migliori possibilità di segnare bene sui puzzle rispetto ai membri del gruppo 2.

Punteggi, ID e gruppo dovrebbero essere tutti effetti casuali e time1 e time2 dovrebbero essere effetti fissi (si noti che potrebbe variare per gli stessi dati se si hanno pensieri diversi nel modello, ad esempio si potrebbe voler considerare il gruppo come fisso a seconda della domanda di ricerca).

Dato che il modello, questa sarebbe la versione più completa del modello, utilizzando aov():

aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df) 

ecco l'output:

> summary(aov(var~time1*time2 + Error(group/id/(time1*time2)),data=df)) 

Error: group 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 1 771.7 771.7    

Error: group:id 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 18 243.8 13.55    

Error: group:id:time1 
      Df Sum Sq Mean Sq F value Pr(>F)  
time1  2 7141 3571 181.6 <2e-16 *** 
Residuals 38 747  20     
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: group:id:time2 
      Df Sum Sq Mean Sq F value Pr(>F)  
time2  2 16353 8176 434.6 <2e-16 *** 
Residuals 38 715  19     
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: group:id:time1:time2 
      Df Sum Sq Mean Sq F value Pr(>F) 
time1:time2 4 214.5 53.63 5.131 0.00103 ** 
Residuals 76 794.3 10.45     
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
Warning message: 
In aov(var ~ time1 * time2 + Error(group/id/(time1 * time2)), data = df) : 
    Error() model is singular 

(Insieme con i link di cui sopra, ecco alcune ulteriori indicazioni su random vs. fixed effects)