paarweise ANOVA der Teilmenge der Daten

Question

Ich muss mehrere paarweise ANOVA in R durchführen, und korrigieren Sie die p-Werte mit Bonferroni. Allerdings muss ich nicht jede KLASSE miteinander vergleichen. Unten ist mein Datenformat und : die Paare, von denen ich das log10relquant kontrastieren muss. Kennt jemand von euch, wie ich das ausführen könnte? Ich benutze die dplyr, lsmeans und broom Pakete.

SEX  EXPERIENCED AGE CLASS compound relquant log10relquant 

1 FEMALE   NO  1D  1F  C14 0.004012910  -2.396541 
2 FEMALE   NO  1D  1F  C14 0.003759812  -2.424834 
3 FEMALE   NO  1D  1F  C14 0.003838553  -2.415832 
4 FEMALE   NO  1D  1F  C14 0.003582754  -2.445783 
5 MALE   NO  1D  1M  C14 0.005099237  -2.292495 
6 MALE   NO  1D  1M  C14 0.005379093  -2.269291 

selcontrasts <- c("1F - 1M", "4F - 4M", "4EF - 4EM", 
        "7F - 7M", "7EF - 7EM", # sex differences 
       "1M - 4M", "4M - 7M", "1M - 7M", "1F - 4F", 
       "4F - 7F", "1F - 7F", # age differences 
       "4M - 4EM", "7M - 7EM", "4F - 4EF", 
       "7F - 7EF" # social experience) 

x=list(selcontrasts) 

Derzeit ist dieses Ich bin mit dem gesamten Datensatz zu paaren (so jede Klasse Compair) anstelle der ausgewählten Kontraste:

pvalsage=data.frame(datagr %>% 
    do(data.frame(summary(contrast(lsmeans(
      aov(log10relquant ~ CLASS, data = .), ~ CLASS),    
      method="pairwise",adjust="none"))))) 

Um nur die ausgewählten Kontraste in der Liste x zu tun, habe ich versucht, :

pvalsage=data.frame(datagr %>% 
    do(data.frame(summary(contrast(lsmeans(
     aov(log10relquant ~ CLASS, data = .),~ CLASS), 
     method = x, adjust="none"))))) 

Aber ich bekomme die Fehlermeldung:

error in contrast.ref.grid(lsmeans(aov(log10relquant ~ CLASS, data = .), : 
Nonconforming number of contrast coefficients 

Answer 1

Wenn ich die Frage richtig verstehe (und ich sehr wohl nicht), sind wirklich drei Faktoren beteiligt: ​​SEX (zwei Ebenen), EXPERIENCED (zwei Ebenen) und AGE (3 Ebenen, nämlich 1, 4 und 7) . Und was erforderlich ist, sind separate Vergleiche der Ebenen jedes Faktors für jede Kombination der anderen beiden.

Wenn dies der Fall ist, dann macht es die Kombination der drei Faktoren in einem namens CLASS nur viel schwieriger, weil es es viel schwieriger macht, die Ebenen der Faktoren getrennt zu verfolgen. Es ist einfacher, ein Modell zu erstellen, das alle drei Faktoren berücksichtigt, die Mittel für jede Faktorkombination zu schätzen und dann die erforderlichen Vergleiche unter Verwendung der Variablen by durchzuführen. So für jeden Datensatz dat, was Sie tun:

require(emmeans) 
mod = aov(log10relquant ~ SEX * EXPERIENCED * AGE, data = dat) 
emm = emmeans(mod, ~ SEX * EXPERIENCED * AGE) 
rbind(pairs(emm, by = c("EXPERIENCED", "AGE")), 
     pairs(emm, by = c("SEX", "EXPERIENCED")), 
     pairs(emm, by = c("SEX", "AGE")), 
     adjust = "bonferroni") 

Ich versuche nicht, diese in dem funktionalen-Programmierparadigma einzuzubetten; Ich überlasse es dem OP, diese Details herauszufinden.

Anmerkung: Das emmeans Paket (geschätztes Randmittel) ist eine Fortsetzung des lsmeans, die in der Zukunft weiterentwickelt werden. Es funktioniert genauso.

PS - Mit Blick auf den Code in der Frage, bin ich besorgt, dass die Endergebnisse nicht zeigen die tatsächlichen Schätzungen verglichen werden (die EMMs), nur die Vergleiche; und die weitere Implikation bei der Benennung, dass wirklich nur P-Werte gesucht werden. Das grabt mich an. Ich mag es nicht, Leuten zuzusehen, die direkt zu statistischen Tests gehen, ohne auch nur die Mengen zu betrachten, die getestet werden.

Answer 2

Sie könnten trotzdem den paarweisen Kontrast machen und dann die Zeilen, die nur Ihre selcontrasts enthalten, in einen neuen Datenrahmen filtern, gefolgt von p.adjust = bonferroni mit nur den Kontrasten, die Sie interessieren.

oder Sie können eine mycontr.lmsc Funktion und definieren selcontrasts schreiben und verwenden, die als Methode =

(Y)