# Permanova with several taxa (loop that analyses the whole data table, generates output table )
# Ralf Schwamborn, 2017
# Data format: First column = Factor, other columns = taxa
# setwd()
# Example: Two Lakes
# Lake A, Lake B, 20 samples each
A <- rep ("A", 20)
B <- rep ("B", 20)
LAKE <- as.factor(c(A,B))
Fish_dens <- rnorm(mean = 20, sd = 2, 40)
Zoop_dens <- rnorm(mean = 200, sd = 2, 40)
Pytop_dens <- c((rnorm(mean = 1000, sd = 2, 20)), (rnorm(mean = 4000, sd = 2, 20)))
Lake_data <- data.frame(LAKE,Fish_dens,Zoop_dens,Pytop_dens)
attach(Lake_data)
## The following objects are masked _by_ .GlobalEnv:
##
## Fish_dens, LAKE, Pytop_dens, Zoop_dens
# Run Permutation test
library(coin)
## Loading required package: survival
independence_test(Fish_dens~ LAKE)
##
## Asymptotic General Independence Test
##
## data: Fish_dens by LAKE (A, B)
## Z = -0.87968, p-value = 0.379
## alternative hypothesis: two.sided
independence_test(Zoop_dens~ LAKE)
##
## Asymptotic General Independence Test
##
## data: Zoop_dens by LAKE (A, B)
## Z = -1.1084, p-value = 0.2677
## alternative hypothesis: two.sided
independence_test(Pytop_dens~ LAKE)
##
## Asymptotic General Independence Test
##
## data: Pytop_dens by LAKE (A, B)
## Z = -6.245, p-value = 4.238e-10
## alternative hypothesis: two.sided
#loop that reads the table performs a Permanova and generates output table with "p" values
input.table <- Lake_data
independence_test(input.table[,1]~ input.table[,2])
##
## Asymptotic General Independence Test
##
## data: input.table[, 1] by input.table[, 2]
## Z = -0.87968, p-value = 0.379
## alternative hypothesis: two.sided
independence_test(input.table[,1]~ input.table[,3])
##
## Asymptotic General Independence Test
##
## data: input.table[, 1] by input.table[, 3]
## Z = -1.1084, p-value = 0.2677
## alternative hypothesis: two.sided
independence_test(input.table[,1]~ input.table[,4])
##
## Asymptotic General Independence Test
##
## data: input.table[, 1] by input.table[, 4]
## Z = -6.245, p-value = 4.238e-10
## alternative hypothesis: two.sided
# Obtain the "p" value (tanks to Marcin Kosinski, see: https://stackoverflow.com/questions/32639460/how-to-extract-value-from-one-s4-class)
res.perm <- independence_test(input.table[,1]~ input.table[,4])
str(res.perm)
## Formal class 'ScalarIndependenceTest' [package "coin"] with 7 slots
## ..@ parameter : chr "mu"
## ..@ nullvalue : num(0)
## ..@ distribution:Formal class 'AsymptNullDistribution' [package "coin"] with 10 slots
## .. .. ..@ seed : int NA
## .. .. ..@ q :function (p)
## .. .. ..@ d :function (x)
## .. .. ..@ support :function ()
## .. .. ..@ parameters : list()
## .. .. ..@ pvalue :function (q)
## .. .. ..@ midpvalue :function (q)
## .. .. ..@ pvalueinterval:function (q)
## .. .. ..@ p :function (q)
## .. .. ..@ name : chr "Univariate Normal Distribution"
## ..@ statistic :Formal class 'ScalarIndependenceTestStatistic' [package "coin"] with 15 slots
## .. .. ..@ alternative : chr "two.sided"
## .. .. ..@ paired : logi FALSE
## .. .. ..@ teststatistic : Named num -6.24
## .. .. .. ..- attr(*, "names")= chr "A"
## .. .. ..@ standardizedlinearstatistic: Named num -6.24
## .. .. .. ..- attr(*, "names")= chr "A"
## .. .. ..@ linearstatistic : num 20006
## .. .. ..@ expectation : Named num 50009
## .. .. .. ..- attr(*, "names")= chr "A"
## .. .. ..@ covariance :Formal class 'Variance' [package "coin"] with 1 slot
## .. .. .. .. ..@ variance: Named num 23081996
## .. .. .. .. .. ..- attr(*, "names")= chr "A"
## .. .. ..@ xtrans : num [1:40, 1] 1001 1001 998 998 999 ...
## .. .. .. ..- attr(*, "assign")= int 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : NULL
## .. .. .. .. ..$ : chr ""
## .. .. ..@ ytrans : num [1:40, 1] 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:40] "1" "2" "3" "4" ...
## .. .. .. .. ..$ : chr "A"
## .. .. .. ..- attr(*, "assign")= int 1
## .. .. ..@ xtrafo :function (data, numeric_trafo = id_trafo, factor_trafo = f_trafo,
## ordered_trafo = of_trafo, surv_trafo = logrank_trafo, var_trafo = NULL,
## block = NULL)
## .. .. ..@ ytrafo :function (data, numeric_trafo = id_trafo, factor_trafo = f_trafo,
## ordered_trafo = of_trafo, surv_trafo = logrank_trafo, var_trafo = NULL,
## block = NULL)
## .. .. ..@ x :'data.frame': 40 obs. of 1 variable:
## .. .. .. ..$ input.table[, 4]: num [1:40] 1001 1001 998 998 999 ...
## .. .. .. ..- attr(*, "terms")=Classes 'terms', 'formula' language ~input.table[, 4]
## .. .. .. .. .. ..- attr(*, "variables")= language list(input.table[, 4])
## .. .. .. .. .. ..- attr(*, "factors")= int [1, 1] 1
## .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. ..$ : chr "input.table[, 4]"
## .. .. .. .. .. .. .. ..$ : chr "input.table[, 4]"
## .. .. .. .. .. ..- attr(*, "term.labels")= chr "input.table[, 4]"
## .. .. .. .. .. ..- attr(*, "order")= int 1
## .. .. .. .. .. ..- attr(*, "intercept")= int 1
## .. .. .. .. .. ..- attr(*, "response")= int 0
## .. .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. .. .. .. ..- attr(*, "predvars")= language list(input.table[, 4])
## .. .. .. .. .. ..- attr(*, "dataClasses")= Named chr "numeric"
## .. .. .. .. .. .. ..- attr(*, "names")= chr "input.table[, 4]"
## .. .. ..@ y :'data.frame': 40 obs. of 1 variable:
## .. .. .. ..$ input.table[, 1]: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..- attr(*, "terms")=Classes 'terms', 'formula' language ~input.table[, 1]
## .. .. .. .. .. ..- attr(*, "variables")= language list(input.table[, 1])
## .. .. .. .. .. ..- attr(*, "factors")= int [1, 1] 1
## .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. ..$ : chr "input.table[, 1]"
## .. .. .. .. .. .. .. ..$ : chr "input.table[, 1]"
## .. .. .. .. .. ..- attr(*, "term.labels")= chr "input.table[, 1]"
## .. .. .. .. .. ..- attr(*, "order")= int 1
## .. .. .. .. .. ..- attr(*, "intercept")= int 1
## .. .. .. .. .. ..- attr(*, "response")= int 0
## .. .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. .. .. .. ..- attr(*, "predvars")= language list(input.table[, 1])
## .. .. .. .. .. ..- attr(*, "dataClasses")= Named chr "factor"
## .. .. .. .. .. .. ..- attr(*, "names")= chr "input.table[, 1]"
## .. .. ..@ block : Factor w/ 1 level "0": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. ..@ weights : num [1:40] 1 1 1 1 1 1 1 1 1 1 ...
## ..@ estimates : list()
## ..@ method : chr "General Independence Test"
## ..@ call : language independence_test.IndependenceProblem(object = <S4 object of class "IndependenceProblem">)
res.perm@distribution@pvalue(res.perm@statistic@teststatistic)
## [1] 4.238155e-10
p_value <- res.perm@distribution@pvalue(res.perm@statistic@teststatistic)
# As Loop
taxa.names <- c(names(input.table[2:(ncol(input.table))]))
outp.table <- data.frame(taxa.names, (rep(0, length(taxa.names)) ))
names(outp.table) <- c("Taxon", "p")
for (i in 2 : (length(taxa.names)+1) ) {
res.perm <- independence_test(input.table[,1]~ input.table[,i])
p_value <- res.perm@distribution@pvalue(res.perm@statistic@teststatistic)
p_value <- round(p_value, 6)
outp.table[(i-1),2] <- p_value
}
outp.table
## Taxon p
## 1 Fish_dens 0.379035
## 2 Zoop_dens 0.267700
## 3 Pytop_dens 0.000000
# As function (with loop inside)
Permanova_table.fun1 <- function(input.tab) {
input.table <- input.tab
taxa.names <- c(names(input.table[2:(ncol(input.table))]))
outp.table <- data.frame(taxa.names, (rep(0, length(taxa.names)) ))
names(outp.table) <- c("Taxon", "p")
for (i in 2 : (length(taxa.names)+1) ) {
res.perm <- independence_test(input.table[,1]~ input.table[,i])
p_value <- res.perm@distribution@pvalue(res.perm@statistic@teststatistic)
p_value <- round(p_value, 6)
outp.table[(i-1),2] <- p_value
}
; outp.table
}
# Test the function
Permanova_table.fun1(Lake_data)
## Taxon p
## 1 Fish_dens 0.379035
## 2 Zoop_dens 0.267700
## 3 Pytop_dens 0.000000