Hello, I post this question in geomorph forum as well ( https://groups.google.com/forum/#!topic/geomorph-r-package/K_fggqrh200) because this question is a direct extension of several geomorph-related questions I asked there earlier. I post here since this question essentially goes beyond geomorph and is about coding in R to derive individual magnitude of fluctuating asymmetry (FA).
A recent paper by Gunz (Sci Adv. 2020;6(7):eaax9935) suggested that magnitude of FA can be quantified as the Procrustes distance between each individual and its relabled reflection after directional asymmetry (DA) has been removed. My purpose is to translate this sentence into R code. Gunz cited Mardia and Bookstein's paper "Statistical assessment of bilateral symmetry of shapes", where Equation 3.2 below is about decomposition of total asymmetry into DA and FA: [image: Snipaste_2020-04-18_18-11-26.png] In this formula, X^P and Y^P are Procrustes shape coordinates for original and relabled reflection of all specimens, respectively. X_bar and Y_bar are the mean shape of all X^P and Y^P, respectively. D^2 is squared Procrustes distance between two configurations, which is calculated in squared Euclidean distance. It can be calculated using the *dist* function. Using the scallop data from geomorph (5 specimens in total which are of object symmetry), I tried to calculate magnitude of FA for the 1st specimen in the dataset following Mardia's equation. You may directly copy and paste the code to replicate my results. *library(Morpho)library(geomorph)data(scallops)############################### Code adapted from bilat.symmetry source code############################A <- gpagen(scallops$coorddata)size <- A$CsizeA <- A$coordsn <- dim(A)[3]k <- dim(A)[2]p <- dim(A)[1]ind <- scallops$indnind <- nlevels(ind)land.pairs <- scallops$land.pairsnpairs <- nrow(land.pairs)A2 <- AA2[land.pairs[, 1], , ] <- A[land.pairs[, 2], , ]A2[land.pairs[, 2], , ] <- A[land.pairs[, 1], , ]A2[, 1, ] <- A2[, 1, ] * (-1)A_long <- array(c(A, A2), c(p, k, 2 * n))ind <- factor(rep(ind, 2))side <- gl(2, n)gpa.res <- gpagen(A_long)# A_long below contains Procrustes shape coordinates for all original specimens and # their relabled reflections. Ordering of specimens is listed belowA_long <- gpa.res$coords* *# Ordering of specimens in A_long* *# row A_long * *# 1 ind1 * *# 2 ind2 * *# 3 ind3 * *# 4 ind4 * *# 5 ind5 * *# 6 ind1RelabledReflection* *# 7 ind2RelabledReflection* *# 8 ind3RelabledReflection* *# 9 ind4RelabledReflection* *# 10 ind5RelabledReflection* *############################ Total asymmetry for the 1st specimen ####################################TotAsym_score_s1 <- dist(two.d.array(A_long[,,c(1, 6)]))############################ Magnitude of DA for the 1st specimen ####################################mean_Original <- mshape(A_long[,,c(1:5)])mean_RelabRefl <- mshape(A_long[,,c(6:10)])OriginalRelabRefl <- bindArr(mean_Original, mean_RelabRefl, along = 3)DA_score <- dist(two.d.array(OriginalRelabRefl))############################ FA score calculation using 2 methods ############################**########* *# Method 1 for calculating FA score: total asymmetry minus DAFA_score1a_s1 <- sqrt(TotAsym_score_s1^2 - DA_score^2)FA_score1b_s1 <- TotAsym_score_s1 - DA_score# Method 2 for calculating FA score: by Madria's algorithm for FAdiff_Orginal_meanOriginal_s1 <- A_long[,,1] - mean_Originaldiff_RelabRefl_meanRelabRefl_s1 <- A_long[,,6] - mean_RelabRefldiff_both_s1 <- bindArr(diff_Orginal_meanOriginal_s1, diff_RelabRefl_meanRelabRefl_s1, along = 3)FA_score2_s1 <- dist(two.d.array(diff_both_s1))* Method 1 obtained FA by total asymmetry minus DA. Two FA scores are listed since I am not sure whether squared distance or distance should be used during calculation. Method 2 obtained FA according to the formula for FA shown in the figure above. Method 1 and 2 should certainly give identical FA value if calculated correctly. However, *FA_score1a_s1, FA_score1b_s1, and FA_score2_s1 are different*.* I must be wrong somewhere but I cannot figure out*. I know MorphoJ could do the calculation automatically. But I wish to know the steps behind. Coding by myself helps. I guess I am close to the answer. It would be very much appreciated if anyone could comment on the code, or directly answer by coding. Thank you! Best regards, YF Wen -- 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/da599158-d388-4e66-845e-daa756c7f486%40googlegroups.com.
