1

From what I've read, the Monte Carlo correction is often used in pairwise PERMANOVA comparisons to improve the accuracy of the p-value when given a low number of sample replicates and limited number of possible permutations. This is easily implemented in the PRIMER software, but I can't find an existing function to do this in R.

In the basic p.adjust() function in R, the possible methods are c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

So I'm wondering if, firstly, anyone knows of an existing package or function that allows for monte carlo p adjustments, or, failing that, if someone could please describe the method to me in fairly simple terms so that I may cobble together a function myself?

I'll include the pairwise permanova function I am working with (this function was obtained from https://www.researchgate.net/post/How_can_I_do_PerMANOVA_pairwise_contrasts_in_R)

data(iris)

library(vegan)

pairwise.adonis <- function(x,factors, sim.function = 'vegdist', sim.method = 'bray', p.adjust.m ='bonferroni')
{
  library(vegan)

  co = combn(unique(as.character(factors)),2)
  pairs = c()
  F.Model =c()
  R2 = c()
  p.value = c()


  for(elem in 1:ncol(co)){
    if(sim.function == 'daisy'){
      library(cluster); x1 = daisy(x[factors %in% c(co[1,elem],co[2,elem]),],metric=sim.method)
    } else{x1 = vegdist(x[factors %in% c(co[1,elem],co[2,elem]),],method=sim.method)}

    ad = adonis(x1 ~ factors[factors %in% c(co[1,elem],co[2,elem])] );
    pairs = c(pairs,paste(co[1,elem],'vs',co[2,elem]));
    F.Model =c(F.Model,ad$aov.tab[1,4]);
    R2 = c(R2,ad$aov.tab[1,5]);
    p.value = c(p.value,ad$aov.tab[1,6])
  }
  p.adjusted = p.adjust(p.value,method=p.adjust.m)
  sig = c(rep('',length(p.adjusted)))
  sig[p.adjusted <= 0.05] <-'.'
  sig[p.adjusted <= 0.01] <-'*'
  sig[p.adjusted <= 0.001] <-'**'
  sig[p.adjusted <= 0.0001] <-'***'

  pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted,sig)
  print("Signif. codes:  0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1")
  return(pairw.res)

} 


pairwise.adonis(iris[,1:4],iris[,5])
J.Con
  • 4,101
  • 4
  • 36
  • 64
  • This is not really a SO question so much as one for [CrossValidated](https://stats.stackexchange.com/), I suggest you ask there. – r2evans Sep 29 '17 at 22:31
  • I did and it got me the Tumbleweed badge! – J.Con Sep 29 '17 at 23:50
  • Oh, bummer ... I got that badge [some time ago](https://stackoverflow.com/questions/22994168/taskhandler-that-executes-before-the-expr). I still think it's more appropriate there (since it's a bit more the theory than about programming itself), but you apparently did not find any takers. Not much to offer at the moment unless you're willing to contract somebody code it for you :-( – r2evans Sep 29 '17 at 23:57
  • This can be regarded as a technical question of implementation that suits SO. However, the code would really need some exercise in **R** style. – Jari Oksanen Oct 02 '17 at 08:05
  • Hi @JariOksanen, thank you for your comment. I see you faithfully comment on most questions tagged with vegan. Can you please elaborate on what you mean by 'the code would really need some exercise in R style'? Thank you – J.Con Oct 03 '17 at 10:13
  • Most obvious thing is that in **R** calculation should be separated from presentation (printing): instead of hard-coding formats among calculations, you should use `printCoefmat` to get the stars etc. Then `library()` should not be within `for` loop. Also, `subset` command for `data` could be handy in delimiting analysis to a certain pair of variables (`adonis2` has no `subset=` argument -- perhaps it should -- but you can use `subset` within `data=`). – Jari Oksanen Oct 03 '17 at 13:26
  • That is good advice, thank you very much. Do you have any comments on monte carlo corrections, per chance? Thank you – J.Con Oct 04 '17 at 23:06
  • My only comment is a non-technical one: I don't use these corrections. There is a lot of literature that warns against their use (providing sound procedures are used in reporting the results). – Jari Oksanen Oct 10 '17 at 11:51
  • Fair enough! I appreciate your advice. – J.Con Oct 11 '17 at 21:18
  • Hi, could you provide references agaisnt the use of Montecarlo in Permanova designs?. Why the permutation test in PERMANOVA could yield significant different results when unique permutations are low (e.g. < 10) compared to Monte Carlo?. – AnastD Oct 03 '19 at 14:17

0 Answers0