vegandevs / vegan

R package for community ecologists: popular ordination methods, ecological null models & diversity analysis
https://vegandevs.github.io/vegan/
GNU General Public License v2.0
453 stars 98 forks source link

Procrustes analysis #408

Open Olga2419 opened 3 years ago

Olga2419 commented 3 years ago

I am trying to find out if there is a correlation between two data sets (if one of them, which is called "varechem", depends on the other, which is called "varespec"). I ran the following code and this is what I got. I have never done such an analysis, so I would be grateful if someone could confirm that the analysis I have done is correct and can help me with its interpretation.

library(vegan) varespec <- read.csv(file = "antib_av.csv", row.names=1, header=TRUE) varespec

   ENR    LVX    CPR    ENX   FMQ    NFX  MTZ   OTC    AZM    CLR    ERY   JOS    

Jan 4.260 157.47 64.40 3.990 8.33 21.370 1.45 0.000 731.43 31.670 5.180 1.46 Feb 3.370 222.36 57.04 1.920 13.89 17.090 2.57 0.000 0.00 98.650 22.180 5.19 Mar 2.280 226.86 46.23 1.550 19.32 15.980 0.71 0.000 0.00 51.580 14.170 2.21 Jul 0.094 295.76 99.93 0.264 5.20 26.010 0.42 5.950 251.94 12.200 9.340 0.44 Aug 3.810 266.13 54.44 0.039 2.78 18.400 3.28 0.128 324.91 3.710 0.050 2.49 Sep 0.167 240.51 55.39 0.236 0.00 16.400 3.43 6.280 331.28 0.062 0.028 4.10 Nov 1.910 156.19 38.58 0.322 3.19 10.390 1.98 0.249 244.62 5.220 0.191 3.10 Oct 14.670 376.19 288.63 26.050 0.00 54.460 2.05 0.000 659.38 15.250 8.450 3.55 Dec 0.000 174.27 40.87 0.000 0.00 0.518 6.84 0.000 526.03 0.273 0.187 10.66

>vare.mds <- metaMDS(varespec, trace = FALSE)
>vare.mds

Call: metaMDS(comm = varespec, trace = FALSE)

global Multidimensional Scaling using monoMDS Data: wisconsin(sqrt(varespec)) Distance: bray

Dimensions: 2 Stress: 0.05444056 Stress type 1, weak ties Two convergent solutions found after 20 tries Scaling: centring, PC rotation, halfchange scaling Species: expanded scores based on ‘wisconsin(sqrt(varespec))’

>varechem<- read.csv(file = "int_av.csv", row.names=1, header=TRUE)
>varechem

     Intl1     Intl3

Jan 149679690 0 Feb 1095036490 54966409 Mar 615081914 29106165 Jul 933975403 165422426 Aug 1186102453 204000819 Sep 2687029085 539090113 Nov 1295730991 271652788 Oct 320871912 115514535 Dec 4031716781 430955068

rankindex(scale(varechem), varespec, c("euc","man","bray","jac","kul"))

    euc         man        bray         jac         kul 

0.01261261 -0.02265122 0.02599743 0.02599743 0.07078507

pc <- prcomp(varechem, scale = TRUE) ###PCA was performed with data (int) pc<- scores(pc, display = "sites", choices = 1:2) edis <- vegdist(pc, method = "bray") vare.dis <- vegdist(wisconsin(sqrt(varespec))) mantel(vare.dis, edis)

Mantel statistic based on Pearson's product-moment correlation

Call: mantel(xdis = vare.dis, ydis = edis)

Mantel statistic r: 0.4551 Significance: 0.007

Upper quantiles of permutations (null model): 90% 95% 97.5% 99% 0.288 0.340 0.383 0.428 Permutation: free Number of permutations: 999

> pc <- scores(pc, choices = 1:2) 
> pro <- protest(vare.mds, pc, scores = "sites", permutations = 999)
>plot(pro)

pro

Call: protest(X = vare.mds, Y = pc, scores = "sites", permutations = 999)

Procrustes Sum of Squares (m12 squared): 0.597 Correlation in a symmetric Procrustes rotation: 0.6349 Significance: 0.034 Permutation: free Number of permutations: 999

summary(pro, digits = getOption("digits"))

Call: protest(X = vare.mds, Y = pc, scores = "sites", permutations = 999)

Number of objects: 9 Number of dimensions: 2

Procrustes sum of squares:
0.5969577 Procrustes root mean squared error: 0.2575435 Quantiles of Procrustes errors: Min 1Q Median 3Q Max 0.09225553 0.18457338 0.25827092 0.28743807 0.41773741 Rotation matrix: [,1] [,2] [1,] 0.9678393 0.2515694 [2,] -0.2515694 0.9678393

Translation of averages: [,1] [,2] [1,] 2.187383e-18 4.797596e-18

Scaling of target: [1] 0.6348561

jarioksa commented 3 years ago

I haven't written Procrustes analysis for such an analysis, but I know people use it for the purpose. I am not qualified to say if that makes sense. Please consult the paper on Protest cited in the ?protests help for an authoritative answer.