Answers.Body: <p>I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.</p>
  <p>The number of constrained axes is a parameter chosen as an input to the <a href="https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/0012-9658%282003%29084%5B0511%3ACAOPCA%5D2.0.CO%3B2" rel="nofollow noreferrer">CAPscale method</a>:</p>
  <blockquote>
  <p>An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for <span class="math-container">$m$</span>? One might be concerned that the method could suffer from a poor choice for <span class="math-container">$m$</span>. If <span class="math-container">$m$</span> is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if <span class="math-container">$m$</span> is too large relative to <span class="math-container">$N$</span>, then a misleading canonical plot could result. For example, we cannot set <span class="math-container">$m = N - 1$</span>, for then the squared canonical correlations would all be equal to 1.0, because we have <span class="math-container">$N$</span> points!</p>
  </blockquote>
  <p>So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.</p>
  <p>I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute  p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:</p>
  <blockquote>
  <p>As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).</p>
  </blockquote>
  <p>This seems to be <a href="https://rdrr.io/rforge/vegan/man/anova.cca.html" rel="nofollow noreferrer">directly implemented as the <code>anova.cca()</code> function</a> in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.</p>
  <p>For how others have approached similar questions you can see <a href="https://www.researchgate.net/post/How-to-analyse-the-effect-of-a-factor-on-species-distribution-and-a-species-functional-trait" rel="nofollow noreferrer">this discussion</a>. For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). <a href="https://rdrr.io/rforge/vegan/man/capscale.html" rel="nofollow noreferrer">Here is the CAPscale function</a> in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so <a href="https://jacobrprice.github.io/2017/08/26/phyloseq-to-vegan-and-back.html" rel="nofollow noreferrer">I think that you can pass phyloseq back into vegan functions</a> and e.g. use those analyses directly on your <code>cap_ord</code> object. But see that link for more info.</p>
  <p>For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see <a href="https://si.biostat.washington.edu/sites/default/files/modules/Lecture6_0.pdf" rel="nofollow noreferrer">here</a>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-20T21:48:47.527Z
