# The following is example R-code that can be used to run the analysis # we performed on the salamander data in Appendix A. Modifications are # likely necessary for other data. Nevertheless, computational steps # should be obvious # First step, import and adjust the necessary matrices # Matrices for our example are provided in Ecological Archives y.mat<-read.csv("file_name",header=TRUE) # A file name and location needs to be specified y.mat<-as.matrix(y.mat) dim(y.mat)<-c(336,3) x.mat.full<-read.csv("file_name",header=TRUE) # A file name and location needs to be specified x.mat.full<-as.matrix(x.mat.full) dim(x.mat.full)<-c(336,4) x.mat.red<-x.mat.full[,-4] # This removes the coding for the interaction effect # Second step, estimate parameters for full and reduced models b.mat.full<-solve((t(x.mat.full)%*%x.mat.full))%*%(t(x.mat.full)%*%y.mat) b.mat.red<-solve((t(x.mat.red)%*%x.mat.red))%*%(t(x.mat.red)%*%y.mat) # Third step, Estimate LS means # In this case, but not necessarily always, a = species and b = locality # Vectors will be caculated for 2 species, a1 (P. jordani) and a2 (P. teyahalee) # In our example the following dummy variables were used: # P. jordani = 1, P. teyahalee = -1, sympatry = 1, allopatry = -1 # Thus, LS means could be calculated as follows: a1b1<-cbind(1,1,1,1) a1b2<-cbind(1,1,-1,-1) a2b1<-cbind(1,-1,1,-1) a2b2<-cbind(1,-1,-1,1) x.ls.full<-rbind(a1b1,a1b2,a2b1,a2b2) x.ls.red<-x.ls.full[,-4] obs.ls.full<-x.ls.full%*%b.mat.full # Observed ls means (full model) obs.ls.red<-x.ls.red%*%b.mat.red # Observed ls means (reduced model) # Fourth Step, vector and statistics calculations obs.a1.vect<-obs.ls.full[1,]-obs.ls.full[2,] # These are the phenotypic change vectors obs.a2.vect<-obs.ls.full[3,]-obs.ls.full[4,] obs.d.a1<-sqrt(t(obs.a1.vect)%*%obs.a1.vect) # These are lengths of vectors obs.d.a2<-sqrt(t(obs.a2.vect)%*%obs.a2.vect) obs.contrast<-abs(obs.d.a1-obs.d.a2) obs.angle<-acos(t((obs.a1.vect)/obs.d.a1)%*%((obs.a2.vect)/obs.d.a2)) obs.angle<-obs.angle*180/pi # This step is only necessary to convert radians to degrees # Fifth Step, set-up premutation procedure y.hat<-x.mat.red%*%b.mat.red # Predicted values from reduced model y.res<-y.mat-y.hat # Resdiuals of reduced mode (these are the permuted units) # PERMUTATION PROCEDURE permute<-999 # This is the number of random iterations to occur (can be varied) # Need to set-up distributions to be generated dist.d1<-NULL dist.d1<-rbind(dist.d1,obs.d.a1) # Observed value is first random value dist.d2<-NULL dist.d2<-rbind(dist.d2,obs.d.a2) # Observed value is first random value dist.contrast<-NULL dist.contrast<-rbind(dist.contrast,obs.contrast) # Observed value is first random value dist.angle<-NULL dist.angle<-rbind(dist.angle,obs.angle) # Observed value is first random value # In addition to saving random values, it is wise to save the outcome of comparisons # of observed and random values. # This can be done with logical statements (below). # Separate distributions are created for these comparisons. # The 'p' indicates that these distributions will be used # to calculate empirical probabilities. pdist.contrast<-NULL pdist.contrast<-rbind(pdist.contrast,1) pdist.angle<-NULL pdist.angle<-rbind(pdist.angle,1) # Create an array from 1 to number of object in data set # This will be randomized later line<-array(1:(length(x.mat.full[,1])),dim=c(length(x.mat.full[,1]))) for(i in 1:permute){ line.rand<-sample(line,replace=FALSE) y.res.temp<-cbind(line.rand,y.res) z<-(order(line.rand)) y.res.temp2<-as.matrix(y.res.temp[z,]) y.res.rand<-y.res.temp2[,-1] # Rows of residuals are now randomized # Create random values y.rand<-y.hat+y.res.rand # Estimate parameters b.mat.rand<-solve((t(x.mat.full)%*%x.mat.full))%*%(t(x.mat.full)%*%y.rand) # Calculate LS means rand.ls.full<-x.ls.full%*%b.mat.rand # Repeat fourth step for random data! rand.a1.vect<-rand.ls.full[1,]-rand.ls.full[2,] # These are the phenotypic change vectors rand.a2.vect<-rand.ls.full[3,]-rand.ls.full[4,] rand.d.a1<-sqrt(t(rand.a1.vect)%*%rand.a1.vect) # These are lengths of vectors rand.d.a2<-sqrt(t(rand.a2.vect)%*%rand.a2.vect) rand.contrast<-abs(rand.d.a1-rand.d.a2) rand.angle<-acos(t((rand.a1.vect)/rand.d.a1)%*%((rand.a2.vect)/rand.d.a2)) rand.angle<-rand.angle*180/pi # This step is only necessary to convert radians to degrees # Append distributions dist.d1<-rbind(dist.d1,rand.d.a1) dist.d2<-rbind(dist.d2,rand.d.a2) dist.contrast<-rbind(dist.contrast,rand.contrast) dist.angle<-rbind(dist.angle,rand.angle) aa<-ifelse(rand.contrast>=obs.contrast,1,0) bb<-ifelse(rand.angle>=obs.angle,1,0) pdist.contrast<-rbind(pdist.contrast,aa) pdist.angle<-rbind(pdist.angle,bb) } # Empirical probailities are calculated as follows p.contrast<-sum(pdist.contrast)/(permute+1) p.angle<-sum(pdist.angle)/(permute+1) # OUTPUT out.head<-cbind("iteration","d1","d2","abs_diff","angle") iter<-array(1:permute) iter.lab<-"observed" iter.lab<-append(iter.lab,iter) out.mat<-cbind(iter.lab,dist.d1,dist.d2,dist.contrast,dist.angle) out.mat<-rbind(out.head,out.mat) write.csv(out.mat,file="path/output.csv",row.names=FALSE) # A path needs to be specified # The previous step should produce a table with every random distance, contrast, and angle. # Probabilities can be calculated manually from this table, or as above in R. # write.table can be used in place of write.csv. # The following step provides the results of analysis in text form result1<-"The length of vector 1 is " result1<-append(result1,obs.d.a1) result2<-"The length of vector 2 is " result2<-append(result2,obs.d.a2) result3<-"The absolute difference between these lengths is, |d1 - d2|= " result3<-append(result3,obs.contrast) result4<-"Prand = " result4<-append(result4,p.contrast) result5<-"The angle between vectors is " result5<-append(result5,obs.angle) result6<-"Prand = " result6<-append(result6,p.angle) result<-rbind(result1,result2,result3,result4,result5,result6) write.table(result,file="path/result.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)