############################################################
#Functions for generating different correlation structures
#Routines from HARDIN et al., 2011
###########################################################
simcor = function (k = 6, size = c(10, 10, 10, 10, 10, 10), rho = c(0.7,0.7, 0.5, 0.9, 0.85, 0.4), delta = 0.39, epsilon = 0.99 -max(rho), eidim = 2) 
{
  ndim <- sum(size)
  bigcor <- matrix(rep(delta, ndim * ndim), ncol = ndim)
  for (i in 1:k) {
    cor <- matrix(rep(rho[i], size[i] * size[i]), ncol = size[i])
    if (i == 1) {bigcor[1:size[1], 1:size[1]] <- cor}
    if (i != 1) {bigcor[(sum(size[1:(i - 1)]) + 1):sum(size[1:i]), 
                        (sum(size[1:(i - 1)]) + 1):sum(size[1:i])] <- cor}
  }
  diag(bigcor) <- 1 - epsilon
  
  eivect <- c()
  for (i in 1:ndim) {
    ei <- runif(eidim, -0.01, 0.01)
    eivect <- cbind(eivect, sqrt(epsilon) * ei/sqrt(sum(ei^2)))
  }
  bigE <- t(eivect) %*% eivect
  cor.nz <- bigcor + bigE
  cor.nz
}

################################
##  Simulating the Toeplitz Matrix
################################

# this function calculates an AR(1) Toeplitz matrix
# with a block diagonal structure.  The size and base correlation for each
# block is user specified


# k is the number of groups
# size is a vector of length k specifying the size of each group 
# rho is a vector of length k specifying base correlation values
# epsilon <- (1-max(rho))/(1+max(rho)) - .01 
# eidim is the space from which the noise is generated, the smaller the more noise

simcorTop <- function(k=6,size=c(10,5,8,2,15,50),rho=c(.7,.7,.5,.9,.85,.4),epsilon= .01, eidim=2) {
  
  
  ndim <- sum(size)# dim of correlation matrix
  bigcor<- matrix(rep(0, ndim*ndim), ncol=ndim)
  
  ### generating the basic correlation matrix
  
  
  for (i in 1:k){
    
    top <- c(1,rho[i]^(seq(1:(size[i]-1))))
    cor <- toeplitz(top)
    
    if (i==1){bigcor[1:size[1], 1:size[1]] <- cor}
    if (i!=1){bigcor[(sum(size[1:(i-1)]) + 1):sum(size[1:i]),
                     (sum(size[1:(i-1)]) + 1):sum(size[1:i])] <- cor}
  }
  
  diag(bigcor) <- 1 - epsilon
  
  
  ### adding noise to the correlation matrix
  
  eivect <- c( )
  for (i in 1:ndim) {
    ei <- runif(eidim, -1, 1)
    eivect <- cbind(eivect, sqrt(epsilon) * ei / sqrt(sum(ei^2) ) )
  }
  
  
  bigE <- t(eivect) %*% eivect
  cor.nz <- bigcor + bigE
  cor.nz
  
}
####################################

simcor.H <- function(k=6, size=c(10,5,8,7,15,50), 
                     rho=rbind(c(.9,.7), c(.7,.7), c(.7,.2), c(.5,.3), c(.9,.85), c(.3,.2)), power=1,
                     epsilon=.08, eidim=2){
  
  
  ndim <- sum(size)# dim of correlation matrix
  bigcor<- matrix(rep(0, ndim*ndim), ncol=ndim)
  
  ### generating the basic correlation matrix
  
  
  for (i in 1:(k) ){
    
    cor <- toeplitz(rho.func(rho[i,1],rho[i,2],power,size[i]) )
    
    if (i==1){bigcor[1:size[1], 1:size[1]] <- cor}
    if (i!=1){bigcor[(sum(size[1:(i-1)]) + 1):sum(size[1:i]),
                     (sum(size[1:(i-1)]) + 1):sum(size[1:i])] <- cor}
  }
  diag(bigcor) <- 1 - epsilon
  
  
  ### adding noise to the correlation matrix
  
  eivect <- c( )
  for (i in 1:ndim) {
    ei <- runif(eidim, -1, 1)
    eivect <- cbind(eivect, sqrt(epsilon) * ei / sqrt(sum(ei^2) ) )
  }
  
  
  bigE <- t(eivect) %*% eivect
  cor.nz <- bigcor + bigE
  cor.nz
}

rho.func <- function(r.max, r.min, power,p){
  rhovec <-c()
  
  rhovec[1] <- 1
  for(i in 2:p){
    rhovec[i] <- r.max - ((i-2)/(p-2))^power*(r.max-r.min)
  }
  rhovec}


################################################
###############################################
library(corrplot)
library(vegan)
library(ggplot2)
library(sciplot)
library(agricolae)
library(vegan)
library(gtable)
library(gridExtra)


#Loading the data tables
CD<-read.table(file="CDmolepercent.txt",header=TRUE,row.names=1) # PLFA and Soil fertility used in Lisboa et al 2014b

plfa<-CD[,1:20] # Microbial structure (PLFA). The Y data table

fert<-CD[,31:45] # Soil fertility data table. The X data table

fert.m<-as.matrix(fert)

