#### Adams and Collyer Example Data #########
#  Phenotypic Trajectory Analysis	

# CODE ASSUMES factor 1 = lineage, factor 2=evol. levels
rm(list=ls())
mydata<-read.csv("example1.data.csv",header=T)
taxa<-as.factor(mydata[,1])
level<-as.factor(mydata[,2])
Y<-as.matrix(mydata[,-(1:2)])	
taxaBylevel<-as.factor(paste(taxa,level))  #for LS means at each level

Y<-prcomp(Y)$x	#TO USE IF > 2 variables: for plotting PC axes at the end
######################
##### Basic calculations

#data parameters (n=#spec, p=#ldmk, k=#dim)
n<-length(levels(taxa)) 
p<-length(levels(level))  
k<-ncol(Y) #dimensions of phenotypic data

#Full model
lm.full<-lm(Y~taxa*level,model=T,x=T,y=T,qr=T) 
yhat.full<-predict(lm.full)
summary(manova(lm.full))

#Reduced model
lm.red<-lm(Y ~ taxa + level, model=T,x=T,y=T,qr=T)   #just the covariate
yhat.red<-predict(lm.red)
res.red<-resid(lm.red)

#LS means
lsmeans.obs<-NULL
for (i in 1:k){
  mean.temp<-tapply(yhat.full[,i],taxaBylevel,mean)
  lsmeans.obs<-cbind(lsmeans.obs,mean.temp)
}

#Ls means organized as 3-d arrays
arrayspecs<-function(A){
n.m<-NULL # n.m stands for "new" means
for(i in 1:n){
  temp<-as.matrix(A[((1+(i-1)*p):(i*p)),1:k])
  n.m<-cbind(n.m,temp)}
trajectories<-array(n.m,dim=c(p,k,n))
}

###########################
####Functions calculating Trajectory attributes

# Pathlength distance
pathdist<-function(M) {as.matrix(dist(M))}	#create distance matrix
trajsize<-function(M){
traj.pathdist<-array(0,dim=c(n,1))   		#loop across trajectories
for (i in 1:n){
  temp<-pathdist(M[,,i])
  for (j in 1:(p-1)){
    traj.pathdist[i]<-traj.pathdist[i]+temp[j,j+1]
  }
}
traj.size.dist<-as.matrix(dist(traj.pathdist))		#TrajSizeDiff
}

#trajectory direction 
orient<-function(M) {(svd(var(M))$v[1:k,1])} 	#find orientation
trajorient<-function(M){
traj.orient<-array(NA,dim=c(n,k))   		#loop across trajectories
check.1<-array(NA,dim=c(n))
for (i in 1:n){
  temp<-orient(M[,,i])
  traj.orient[i,]<-temp
  check.1[i]<-M[1,,i]%*%traj.orient[i,]  #check startingpoint location
  check.1[i]<-check.1[i]/abs(check.1[i])
  if(check.1[i]==-1) traj.orient[i,]<--1*traj.orient[i,]
}
options(warn=-1)				#b/c acos of 1 (e.g, on diagonal) yields warning
traj.ang.diff<-(180/pi)*acos(traj.orient%*%t(traj.orient))
#diag(traj.ang.diff)<-0
}

#trajectory shape
### GPA: following J. Claude 2008: Morphometrics in R
trans<-function(A){scale(A,scale=F)} 	##TRANSLATION
csize<-function(A)				##CSIZE
  {p<-dim(A)[1]
  size<-sqrt(sum(apply(A,2,var))*(p-1))
  list("centroid_size"=size,"scaled"=A/size)}
mshape<-function(A){apply(A,c(1,2),mean)}	#meanshape	

pPsup<-function(M1,M2){				## OPA rotation 1-->2
  k<-ncol(M1)
  Z1<-trans(csize(M1)[[2]])
  Z2<-trans(csize(M2)[[2]])
  sv<-svd(t(Z2)%*%Z1)
  U<-sv$v; V<-sv$u; Delt<-sv$d
  sig<-sign(det(t(Z1)%*%Z2))
  Delt[k]<-sig*abs(Delt[k]); V[,k]<-sig*V[,k]
  Gam<-U%*%t(V)
  beta<-sum(Delt)
  list(Mp1=beta*Z1%*%Gam,Mp2=Z2,rotation=Gam,scale=beta,
       df=sqrt(1-beta^2))}

