post hoc pairwise phylogenetic comparisons

[Daniel Fernandes da Silva]

Dear all,

I performed an ANOVA from OLS (function procD.lm) and also a PGLS (procD.pgls) with the same database and in both cases I got a significant result for shape variation among groups with distinct habitat use. Subsequently, I chose the most appropriate model from PGLS results (in this case a "fit.common") and performed paired comparisons (function pairwise) between the groups to find out which one actually differed from each other. However, researching on the application of these analysis, I recovered a discussion in this group (back in 2016):  https://groups.google.com/g/geomorph-r-package/c/wLLhEgjEPrE?pli=1

stating that, in that time, there was no function in Geomorph to do this for phylogenetic data because the statistical properties of phylogenetic comparative methods applied to pairwise tests were not evaluated. As far as I could search, I could not find any updates on this, including the manuals of recent versions of RRPP and Geomorph packages or some published work. Based on this, I followed Adam's suggestion in that 2016 discussion and performed paired comparisons using the proper model of the ANOVA from OLS and compared it with the results of PGLS paired comparisons. The significant ANOVA from PGLS recovered non significant paired comparisons FOR ALL "test.types" ("dist", "VC", "DL"), which is strange to me since PGLS was significant for habitat use. On the other hand, the significant ANOVA from OLS led to up to 3 out of 6 significant pairwise comparisons depending on the "test.type". 

Sorry for the huge text, but I would like to know if: (1) is there any methodological approach to lead with these questions that I missed? (2) if my approach is correct, should I consider the results of OLS or PGLS paired comparisons? 

I would really appreciate it if anyone could help me with these issues.

Sincerely,

Daniel

[Mike Collyer]

Dear Daniel,

Wow, 8 years passes so fast!  Yes, functions have been updated and yes, you can use pairwise with a PGLS fit.  There is nothing wrong with doing this.  But to answer the other question, no, I do not think theoretical advances have been made in the area of pairwise comparisons for GLS models.  The reason why this is challenging, but isn’t for OLS models, is that although it is easy to make predictions with a GLS model (estimate GLS means), and compare those means, it is not possible to do this with a good sampling distribution with respect to appropriate standard error of the estimates.

Hope fully I can explain this.  The residual covariance for any linear model can be calculated as V = R^T C^-1 R, where ^t means matrix transpose, ^-1 means matrix inverse, R is the matrix of residuals, and C is the matrix of phylogenetic covariances.  R is found from GLS estimation.  C is a N x N matrix for N species.  But if we let C = I, an identity matrix (independence), R becomes OLS residuals instead of GLS residuals.  Algebra still works the same.  We can estimate the standard error of the model fit from V.  However, if we make estimates for g groups rather than N species, there is no way to include C in the calculation, because g != N.  This is of little concern for OLS estimation, because a g x g identity matrix can be used instead of an N x N identity matrix.

With pairwise comparisons, the RRPP distribution of pairwise mean differences is consistent with the distribution of F statistics in ANOVA.  I can show this below.  Notice the bolded P-values.

> library(phytools)

> library(RRPP)

> tree <- pbtree(n= 20)

> y <- fastBM(tree)

> g <- factor(rep(1:2, each = 10))

> names(g) <- names(y)

> y[1:10] <- y[1:10] + 2

> fit.ols <- lm.rrpp(y ~ g)

> fit.gls <- lm.rrpp(y ~ g, Cov = vcv(tree))

> 

> pw.ols <- pairwise(fit.ols, groups = g)

> pw.gls <- pairwise(fit.gls, groups = g)

> 

> anova(fit.ols)

Analysis of Variance, using Residual Randomization

Permutation procedure: Randomization of null model residuals 

Number of permutations: 1000 

Estimation method: Ordinary Least Squares 

Sums of Squares and Cross-products: Type I 

Effect sizes (Z) based on F distributions

          Df      SS     MS     Rsq      F      Z Pr(>F)  

g          1  9.3241 9.3241 0.31269 8.1891 2.0588  0.016 *

Residuals 18 20.4948 1.1386 0.68731                       

Total     19 29.8190                                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Call: lm.rrpp(f1 = y ~ g)

> summary(pw.ols)

Pairwise comparisons

Groups: 1 2 

RRPP: 1000 permutations

LS means:

Vectors hidden (use show.vectors = TRUE to view)

Pairwise distances between means, plus statistics

           d UCL (95%)        Z Pr > d

