提取 `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
我想从 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