partizione di anova e confronti (ortogonale singolo df) in r

Question

Voglio fare singolo contrasto df ortogonale in anova (modello fisso o misto). Qui è solo esempio:

require(nlme) 
data (Alfalfa) 
    Variety: a factor with levels Cossack, Ladak, and Ranger 
    Date : a factor with levels None S1 S20 O7 
    Block: a factor with levels 1 2 3 4 5 6 
    Yield : a numeric vector 

Questi dati sono descritti in Snedecor e Cochran (1980) come esempio di un disegno split-plot. La struttura di trattamento utilizzata nell'esperimento era un fattoriale completo di 3 \ times4, con tre varietà di erba medica e quattro date del terzo taglio nel 1943. Le unità sperimentali erano organizzate in in sei blocchi, ciascuno suddiviso in quattro lotti. Le varietà di erba medica (Cossac, Ladak e Ranger) sono state assegnate casualmente ai blocchi e le date del terzo taglio (Nessuna, S1-1 settembre, S20-settembre 20, e O7-7 ottobre) sono state assegnate in modo casuale a le trame. Tutte e quattro le date sono state utilizzate su ciascun blocco.

model<-with (Alfalfa, aov(Yield~Variety*Date +Error(Block/Date/Variety))) 

    > summary(model) 

Error: Block 
      Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

voglio eseguire alcune confronto (contrasti ortogonali all'interno di un gruppo), ad esempio per data, due contrasti:

(a) S1 vs others (S20 O7) 
    (b) S20 vs 07, 

per il fattore varietà due contrasti:

(c) Cossack vs others (Ladak and Ranger) 
    (d) Ladak vs Ranger 

Così l'output di anova sarebbe:

Error: Block 
       Df Sum Sq Mean Sq F value Pr(>F) 
Residuals 5 4.15 0.83 

Error: Block:Date 
      Df Sum Sq Mean Sq F value Pr(>F) 
Date  3 1.9625 0.6542 17.84 3.29e-05 *** 
     (a) S1 vs others ?  ? 
     (b) S20 vs 07  ?  ? 
Residuals 15 0.5501 0.0367 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Block:Date:Variety 
      Df Sum Sq Mean Sq F value Pr(>F) 
Variety  2 0.1780 0.08901 1.719 0.192 
    (c) Cossack vs others ?  ? ? 
    (d) Ladak vs Ranger ?  ?  ? 
Variety:Date 6 0.2106 0.03509 0.678 0.668 
Residuals 40 2.0708 0.05177 

Come posso eseguire questo? ....................

Answer 1

Prima di tutto, perché utilizzare ANOVA? È possibile utilizzare lme dal pacchetto nlme e, oltre ai test di ipotesi aov, si ottengono anche stime interpretabili delle dimensioni dell'effetto e le direzioni degli effetti. In ogni caso, mi vengono in mente due approcci:

• Specificare i contrasti sulle variabili manualmente, come spiegato here.
• Installare il pacchetto multcomp e utilizzare glht.

glht è un po 'supponente sui modelli che sono multivariati nei loro predittori. Per farla breve, però, se dovessi creare una matrice diagonale cm0 con le stesse dimen- sioni e dimensioni del modello vcov (supponiamo che sia una misura lme chiamata model0), allora summary(glht(model0,linfct=cm0)) dovrebbe fornire le stesse stime, SE e test statistiche come summary(model0)$tTable (ma errati valori di p). Ora, se si scherza con combinazioni lineari di righe da cm0 e si creano nuove matrici con lo stesso numero di colonne di cm0 ma queste combinazioni lineari come righe, alla fine si scoprirà il modello per creare una matrice che fornirà l'intercetta stima per ogni cella (confrontarla con predict(model0,level=0)). Ora, un'altra matrice con differenze tra varie righe di questa matrice ti darà le differenze tra i gruppi corrispondenti. Lo stesso approccio, ma con valori numerici impostati su 1 anziché su 0, può essere utilizzato per ottenere le stime dell'inclinazione per ogni cella. Quindi le differenze tra queste stime di inclinazione possono essere utilizzate per ottenere differenze di inclinazione tra i gruppi.

Tre cose da tenere a mente:

• Come ho già detto i valori di p stanno per essere sbagliato per modelli diversi lm, (forse, non hanno provato) aov, e alcuni modelli di sopravvivenza. Questo perché glht presuppone una distribuzione z anziché una distribuzione t per impostazione predefinita (eccetto per lm). Per ottenere valori di p corretti, prendere la statistica test glht calcola e fare manualmente 2*pt(abs(STAT),df=DF,lower=F) per ottenere il p-value a due code dove STAT è la statistica test restituito da glht e DF è il df dal tipo corrispondente di contrasto predefinito in summary(model0)$tTable.
• I contrasti probabilmente non testano più ipotesi indipendenti e è necessaria una correzione di test multipla, se non lo fosse già. Esegui i p-value tramite p.adjust.
• Questa è la mia personale distillazione di molti handwaving da parte di professori e colleghi e molta lettura di Crossvalidated e Stackoverflow su argomenti correlati. Potrei sbagliarmi in molti modi, e se lo sono, spero che qualcuno più ben informato ci corregga entrambi.