1:2 1.365587  1.148485 1.957779  0.016

The P-values are exactly the same when performing pairwise comparisons between the only two groups as with performing ANOVA.  But if we do the same thing with the PGLS model fit:

 

> anova(fit.gls)

Analysis of Variance, using Residual Randomization

Permutation procedure: Randomization of null model residuals 

Number of permutations: 1000 

Estimation method: Generalized Least-Squares (via OLS projection) 

Sums of Squares and Cross-products: Type I 

Effect sizes (Z) based on F distributions

          Df     SS     MS  Rsq      F      Z Pr(>F)  

g          1  5.104 5.1040 0.19 4.2223 1.6193  0.045 *

Residuals 18 21.759 1.2088 0.81                       

Total     19 26.863                                   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Call: lm.rrpp(f1 = y ~ g, Cov = vcv(tree))

> summary(pw.gls)

Pairwise comparisons

Groups: 1 2 

RRPP: 1000 permutations

LS means:

Vectors hidden (use show.vectors = TRUE to view)

Pairwise distances between means, plus statistics

          d UCL (95%)        Z Pr > d

1:2 1.57239  1.964684 1.202466  0.117

Why did this happen?  There is no apparent way to transform the distances by the phylogenetic covariance matrix, C, for the two estimates made.  The distance between PGLS estimates is not significant even though the between-group variance is.  In ANOVA, variances are appropriate because the residual covariance matrix is the GLS covariance matrix.  In pairwise comparisons, not as much.

I believe the solution is probably to use instead of d, the statistic d/se, where se is the estimated standard error in every random permutation.  This would accomplish the same thing, in principle, as a parametric t-test.  Chances are the standard error would tend to be higher in random permutations in the above example, so the observed d would be relatively larger, if divided by se.

This is the research that still has to be done: find a procedure that produces consistent distributions, despite the inability to properly transform statistics.  If the example above were real data and I found conflicting results, I would interpret the results as the effect of group was slightly significant but evidence was lacking to suggest their means are distinct.  In other words, i would not have sufficient statistical power to detect a meaningful difference between means.

Hope that helps!

Mike

--
You received this message because you are subscribed to the Google Groups "geomorph R package" group.
To unsubscribe from this group and stop receiving emails from it, send an email to geomorph-r-pack...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/b8eb2a85-0a92-4bfc-a4e2-0f047465e236n%40googlegroups.com.

[Mike Collyer]

Dear Colleagues,

This email is long.  The beginning is important and rest is intended for people who might find interest in the details.

I wish to alert you about an important update to the RRPP::pairwise function, plus make you aware of some other updates.  The important update is because I discovered a coding error, which results in ignoring input of a null model in the pairwise function, if the full model fit was made with verbose = TRUE.  Although pairwise would calculate new coefficients across random permutations, it was failing to replace them in the model if coefficients were already there, which would be the case with verbose results.  

To update RRPP, please use the Github version or RRPP,

devtools::install_github("mlcollyer/RRPP", build_vignettes = TRUE)

I found this error and fixed it while I was updating the pairwise function, related to the thread (below) started by Daniel last week.  I speculated that a standardized distance, d / pooled se, might be better suited for PGLS models.  (This is because the total sums of squares [TSS] varies across RRPP permutations for a single-factor model, unlike OLS models.)  This was easy to program and I found, using the example I used before in the thread, that P-values were the same between ANOVA and the standardized pairwise distance, for a two-factor model, unlike the demonstration with Euclidean distances.

I also realized that this was one way to generalize a standardized distance.  Another way is to use the Mahalanobis distance.  Interestingly, Mahalanobis distance tracks with MANOVA results for two sample tests, whether OLS or PGLS, the way Euclidean distances track with ANOVA for OLS models, and standardized distances track with ANOVA results for PGLS models.

Therefore, the pairwise function now has two additional test types, like 'dist': ‘stdist' and ‘mdist’.  In order to have a test for Mahalanobis distances, it was essential to have a function to calculate pairwise Mahalanobis distances.  Therefore, I made this function available, as well: RRPP::mahal_dist, which functions the same as the stats::dist function, but allows a covariance matrix to be input (see the help file for an example).  This function and the pairwise test require more observations than variables, a limitation that cannot be avoided.

A few notes I believe are important with these updates:

• The distance between LS means, d, has morphological meaning.  For Procrustes coordinates, d is the Procrustes distance; for Procrustes coordinates projected into a space tangent to shape space (typical), d is Euclidean distance that approximates Procrustes distance.  
• Either standardized or Mahalanobis distances are statistics, and should not be interpreted as the amount of morphological change.  Even if using them for tests, Euclidean distances should be reported as descriptive measures of morphological change.  (Just like one would never call a t-statistic the mean difference.)
• For OLS models, d and standardized d provide the same results for two-sample models, as the latter is basically a linear transformation of the former.  Only with PGLS models would the standardized d make more sense for tests.
• With 3+ groups and OLS models, more pairwise comparisons with standardized d might be deemed “significant” because the pooled SE tends to be larger in random permutations, not because distances are necessarily smaller.  This statistical artifact is probably not good justification for biological interpretations.
• P-values (as well as effect sizes) tended be highly correlated among distance types, so using Euclidean distance is still a good way to go.  Daniel noted that a significant ANOVA for a PGLS model might not result in significant pairwise differences.  This could be because of using Euclidean distances rather than standardized distances.  However, if a distance switches from non-significant to significant by updating from Euclidean to standardized methods, the shape difference is probably not a large effect.  Several borderline cases can make an ANOVA effect seem significant.  It is always smart to focus on effect sizes.

Sorry for the long email.  There are more updates coming down the pike.  There are always more improvements to be made.

Best,

Mike

[Juan Esteban Vrdoljak]

Hi Mike, thanks for your very detailed explanations. I am left with a doubt, in which cases would you prefer to use mahalanobis distances instead of the other options?
Regards
J

[Daniel Fernandes da Silva]

Dear Mike,

Once again thank you very much for your prompt answer and patience to explain and demonstrate in detail the impossibility to transform the distances by the phylogenetic covariance matrix. Now I think I understand why it is possible to have a significant effect considering a PGLS model fit and yet, when performing pairwise comparisons, obtain no significant results. In fact, considering my data, I believe your example is exactly what happened: the effect of group has no sufficient statistical power to detect a meaningful difference between means, so I found conflicting results and must interpret that there is no significant group effect based on the results of PGLS pairwise comparisons.  

Thank you so very much.

Daniel

You received this message because you are subscribed to a topic in the Google Groups "geomorph R package" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/geomorph-r-package/rXtKYXe5AIc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to geomorph-r-pack...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/7536C5B3-A672-49B3-AD5B-DD57D90C0053%40gmail.com.

--

Daniel S. Fernandes
Universidade Federal do Rio de Janeiro
Instituto de Biologia - Dep. de Zoologia

Laboratório TaxoN (Laboratório de Répteis)

Av. Carlos Chagas Filho, 791

Antigo Polo Biorio - Prédio Compartilhado

Ilha do Fundão - Rio de Janeiro - RJ - BRASIL

CEP: 21941-599

e-mail: danfe...@gmail.com
Tel: (21)3938-6366

Currículo Lattes: http://lattes.cnpq.br/1940653720704787 

Orcid: https://orcid.org/0000-0001-6155-6246

Site: http://www.herpetologiamuseunacional.com.br/equipe_daniel_fernandes.html

[Mike Collyer]

Hi Juan,

If you are using MANOVA instead of ANOVA, then it’s better to use Mahalanobis distance. The reason these are consistent is both MANOVA and Mahalanobis distance use the inverse of the residual sums of squares and cross-products matrix in their calculations. The former for standardizing an effect sums of squares and cross-products matrix and the latter for the distance between vectors (also scaled by degrees of freedom). 

ANOVA divides the trace of the an effect covariance matrix by the trace of the residual covariance matrix. The square root of the same matrix trace is used to calculate the pooled standard error used for standardized distance. 

To summarize, MANOVA and Mahalanobis distance involve multiplication of inverse residual covariance matrices or their transformations, and ANOVA and standardized distance involve dividing by traces of residual covariance matrices or their transformations. Thus, there is an inherent matching of types. 

Hope that helps!

Mike

Sent from my iPhone

On Jul 2, 2024, at 9:28 AM, Juan Esteban Vrdoljak <juan.v...@gmail.com> wrote:

Hi Mike, thanks for your very detailed explanations. I am left with a doubt, in which cases would you prefer to use mahalanobis distances instead of the other options?

To view this discussion on the web, visit https://groups.google.com/d/msgid/geomorph-r-package/06d31149-0118-44c4-85de-45a47b59483fn%40googlegroups.com.