#                                      Appendix S5

#                           Test function tpaired.krandtest.R

#' Paired t-tests of differences between T1 and T2 for each species
#' 
#' This function computes paired t-tests for each species, for abundances observed 
#' at time 1 (T1) and at time 2 (T2). The test is one-tailed in the direction of the sign 
#' (+ or –) of the t statistic.
#' 
###' The required function "tpaired.randtest.R" is listed in the lower part of this file.
###' It is also available in the R package "adespatial" on CRAN.
#' 
#' Arguments –
#' mat1: site-by-species data at time T1 (data.frame or matrix).
#' mat2: site-by-species data at time T2 (data.frame or matrix).
#' nperm: Number of permutations. Use 999, 9999, or more, to allow for correction of 
#'    p-values for multiple tests.
#' list.all: if FALSE, output only lists t.test results for species with t.stat not 0;
#'    if TRUE, the output lists t.test results for all species; when the t-statistic 
#'    is 0, the p-values in the output table receive the code -999 and Sign(T1-T2) is 0.
#' 
#' Details –
#' The species that do not vary in either data set are discarded before the calculations 
#' of the paired t-tests begin.
#' p-values should be corrected for multiple testing. Use function p.adjust() of {stats}:
#'    p.adjust(res$t.test$p.param)  or  p.adjust(res$t.test$p.perm)
#' Correction methods "holm" (default) and "hochberg" are fine for this type of analysis.
#' 
#' Value –
#' 1. A table with species in rows and 6 columns: 
#'    "mean(T1-T2)","t.stat","p.param","p.perm","p<=0.05","Sign(T1-T2)"
#'    The parametric and permutational p-values are not corrected for multiple tests.
#'    A star is shown in column "p<=0.05" if the parametric p-value is ≤ 0.05.
#' 2. A list of names of the species tested; their t statistics were not 0.
#' 3. A list of names of the species not tested because their t-statistics were 0.
#' 
#' Author – Pierre Legendre \email{pierre.legendre@@umontreal.ca}
#' 
#' References –
#' 
#' Legendre, P. 2019. A temporal beta-diversity index to identify sites that have changed 
#' in exceptional ways in space-time surveys. Ecology and Evolution 2019; 00:1–15.
#' 
#' van den Brink, P. J. & C. J. F. ter Braak. 1999. Principal response curves: analysis of 
#' time-dependent multivariate responses of biological community to stress. Environmental 
#' Toxicology and Chemistry 18: 138–148.
#' 
#' Example –
#' 
#' ### Invertebrate communities subjected to insecticide treatment.
#' 
#' ## As an example in their paper on Principal Response Curves (PRC), van den Brink & ter 
#' Braak (1999) used observations on the abundances of 178 invertebrate species 
#' (macroinvertebrates and zooplankton) subjected to treatments in 12 mesocosms by the 
#' insecticide chlorpyrifos. The mesocosms were sampled at 11 occasions. The data, 
#' available in the {vegan} package, are log-transformed species abundances, 
#' y.tranformed = loge(10*y+1).
#' 
#' ## The data of survey #4 will be compared to those of survey #11 in this example. 
#' Survey #4 was carried out one week after the insecticide treatment, whereas the 
#' fauna of the mesocosms was considered to have fully recovered from the treatment 
#' at the time of survey #11.
#' 
#' require(vegan)
#' data(pyrifos)
#' 
#' ## The mesocosms had originally been attributed at random to the treatments. However,  
#' to facilitate presentation of the results, they will be listed here in order of 
#' increased insecticide doses: {0, 0, 0, 0, 0.1, 0.1, 0.9, 0.9, 6, 6, 44, 44} μg/L.
#' 
#' survey4.order = c(38,39,41,47,37,44,40,46,43,48,42,45)
#' survey11.order = c(122,123,125,131,121,128,124,130,127,132,126,129)
#' 
#' ## Paired t-tests of differences between survey.4 and survey.11 for the p species
#' 
#' res <- tpaired.krandtest(pyrifos[survey4.order,],pyrifos[survey11.order,])
#' 
#' ===============
#' 
tpaired.krandtest <- function(mat1, mat2, nperm=999, list.all=FALSE)
{
epsilon <- .Machine$double.eps
n = nrow(mat1)
p = ncol(mat1)
if(nrow(mat2) != n) stop("Unequal number of rows")
if(ncol(mat2) != p) stop("Unequal number of species")

# Check species names in files mat1 and mat2
sp.diff <- which(colnames(mat1) != colnames(mat2))
if(length(sp.diff)>0) 
cat("Warning – The following species names differ between mat1 & mat2:\n",sp.diff,"\n\n")

# Select the species that have variances > epsilon. Discard the other species.
comm.gr = rbind(mat1, mat2)
tmp = apply(comm.gr, 2, var)
sel.sp = which(tmp > epsilon)
pp <- length(sel.sp)
cat((p-pp),"species were eliminated because they did not vary in the combined data set\n")
cat(pp,"species retained:\n")
mat1.sel = mat1[,sel.sp]
mat2.sel = mat2[,sel.sp]

# 
res = matrix(0,pp,4)
sp.names <- rownames(res) <- colnames(mat1.sel)

#
Star <- NA
Keep <- NA
Reject <- NA
for(j in 1:pp) {
  vec1 = mat1[,j]
  vec2 = mat2[,j]
  	# if(any(abs(vec2-vec1) > 0)) Keep = c(Keep,j) else Reject = c(Reject,j)
	cat(j," ")      # Print ID of species being processed
	tmp <- t.test(vec2, vec1, paired=TRUE)
	if(tmp$estimate == 0) {
		Reject = c(Reject,j) 
		res[j,3:4] = c(-999,-999)
		Star <- c(Star, " ")
		} else {
		Keep = c(Keep,j) 
		if(sign(tmp$estimate)<0) tail="less" else if(sign(tmp$estimate)>=0) tail="greater"
		t.res <- tpaired.randtest(vec2, vec1, nperm=nperm, alternative=tail, silent=TRUE)
		res[j,1:4] <- c(t.res$estim, t.res$t.ref, t.res$p.param, t.res$p.perm)
		Star <- c(Star, ifelse(t.res$p.perm>0.05, " ","*"))
		}
  }
cat("\n\n")
Keep <- Keep[-1]
Reject <- Reject[-1]
tmp <- length(Reject)
if(tmp>0) cat(tmp,"species not tested because t.stat = 0. See 'No_test' output list\n\n")
res = data.frame(res,Star[-1],sign(res[,1]))
colnames(res) <- c("mean(T2-T1)","t.stat","p.param","p.perm","p<=0.05","Sign(T2-T1)")
if(list.all) 
	out <- list(t.tests=res, Tested=sp.names[Keep], No_test=sp.names[Reject])
	else out <- list(t.tests=res[Keep,],Tested=sp.names[Keep],No_test=sp.names[Reject])
out
}