pgpa<-function(A)
  {p<-dim(A)[1]; k<-dim(A)[2]; n<-dim(A)[3]  
  temp2<-temp1<-array(NA,dim=c(p,k,n)); Siz<-numeric(n)#; Qm2<-numeric(n)
  for (i in 1:n)
    {Acs<-csize(A[,,i])
    Siz[i]<-Acs[[1]]
    temp1[,,i]<-trans(Acs[[2]])}
  Qm1<-dist(t(matrix(temp1,k*p,n)))
  Q<-sum(Qm1); iter<-0
  while (abs(Q)> 0.00001)
    {for(i in 1:n){
      M<-mshape(temp1[,,-i])
      temp2[,,i]<-pPsup(temp1[,,i],M)[[1]]}
    Qm2<-dist(t(matrix(temp2,k*p,n)))
    Q<-sum(Qm1)-sum(Qm2)
    Qm1<-Qm2
    iter=iter+1
    temp1<-temp2}
  list("rotated"=temp2,"it.number"=iter,"Q"=Q,"intereucl.dist"=Qm2,"mshape"=
	csize(mshape(temp2))[[2]],"cent.size"=Siz)
}

## loop for GPA and shape distances
trajshape<-function(M){
  x<-pgpa(M)
  traj.shape.dist<-as.matrix(x$intereucl.dist) 
}

## TrajectorySummaryStat
sumstat<-function(M){
M<-as.dist(M)
x<-if(length(M)>1)(x=var(M)) else 0
}


##################MAIN LOOP################
traj.specs.obs<-arrayspecs(lsmeans.obs)
trajsize.obs<-trajsize(traj.specs.obs)
trajdir.obs<-trajorient(traj.specs.obs)
diag(trajdir.obs)<-0 #b/c some NA/N values on diagonal)
trajshape.obs<-trajshape(traj.specs.obs)
sumstatsize.obs<-sumstat(trajsize.obs)
sumstatdir.obs<-sumstat(trajdir.obs)
sumstatshape.obs<-sumstat(trajshape.obs)

### PERMUTATION PROCEDURE
permute<-999
PSize<-POrient<-PShape<-array(1,dim=c(n,n))
PSumSize<-PSumOrient<-PSumShape<-1
for(i in 1:permute){
	res.rand<-res.red[sample(nrow(res.red)),]	#Resample residuals 
	y.rand<-yhat.red+res.rand			#Residuals + predicted values
	yhat.rand<-predict(lm(y.rand~taxa*level,model=T,x=T,y=T,qr=T))
	lsmeans.rand<-NULL
	for (j in 1:k){
	  mean.temp<-tapply(yhat.rand[,j],taxaBylevel,mean)
	  lsmeans.rand<-cbind(lsmeans.rand,mean.temp)
	}
	traj.specs.rand<-arrayspecs(lsmeans.rand)
	trajsize.rand<-trajsize(traj.specs.rand)
	trajdir.rand<-trajorient(traj.specs.rand)
	diag(trajdir.rand)<-0 #b/c some NA/N values on diagonal)
	trajshape.rand<-trajshape(traj.specs.rand)
	sumstatsize.rand<-sumstat(trajsize.rand)
	sumstatdir.rand<-sumstat(trajdir.rand)
	sumstatshape.rand<-sumstat(trajshape.rand)

	PSize<-ifelse(trajsize.rand>=trajsize.obs, PSize+1,PSize) 
	POrient<-ifelse(trajdir.rand>=trajdir.obs,POrient+1,POrient) 
	PShape<-ifelse(trajshape.rand>=trajshape.obs,PShape+1,PShape) 

	PSumSize<-if(sumstatsize.rand>=sumstatsize.obs) (PSumSize+1) else PSumSize
	PSumSOrient<-if(sumstatdir.rand>=sumstatdir.obs) (PSumOrient+1) else PSumOrient
	PSumShape<-if(sumstatshape.rand>=sumstatshape.obs) (PSumShape+1) else PSumShape
}  #end permute

PSize<-PSize/(permute+1)
POrient<-POrient/(permute+1)
PShape<-PShape/(permute+1)

PSumSize<-PSumSize/(permute+1)
PSumOrient<-PSumOrient/(permute+1)
PSumShape<-PSumShape/(permute+1)

trajsize.obs
PSize
trajdir.obs
POrient
trajshape.obs
PShape

######NOTE: summary statistics only valid for 2+ trajectories!!
sumstatsize.obs
PSumSize
sumstatdir.obs
PSumOrient
sumstatshape.obs
PSumShape

#########Plot Data and LS means 
plot(Y[,1:2],type="n",xlab="PC I", ylab="PC II",asp=1)
points(Y[,1:2])
for (i in 1:n){		 	
  for (j in 1:(p-1)){		
    points(traj.specs.obs[(j:(j+1)),1,i],traj.specs.obs[(j:(j+1)),2,i],type="l",pch=21,col="black")    
  }
}
for (i in 1:n){		 	
  for (j in 2:(p-1)){		
    points(traj.specs.obs[j,1,i],traj.specs.obs[j,2,i],pch=21,bg="gray",col="black",cex=1.5)
  }
}
for (i in 1:n){
  points(traj.specs.obs[1,1,i],traj.specs.obs[1,2,i],pch=21,bg="white",col="black",cex=1.5)
}
for (i in 1:n){
  points(traj.specs.obs[p,1,i],traj.specs.obs[p,2,i],pch=21,bg="black",col="black",cex=1.5)
}