0

I work on a study comparing the amino acid profile of 5 sample groups with 3 replicates in each group. This is how the data looks like except that there are 17 columns (one for each amino acid) in total. The values represent the measured level (in micromole per gram sample, not percentage of total) of individual amino acid in each sample. 

Group Asp   Asn   Glu   Gln   His  ...
  A    13    8     9    15    12
  A    5     7     6    11     1
  A    10   15     1    14    13
  B    13   10    11     5     7
  B    12    5    10     9     4
  B    6     2    13     1     5
  C    9     4     3    15     1 
  C    2     4    10     1     5
  C    5     9     8     7    14
  D    5    13     2    10     6
  D    7    13     9     1     5
  D    13   12     5     8    15
  E    15    7     4     2     1
  E    1     2    15    14     8
  E    1     3     2    13     5

I used an ANOVA to compare the mean sum of total amino acids between groups, then I ran a PCA to visualize differences in profile. I was asked by a reviewer to also analyze for significant differences in individual amino acids.

I generally find PCAs sufficient for this type of analysis. Besides, the number of samples per group is small (n = 3) and there may be correlation between the response variables (i.e. individual amino acid content) meaning the results of multiple univariate ANOVAs (for each amino acid) may be misleading (inflating type I error).

Correlation seems to be an issue for multivariate ANOVA (MANOVA) as well. In this specific case a MANOVA would be rank defficient anyway.

But there are examples of comparable data published in the scientific literature where multiple univariate ANOVAs are used to detect difference in individual amino acids between groups.

Does anyone has good arguments for or against using multiple univariate ANOVAs or other quantitative methods in such cases?

Drosof
  • 65
  • 7
  • Are the values in the table in units like micromoles per sample, or are they percentages of each individual amino acid with respect to the total of all amino acids? Also, out of curiosity, why are there only 17 columns? Particularly interesting as you have separate values for Asp and Asn and for Glu and Gln, the pairs that some analytical methods can't distinguish. – EdM Jan 30 '20 at 19:37
  • Yes the values are given as individual content (like I micromole per gram sample) and not in percentage of the total. I edited the question to clarify. Each column correspond to an amino acid detected in at least one of the samples. There are some HPLC methods to distinguish Asp/Glu from Asn/Gln. – Drosof Jan 30 '20 at 22:49

1 Answers1

1

Your approach, starting with ANOVA to compare total amino acid content among groups and then using PCA to illustrate underlying differences in composition, seems sound. Your reluctance to just go forward with 17 univariate ANOVAs is well founded. Besides not accounting for the possibility of correlations among the response variables, the simple fact of performing 17 univariate ANOVAs will require substantial correction for multiple comparisons.

The simplest way forward might be to just do the requested 17 univariate ANOVAs anyway (after all, if the reviewer wants to know the results then other readers probably would, too). You do this, however, while highlighting the difficulties involved in interpreting the results. That would get you quickly beyond the issue raised by the reviewer, and could even give you the opportunity to demonstrate the superiority of your PCA-based approach.

For MANOVA, it's not immediately clear to me why you think that it would be rank deficient, unless the analyses were done on percentages rather than raw values. (Maybe I just haven't thought that through enough.) The literature on handling necessarily rank-deficient compositional data and associated threads on CrossValidated could provide guidance.

You could consider multivariate partial least squares (or "projection to latent structures"), a technique frequently used in chemometrics, which I think in your case would be closely related to your PCA approach.

Note that the problem of having multiple correlated response variables is even more dramatic in gene-expression studies, which can have on the order of 20,000 response variables with only a handful of different conditions and replicates. Although it might be overkill for your application, you might consider exploring how information is pooled among observations and how the multiple comparisons problem is dealt with by packages like limma for differential expression analysis.

EdM
  • 57,766
  • 7
  • 66
  • 187
  • I guess running `manova()` on 17 response variables with only 3 values per group is the reason why the model is rank deficient. I also ran a PLSR on this data. – Drosof Feb 02 '20 at 19:18
  • In the case I would perform the multiple univariate ANOVAs, I assume the threshold _p_-value from the Bonferroni correction would be 0.05/17 = 0,00294. Is it the treshold one would look at when running the multiple comparisons using multiple `pairwise.t.test(aa, Group, p.adjust.method="bonferroni")`? As you mention it still doesn't account for possible correlations so perhaps I should just argue for not using this method. – Drosof Feb 02 '20 at 19:54
  • 1
    @Drosof your Bonferroni p value is OK, but that correction is extremely conservative. The Holm method would be better for controlling family-wise error rate. You could consider controlling false-discovery rate instead (Benjamini-Hochberg), as used in large-scale gene-expression studies. For dealing with your reviewer, I would recommend placing the results of the 17 individual ANOVAs and discussion of the difficulties in interpreting the results in the paper itself or in a supplement; that way you address the reviewer's concern without a back-and-forth argument. – EdM Feb 03 '20 at 15:09
  • Thanks for your thourough answer and comments. What about the [Benjamini-Yekutieli](https://projecteuclid.org/download/pdf_1/euclid.aos/1013699998) correction ? It accounts for correlations among response variables. Following [that post](https://stats.stackexchange.com/questions/63441/what-are-the-practical-differences-between-the-benjamini-hochberg-1995-and-t), I'm not sure `"BY"` method is appropriate in this case though. – Drosof Feb 04 '20 at 08:51
  • @Drosof I don't have any experience with the B-Y correction. My understanding is that gene-expression analysis uses B-H despite the inherent correlations among expression of 10,000 to 20,000 genes. The R `p.adjust()` function allows for either. – EdM Feb 04 '20 at 15:16