#'                              Permutational paired t-test
#' 
#' This function computes a permutation test of comparison of the means
#' of two paired vectors (related samples). For each object, permutations
#' are restricted to the two related observations.
#'
#' Arguments of the function –
#'    vec1, vec2: the two vectors to be compared
#'    nperm = number of permutations (default value: 999)
#'    alternative = c("two.sided", "less", "greater"). Default value: "two.sided"
#'    silent = FALSE: calculation results are printed to the R console.
#'           = TRUE : calculation results are not printed to R console (for simulations).
#'
#' Values returned – A list containing the following results:
#'    estim : mean of the differences
#'    t.ref : reference value of the t-statistic
#'    p.param : parametric p-value
#'    p.perm : permutational p-value
#'    nperm : number of permutations
#'
#' Authors – Pierre Legendre \email{pierre.legendre@@umontreal.ca}
#'   Permutation code improved by Guillaume Blanchet. 
#' 
#' Reference –
#' Zar, J. H. 1999. Biostatistical analysis. 4th edition. Prentice Hall, New Jersey.
#'
#' Example –
#' ## Deer leg length data from Zar (1999, p. 162)
#'
#' deer = matrix( c(142,140,144,144,142,146,149,150,142,148,138,136,147,139,143,141,143, 
#'    145,136,146),10,2)
#' rownames(deer) = paste("Deer",1:10,sep=".")
#' colnames(deer) = c('Hind.leg', 'Fore.leg')
#'
#' res = tpaired.randtest(deer[,1], deer[,2])   # Two-tailed test by default
#'
#' ## Compare the results to:  res2 = t.test(deer[,1], deer[,2], paired=TRUE) 
#'
#' ===============
#'
tpaired.randtest <- function(vec1, vec2, nperm=999, alternative="two.sided", silent=FALSE)
{
n1 <- length(vec1)
n2 <- length(vec2)
if(n1 != n2) stop("The two vectors have different lengths. They cannot be paired.")

alt <- match.arg(alternative, c("two.sided", "less", "greater"))

res = t.test(vec1, vec2, paired=TRUE, alternative=alt)
t.ref =  res$statistic

# Print these first results
if(!silent) { 
	cat('\nt-test comparing the means of two related samples','\n','\n')
	cat('Number of objects:',n1,'\n')
	cat('Mean of the differences:',res$estimate,'\n')
	cat('t statistic (paired observations):',t.ref,'\n')
	cat('95 percent confidence interval of t:',res$conf.int,'\n')
	cat('Degrees of freedom:',res$parameter,'\n')
	cat('Alternative hypothesis:',alt,'\n')
	cat('Prob (parametric):',res$p.value,'\n')
	}

# Perform the permutation test
# Permutations are restricted to the two related observations for each object.
nPGE <- 1

for(i in 1:nperm)
	{
	mat <- cbind(vec1,vec2)
	topermute <- rbinom(n1,1,0.5)
	mat[topermute==1,] <- mat[topermute==1,2:1]
	
	res.perm = t.test(mat[,1], mat[,2], paired=TRUE, alternative=alt)
	t.perm = res.perm$statistic
	
	if(alt == "two.sided") if( abs(t.perm) >= abs(t.ref) ) nPGE <- nPGE+1
	if(alt == "less")      if(t.perm <= t.ref) nPGE <- nPGE+1
	if(alt == "greater")   if(t.perm >= t.ref) nPGE <- nPGE+1
	}

# Compute and print the permutational p-value
P <- nPGE/(nperm+1)
if(!silent) cat('Prob (',nperm,'permutations):', formatC(P,digits=5,width=7,format="f"),'\n')
#
return(list(estim=res$est[[1]], t.ref=t.ref, p.param=res$p.value, p.perm=P, nperm=nperm))
}