plfa.m<-as.matrix(plfa)

LU<-CD$LandUse #Factor (4 levels: pasture, forest fragment, soybean crop, integrated ecosystem)  that is going to be used for evaluating the effects of the increasing correlation on ANOVA results using PAM.

p1<-corrplot(cor(fert.m))
p2<-corrplot(cor(plfa.m))


#Simulating soil fertility data tables with increasing number of columns covarying at different correlation levels (0.9, 0.7, 0.5, and 0.2)
#*******************************************************************************************************************************


#NOTICE! Here we just show the 13 Pre-Transformations 
# applied to X (Soil fertility) and Y (PLFA) of data tables before the Procrustean analysis.
# This is because for each number of columns covarying in the Soil fertility data table (06, 09, 12, and 15)
# four correlation levels were explored (0.9, 0.7, 0.5, and 0.2). 
# Thus:

# 06 columns covarying in the Soil fertility data table
#*******************************************************
# 13 Procrustes relationships for correlation of 0.9
# 13 Procrustes relationships for correlation of 0.7
# 13 Procrustes relationships for correlation of 0.5
# 13 Procrustes relationships for correlation of 0.2

# 09 columns covarying in the Soil fertility data table
#*******************************************************
# 13 Procrustes relationships for correlation of 0.9
# 13 Procrustes relationships for correlation of 0.7
# 13 Procrustes relationships for correlation of 0.5
# 13 Procrustes relationships for correlation of 0.2

# 12 columns covarying in the Soil fertility data table
#*******************************************************
# 13 Procrustes relationships for correlation of 0.9
# 13 Procrustes relationships for correlation of 0.7
# 13 Procrustes relationships for correlation of 0.5
# 13 Procrustes relationships for correlation of 0.2

# 15 columns covarying in the Soil fertility data table
#*******************************************************
# 13 Procrustes relationships for correlation of 0.9
# 13 Procrustes relationships for correlation of 0.7
# 13 Procrustes relationships for correlation of 0.5
# 13 Procrustes relationships for correlation of 0.2


#setting simulation objects:

procr1<-array(NA,dim=c(100,1)) # Generates an array with 100 simulated Procrustes relationships
pprocr1<-array(NA,dim=c(100,1)) # generates ans array with 100 simulated Procrustes significance tests (permutations-999)
p1.PAM<-rep(NA,100) # generates 100 simulates Procrustean association metrics (PAMs), i.e, residual vectors

procr2<-array(NA,dim=c(100,1))
pprocr2<-array(NA,dim=c(100,1))
p2.PAM<-rep(NA,100)

procr3<-array(NA,dim=c(100,1))
pprocr3<-array(NA,dim=c(100,1))
p3.PAM<-rep(NA,100) 


procr4<-array(NA,dim=c(100,1))
pprocr4<-array(NA,dim=c(100,1))
p4.PAM<-rep(NA,100)

procr5<-array(NA,dim=c(100,1))
pprocr5<-array(NA,dim=c(100,1))
p5.PAM<-rep(NA,100)


procr6<-array(NA,dim=c(100,1))
pprocr6<-array(NA,dim=c(100,1))
p6.PAM<-rep(NA,100)

procr7<-array(NA,dim=c(100,1))
pprocr7<-array(NA,dim=c(100,1))
p7.PAM<-rep(NA,100)

procr8<-array(NA,dim=c(100,1))
pprocr8<-array(NA,dim=c(100,1))
p8.PAM<-rep(NA,100)

procr9<-array(NA,dim=c(100,1))
pprocr9<-array(NA,dim=c(100,1))
p9.PAM<-rep(NA,100)


procr10<-array(NA,dim=c(100,1))
pprocr10<-array(NA,dim=c(100,1))
p10.PAM<-rep(NA,100)

procr11<-array(NA,dim=c(100,1))
pprocr11<-array(NA,dim=c(100,1))
p11.PAM<-rep(NA,100)


procr12<-array(NA,dim=c(100,1))
pprocr12<-array(NA,dim=c(100,1))
p12.PAM<-rep(NA,100)


procr13<-array(NA,dim=c(100,1))
pprocr13<-array(NA,dim=c(100,1))
p13.PAM<-rep(NA,100)

#SF09

procr14<-array(NA,dim=c(100,1)) 
pprocr14<-array(NA,dim=c(100,1)) 
p14.PAM<-rep(NA,100) 

procr15<-array(NA,dim=c(100,1))
pprocr15<-array(NA,dim=c(100,1))
p15.PAM<-rep(NA,100)

procr16<-array(NA,dim=c(100,1))
pprocr16<-array(NA,dim=c(100,1))
p16.PAM<-rep(NA,100) 


procr17<-array(NA,dim=c(100,1))
pprocr17<-array(NA,dim=c(100,1))
p17.PAM<-rep(NA,100)

procr18<-array(NA,dim=c(100,1))
pprocr18<-array(NA,dim=c(100,1))
p18.PAM<-rep(NA,100)


procr19<-array(NA,dim=c(100,1))
pprocr19<-array(NA,dim=c(100,1))
p19.PAM<-rep(NA,100)

procr20<-array(NA,dim=c(100,1))
pprocr20<-array(NA,dim=c(100,1))
p20.PAM<-rep(NA,100)

