#This file is showing how to simulate Procrustes analysis #in order to verify if PAM is able to come out the multivariate correlation between two data tables for a specific treatment # Considering two simulated data tables X and Y varying in their dimensions ( n=40 and p = 25, 45, and 80) # We applied between the firt 10 rows of both X and Y data tables four correlations levels of 90%, 70%, and # These first 10 rows in X and B correponding to a hypothetical treatment A. # The rest of the line represents 3 others treatment (B, C, and D) which have correlation not higher than 10% #The following sequence shows what will be done # X (n=40, p = 25) Vs Y (n = 40, p = 25) #******************************************************* # 13 Procrustes relationships in 90% of the correlation for the treatment A # 13 Procrustes relationships in 70% of the correlation for the treatment A # 13 Procrustes relationships in 50% of the correlation for the treatment A # 13 Procrustes relationships in 20% of the correlation for the treatment A # X (n=40, p = 25) Vs Y (n = 40, p = 45) #******************************************************* # 13 Procrustes relationships in 90% of the correlation for the treatment A # 13 Procrustes relationships in 70% of the correlation for the treatment A # 13 Procrustes relationships in 50% of the correlation for the treatment A # 13 Procrustes relationships in 20% of the correlation for the treatment A # X (n=40, p = 25) Vs Y (n = 40, p = 80) #******************************************************* # 13 Procrustes relationships in 90% of the correlation for the treatment A # 13 Procrustes relationships in 70% of the correlation for the treatment A # 13 Procrustes relationships in 50% of the correlation for the treatment A # 13 Procrustes relationships in 20% of the correlation for the treatment A ############################################################ #Functions for generating differents covariance structure #Routnes 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(corrplot) library(vegan) library(car) library(car) library(agricolae) #Creating objects for simulations #*********************************************************************** procr1<-array(NA,dim=c(100,1)) pprocr1<-array(NA,dim=c(100,1)) p1.PAM<-rep(NA,100) 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) ##################################################################################### for(times in 1:100){ #Starting the simulations # 2) transfering the covariance structure (90%)for the "BIG" data frame: M25col<-simcor(k=1,size=c(160),rho=c(0.91),epsilon=0.91-0.90) # 90% covariance Total90<-t(chol(M25col))%*% matrix(rnorm(160*10), nrow=160, ncol=10) M90.df<-as.data.frame(t(Total90)) M90.df #corrplot(cor(M90.df)) # 3) dividing the "BIG" dataframe into two matrices wth equak number of columns and lines ( n = 10, p=45, cov= 90%) AX25.df<-M90.df[,1:80] AY25.df<-M90.df[,81:160] #corrplot(cor(AX45.df)) #Checking the covariance incorporated to the AX90 and #corrplot(cor(AY45.df)) # 4) Now we are going to make a matrix with the rest of treatments (B, C, D). #These matrix will have dimenison (n = 30, p = 45, COV X and Y = < 10%) RM90col<-simcor(k=1,size=c(160),rho=c(0.05),epsilon=0.05-0.04) #< 10% covariance # 4) Incorporating the <10% covariance into the data frame (n=30,p= 45 ) RZ<-matrix(rnorm(160*30), nrow=30, ncol=160) RZ.m<-as.matrix(RZ) RTotal90<-t(chol(RM90col))%*% matrix(rnorm(160*30), nrow=160, ncol=30) RM90.df<-as.data.frame(t(RTotal90)) #RM90.df #corrplot(cor(RM90.df)) # 5) Spliting the "BIG" low covariance matrix (< 10%) into two part (X and Y) RX25.df<-RM90.df[,1:80] # creates the X part of the data frame with <10% between B, c, and D categorical levels RY25.df<-RM90.df[,81:160] # creates the Y part of the data frame with <10% between B, c, and D categorical levels #corrplot(cor(RX45.df)) # checking the covariance incorporated #corrplot(cor(RY45.df)) # checking the covariance incorporated #6) Linking the (X90 to X0) in order to build a matrix ( n = 40, p = 45, Cov for A = 90% and Cov for B, C, abd D < 10%) colnames(AX25.df)<-paste("var",seq(1,80),sep="") row.names(AX25.df)<-paste("site",seq(1,10),sep="") row.names(RX25.df)<-paste("site",seq(11,40),sep="") AX25.m<-as.matrix(AX25.df) # Categorical A part of the X matrix ( Cov 90%) RX25.m<-as.matrix(RX25.df) # Other categorical level part of the matrix (B, C, and D, Cov < 10%) X_A90.p80<-rbind(AX25.m,RX25.m) #Linking # 7) Creating Y_A90 colnames(AY25.df)<-paste("var",seq(1,80),sep="") row.names(AY25.df)<-paste("site",seq(1,10),sep="") row.names(RY25.df)<-paste("site",seq(11,40),sep="") AY25.m<-as.matrix(AY25.df) # Categorical A part of the X matrix ( Cov 90%) RY25.m<-as.matrix(RY25.df) # Other categorical level part of the matrix (B, C, and D, Cov < 10%) Y_A90.p80<-rbind(AY25.m,RY25.m) # linking categorical A part of Y (cov 90%) to others categorical levels (B, C, and D; Cov < 10%) #Y_A90.p45 #In the total #corrplot(cor(Y_A90.p45[1:10,])) #Em A #corrplot(cor(Y_A90.p45[11:40,])) ###################### #Finally we have created two simulated matrices with a given between covariance level #for the categorical level A and other between covariance level for the others categorical #levels (B, C, and D). #Set aside Y_A90 and X_A90. These will be related by Procrustes analysis # using Raw, distance matrices, and ordination axes as describe in the main text of the manuscript XA90.p80<-X_A90.p80 YA90.p80<-Y_A90.p80 ############ #Initiating the simulations the Procrustes relationshis ######################################################################################## # here we are just showing the 13 Pre-Transformations before the Procrustean # relationship between X (n=40, p=80) Vs Y (n=40, p=80) Cov for treatment A = 90% ##### # For doing the rest of simulations it is enough changes the settings above ############################################### procr1[times,1]<-protest(XA90.p80^2,YA90.p80^2)$sig pprocr1[times,1]<-protest(XA90.p80^2,YA90.p80^2)$t0 p1.PAM[times]<-list(residuals(XA90.p80^2,YA90.p80^2))) ########### procr2[times,1]<-protest(vegdist(sqrt(XA90.p80^2),"bray"),vegdist(sqrt(YA90.p80^2),"bray"))$sig pprocr2[times,1]<-protest(vegdist(sqrt(XA90.p80^2),"bray"),vegdist(sqrt(YA90.p80^2),"bray"))$t0 p2.PAM[times]<-list(residuals(protest(vegdist(sqrt(XA90.p80^2),"bray"),vegdist(sqrt(YA90.p80^2),"bray")))) ############ procr3[times,1]<-protest(vegdist(decostand((XA90.p80^2),"hel"),"bray"),vegdist(decostand(YA90.p80^2,"hel"),"bray"))$sig pprocr3[times,1]<-protest(vegdist(decostand((XA90.p80^2),"hel"),"bray"),vegdist(decostand(YA90.p80^2,"hel"),"bray"))$t0 p3.PAM[times]<-list(residuals(protest(vegdist(decostand((XA72^2),"hel"),"bray"),vegdist(decostand(YA72^2,"hel"),"bray")))) ########### procr4[times,1]<-protest(vegdist(log((XA90.p80^2)+1),"euc"),vegdist(log((YA90.p80^2)+1),"euc"))$sig pprocr4[times,1]<-protest(vegdist(log((XA90.p80^2)+1),"euc"),vegdist(log((YA90.p80^2)+1),"euc"))$t0 p4.PAM[times]<-list(residuals(protest(vegdist(log((XA90.p80^2)+1),"euc"),vegdist(log((YA90.p80^2)+1),"euc")))) ########## procr5[times,1]<-protest(vegdist(decostand((XA90.p80^2),"hel"),"euc"),vegdist(decostand((YA90.p80^2),"hel"),"euc"))$sig pprocr5[times,1]<-protest(vegdist(decostand((XA90.p80^2),"hel"),"euc"),vegdist(decostand((YA90.p80^2),"hel"),"euc"))$t0 p5.PAM[times]<-list(residuals(protest(vegdist(decostand((XA90.p80^2),"hel"),"euc"),vegdist(decostand((YA90.p80^2),"hel"),"euc")))) ########## procr6[times,1]<-protest(scores(metaMDS(vegdist(sqrt((XA90.p80^2)^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig pprocr6[times,1]<-protest(scores(metaMDS(vegdist(sqrt(XA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$t0 p6.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(XA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites"))))) ######## procr7[times,1]<-protest(scores(metaMDS(vegdist(decostand((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")))$sig pprocr7[times,1]<-protest(scores(metaMDS(vegdist(decostand((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"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((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2,3),display=c("sites"))))) ########### procr8[times,1]<-protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2,3),display=c("sites")), scores(rda(log((YA90.p80^2)+1)),choice=c(1,2,3),display=c("sites")))$sig pprocr8[times,1]<-protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2,3),display=c("sites")), scores(rda(log((YA90.p80^2))),choice=c(1,2,3),display=c("sites")))$t0 p8.PAM[times]<-list(residuals(protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2,3),display=c("sites")), scores(rda(log((YA90.p80^2)+1)),choice=c(1,2,3),display=c("sites"))))) ############## procr9[times,1]<-protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites")))$sig pprocr9[times,1]<-protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites")))$t0 p9.PAM[times]<-list(residuals(protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2,3),display=c("sites"))))) ############# procr10[times,1]<-protest(scores(metaMDS(vegdist(sqrt(XA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig pprocr10[times,1]<-protest(scores(metaMDS(vegdist(sqrt(XA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0 p10.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(sqrt(XA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(sqrt(YA90.p80^2),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites"))))) ########## procr11[times,1]<-protest(scores(metaMDS(vegdist(decostand((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$sig pprocr11[times,1]<-protest(scores(metaMDS(vegdist(decostand((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")))$t0 p11.PAM[times]<-list(residuals(protest(scores(metaMDS(vegdist(decostand((XA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites")), scores(metaMDS(vegdist(decostand((YA90.p80^2),"hel"),"bray"),k=3,trymax=50),choice=c(1,2),display=c("sites"))))) ################ procr12[times,1]<-protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2),display=c("sites")), scores(rda(log((YA90.p80^2)+1)),choice=c(1,2),display=c("sites")))$sig pprocr12[times,1]<-protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2),display=c("sites")), scores(rda(log((YA90.p80^2)+1)),choice=c(1,2),display=c("sites")))$t0 p12.PAM[times]<-list(residuals(protest(scores(rda(log((XA90.p80^2)+1)),choice=c(1,2),display=c("sites")), scores(rda(log((YA90.p80^2)+1)),choice=c(1,2),display=c("sites"))))) ############### procr13[times,1]<-protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2),display=c("sites")))$sig pprocr13[times,1]<-protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2),display=c("sites")))$t0 p13.PAM[times]<-list(residuals(protest(scores(rda(decostand((XA90.p80^2),"hel")),choice=c(1,2),display=c("sites")), scores(rda(decostand((YA90.p80^2),"hel")),choice=c(1,2),display=c("sites"))))) } #Final of the simulations!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #PAM #************************************************************************************************************************* MA1<-mapply(as.matrix,p1.PAM) p1<-apply(MA1,1,mean) 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) ############################################################### #After getting the 52 PAMs (13 Pre-Transformation for each level ov covariance between X and Y for the treatment A) #these are going to set aside for further analysis # Obtaining different unformation from the results #P value (significance) of correlation (put it into Excel spreedshet!) #*************************** length(which(procr1[,1]<0.05)) #Within 100 simulations for THE SAME PRE-TRANSFORMATION (in this case the number 1) #Calculate the mumber of time in wich the Procrustes corrlation statistic (R) #was considered significat (P<0.05) #Procrustes correlation mean(pprocr40[,1]) #Calculate the mean of the Procrustean correlation statistic (R) based on 100 simulations for THE SAME PRE-TRANSFORMATION (in this case the number 1) range(pprocr40[,1]) # calculate the range of the Procrustean correlation statistic (R) based on 100 simulations for THE SAME PRE-TRANSFORMATION (in this case the number 1) ########################## max(p1)-min(p1) # 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) #################################### #Rest of analysis (ANOVA and post hoc) ######################################### #ANOVA #looping for runing an ANOVA for each of the 100 simulated PAMs for THE SAME PRE-TRANSFORMATION ( in this case the number 1) MA<-mapply(as.matrix,p1.PAM) #Creates the matix 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 a 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 significan ########################################## #Creating a PAM matrix that is going to be used in ordination. BETWEEN.90<-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(BETWEEN.90) write.table(BETWEEN.70, file="BETPAM90.txt", row.names=F, col.names=F)