pmlMix {phangorn} | R Documentation |
Phylogenetic mixture model.
pmlMix(formula, fit, m = 2, omega = rep(1/m, m), ...)
formula |
a formula object (see details). |
fit |
an object of class pml . |
m |
number of mixtures. |
omega |
mixing weights. |
... |
Further arguments passed to or from other methods. |
The formula
object allows to specify which parameter get
optimized. The formula is generally of the form edge + bf + Q
~ rate + shape + ...
, on the left side are the parameters which
get optimized over all mixtures, on the right the parameter which
are optimized specific to each mixture. The parameters available
are "nni", "bf", "Q", "inv", "shape", "edge", "rate"
.
Each parameters can be used only once in the formula.
"rate"
and "nni"
are only available for the right side
of the formula. On the other hand parameters for invariable sites
are only allowed on the left-hand side.
The convergence of the algorithm is very slow and is likely that the
algorithm can get stuck in local optima.
pmlMix
returns a list with elements
logLik |
log-likelihood of the fit |
omega |
mixing weights. |
fits |
fits for the final mixtures. |
Klaus Schliep klaus.schliep@gmail.com
X <- allSitePattern(5) tree <- read.tree(text = "((t1:0.3,t2:0.3):0.1,(t3:0.3,t4:0.3):0.1,t5:0.5);") fit <- pml(tree,X, k=4) weights <- 1000*exp(fit$site) attr(X, "weight") <- weights fit1 <- update(fit, data=X, k=1) fit2 <- update(fit, data=X) (fitMixture <- pmlMix(edge~rate, fit1 , m=4)) (fit2 <- optim.pml(fit2, optGamma=TRUE)) ## Not run: data(Laurasiatherian) dm <- dist.logDet(Laurasiatherian) tree <- NJ(dm) fit=pml(tree, Laurasiatherian) fit = optim.pml(fit) fit2 <- update(fit, k=4) fit2 <- optim.pml(fit2, optGamma=TRUE) fitMix = pmlMix(edge ~ rate, fit, m=4) fitMix # # simulation of mixture models # X <- allSitePattern(5) tree1 <- read.tree(text = "((t1:0.1,t2:0.5):0.1,(t3:0.1,t4:0.5):0.1,t5:0.5);") tree2 <- read.tree(text = "((t1:0.5,t2:0.1):0.1,(t3:0.5,t4:0.1):0.1,t5:0.5);") tree1 <- unroot(tree1) tree2 <- unroot(tree2) fit1 <- pml(tree1,X) fit2 <- pml(tree2,X) weights <- 2000*exp(fit1$site) + 1000*exp(fit2$site) attr(X, "weight") <- weights fit1 <- pml(tree1, X) fit2 <- optim.pml(fit1) logLik(fit2) AIC(fit2, k=log(3000)) fitMixEdge = pmlMix( ~ edge, fit1, m=2) logLik(fitMixEdge) AIC(fitMixEdge, k=log(3000)) fit.p <- pmlPen(fitMixEdge, .25) logLik(fit.p) AIC(fit.p, k=log(3000)) ## End(Not run)