33# # adjusting pseudo-F values. Also calculates multivariate dispersion for
44# # categorical variables.
55
6+ require(vegan )
7+ require(tibble )
8+
69# First, a quick function to calculate omega squared based on:
710 # https://academic.oup.com/bioinformatics/article/31/15/2461/188732
811 # http://www.real-statistics.com/multiple-regression/other-measures-effect-size-anova/
@@ -42,7 +45,7 @@ PERMANOVA.omega2 <- function(adonis2.object, num.control.vars) {
4245 # # by.adonsi2 is to pass through for the by = argument in adonis2
4346
4447
45- group.PERMANOVA <- function (var.names , var.table , var.table.c , control.vars = " " , species.table , species.table.c , num.control.vars = 0 , by.adonis2 = " terms" , perms = 99999 , method = " bray" ) {
48+ group.PERMANOVA <- function (var.names , var.table , var.table.c , control.vars = " " , species.table , species.table.c , num.control.vars = 0 , by.adonis2 = " terms" , AIC.type = " AICc " , perms = 99999 , method = " bray" ) {
4649
4750 # # need to order columns otherwise the order of the columns and the order of col.numbers is incorrect.
4851
@@ -51,14 +54,16 @@ group.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars = ""
5154
5255
5356
54- output <- data.frame ( row .names = var.names ,
57+ output <- tibble( var .names = var.names ,
5558 var.explnd = rep(NA , times = length(var.names )),
5659 avg.var.explnd = rep(NA ),
5760 pseudo.F = rep(NA ),
5861 F .pval = rep(NA ),
5962 adj.F.pval = rep(NA ),
6063 disp.F.pval = rep(NA ),
61- adj.disp.F.pval = rep(NA )
64+ adj.disp.F.pval = rep(NA ),
65+ omega2 = rep(NA ),
66+ AIC.stat = rep(NA )
6267 )
6368 col.numbers <- which(colnames(var.table ) %in% var.names )
6469
@@ -93,7 +98,8 @@ group.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars = ""
9398 output $ omega2 [i ] <- PERMANOVA.omega2(temp , num.control.vars )
9499 output $ pseudo.F [i ] <- temp $ F [num.control.vars + 1 ]
95100 output $ F .pval [i ] <- temp $ `Pr(>F)` [num.control.vars + 1 ]
96-
101+ output $ AIC.stat [i ] <- AICc.PERMANOVA2(temp )[[AIC.type ]]
102+
97103
98104
99105 # the ordering is so that col.names[i] works!
@@ -128,6 +134,8 @@ group.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars = ""
128134 p.adjust(p = output $ F .pval , method = " holm" )
129135 output $ adj.disp.F.pval <-
130136 p.adjust(p = output $ disp.F.pval , method = " holm" )
137+ output $ AIC.delta <-
138+ output $ AIC.stat - min(output $ AIC.stat )
131139
132140return (output )
133141
@@ -136,7 +144,7 @@ return(output)
136144}
137145
138146
139- group.univ.PERMANOVA <- function (var.names , var.table , var.table.c , control.vars = " " , species.vector , species.vector.c , num.control.vars = 0 , by.adonis2 = " terms" , perms = 99999 , method = " euclid" ) {
147+ group.univ.PERMANOVA <- function (var.names , var.table , var.table.c , control.vars = " " , species.vector , species.vector.c , num.control.vars = 0 , by.adonis2 = " terms" , perms = 99999 , method = " euclid" , AIC.type = " AICc " ) {
140148
141149 # # need to order columns otherwise the order of the columns and the order of col.numbers is incorrect.
142150
@@ -145,13 +153,14 @@ group.univ.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars
145153
146154
147155
148- output <- data.frame ( row .names = var.names ,
156+ output <- tibble( var .names = var.names ,
149157 var.explnd = rep(NA , times = length(var.names )),
150158 avg.var.explnd = rep(NA ),
151159 pseudo.F = rep(NA ),
152160 F .pval = rep(NA ),
153161 adj.F.pval = rep(NA ),
154- AICc = rep(NA )
162+ omega2 = rep(NA ),
163+ AIC.stat = rep(NA )
155164 )
156165 col.numbers <- which(colnames(var.table ) %in% var.names )
157166
@@ -188,7 +197,7 @@ group.univ.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars
188197 output $ omega2 [i ] <- PERMANOVA.omega2(temp , num.control.vars )
189198 output $ pseudo.F [i ] <- temp $ F [num.control.vars + 1 ]
190199 output $ F .pval [i ] <- temp $ `Pr(>F)` [num.control.vars + 1 ]
191- output $ AICc [i ] <- AICc.PERMANOVA2(temp )[3 ]
200+ output $ AIC.stat [i ] <- AICc.PERMANOVA2(temp )[[ AIC.type ] ]
192201
193202
194203
@@ -198,6 +207,8 @@ group.univ.PERMANOVA <- function(var.names, var.table, var.table.c, control.vars
198207 }
199208
200209 output $ adj.F.pval <- p.adjust(p = output $ F .pval , method = " holm" )
210+ output $ delta.aic <- output $ AIC.stat - min(output $ AIC.stat )
211+
201212
202213 return (output )
203214
0 commit comments