Answers.Text: Thanks, I will try to understand these all.
Answers.UserDisplayName: bioinfonext
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-20T20:05:49.670Z
Answers.LastActivityDate: 2021-01-20T20:05:49.670Z
Answers.UserDisplayName: Maximilian Press
Body: <p>I used the below code to do CAP analysis to link microbial data with elements, but not sure how to do statistical analysis to show which element variation is responsible most for the difference in microbial composition. I can only able to see two axis in figure so how I can get complete information about the variation based on these mineral element and individually like which mineral responsible most for microbial diversity?</p>
  <pre><code>    # CAP ordinate
      cap_ord &lt;- ordinate(
          physeq = phyloseq_prop, 
          method = &quot;CAP&quot;,
          distance = bray_not_na,
          formula = ~ ~ P + Mn + Cu +Zn + Rb + Mg )
      
          
          # CAP plot
      cap_plot &lt;- plot_ordination(
        physeq =phyloseq_prop , 
        ordination = cap_ord, 
          color = &quot;SampleType&quot;, 
          axes = c(1,2)
      ) + 
          aes(shape = Region) + 
          geom_point(aes(colour = SampleType), alpha = 0.4, size = 4) + 
          geom_point(colour = &quot;grey90&quot;, size = 1.5) + 
          scale_color_manual(values = c(&quot;#CBD588&quot;, &quot;#5F7FC7&quot;, &quot;orange&quot;,&quot;#DA5724&quot;, &quot;#508578&quot;, &quot;#CD9BCD&quot;,
         &quot;#AD6F3B&quot;, &quot;#673770&quot;,&quot;#D14285&quot;,&quot;#652926&quot;, &quot;#C84248&quot;, 
        &quot;#8569D5&quot;)
          )
          
          # Now add the environmental variables as arrows
      arrowmat &lt;- vegan::scores(cap_ord, display = &quot;bp&quot;)
      
      # Add labels, make a data.frame
      arrowdf &lt;- data.frame(labels = rownames(arrowmat), arrowmat)
      
      # Define the arrow aesthetic mapping
      arrow_map &lt;- aes(xend = CAP1, 
          yend = CAP2, 
          x = 0, 
          y = 0, 
          shape = NULL, 
          color = NULL, 
          label = labels)
      
      label_map &lt;- aes(x = 1.3 * CAP1, 
          y = 1.3 * CAP2, 
          shape = NULL, 
          color = NULL, 
          label = labels)
      
      arrowhead = arrow(length = unit(0.02, &quot;npc&quot;))
      
      pdf(&quot;cap2.pdf&quot;, height=8, width=8)
      # Make a new graphic
      cap_plot + 
        geom_segment(
          mapping = arrow_map, 
          size = .5, 
          data = arrowdf, 
          color = &quot;gray&quot;, 
          arrow = arrowhead
        ) + 
        geom_text(
          mapping = label_map, 
          size = 4,  
          data = arrowdf, 
          show.legend = FALSE
        )
        
        dev.off();
        
          
       cap_ord
      
       Call: capscale(formula = distance ~  P + Mn + Cu +Zn + Rb + Mg, data = data)
  
                                    Inertia Proportion Rank
  Total          2.069e+01  1.000e+00
  Constrained    5.240e+00  2.533e-01   10
  Unconstrained  1.545e+01  7.468e-01   97
  Imaginary     -1.090e-03 -5.268e-05    1
  Inertia is squared Bray distance
  
  Eigenvalues for constrained axes:
    CAP1   CAP2   CAP3   CAP4   CAP5   CAP6   CAP7   CAP8   CAP9  CAP10
  2.5733 0.8789 0.5647 0.3797 0.3037 0.1873 0.1243 0.1030 0.0800 0.0453
  
  Eigenvalues for unconstrained axes:
   MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8
  4.790 1.343 1.138 0.916 0.784 0.649 0.518 0.501
  (Showing 8 of 97 unconstrained eigenvalues)
    
  </code></pre>
  <p>Many thanks</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-19T20:47:38.650Z
Comments.Text: without some information it may be hard to answer. For example, you mention CCA but in your function call to `ordinate()` you are actually using a different technique `method="CAP"`. It might help to actually see the plotted output, to understand whether the issue is in data or a difficulty of the tool. It looks like you are using phyloseq, have you looked at the documentation? https://www.rdocumentation.org/packages/phyloseq/versions/1.16.2/topics/ordinate
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-19T20:49:27.523Z
Comments.Text: If you are just interested in something like a hypothesis test, you might consider using the analysis of similarity (ANOSIM): https://en.wikipedia.org/wiki/Analysis_of_similarities. You don't really need to do ordination for this, it can just take the dissimilarity matrix.
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-19T22:06:01.517Z
Comments.Text: Thanks, sorry for the confusion, I wanted to ask about the CAP analysis. I just got two axis in figure explaining the total 17% variability. I tried the to capture the output of cap_ord as updated in the script but it is also not showing about complete variability explained by these all mineral element
Comments.UserDisplayName: bioinfonext
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-20T14:27:33.983Z
Comments.Text: Can I say around 25% variability explained by these variables,  from above mention data;Constrained 5.240e+00 2.533e-01
Comments.UserDisplayName: bioinfonext
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-20T16:06:26.390Z
Comments.Text: having the plot would help us to better understand your question and guide you better. Also if you could expand the data returned. Does it return any measure of variability or weight for each of your variables? (This is not shown on the output you pasted, but not sure if there is something else it isn't printed or simply it doesn't return that information). From what package is the ordinate function?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-20T19:09:30.793Z
Comments.Text: Thanks, I have updated the result
Comments.UserDisplayName: bioinfonext
Tags: <r><microbial-genomics><phyloseq>
Title: Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements
ThreadText: # Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements
  
  (Asked by: username_0 on 2021-01-19)
  
  I used the below code to do CAP analysis to link microbial data with elements, but not sure how to do statistical analysis to show which element variation is responsible most for the difference in microbial composition. I can only able to see two axis in figure so how I can get complete information about the variation based on these mineral element and individually like which mineral responsible most for microbial diversity?
  
  ```
      # CAP ordinate
      cap_ord <- ordinate(
          physeq = phyloseq_prop, 
          method = "CAP",
          distance = bray_not_na,
          formula = ~ ~ P + Mn + Cu +Zn + Rb + Mg )
      
          
          # CAP plot
      cap_plot <- plot_ordination(
        physeq =phyloseq_prop , 
        ordination = cap_ord, 
          color = "SampleType", 
          axes = c(1,2)
      ) + 
          aes(shape = Region) + 
          geom_point(aes(colour = SampleType), alpha = 0.4, size = 4) + 
          geom_point(colour = "grey90", size = 1.5) + 
          scale_color_manual(values = c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
         "#AD6F3B", "#673770","#D14285","#652926", "#C84248", 
        "#8569D5")
          )
          
          # Now add the environmental variables as arrows
      arrowmat <- vegan::scores(cap_ord, display = "bp")
      
      # Add labels, make a data.frame
      arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
      
      # Define the arrow aesthetic mapping
      arrow_map <- aes(xend = CAP1, 
          yend = CAP2, 
          x = 0, 
          y = 0, 
          shape = NULL, 
          color = NULL, 
          label = labels)
      
      label_map <- aes(x = 1.3 * CAP1, 
          y = 1.3 * CAP2, 
          shape = NULL, 
          color = NULL, 
          label = labels)
      
      arrowhead = arrow(length = unit(0.02, "npc"))
      
      pdf("cap2.pdf", height=8, width=8)
      # Make a new graphic
      cap_plot + 
        geom_segment(
          mapping = arrow_map, 
          size = .5, 
          data = arrowdf, 
          color = "gray", 
          arrow = arrowhead
        ) + 
        geom_text(
          mapping = label_map, 
          size = 4,  
          data = arrowdf, 
          show.legend = FALSE
        )
        
        dev.off();
        
          
       cap_ord
      
       Call: capscale(formula = distance ~  P + Mn + Cu +Zn + Rb + Mg, data = data)
  
                                    Inertia Proportion Rank
  Total          2.069e+01  1.000e+00
  Constrained    5.240e+00  2.533e-01   10
  Unconstrained  1.545e+01  7.468e-01   97
  Imaginary     -1.090e-03 -5.268e-05    1
  Inertia is squared Bray distance
  
  Eigenvalues for constrained axes:
    CAP1   CAP2   CAP3   CAP4   CAP5   CAP6   CAP7   CAP8   CAP9  CAP10
  2.5733 0.8789 0.5647 0.3797 0.3037 0.1873 0.1243 0.1030 0.0800 0.0453
  
  Eigenvalues for unconstrained axes:
   MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8
  4.790 1.343 1.138 0.916 0.784 0.649 0.518 0.501
  (Showing 8 of 97 unconstrained eigenvalues)
    
  
  ```
  
  Many thanks
  
  **Tags:** r, microbial-genomics, phyloseq
  
  > **Comments:**
  > **username_1**: without some information it may be hard to answer. For example, you mention CCA but in your function call to `ordinate()` you are actually using a different technique `method="CAP"`. It might help to actually see the plotted output, to understand whether the issue is in data or a difficulty of the tool. It looks like you are using phyloseq, have you looked at the documentation? https://www.rdocumentation.org/packages/phyloseq/versions/1.16.2/topics/ordinate  
  > **username_1**: If you are just interested in something like a hypothesis test, you might consider using the analysis of similarity (ANOSIM): https://en.wikipedia.org/wiki/Analysis_of_similarities. You don't really need to do ordination for this, it can just take the dissimilarity matrix.  
  > **username_0**: Thanks, sorry for the confusion, I wanted to ask about the CAP analysis. I just got two axis in figure explaining the total 17% variability. I tried the to capture the output of cap_ord as updated in the script but it is also not showing about complete variability explained by these all mineral element  
  > **username_0**: Can I say around 25% variability explained by these variables,  from above mention data;Constrained 5.240e+00 2.533e-01  
  > **username_2**: having the plot would help us to better understand your question and guide you better. Also if you could expand the data returned. Does it return any measure of variability or weight for each of your variables? (This is not shown on the output you pasted, but not sure if there is something else it isn't printed or simply it doesn't return that information). From what package is the ordinate function?  
  > **username_0**: Thanks, I have updated the result  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-01-20. Score: 1 | Accepted answer
  
  I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.
  
  The number of constrained axes is a parameter chosen as an input to the [CAPscale method](https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/0012-9658%282003%29084%5B0511%3ACAOPCA%5D2.0.CO%3B2):
  
  > An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for $m$? One might be concerned that the method could suffer from a poor choice for $m$. If $m$ is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if $m$ is too large relative to $N$, then a misleading canonical plot could result. For example, we cannot set $m = N - 1$, for then the squared canonical correlations would all be equal to 1.0, because we have $N$ points!
  
  So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.
  
  I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:
  
  > As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).
  
  This seems to be [directly implemented as the `anova.cca()` function](https://rdrr.io/rforge/vegan/man/anova.cca.html) in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.
  
  For how others have approached similar questions you can see [this discussion](https://www.researchgate.net/post/How-to-analyse-the-effect-of-a-factor-on-species-distribution-and-a-species-functional-trait). For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). [Here is the CAPscale function](https://rdrr.io/rforge/vegan/man/capscale.html) in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so [I think that you can pass phyloseq back into vegan functions](https://jacobrprice.github.io/2017/08/26/phyloseq-to-vegan-and-back.html) and e.g. use those analyses directly on your `cap_ord` object. But see that link for more info.
  
  For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see [here](https://si.biostat.washington.edu/sites/default/files/modules/Lecture6_0.pdf).
  
  > **Comments:**
  > **username_0**: Thanks, I will try to understand these all.  
  
  ---
  
  