procr21<-array(NA,dim=c(100,1))
pprocr21<-array(NA,dim=c(100,1))
p21.PAM<-rep(NA,100)

procr22<-array(NA,dim=c(100,1))
pprocr22<-array(NA,dim=c(100,1))
p22.PAM<-rep(NA,100)


procr23<-array(NA,dim=c(100,1))
pprocr23<-array(NA,dim=c(100,1))
p23.PAM<-rep(NA,100)

procr24<-array(NA,dim=c(100,1))
pprocr24<-array(NA,dim=c(100,1))
p24.PAM<-rep(NA,100)


procr25<-array(NA,dim=c(100,1))
pprocr25<-array(NA,dim=c(100,1))
p25.PAM<-rep(NA,100)


procr26<-array(NA,dim=c(100,1))
pprocr26<-array(NA,dim=c(100,1))
p26.PAM<-rep(NA,100)

#SF12
procr27<-array(NA,dim=c(100,1)) # Generates an array with 100 simulated Procrustes relationships
pprocr27<-array(NA,dim=c(100,1)) # generates ans array with 100 simulated Procrustes significance tests (permutations-999)
p27.PAM<-rep(NA,100) # generates 100 simulates Procrustean association metrics (PAMs), i.e, residual vectors

procr28<-array(NA,dim=c(100,1))
pprocr28<-array(NA,dim=c(100,1))
p28.PAM<-rep(NA,100)

procr29<-array(NA,dim=c(100,1))
pprocr29<-array(NA,dim=c(100,1))
p29.PAM<-rep(NA,100) 


procr30<-array(NA,dim=c(100,1))
pprocr30<-array(NA,dim=c(100,1))
p30.PAM<-rep(NA,100)

procr31<-array(NA,dim=c(100,1))
pprocr31<-array(NA,dim=c(100,1))
p31.PAM<-rep(NA,100)


procr32<-array(NA,dim=c(100,1))
pprocr32<-array(NA,dim=c(100,1))
p32.PAM<-rep(NA,100)

procr33<-array(NA,dim=c(100,1))
pprocr33<-array(NA,dim=c(100,1))
p33.PAM<-rep(NA,100)

procr34<-array(NA,dim=c(100,1))
pprocr34<-array(NA,dim=c(100,1))
p34.PAM<-rep(NA,100)

procr35<-array(NA,dim=c(100,1))
pprocr35<-array(NA,dim=c(100,1))
p35.PAM<-rep(NA,100)


procr36<-array(NA,dim=c(100,1))
pprocr36<-array(NA,dim=c(100,1))
p36.PAM<-rep(NA,100)

procr37<-array(NA,dim=c(100,1))
pprocr37<-array(NA,dim=c(100,1))
p37.PAM<-rep(NA,100)


procr38<-array(NA,dim=c(100,1))
pprocr38<-array(NA,dim=c(100,1))
p38.PAM<-rep(NA,100)


procr39<-array(NA,dim=c(100,1))
pprocr39<-array(NA,dim=c(100,1))
p39.PAM<-rep(NA,100)

procr40<-array(NA,dim=c(100,1)) # Generates an array with 100 simulated Procrustes relationships
pprocr40<-array(NA,dim=c(100,1)) # generates ans array with 100 simulated Procrustes significance tests (permutations-999)
p40.PAM<-rep(NA,100) # generates 100 simulates Procrustean association metrics (PAMs), i.e, residual vectors

procr41<-array(NA,dim=c(100,1))
pprocr41<-array(NA,dim=c(100,1))
p41.PAM<-rep(NA,100)

procr42<-array(NA,dim=c(100,1))
pprocr42<-array(NA,dim=c(100,1))
p42.PAM<-rep(NA,100) 


procr43<-array(NA,dim=c(100,1))
pprocr43<-array(NA,dim=c(100,1))
p43.PAM<-rep(NA,100)

procr44<-array(NA,dim=c(100,1))
pprocr44<-array(NA,dim=c(100,1))
p44.PAM<-rep(NA,100)


procr45<-array(NA,dim=c(100,1))
pprocr45<-array(NA,dim=c(100,1))
p45.PAM<-rep(NA,100)

procr46<-array(NA,dim=c(100,1))
pprocr46<-array(NA,dim=c(100,1))
p46.PAM<-rep(NA,100)

procr47<-array(NA,dim=c(100,1))
pprocr47<-array(NA,dim=c(100,1))
p47.PAM<-rep(NA,100)

procr48<-array(NA,dim=c(100,1))
pprocr48<-array(NA,dim=c(100,1))
p48.PAM<-rep(NA,100)


procr49<-array(NA,dim=c(100,1))
pprocr49<-array(NA,dim=c(100,1))
p49.PAM<-rep(NA,100)

procr50<-array(NA,dim=c(100,1))
pprocr50<-array(NA,dim=c(100,1))
p50.PAM<-rep(NA,100)


procr51<-array(NA,dim=c(100,1))
pprocr51<-array(NA,dim=c(100,1))
p51.PAM<-rep(NA,100)


procr52<-array(NA,dim=c(100,1))
pprocr52<-array(NA,dim=c(100,1))
p52.PAM<-rep(NA,100)




