For what it’s worth… I have neither tried to fully penetrate Paolo’s code nor the technical details in the paper he sent, but I did notice his code used the default of gpagen to project Procrustes residuals into tangent space. Therefore, the density plot results appear to confound methodological differences and whether configurations were projected into tangent space. Updating that portion of the code (not projecting Procrustes residuals) did not seem to resolve anything, and I would not expect it, since the Procrustes distances were measured among specimens after alignment to the mean configuration via GPA, not alignment between pairs of configurations. However, an interesting observation. The “distortion” between distances from the pairwise approach Paolo provided code for and the distances from the GPA-aligned configurations is on par with the distortion one would expect from projection into tangent space; i.e., so small in this empirical example as to be rounding error if using more than 6 decimal places.
It is an interesting debatable point. Calculating Procrustes distances by pairs means that the same result will be found between pairs, irrespective of whether a full data set is expanded or subsetted (in which case the means could change, and thus, the Procrustes residuals could also change). Alternatively, without the mean, without projection into a tangent space, and without visualization of points in the space, perhaps with a mapping of shape change (transformation grids), I wonder how one sees the forest from the trees? Anyway, to confirm, yes, gpagen uses an “Euclidean fashion” to calculate distances based on Procrustes residuals that have been obtained via GPA and does not subset the data to all possible n(n - 1) / 2 pairs for n specimens, with an extra rotation step in each case. (I did confirm that after performing GPA with gpagen, and then performing it again on the Procrustes residuals obtained only for, e.g., specimens 1 and 2, specimens 1 and 3, …, I get the exact same results as obtained from Paolo’s code, meaning the pairwise approach is GPA performed on Procrustes residuals from GPA, by pairs.) So two important points from this. First, in gpagen, I should have been clear before that the Procrustes distances if calculated on projected residuals are actually Euclidean distances in the tangent space. Do not project if you wish to have Procrustes distances. Second, the method Paolo is outlining is not Procrustes distance from GPA. It is Procrustes distance based only on two shapes, which is tantamount to performing GPA, and subsequently realigning configurations as if all other n - 2 configurations do not matter. Whether the other n - 2 configurations do not matter might be interesting to discuss. Cheers! Mike > On Oct 17, 2022, at 11:51 AM, Paolo Piras <[email protected]> wrote: > > It is very probable [correct me if I'm wrong, but see the last line of code > below] that gpagen(), and other functions from other packages, computes > pairwise distances in the Euclidean fashion after gpa alignment. While there > is nothing wrong in the Euclidean fashion (in most cases) I think that these > distances still contain rotations, albeit very small. In fact...aligned > specimens are aligned with the consensus but not among them. At any pairwise > comparison one should align each pair. In fact, the concept of "alignment" > pertains to single pairs. This will have ***very little*** impact on the > majority of cases and my message is just a technicality. However, if one > wants to perform some kind of parallel transport in the shape space or size > and shape space he/she should be aware about this. In particular, see section > 3 in the attached paper. > Let's try this with the help of an ad hoc written function (that allows both > Euclidean/Kendall option as well as shape-space/size-and-shape-space option) > that is a little bit slow but works... > There are also some ancillary functions... > All the best > Paolo > > > Just select all and run > > library(geomorph) > library(Morpho) > library(shapes) > library(Matrix) > data(plethodon) > > array2mat<-function(array,inds,vars){ > > if(class(array)=="matrix"){array<-array(array,dim=c(nrow(array),ncol(array),1))} > X1 <-aperm(array,c(3,2,1)) > dim(X1)<- c(inds, vars) > > if(!is.null(dimnames(array)[3])){rownames(X1)<-unlist(dimnames(array)[3])}else{rownames(X1)<-c(1:nrow(X1))} > return(X1) > } > scaleshapes<-function(array){ > k<-dim(array)[1] > m<-dim(array)[2] > n<-dim(array)[3] > library(abind) > library(shapes) > library(Morpho) > divmat<-function(mat){ > mat2<-mat/cSize(mat) > mat2 > } > > if(is.matrix(array)==T){array<-array(array,dim=c(nrow(array),ncol(array),1))} > scaled<-vecx(t(apply(array,3,divmat)),revert=T,lmdim=m) > return(scaled) > } > > centroids<-function(array){ > meanmat<-function(mat){ > mean<-apply(mat,2,mean) > mean > } > > if(class(array)[1]=="matrix"){array<-array(array,dim=c(dim(array)[1],dim(array)[2],1))} > centroidi<-t(apply(array,3,meanmat)) > return(centroidi) > } > > centershapes<-function(array,vectors=NULL){ > > if(is.matrix(array)==T){array<-array(array,dim=c(nrow(array),ncol(array),1))} > if(is.null(vectors)==T){centros<-centroids(array)}else{centros=vectors} > k<-dim(array)[1] > m<-dim(array)[2] > n<-dim(array)[3] > prov<-array(rep(t(centros),each=k),dim=c(k,m,n)) > array2<-array-prov > return(array2) > } > eucdist<-function(a,b,scale=T,center=T,doopa=T){ > if(center==T){ > a<-centershapes(a)[,,1] > b<-centershapes(b)[,,1] > } > > if(doopa==T){ > theopa<-rotonto(a,b,scale=F,reflection=F) > a<-a > b<-theopa$yrot} > > k<-nrow(a) > m<-ncol(a) > > if(scale==T){ > a<-scaleshapes(array(a,dim=c(k,m,1)))[,,1] > b<-scaleshapes(array(b,dim=c(k,m,1)))[,,1] > } > > dist<-dist(rbind(array2mat(array(a,dim=c(k,m,1)),1,k*m),array2mat(array(b,dim=c(k,m,1)),1,k*m))) > dist > } > > paireucdist<-function(array,scale=T,doopa=T,euc=T,noneuc=c("kendall","ssriemdist")){ > confs<-combn(dim(array)[3],2) > if(!is.null(dimnames(array)[[3]])){ > contrlist<-NULL > for(i in 1:ncol(confs)){ > > contrveci<-paste(dimnames(array)[[3]][confs[,i][1]],dimnames(array)[[3]][confs[,i][2]],sep="/") > contrlist<-c(contrlist,list(contrveci)) > }} > res<-NULL > for(i in 1:ncol(confs)){ > > if(euc==T){resi<-eucdist(array[,,confs[,i][1]],array[,,confs[,i][2]],scale=scale,doopa=doopa)}else{ > > if(noneuc[1]=="kendall"){resi<-kendalldist(array[,,confs[,i][1]],array[,,confs[,i][2]])} > > if(noneuc[1]=="ssriemdist"){resi<-ssriemdist(array[,,confs[,i][1]],array[,,confs[,i][2]])} > } > res<-c(res,resi) > print(paste(i,"of",length(confs)/2,sep=" ")) > } > if(is.null(dimnames(array)[[3]])){cbind(t(combn(dim(array)[3],2)),res)}else{ > data.frame(unlist(contrlist),res,apply(t(confs),2,paste))} > } > Y.gpa <-gpagen(plethodon$land, PrinAxes = FALSE,verbose = T) > mat<-Y.gpa$procD > > ##Euclidean fashion > pairwisedist<-paireucdist(plethodon$land) > plot(pairwisedist[,3],mat,asp=1) ##### of course they are very close in this > case > plot(density(pairwisedist[,3]/mat)) #### not identical > ##Kendall fashion > pairwisedist2<-paireucdist(plethodon$land,euc=F,noneuc="kendall") > plot(pairwisedist2[,3],mat,asp=1) ##### of course they are very close in this > case > plot(density(pairwisedist2[,3]/mat)) #### not identical > > ### scholastic > plot(density(pairwisedist2[,3]/pairwisedist[,3])) #### not identical > > > ### > mat2<-dist(array2mat(Y.gpa$coords,40,24)) > range(mat2-mat) #### this confirms that "pairwise" distances in gpagen() are > computed without re-aligning each pair. > > Il giorno lun 17 ott 2022 alle ore 16:49 Dagmawit Getahun <[email protected] > <mailto:[email protected]>> ha scritto: > Hi Mike, > > Thank you so much! This was very helpful and I was able to get the distances. > > Best, > Dag > > On Sun, Oct 16, 2022 at 9:40 PM Mike Collyer <[email protected] > <mailto:[email protected]>> wrote: > Dear Dag, > > When running the GPA with the gpagen function in geomorph, use the argument > verbose = TRUE and it will include the Procrustes distances among specimens > as output. This can take a while with large data sets, which is why it is not > set as a default. > > Cheers! > Mike > > Sent from my iPhone > >> On Oct 16, 2022, at 9:33 PM, Dagmawit Getahun <[email protected] >> <mailto:[email protected]>> wrote: >> >> Hi all, >> >> I have a dataset of 3D cranial landmarks in the Morphologika file format as >> a text file. I was able to upload it and run a GPA as well as PCA in R using >> geomorph and I wanted to know if there was an R code to analyze the actual >> Procrustes distances between the specimens in my dataset and whether I can >> do that in the format that my data set is in? >> >> Thanks in advance for the help! >> >> Best, >> Dag >> >> -- >> You received this message because you are subscribed to the Google Groups >> "Morphmet" group. >> To unsubscribe from this group and stop receiving emails from it, send an >> email to [email protected] >> <mailto:[email protected]>. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/morphmet2/a9136b76-7691-4d80-9096-6c59afdfe8a4n%40googlegroups.com >> >> <https://groups.google.com/d/msgid/morphmet2/a9136b76-7691-4d80-9096-6c59afdfe8a4n%40googlegroups.com?utm_medium=email&utm_source=footer>. > > > -- > Dagmawit Abebe Getahun > Doctoral Candidate > Physical Anthropology > The Graduate Center , CUNY > New York Consortium in Evolutionary Primatology > > Adjunct Lecturer, Anthropology Department > Lehman College, CUNY > > > -- > You received this message because you are subscribed to the Google Groups > "Morphmet" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected] > <mailto:[email protected]>. > To view this discussion on the web visit > https://groups.google.com/d/msgid/morphmet2/CALLJY8%2BejP-9F0-V-%3Dq4nSFADtBcbxb%3DOFVKhdLo2Y%3DE%2B7UJzA%40mail.gmail.com > > <https://groups.google.com/d/msgid/morphmet2/CALLJY8%2BejP-9F0-V-%3Dq4nSFADtBcbxb%3DOFVKhdLo2Y%3DE%2B7UJzA%40mail.gmail.com?utm_medium=email&utm_source=footer>. > <Varano et al BIOMAT.pdf> -- You received this message because you are subscribed to the Google Groups "Morphmet" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/morphmet2/2F383CC2-0B2D-4F5E-9926-586726A4AFE3%40gmail.com.
