提取 `mvabund::anova.manyglm` 的 p 值并将测试作为一个 table

Extract `mvabund::anova.manyglm`'s p-values and tests as one table

我想从 anova.manyglm 的结果中获取单变量测试部分,但默认属性给出了两个单独的 table,所以我试图合并它们,以便测试和同一变量的 p 值列彼此相邻。

library(mvabund)
data("Tasmania")
attach(Tasmania)
tasmvabund=mvabund(copepods)
tas_mod = manyglm(tasmvabund~ block*treatment)
tas_sum <- summary.manyglm(tas_mod)
tas_aov <- anova.manyglm(tas_mod, resamp = "perm.resid",  p.uni= "adjusted", nBoot = 50)



  > str(tas_aov)
List of 11
 $ family      : chr "negative.binomial"
 $ p.uni       : chr "adjusted"
 $ test        : chr "Dev"
 $ cor.type    : chr "I"
 $ resamp      : chr "perm.resid"
 $ nBoot       : num 50
 $ shrink.param: num [1:4] 0 0 0 0
 $ n.bootsdone : num 50
 $ table       :'data.frame':   4 obs. of  4 variables:
  ..$ Res.Df  : int [1:4] 15 12 11 8
  ..$ Df.diff : num [1:4] NA 3 1 3
  ..$ Dev     : num [1:4] NA 117.5 66.9 37.4
  ..$ Pr(>Dev): num [1:4] NA 0.0196 0.0196 0.0588
  ..- attr(*, "heading")= chr [1:2] "Analysis of Deviance Table\n" "Model: tasmvabund ~ block * treatment"
  ..- attr(*, "title")= chr "\nMultivariate test:\n"
 $ uni.p       : num [1:4, 1:12] NA 0.4706 0.0196 0.451 NA ...
  ..- attr(*, "title")= chr "Univariate Tests:"
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "(Intercept)" "block" "treatment" "block:treatment"
  .. ..$ : chr [1:12] "Ameira" "Adopsyllus" "Ectinosoma" "Ectinosomat" ...
 $ uni.test    : num [1:4, 1:12] NA 4.93 13.94 6.35 NA ...
  ..- attr(*, "title")= chr "Univariate Tests:"
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "(Intercept)" "block" "treatment" "block:treatment"
  .. ..$ : chr [1:12] "Ameira" "Adopsyllus" "Ectinosoma" "Ectinosomat" ...
 - attr(*, "class")= chr "anova.manyglm"

这仅打印多变量部分。

tas_aov$table
'data.frame':   4 obs. of  4 variables:
 $ Res.Df  : int  15 12 11 8
 $ Df.diff : num  NA 3 1 3
 $ Dev     : num  NA 117.5 66.9 37.4
 $ Pr(>Dev): num  NA 0.0196 0.0196 0.0588
 - attr(*, "heading")= chr [1:2] "Analysis of Deviance Table\n" "Model: tasmvabund ~ block * treatment"
 - attr(*, "title")= chr "\nMultivariate test:\n"

这仅打印单变量检验的 p 值。

tas_aov$uni.p
 str(tas_aov$uni.p)
 num [1:4, 1:12] NA 0.4706 0.0196 0.451 NA ...
 - attr(*, "title")= chr "Univariate Tests:"
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:4] "(Intercept)" "block" "treatment" "block:treatment"
  ..$ : chr [1:12] "Ameira" "Adopsyllus" "Ectinosoma" "Ectinosomat" ...

这只打印测试 table。

tas_aov$uni.test
 str(tas_aov$uni.test)
 num [1:4, 1:12] NA 4.93 13.94 6.35 NA ...
 - attr(*, "title")= chr "Univariate Tests:"
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:4] "(Intercept)" "block" "treatment" "block:treatment"
  ..$ : chr [1:12] "Ameira" "Adopsyllus" "Ectinosoma" "Ectinosomat" ...

所以我试图合并它们,然后按名称排列列,但我还没有成功。

tab <- cbind(tas_aov$uni.p, tas_aov$uni.test)

tab[ , order(names(tab))]
Error in order(names(tab)) : argument 1 is not a vector

理想情况下,我只想要原始的 Univariate 部分,因为它整齐地打印了列名,即 p 值旁边的测试,但我很高兴只合并 tables。非常感谢任何提示。

一种方法是将两个矩阵旋转更长的时间,然后对它们进行内部连接:

library(tidyverse)
tas_aov$uni.test %>%
  as.data.frame %>%
  rownames_to_column("variable") %>%
  pivot_longer(-variable,names_to = "species", values_to = "test") %>%
inner_join({
tas_aov$uni.p %>%
    as.data.frame %>%
    rownames_to_column("variable") %>%
    pivot_longer(-variable,names_to = "species", values_to = "pval")}) %>%
dplyr::filter(!is.na(test))
# A tibble: 36 x 4
   variable species      test   pval
   <chr>    <chr>       <dbl>  <dbl>
 1 block    Ameira       4.93 0.549 
 2 block    Adopsyllus  17.6  0.0392
 3 block    Ectinosoma   7.71 0.549 
 4 block    Ectinosomat 10.9  0.235 
 5 block    Haloschizo   3.13 0.725 
 6 block    Lepta.A     20.9  0.0196
 7 block    Lepta.B      8.02 0.549 
 8 block    Lepta.C     13.8  0.0784
 9 block    Mictyricola  6.77 0.549 
10 block    Parevansula  6.10 0.549 
# … with 26 more rows

如果你愿意,你可以转向宽屏:

tas_aov$uni.test %>%
  as.data.frame %>%
  rownames_to_column("variable") %>%
  pivot_longer(-variable,names_to = "species", values_to = "test") %>%
inner_join({
tas_aov$uni.p %>%
    as.data.frame %>%
    rownames_to_column("variable") %>%
    pivot_longer(-variable,names_to = "species", values_to = "pval")}) %>%
dplyr::filter(!is.na(test)) %>%
pivot_wider(id_cols = "species", names_from = "variable", 
            values_from = c("test","pval"),names_glue = "{variable}_{.value}") %>%
relocate("species", sort(names(.)))
# A tibble: 12 x 7
   species     block_pval block_test `block:treatment_pval` `block:treatment_test` treatment_pval treatment_test
   <chr>            <dbl>      <dbl>                  <dbl>                  <dbl>          <dbl>          <dbl>
 1 Ameira          0.549        4.93                  0.510             6.35               0.0196       1.39e+ 1
 2 Adopsyllus      0.0392      17.6                   0.765             1.30               0.961        9.10e- 2
 3 Ectinosoma      0.549        7.71                  0.784             0.836              0.0196       1.47e+ 1
 4 Ectinosomat     0.235       10.9                   0.784             1.12               0.922        4.09e- 1
 5 Haloschizo      0.725        3.13                  0.922             0.000168           0.490        2.10e+ 0
 6 Lepta.A         0.0196      20.9                   0.196            13.1                0.961        1.46e- 1
 7 Lepta.B         0.549        8.02                  0.471             7.40               0.725        1.36e+ 0
 8 Lepta.C         0.0784      13.8                   0.784             0.340              0.255        5.80e+ 0
 9 Mictyricola     0.549        6.77                  0.510             5.27               0.0392       1.17e+ 1
10 Parevansula     0.549        6.10                  0.706             1.73               1            1.78e-15
11 Quin            0.549        6.02                  0.922             0.000128           0.0784       8.81e+ 0
12 Rhizothrix      0.216       11.6                   0.922             0.00000308         0.137        7.83e+ 0