null p-values for logistic phyloglm seem biased?

[pbradz]

Hi folks,

Thanks so much for writing and maintaining such a useful package!

I noticed something odd when using phyloglm to associate two binary traits. Under the null, when each trait is being generated independently along the same tree, I would expect to see a more-or-less "flat" p-value distribution. Instead, the distributions I am seeing seem pretty sharply tilted towards the right (i.e., conservative). How much they are skewed also seems to depend on the tree: random coalescent trees seem to consistently give more conservative results than random bifurcating trees made ultrametric with chronos(). Turning the bootstrap on doesn't seem to make an appreciable difference.

This is the code I'm using:

library(ape)
library(phylolm)
library(tidyverse)

get_null_pvals <- function(test_tree, N=200, bstrap=0, alpha=2, beta=1) {
  force(beta)
  force(alpha)
  map_dbl(1:N, function(.) {
    traits <- rbinTrait(n=2, test_tree, beta=beta, alpha=alpha)
    result <- tryCatch(
      phyloglm(traits[, 2] ~ traits[, 1], phy=test_tree, boot=bstrap),
      error = function(e) {
        return(NA)
      })
    if (class(result)=="phyloglm") {
      pv <- coefficients(summary(result))[2, "p.value"]
      return(pv)
    } else if (class(result) == "logical") {
      return(NA)
    }
  }, .progress=TRUE)
}
rc_pv <- get_null_pvals(rcoal(n=100), alpha=2, beta=1)
rt_pv <- get_null_pvals(chronos(rtree(n=100)), alpha=2, beta=1)
hist(rc_pv)
hist(rt_pv)

Curious if anyone has thoughts as to why this would be happening? Am I missing something obvious?

Patrick

P.S. I have not seen this behavior with phylolm and rTrait, regardless of the tree I use, e.g.:

[lamho86]

Hi Patrick,

This is really odd. I am not sure what is happening. Thank you for pointing this out. I will need to think more about it.
Lam