Add support for phylolm to DHARMa

[florianhartig]

A user requests support for https://github.com/lamho86/phylolm

Perspective for this request: Currently, the package does not allow to simulate from the fitted model, which prevents an easy integration into DHARMa. Support could be added if the package developers implement a simulate function, see here

Interim solution for users that, for some reason, have to urgently use DHARMa with this package: take your fitted model, create a simulate function for this model structure yourself, and then use createDHARMa (see help), this will allow most options of the package to be run. See further comments on support of new packages in here. Note that the vignette has some further comments / examples on creating custom simulation functions and reading them into DHARMa

See also comments here: lamho86/phylolm#27

[florianhartig]

Started to address this in https://github.com/florianhartig/DHARMa/tree/129-phylolm

[albanmathieu-pro]

Hi @florianhartig ,
I want to do exactly what you are working on this branch,
I've seen your discussion here, lamho86/phylolm#27

I tried to install your dev branch 129-phylolm

devtools::install_github(repo = "florianhartig/DHARMa", subdir = "DHARMa", ref = "129-phylolm",
 dependencies = T, build_vignettes = T, force = TRUE)

but i can not run your test

set.seed(123456)
tre = rtree(50)
x = rTrait(n=1,phy=tre)
X = cbind(rep(1,50),x)
y = rbinTrait(n=1,phy=tre, beta=c(-1,0.5), alpha=1 ,X=X)
dat = data.frame(trait01 = y, predictor = x)
fit = phyloglm(trait01~predictor,phy=tre,data=dat,boot=100)

summary(fit)
coef(fit)
vcov(fit)

DHARMa::simulateResiduals(fit, plot = T)

here is the error:

Error in nobs.default(fittedModel) : no 'nobs' method is available
In addition: Warning message:
In checkModel(fittedModel) :
  DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work!

Can you help me figure out please ?

[MorceletL]

Hi,
I'm currently running phylolm models and I would like to first check if the model actually managed to correct for the phylogenetic signal, and then do inferences.
I've been reading the different comments you wrote here : https://github.com/lamho86/phylolm/issues/27, and tried to generate residuals from a model runned with bootstrap:
model <- phylolm(log(Time_RRH) ~ PC2, tree, data = breeding_RRH, REML =F, boot = 100, save = T)
simulateResiduals(model, rotation = T)
Error in Method("family") : no method for 'family' applicable for a class "phylolm" object
Warning message: In checkModel(fittedModel) : DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work!

I wanted to know if you managed to build a solution to interpret and use phylolm model.
Thanks in advance
Lucile

[florianhartig]

Hi all,

first of all, to run this, you need to also install an updated version of the phylolm package, this is not yet on CRAN.

devtools::install_github("lamho86/phylolm")
devtools::install_github(repo = "florianhartig/DHARMa", subdir = "DHARMa", ref = "129-phylolm",
 dependencies = T, build_vignettes = T, force = TRUE)

I have requested that the phylolm developers push an update to CRAN so that this functionality becomes available directly, see here lamho86/phylolm#27

Second, so far I had implemented phyloglm, but not phylolm, but I have added this now, so now it should work for both models.

If everything is installed, it should work now

set.seed(123456)
tre = rcoal(60)
taxa = sort(tre$tip.label)
b0=0; b1=1;
x <- rTrait(n=1, phy=tre,model="BM",
            parameters=list(ancestral.state=0,sigma2=10))
y <- b0 + b1*x + 
  rTrait(n=1,phy=tre,model="lambda",parameters=list(
    ancestral.state=0,sigma2=1,lambda=0.5))
dat = data.frame(trait=y[taxa],pred=x[taxa])
fit = phylolm(trait~pred,data=dat,phy=tre,model="lambda")
summary(fit)

res = DHARMa::simulateResiduals(fit, plot = T)

res = DHARMa::simulateResiduals(fit, plot = T, rotation = "estimated")

@MorceletL - the rotation argument should in principle allow you to test if the phylogenetic signal is taken out but also this is highly experimental and I also haven't checked in detail that simulations in phylolm are done properly! Proceed with care!

