# 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)