for(times in 1:100){
  
  # this part of simulation is going to create 100 different data frames with a correlation level almost zero.
  # which are going to be used latter
  #******************************************
  fert.unc<-t((solve(t(chol(cov(fert.m)))))%*%t(fert.m))
  fert.perm<-fert.unc[,sample(15,replace=F)]
  
  Mcol<-simcor(k=1,size=c(15),rho=c(0.05),epsilon=0.05-0.04) #   Correlation < 0.1
  Total<-t(chol(Mcol))%*% t(fert.perm)
  M.df<-as.data.frame(t(Total))
  M.df
  
  
 
  # 100 simulations incorporating a specified correlation level (0.9, 0.7, 0.5, and 0.2)
  # within the whole Soil fertility data table (n=53, p=15)
  
  SF.cor<-simcor(k=1,size=c(15),rho=c(0.20),epsilon=0.21-0.20) # Here a correlation of 0.9 has been  incorporated "rho=c(0.90),epsilon=0.91-0.90"
  SF.sim<-t(chol(SF.cor))%*% t(fert.perm)
  SF.t<-t(SF.sim)
  SF.df<-data.frame(SF.t)
  SF.df
  
  #Increasing the correlation level within the Soil fertility data table by binding at a specific level
  #columns from the data frame with correlation of almost 0 and columns from the soil fertility data frame with a specified correlation level
  
  
  # Creating a Soil Fertility (SF) data table having only 06 columns covarying at a specific level
  SF6cor<-cbind(SF.df[,1:6],M.df[,1:9]) # Only 6 columns in Soil Fertility matrix are going to be related at the correlation level specified previously
  
  # Creating a Soil Fertility (SF) data table having only 09 columns covarying at a specific level
  SF9cor<-cbind(SF.df[,1:9],M.df[,1:6]) # Only 9 columns in Soil Fertility matrix are going to be related at the correlation level specified previously
  
  # Creating a Soil Fertility (SF) data table having only 12 columns (15) covarying at a specific level
  SF12cor<-cbind(SF.df[,1:12],M.df[,1:3]) # Only 9 columns in Soil Fertility matrix are going to be related at the correlation level specified previously
  
  # Creating a Soil Fertility (SF) data table having all  columns covarying
  SF15cor<-SF.df # All columns in Soil Fertility matrix are going to be related at the correlation level specified previously
  
 
  
  #Procrustes analysis
  fert.m<-as.matrix(fert) # Soil fertility data table (X)
  plfa.m<-as.matrix(plfa) # Soil microbial structure (PLFA) data table (Y)
  
  #1th Pre-transformation:
  #simulatiing with the original covariance structure for both SF (fert.m) and FA (plfa.m)
  
  procr1[times,1]<-protest((SF6cor),plfa.m)$sig
  pprocr1[times,1]<-protest((SF6cor),plfa.m)$t0
  p1.PAM[times]<-list(residuals(protest((SF6cor^2),plfa.m)))
  
  #2th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a distance matrix formed
  # based on Bray-Curtis distance of square-rooted data
  
  procr2[times,1]<-protest(vegdist(sqrt(SF6cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$sig
  pprocr2[times,1]<-protest(vegdist(sqrt(SF6cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$t0
  p2.PAM[times]<-list(residuals(protest(vegdist(sqrt(SF6cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))))
  
  #3th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed  and a distance matrix formed
  # based on Bray-Curtis distance of Hellinger-transformed data
  
  procr3[times,1]<-protest(vegdist(decostand((SF6cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$sig
  pprocr3[times,1]<-protest(vegdist(decostand((SF6cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$t0
  p3.PAM[times]<-list(residuals(protest(vegdist(decostand((SF6cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))))
  
  
  #4th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a  distance matrix formed
  # based on Euclidean distance of log(X+1) transformed data
  
  procr4[times,1]<-protest(vegdist(log((SF6cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$sig
  pprocr4[times,1]<-protest(vegdist(log((SF6cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$t0
  p4.PAM[times]<-list(residuals(protest(vegdist(log((SF6cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))))
  
  
  
  #5th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) after both data tables  transformed and a  distance matrix formed
  # based on Euclidean distance of Hellinger - transformed data
  
  procr5[times,1]<-protest(vegdist(decostand((SF6cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$sig
  pprocr5[times,1]<-protest(vegdist(decostand((SF6cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$t0
  p5.PAM[times]<-list(residuals(protest(vegdist(decostand((SF6cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))))
  
  
  #6th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using square-root transformed data
  
  
  
  procr6[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                           scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig
  
  pprocr6[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  p6.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                        scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  #7th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using Hellinger transformed data
  
  procr7[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                           scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig  
  
  
  pprocr7[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p7.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                        scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  
  #8th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from log(x+1) transformed data
  
  procr8[times,1]<-protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                           scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr8[times,1]<-protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                            scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p8.PAM[times]<-list(residuals(protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                                        scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))))
  
  
  
  #9th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from Hellinger transformed data
  
  procr9[times,1]<-protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                           scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr9[times,1]<-protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p9.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                                        scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))))
  
  #10th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using square root transformed data
  
  procr10[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr10[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p10.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF6cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #11th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using Hellinger transformed data
  
  
  procr11[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr11[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p11.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF6cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from log(x+1) transformed data
  
  
  procr12[times,1]<-protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2),display=c("sites")),
                            scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$sig
  
  pprocr12[times,1]<-protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2),display=c("sites")),
                             scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$t0
  
  p12.PAM[times]<-list(residuals(protest(scores(rda(log((SF6cor^2)+1)),choice=c(1,2),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from Hellinger transformed data
  
  
  procr13[times,1]<-protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$sig
  
  
  pprocr13[times,1]<-protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$t0
  
  
  p13.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF6cor^2),"hel")),choice=c(1,2),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))))
  
  ####################################################################################################################
  #1th Pre-transformation:
  #simulatiing with the original covariance structure for both SF (fert.m) and FA (plfa.m)
  
  procr14[times,1]<-protest((SF9cor^2),plfa.m)$sig
  pprocr14[times,1]<-protest((SF9cor^2),plfa.m)$t0
  p14.PAM[times]<-list(residuals(protest((SF9cor^2),plfa.m)))
  
  #2th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a distance matrix formed
  # based on Bray-Curtis distance of square-rooted data
  
  procr15[times,1]<-protest(vegdist(sqrt(SF9cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$sig
  pprocr15[times,1]<-protest(vegdist(sqrt(SF9cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$t0
  p15.PAM[times]<-list(residuals(protest(vegdist(sqrt(SF9cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))))
  
  #3th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed  and a distance matrix formed
  # based on Bray-Curtis distance of Hellinger-transformed data
  
  procr16[times,1]<-protest(vegdist(decostand((SF9cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$sig
  pprocr16[times,1]<-protest(vegdist(decostand((SF9cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$t0
  p16.PAM[times]<-list(residuals(protest(vegdist(decostand((SF9cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))))
  
  
  #4th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a  distance matrix formed
  # based on Euclidean distance of log(X+1) transformed data
  
  procr17[times,1]<-protest(vegdist(log((SF9cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$sig
  pprocr17[times,1]<-protest(vegdist(log((SF9cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$t0
  p17.PAM[times]<-list(residuals(protest(vegdist(log((SF9cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))))
  
  
  
  #5th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) after both data tables  transformed and a  distance matrix formed
  # based on Euclidean distance of Hellinger - transformed data
  
  procr18[times,1]<-protest(vegdist(decostand((SF9cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$sig
  pprocr18[times,1]<-protest(vegdist(decostand((SF9cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$t0
  p18.PAM[times]<-list(residuals(protest(vegdist(decostand((SF9cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))))
  
  
  #6th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using square-root transformed data
  
  
  
  procr19[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig
  
  pprocr19[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  p19.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  #7th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using Hellinger transformed data
  
  procr20[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig  
  
  
  pprocr20[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p20.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  
  #8th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from log(x+1) transformed data
  
  procr21[times,1]<-protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                            scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr21[times,1]<-protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                             scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p21.PAM[times]<-list(residuals(protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))))
  
  
  
  #9th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from Hellinger transformed data
  
  procr22[times,1]<-protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr22[times,1]<-protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p22.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))))
  
  #10th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using square root transformed data
  
  procr23[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr23[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p23.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF9cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #11th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using Hellinger transformed data
  
  
  procr24[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr24[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p24.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF9cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from log(x+1) transformed data
  
  
  procr25[times,1]<-protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2),display=c("sites")),
                            scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$sig
  
  pprocr25[times,1]<-protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2),display=c("sites")),
                             scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$t0
  
  p25.PAM[times]<-list(residuals(protest(scores(rda(log((SF9cor^2)+1)),choice=c(1,2),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from Hellinger transformed data
  
  
  procr26[times,1]<-protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$sig
  
  
  pprocr26[times,1]<-protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$t0
  
  
  p26.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF9cor^2),"hel")),choice=c(1,2),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))))
  
  ################################################################################################################################
  #1th Pre-transformation:
  #simulatiing with the original covariance structure for both SF (fert.m) and FA (plfa.m)
  
  procr27[times,1]<-protest((SF12cor^2),plfa.m)$sig
  pprocr27[times,1]<-protest((SF12cor^2),plfa.m)$t0
  p27.PAM[times]<-list(residuals(protest((SF12cor^2),plfa.m)))
  
  #2th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a distance matrix formed
  # based on Bray-Curtis distance of square-rooted data
  
  procr28[times,1]<-protest(vegdist(sqrt(SF12cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$sig
  pprocr28[times,1]<-protest(vegdist(sqrt(SF12cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$t0
  p28.PAM[times]<-list(residuals(protest(vegdist(sqrt(SF12cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))))
  
  #3th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed  and a distance matrix formed
  # based on Bray-Curtis distance of Hellinger-transformed data
  
  procr29[times,1]<-protest(vegdist(decostand((SF12cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$sig
  pprocr29[times,1]<-protest(vegdist(decostand((SF12cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$t0
  p29.PAM[times]<-list(residuals(protest(vegdist(decostand((SF12cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))))
  
  
  #4th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a  distance matrix formed
  # based on Euclidean distance of log(X+1) transformed data
  
  procr30[times,1]<-protest(vegdist(log((SF12cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$sig
  pprocr30[times,1]<-protest(vegdist(log((SF12cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$t0
  p30.PAM[times]<-list(residuals(protest(vegdist(log((SF12cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))))
  
  
  
  #5th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) after both data tables  transformed and a  distance matrix formed
  # based on Euclidean distance of Hellinger - transformed data
  
  procr31[times,1]<-protest(vegdist(decostand((SF12cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$sig
  pprocr31[times,1]<-protest(vegdist(decostand((SF12cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$t0
  p31.PAM[times]<-list(residuals(protest(vegdist(decostand((SF12cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))))
  
  
  #6th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using square-root transformed data
  
  
  
  procr32[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig
  
  pprocr32[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  p32.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  #7th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using Hellinger transformed data
  
  procr33[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig  
  
  
  pprocr33[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p33.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  
  #8th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from log(x+1) transformed data
  
  procr34[times,1]<-protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                            scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr34[times,1]<-protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                             scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p34.PAM[times]<-list(residuals(protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))))
  
  
  
  #9th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from Hellinger transformed data
  
  procr35[times,1]<-protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr35[times,1]<-protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p35.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))))
  
  #10th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using square root transformed data
  
  procr36[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr36[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p36.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF12cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #11th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using Hellinger transformed data
  
  
  procr37[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr37[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p37.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF12cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from log(x+1) transformed data
  
  
  procr38[times,1]<-protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2),display=c("sites")),
                            scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$sig
  
  pprocr38[times,1]<-protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2),display=c("sites")),
                             scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$t0
  
  p38.PAM[times]<-list(residuals(protest(scores(rda(log((SF12cor^2)+1)),choice=c(1,2),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from Hellinger transformed data
  
  
  procr39[times,1]<-protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$sig
  
  
  pprocr39[times,1]<-protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$t0
  
  
  p39.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF12cor^2),"hel")),choice=c(1,2),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))))
  ################################################################################################################################
  #1th Pre-transformation:
  #simulatiing with the original covariance structure for both SF (fert.m) and FA (plfa.m)
  
  procr40[times,1]<-protest((SF15cor^2),plfa.m)$sig
  pprocr40[times,1]<-protest((SF15cor^2),plfa.m)$t0
  p40.PAM[times]<-list(residuals(protest((SF15cor^2),plfa.m)))
  
  #2th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a distance matrix formed
  # based on Bray-Curtis distance of square-rooted data
  
  procr41[times,1]<-protest(vegdist(sqrt(SF15cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$sig
  pprocr41[times,1]<-protest(vegdist(sqrt(SF15cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))$t0
  p41.PAM[times]<-list(residuals(protest(vegdist(sqrt(SF15cor^2),"bray"),vegdist(sqrt(plfa.m),"bray"))))
  
  #3th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed  and a distance matrix formed
  # based on Bray-Curtis distance of Hellinger-transformed data
  
  procr42[times,1]<-protest(vegdist(decostand((SF15cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$sig
  pprocr42[times,1]<-protest(vegdist(decostand((SF15cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))$t0
  p42.PAM[times]<-list(residuals(protest(vegdist(decostand((SF15cor^2),"hel"),"bray"),vegdist(decostand(plfa.m,"hel"),"bray"))))
  
  
  #4th Pre-transformation:
  #simulating Procrustes relationship between (fert.m) and FA (plfa.m) after both data tables transformed and a  distance matrix formed
  # based on Euclidean distance of log(X+1) transformed data
  
  procr43[times,1]<-protest(vegdist(log((SF15cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$sig
  pprocr43[times,1]<-protest(vegdist(log((SF15cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))$t0
  p43.PAM[times]<-list(residuals(protest(vegdist(log((SF15cor^2)+1),"euc"),vegdist(log(plfa.m+1),"euc"))))
  
  
  
  #5th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) after both data tables  transformed and a  distance matrix formed
  # based on Euclidean distance of Hellinger - transformed data
  
  procr44[times,1]<-protest(vegdist(decostand((SF15cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$sig
  pprocr44[times,1]<-protest(vegdist(decostand((SF15cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))$t0
  p44.PAM[times]<-list(residuals(protest(vegdist(decostand((SF15cor^2),"hel"),"euc"),vegdist(decostand(plfa.m,"hel"),"euc"))))
  
  
  #6th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using square-root transformed data
  
  
  
  procr45[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig
  
  pprocr45[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  p45.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  #7th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (3 axes) from Bray-Curtis matrix distance using Hellinger transformed data
  
  procr46[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig  
  
  
  pprocr46[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p46.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))))
  
  
  #8th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from log(x+1) transformed data
  
  procr47[times,1]<-protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                            scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr47[times,1]<-protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                             scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p47.PAM[times]<-list(residuals(protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2,3),display=c("sites")))))
  
  
  
  #9th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (3 axes) from Hellinger transformed data
  
  procr48[times,1]<-protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$sig
  
  
  pprocr48[times,1]<-protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))$t0
  
  
  p48.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2,3),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2,3),display=c("sites")))))
  
  #10th Pre-transformation:
  #simulating Procrustes relationship between soi fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using square root transformed data
  
  procr49[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr49[times,1]<-protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p49.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(SF15cor^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(sqrt(plfa.m),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #11th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on NMS axes (2 axes) from Bray-Curtis distance using Hellinger transformed data
  
  
  procr50[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                            scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig
  
  pprocr50[times,1]<-protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                             scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0
  
  p50.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((SF15cor^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")),
                                         scores(metaMDS(vegdist(decostand(plfa.m,"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from log(x+1) transformed data
  
  
  procr51[times,1]<-protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2),display=c("sites")),
                            scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$sig
  
  pprocr51[times,1]<-protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2),display=c("sites")),
                             scores(rda(log((plfa.m^2)+1)),choice=c(1,2),display=c("sites")))$t0
  
  p51.PAM[times]<-list(residuals(protest(scores(rda(log((SF15cor^2)+1)),choice=c(1,2),display=c("sites")),
                                         scores(rda(log(plfa.m+1)),choice=c(1,2),display=c("sites")))))
  
  
  #12th Pre-transformation:
  #simulating Procrustes relationship between soil fertility (fert.m) and microbial community (plfa.m) 
  #based on PCA axes (2 axes) from Hellinger transformed data
  
  
  procr52[times,1]<-protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2),display=c("sites")),
                            scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$sig
  
  
  pprocr52[times,1]<-protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2),display=c("sites")),
                             scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))$t0
  
  
  p52.PAM[times]<-list(residuals(protest(scores(rda(decostand((SF15cor^2),"hel")),choice=c(1,2),display=c("sites")),
                                         scores(rda(decostand(plfa.m,"hel")),choice=c(1,2),display=c("sites")))))
  
  
} #HERE WE FINSH THE CODE FOR THE 100 SIMULATIONS of each Pre-Transformations for each number of columns correlating in X data table

#Notice that we could have increased the code by considering 13 Pre-Tranformations X 4 covariance levels 
# for each number of columns in the Soil Fertility data table. 


citation(rjags)

# The following lines stand for PAMs derived from 4 sets of the same 13 Pre-Transformations, 
#one for each number of columns covarying (06, 09, 12, 15) at a specific correlation level (in this case 0.9)

#SF6
MA1<-mapply(as.matrix,p1.PAM) # It compile the 100 simulated PAMs of the same Pre-Transformation into a matrix
p1<-apply(MA1,1,mean) #It takes the average of the 100 simulationed PAMs in a specific Pre-Transformation
p1

MA2<-mapply(as.matrix,p2.PAM)
p2<-apply(MA2,1,mean)

MA3<-mapply(as.matrix,p3.PAM)
p3<-apply(MA3,1,mean)

MA4<-mapply(as.matrix,p4.PAM)
p4<-apply(MA4,1,mean)

MA5<-mapply(as.matrix,p5.PAM)
p5<-apply(MA5,1,mean)

MA6<-mapply(as.matrix,p6.PAM)
p6<-apply(MA6,1,mean)

MA7<-mapply(as.matrix,p7.PAM)
p7<-apply(MA7,1,mean)

MA8<-mapply(as.matrix,p8.PAM)
p8<-apply(MA8,1,mean)

MA9<-mapply(as.matrix,p9.PAM)
p9<-apply(MA9,1,mean)

MA10<-mapply(as.matrix,p10.PAM)
p10<-apply(MA10,1,mean)

MA11<-mapply(as.matrix,p11.PAM)
p11<-apply(MA11,1,mean)

MA12<-mapply(as.matrix,p12.PAM)
p12<-apply(MA12,1,mean)

MA13<-mapply(as.matrix,p13.PAM)
p13<-apply(MA13,1,mean)

#SF9
############################
MA14<-mapply(as.matrix,p14.PAM)
p14<-apply(MA14,1,mean)

MA15<-mapply(as.matrix,p15.PAM)
p15<-apply(MA15,1,mean)

MA16<-mapply(as.matrix,p16.PAM)
p16<-apply(MA16,1,mean)

MA17<-mapply(as.matrix,p17.PAM)
p17<-apply(MA17,1,mean)

MA18<-mapply(as.matrix,p18.PAM)
p18<-apply(MA18,1,mean)

MA19<-mapply(as.matrix,p19.PAM)
p19<-apply(MA19,1,mean)

MA20<-mapply(as.matrix,p20.PAM)
p20<-apply(MA20,1,mean)

MA21<-mapply(as.matrix,p21.PAM)
p21<-apply(MA21,1,mean)

MA22<-mapply(as.matrix,p22.PAM)
p22<-apply(MA22,1,mean)

MA23<-mapply(as.matrix,p23.PAM)
p23<-apply(MA23,1,mean)

MA24<-mapply(as.matrix,p24.PAM)
p24<-apply(MA24,1,mean)

MA25<-mapply(as.matrix,p25.PAM)
p25<-apply(MA25,1,mean)

MA26<-mapply(as.matrix,p26.PAM)
p26<-apply(MA26,1,mean)
################################


#SF12
############################
MA27<-mapply(as.matrix,p27.PAM)
p27<-apply(MA27,1,mean)

MA28<-mapply(as.matrix,p28.PAM)
p28<-apply(MA28,1,mean)

MA29<-mapply(as.matrix,p29.PAM)
p29<-apply(MA29,1,mean)

MA30<-mapply(as.matrix,p30.PAM)
p30<-apply(MA30,1,mean)

MA31<-mapply(as.matrix,p31.PAM)
p31<-apply(MA31,1,mean)

MA32<-mapply(as.matrix,p32.PAM)
p32<-apply(MA32,1,mean)

MA33<-mapply(as.matrix,p33.PAM)
p33<-apply(MA33,1,mean)

MA34<-mapply(as.matrix,p34.PAM)
p34<-apply(MA34,1,mean)

MA35<-mapply(as.matrix,p35.PAM)
p35<-apply(MA35,1,mean)

MA36<-mapply(as.matrix,p36.PAM)
p36<-apply(MA36,1,mean)

MA37<-mapply(as.matrix,p37.PAM)
p37<-apply(MA37,1,mean)

MA38<-mapply(as.matrix,p38.PAM)
p38<-apply(MA38,1,mean)

MA39<-mapply(as.matrix,p39.PAM)
p39<-apply(MA39,1,mean)
################################

################################


#SF15
############################
MA40<-mapply(as.matrix,p40.PAM)
p40<-apply(MA40,1,mean)

MA41<-mapply(as.matrix,p41.PAM)
p41<-apply(MA41,1,mean)

MA42<-mapply(as.matrix,p42.PAM)
p42<-apply(MA42,1,mean)

MA43<-mapply(as.matrix,p43.PAM)
p43<-apply(MA43,1,mean)

MA44<-mapply(as.matrix,p44.PAM)
p44<-apply(MA44,1,mean)

MA45<-mapply(as.matrix,p45.PAM)
p45<-apply(MA45,1,mean)

MA46<-mapply(as.matrix,p46.PAM)
p46<-apply(MA46,1,mean)

MA47<-mapply(as.matrix,p47.PAM)
p47<-apply(MA47,1,mean)

MA48<-mapply(as.matrix,p48.PAM)
p48<-apply(MA48,1,mean)

MA49<-mapply(as.matrix,p49.PAM)
p49<-apply(MA49,1,mean)

MA50<-mapply(as.matrix,p50.PAM)
p50<-apply(MA50,1,mean)

MA51<-mapply(as.matrix,p51.PAM)
p51<-apply(MA51,1,mean)

MA52<-mapply(as.matrix,p52.PAM)
p52<-apply(MA52,1,mean)



###############################################################

#After getting the 52 PAMs, these are going to set aside for further analysis


# Obtaining different information from the results



#P value (significance) of correlation (put it into Excel spreedshet!)
#***************************
length(which(procr52[,1]<0.05))
mean(pprocr52[,1])
max(p52)-min(p52) #  Calculates the range of the residuals within the PAM for a specific Pre-Transformation (in this case, the number 1) ( a measure of Procrustean residual variability)

MA<-mapply(as.matrix,p52.PAM) #Creates the matrix of simulated PAMs for the SAME PRE-TRANSFORMATION (in this case the number 1)
LIST<-list()
for(i in 1:ncol(MA)) # For each simulated PAM (100) from the SAME PRE - TRANSFORMATION it applies an ANOVA having Land Use Type as a 4-level factor
{
  mod<-lm(MA[,i] ~ CD$LandUse,data=CD)
  LIST[[i]]<-Anova(mod,type="III")$"Pr(>F)"[2]
  sig.perc<-(length(which(unlist(LIST)<0.05))/100)
}

sig.perc # It returns the number of ANOVAs came out with significant
##########################################


#Creating a PAM matrix that is going to be used in ordination.

WIT20new2<-cbind(as.matrix(p1),as.matrix(p2),as.matrix(p3),as.matrix(p4),as.matrix(p5),
                    as.matrix(p6),as.matrix(p7),as.matrix(p8),as.matrix(p9),as.matrix(p10),
                    as.matrix(p11),as.matrix(p12),as.matrix(p13),as.matrix(p14),as.matrix(p15),
                    as.matrix(p16),as.matrix(p17),as.matrix(p18),as.matrix(p19),as.matrix(p20),
                    as.matrix(p21),as.matrix(p22),as.matrix(p23),as.matrix(p24),as.matrix(p25),
                    as.matrix(p26),as.matrix(p27),as.matrix(p28),as.matrix(p29),as.matrix(p30),
                    as.matrix(p31),as.matrix(p32),as.matrix(p33),as.matrix(p34),as.matrix(p35),
                    as.matrix(p36),as.matrix(p37),as.matrix(p38),as.matrix(p39),as.matrix(p40),
                    as.matrix(p41),as.matrix(p42),as.matrix(p43),as.matrix(p44),as.matrix(p45),
                    as.matrix(p46),as.matrix(p47),as.matrix(p48),as.matrix(p49),as.matrix(p50),
                    as.matrix(p51),as.matrix(p52))

dim(WIT20new2)

write.table(WIT20new2, file="PAMWIT20new2.txt", row.names=F, col.names=F)


#That is the end of the code!!!