0

R で複数のペアワイズ ANOVA を実行し、bonferroni を使用して p 値を修正する必要があります。ただし、すべてのクラスを互いに比較する必要はありません。以下は私のデータ形式とselcontrasts: log10relquant を対比する必要があるペアです。これを実行する方法を知っている人はいますか？dplyr、lsmeansおよびbroomパッケージを使用します。

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)

現在、選択したコントラストの代わりに、これを使用してデータセット全体をペアにしています(すべてのクラスを比較するため):

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

リストxで選択したコントラストのみを実行するために、次のことを試しました：

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

しかし、私はエラーが発生します:

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

2 に答える 2

0

私がこの質問を正しく理解していれば (おそらくそうではないかもしれませんが)、実際にはSEX(2 つのレベル)、EXPERIENCED(2 つのレベル)、およびAGE(3 つのレベル、つまり 1、4、および 7) の 3 つの要因が関係しています。必要なのは、他の 2 つの組み合わせごとに、各因子のレベルを個別に比較することです。

その場合、3 つの因子を名前付きの 1 つに結合するCLASSと、因子のレベルを個別に追跡することがはるかに難しくなるため、非常に難しくなります。より簡単なのは、3 つの要因すべてを説明するモデルを当てはめ、各要因の組み合わせの平均を推定し、by変数を使用して必要な比較を行うことです。したがって、各データセットdatに対して、次のことを行います。

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

これを関数型プログラミングのパラダイムに組み込もうとはしませんでした。これらの詳細を把握するために、OP に任せます。

注: emmeansパッケージ (推定された限界平均) はlsmeansの継続であり、将来的に非推奨になります。同じように機能します。

PS -- 質問のコードを見ると、比較される実際の見積もり (EMM) が最終結果に表示されず、比較のみが表示されることが懸念されます。そして、実際にはP値のみが求められているという命名のさらなる意味。これは私に感謝します。私は、人々がテストされている量を見ずに統計テストに直行するのを見るのは好きではありません.

于 2018-01-14T22:23:14.743 に答える
0

とにかくペアワイズコントラストを実行してから、セルコントラストのみを含む行を新しいデータフレームにフィルタリングし、その後に興味のあるコントラストのみを含む p.adjust= bonferroni を続けることができます。

または、mycontr.lmsc 関数を記述して selcontrasts を定義し、それを method = として使用することもできます

(Y)

于 2018-01-20T20:49:47.050 に答える