### Sinus deformation analysis

####Libraries####
library(RToolsForDeformetrica)
library(ade4)
library(ggplot2)
library(rgl)
library(geometry)
library(MASS)
library(Morpho)
####Loading deformation fields####
NumberOfSubject = 24# number of subjects (pooled)
KernelWidth = 0.1# kernel width of the deformation kernel

MomentaPathname = 'DeterministicAtlas__EstimatedParameters__Momenta.txt'
ControlPointsPathname = 'DeterministicAtlas__EstimatedParameters__ControlPoints.txt'

Momenta = import_momenta(MomentaPathname, NumberOfSubject)
ControlPoints = import_control_points(ControlPointsPathname)
velocities = compute_velocities(Momenta$momenta, ControlPoints, KernelWidth)
ColumnMean = colMeans(velocities$V)
velocitiesC = velocities$V-matrix(rep(ColumnMean,nrow(velocities$V)),nrow=nrow(velocities$V),ncol=ncol(velocities$V),byrow=TRUE)

####Additional data required for shape analysis####
names_data<-read.csv("names_sinus.csv",header=T,sep=",")

#grouping factor
ecology<-as.factor(as.character(names_data$Ecology[c(1:18,20:22)]))

#data for plots

forms<-as.numeric(names_data$Form)
cols<-as.character(names_data$Color)
fills<-as.character(names_data$Fill)

####shape analysis####
#PCA analysis
PCAsinus<-dudi.pca(velocitiesC,scannf=F,scale=F,center=F,nf=Inf)


PCAplotSinus<-data.frame(PCAsinus$li[,1:3],cols,forms,rownames(PCAsinus$li))
colnames(PCAplotSinus)<-c("PC1","PC2","PC3","colo","shapes","labels")
#PC2 vs PC1 plot
Sinus_p<-ggplot(PCAplotSinus,aes(x=PC1,y=PC2,color=colo))
Sinus_p+geom_point(size=5,colour=cols,fill=cols,shape=forms)+
  theme(panel.background = element_blank(),axis.line = element_line(),axis.text = element_text(size=12,colour = "black"))+coord_equal()
#PC3 vs PC1 plot
Sinus_p<-ggplot(PCAplotSinus,aes(x=PC1,y=PC3,color=colo))
Sinus_p+geom_point(size=5,colour=cols,fill=cols,shape=forms)+
  theme(panel.background = element_blank(),axis.line = element_line(),axis.text = element_text(size=12,colour = "black"))+coord_equal()

#3D plot
PC1<-PCAplotSinus[,1]
PC2<-PCAplotSinus[,2]
PC3<-PCAplotSinus[,3]

plot3d(PCAplotSinus[, 1], PCAplotSinus[, 2], PCAplotSinus[,3], type = "s", xlim = c(min(PCAplotSinus[, 1])*1.2,max(PCAplotSinus[, 1])*1.2), ylim = c(min(PCAplotSinus[, 2])*1.2,max(PCAplotSinus[, 2])*1.2), zlim = c(min(PCAplotSinus[, 3])*1.2,max(PCAplotSinus[, 3])*1.2), asp = 1, xlab = paste("PC1"),ylab = paste("PC2"), zlab = paste("PC3"),col = PCAplotSinus$colo, size = 0.7)

#add text to selected individuals
text3d(PCAplotSinus[23:24,1]+0.05,PCAplotSinus[23:24,2]+0.05,PCAplotSinus[23:24,3]+0.05,c("RSS45","FAM97057"),col = "black", cex = 1,alpha=1)

#3D convex hulls
psGroupHunt<-matrix(c(PC1[grepl("Group",names_data$Ecology)],PC2[grepl("Group",names_data$Ecology)],PC3[grepl("Group",names_data$Ecology)]), ncol=3)
psMeso<-matrix(c(PC1[grepl("Meso",names_data$Ecology)],PC2[grepl("Meso",names_data$Ecology)],PC3[grepl("Meso",names_data$Ecology)]), ncol=3)
ts.surfGroupHunt<-t(convhulln(psGroupHunt))
ts.surfMeso<-t(convhulln(psMeso))
rgl.triangles(	psGroupHunt	[	ts.surfGroupHunt	,1],	psGroupHunt	[	ts.surfGroupHunt	,2],	psGroupHunt	[	ts.surfGroupHunt	,3],col= "#8A0000",alpha=0.4)
rgl.triangles(	psMeso	[	ts.surfMeso	,1],	psMeso	[	ts.surfMeso	,2],	psMeso	[	ts.surfMeso	,3],col= "#BDBDBD",alpha=0.4)


#inspection of PCA data to establish the number of PCs to be included in the CVA
## Build dataset:
ncomp <- 23
dat <- cbind(Group = ecology, as.data.frame(PCAsinus$li[c(1:18,20:22), 1:ncomp]))
accuracies <- rep(NA, ncomp - 1)
names(accuracies) <- paste(2:ncomp, "axes")
for (k in 3:ncol(dat)) {
  mod <- lda(Group ~ ., data = dat[, 1:k],
             prior = rep(1/3, 3),
             CV = TRUE)
  accuracies[k - 2] <- sum(mod$class == dat$Group) / nrow(dat)
}
plot(x = 2:ncomp, y = accuracies,
     type = "b", pch = 16,
     xlab = "Number of PCA axes in LDA model",
     ylab = "X-val classification accuracy")
grid()


##PCA data to be used for CVA
priori=PCAsinus$li[c(1:18,20:22),1:7]
posteriori=PCAsinus$li[-c(1:18,20:22),1:7]

#CVA analysis with three groups
cva=CVA(priori, ecology, cv=T, prior=rep(1/3,3), weighting=T)
classify(cva)
#a posteriori individuals scores
x <- as.matrix(posteriori)
xscale <- sweep(x,2,cva$Grandm)
xscores <- xscale %*% cva$CV

##Function to estimate a posteriori probabilities for post hoc individuals
CVAposterior <- function(myCVA,x,prior=NULL) {
  if (length(dim(myCVA$groupmeans)) == 3) {
    myCVA$groupmeans <- vecx(myCVA$groupmeans)
    myCVA$Grandm <- c(myCVA$Grandm)
    x <- vecx(x)
  }
  x <- as.matrix(x)
  xscale <- sweep(x,2,myCVA$Grandm)
  xscores <- xscale %*% myCVA$CV
  
  if (is.null(prior))
    prior <- myCVA$prior
  GmeanCenter <- sweep(myCVA$groupmeans,2,myCVA$Grandm)
  GmeanScores <- GmeanCenter%*%myCVA$CV
  classprobs <- NULL
  classVec <- myCVA$groups
  for (i in 1:nrow(xscores)) {
    tmpdist <- (rowSums(sweep(GmeanScores,2,xscores[i,])^2))
    post <- Morpho:::probpost(tmpdist,prior)
    classVec[i] <- names(tmpdist)[which(post == max(post))]
    classprobs <- rbind(classprobs,post)
  }
  
  return(classprobs)
  
}

###A posteiori probabilities for post hoc individuals
cvapost=round(CVAposterior(cva,posteriori),2)
rownames(cvapost)<-names_data$ID[-c(1:18,20:22)]
colnames(cvapost)<-levels(ecology)
cvapost

###Typicality probabilities based on CVA scores

tmpprob <- typprobClass(x = xscores, ## the individual to consider
                        data = cva$CVscores, ## the reference data (CV scores)
                        groups = ecology ## the groups involved
                        ,method="chisquare", outlier = 0.05)

round(tmpprob$probs,3)