Changing alternative hypotheses for AMOVA randtest

[Sally Chang]

Hi there,

I was wondering if anyone had experience running a permutation test on an AMOVA done in Poppr and had success in changing the alternative hypotheses, or could give me some insight into how randtest() chooses what hypothesis to use for each test. 

I call randtest() as such: 

> One_amova_randtest <-randtest(OneIndiv_amova,nrepet=999,alter=c("two-tailed","two-tailed","two-tailed"))

And get the following results: 

> One_amova_randtest

class: krandtest 

Monte-Carlo tests

Call: randtest.amova(xtest = OneIndiv_amova, nrepet = 999, alter = c("two-tailed", 

    "two-tailed", "two-tailed"))

Number of tests:   4 

Adjustment method for multiple comparisons:   none 

Permutation number:   999 

                         Test       Obs   Std.Obs   Alter Pvalue

1   Variations within samples 16.666663 -7.864566    less  0.001

2  Variations between samples  2.373188  3.340489 greater  0.001

3 Variations between Location  2.959384 13.671117 greater  0.001

4     Variations between Zone 69.896480  3.953857    less  1.000

other elements: adj.method call 

It chooses "less, greater, great, less" for the alternative hypotheses no matter if I try to use the "alter" option in the function call or not. 

For reference, I am using the following versions of things: 

R version 3.3.2 (2016-10-31)

Platform: x86_64-apple-darwin13.4.0 (64-bit)

Running under: OS X El Capitan 10.11.2

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] poppr_2.2.1    adegenet_2.0.1 ade4_1.7-4    

loaded via a namespace (and not attached):

 [1] Rcpp_0.12.7      spdep_0.6-8      plyr_1.8.4       pegas_0.9       

 [5] LearnBayes_2.15  tools_3.3.2      boot_1.3-18      digest_0.6.10   

 [9] tibble_1.2       nlme_3.1-128     gtable_0.2.0     lattice_0.20-34 

[13] mgcv_1.8-15      fastmatch_1.0-4  Matrix_1.2-7.1   igraph_1.0.1    

[17] shiny_0.14.2     DBI_0.5-1        parallel_3.3.2   coda_0.18-1     

[21] cluster_2.0.5    dplyr_0.5.0      stringr_1.1.0    gtools_3.5.0    

[25] grid_3.3.2       R6_2.2.0         phangorn_2.0.4   sp_1.2-3        

[29] gdata_2.17.0     ggplot2_2.1.0    reshape2_1.4.2   seqinr_3.3-3    

[33] deldir_0.1-12    magrittr_1.5     nnls_1.4         scales_0.4.0    

[37] htmltools_0.3.5  MASS_7.3-45      splines_3.3.2    gmodels_2.16.2  

[41] assertthat_0.1   permute_0.9-4    mime_0.5         ape_3.5         

[45] colorspace_1.2-7 xtable_1.8-2     httpuv_1.3.3     quadprog_1.5-5  

[49] stringi_1.1.2    munsell_0.4.3    vegan_2.4-1     

Any insight would be much appreciated into how choice of/ability to change hypotheses works in the context of randtest() for a poppr.amova() result.

Thanks!

-Sally

[Zhian Kamvar]

Hi Sally,

I think you may have found a bug in ade4's implementation of randtest.amova(). You may consider using the implementations found in pegas::amova() or vegan::adonis(). All you need are your strata and a distance matrix.  

Here's an example of how to do that: http://popgen.nescent.org/2015-12-15-microsatellite-differentiation.html#analysis-of-molecular-variance

If you want to stick with the ade4 implementation of AMOVA (which gives you within-individual variance by default), you can use the function as.krandtest() to fix the alternate hypotheses:

with(One_amova_randtest, as.krandtest(sim, obs, alter =c("two-sided", "two-sided", "two-sided", "two-sided"), call = call, names = names)

For context: poppr.amova provides an interface to AMOVA in the packages ade4 (default) and pegas. The function randtest is found in ade4, and it appears that the alternates are hardcoded in the function (here's the code that's currenlty on CRAN: https://github.com/cran/ade4/blob/c28ea90510fdfe980c776c286d3a47256fdc4b77/R/randtest.amova.R).

The reason why the alternate hypotheses are "less", "greater", "greater", "less" is because they are hardcoded in the function as alter=c("less","greater","greater"). This 3-element vector gets incorrectly recycled to the number of hierarchical levels (note vector recycling is default in this language: http://eriqande.github.io/rep-res-web/lectures/vectorization_recycling_and_indexing.html#recycling).

I hope that helps. I'll include a note about this in a future poppr release and contact the maintainer of ade4 about this bug.

Best,

Zhian

P.S. Thank you for posting the results of sessionInfo()!