In emvolz-phylodynamics/mlesky: Maximum likelihood estimation of effective population size through time, growth rates, and phylodynamic regressions

Here we reanalyze the MRSA phylogeny from Volz & Didelot, Systematic Biology 2018.

require(mlesky)
treeMRSA <- ape::read.tree(system.file('mrsa.nwk', package = 'mlesky'))
treeMRSA$root.time=1967.306
fit <- mlskygrid(treeMRSA)
plot (fit, logy = FALSE)

The analysis above assumed the default smoothing parameter tau=1. We can use cross-validation to optimize this smoothing parameter:

fit <- mlskygrid(treeMRSA, tau = NULL, tau_lower = 1, tau_upper = 20, ncpu = 6)
plot(fit,logy=FALSE)

We can use AIC to optimise res and then use the cross-validation to optimise tau:

res=optim_res_aic(treeMRSA,ncpu=6)
print(res)
fit <- mlskygrid(treeMRSA, tau = NULL, tau_lower = 1, tau_upper = 20, ncpu = 6,res=res)
plot(fit,logy=FALSE)

Let's compare the three models:

fit <- mlskygrid(treeMRSA, tau = NULL, tau_lower = 1, tau_upper = 20, ncpu = 6,res=20,model=1)
plot(fit,logy=FALSE)
fit <- mlskygrid(treeMRSA, tau = NULL, tau_lower = 1, tau_upper = 20, ncpu = 6,res=20,model=2)
plot(fit,logy=FALSE)
fit <- mlskygrid(treeMRSA, tau = NULL, tau_lower = 1, tau_upper = 20, ncpu = 6,res=20,model=3)
plot(fit,logy=FALSE)

emvolz-phylodynamics/mlesky documentation built on Feb. 13, 2025, 1:47 p.m.