I want to stress that this is all currently untested, so this comes with zero guarantees - if you want to use this, I would make sure with simulated data that the model behaves as expected!

[florianhartig]

I have added a new phylogenetic autocorrelation test, you can try

# https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00044.x

library(DHARMa)
library(phylolm)

set.seed(123)
tre = rcoal(60)
b0=0; b1=1;
x <- runif(length(tre$tip.label), 0,1)
y <- b0 + b1*x + 
  rTrait(n=1, phy=tre,model="BM",
         parameters=list(ancestral.state=0,sigma2=10))
dat = data.frame(trait=y,pred=x)

fit = lm(trait~pred,data=dat)
res = simulateResiduals(fit, plot = T)

testPhylogeneticAutocorrelation(res, tree = tre)


fit = phylolm(trait~pred,data=dat,phy=tre,model="BM")
summary(fit)

res = DHARMa::simulateResiduals(fit, plot = T)
res = DHARMa::simulateResiduals(fit, plot = T, rotation = "estimated")

[MorceletL]

Dear Florian,

Thank you very much for these update.

I've been reading the article you mentionned here, and applied it earlier to my model using the package phylosignal, currently using the same data doesn't give the same results, which is a bit strange (it might be due to the scaling of the residuals implemented in DHARMa).

lm2p<-lm(log(Time_RRH) ~ PC2,  data = breeding_RRH)
phylo_res <- phylo4d(tree, tip.data = as.data.frame(residuals(lm2p)))
phyloSignal(phylo_res)

$stat
                    Cmean        I       K    K.star    Lambda
residuals.lm2p. 0.7050401 0.265141 0.41895 0.3425516 0.8570531

$pvalue
                Cmean     I     K K.star Lambda
residuals.lm2p. 0.001 0.001 0.001  0.001  0.001

test<-lm(log(Time_RRH) ~ PC2,  data = breeding_RRH)
test_res <- simulateResiduals(test, plot = T)
testPhylogeneticAutocorrelation(test_res,tree = tree)

	DHARMa Moran's I test for phylogenetic autocorrelation

data:  test_res
observed = -0.110066, expected = -0.022222, sd = 0.050719, p-value = 0.08328
alternative hypothesis: Phylogenetic autocorrelation

Secondly, I would like to have your opinion on the residual part of the model. In the case oh phylolm residuals are given raw (i.e. without considering the phylogenetic correlation as explained by @celineane here https://github.com/lamho86/phylolm/issues/3. The best way to have the phylogenetic residuals is to multiply them by the inverse square root of the covariance matrix (function sqrt_OU_covariance().)

Thus in the rotation argument of DHARMa::simulateResiduals do you agree that this matrix should be added to verify residuals normality and homoscedasticity of the model ? , or shoudl it be the raw residuals ? It is unclear to know which residuals should be analysed to verify the model, since by definition phylolm should correct residuals for phylogenetic correlation.

[florianhartig]

Hi @MorceletL,

about the Revell 2010 paper: sorry about this, I just pasted it there as a reminder to look at the methods in this paper, but it has nothing to do with the later code example.

Your second question is probably the relevant one: in DHARMa, there are two options to rotate residuals, which are described in https://rdrr.io/cran/DHARMa/man/getQuantile.html:

1. you provide the covariance matrix, in which case you set rotation = COV, but then you need to get the COV out of the model in the correct way
2. DHARMa can estimate the covariance from the simulations, which is the option correlation = "estimated", which I use in the example above. I have found that 2 often works surprisingly well and this is of course very convenient because you can't do anything wrong, but of course 1 is more exact.

So, in the example above, I contrast raw quantile with rotated quantile residuals, and you see that the rotated ones look much better.

To verify the model, you should definitely use the rotated residuals! phylolm is correcting for phylogenetic correlation when estimating the parameters, but as for all gls, the correlation remains in the residuals unless you explicitly account for it.

[florianhartig]

p.s.: unless you use software that allows to condition on the covariance structure, which is not possible for gls, but for some Bayesian implementations or also glmmTMB. See discussion in https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13971