Permalink
Browse files

nsim>1 for sim.Mk & sim.multiMk

  • Loading branch information...
liamrevell committed Feb 7, 2018
1 parent e8985ae commit 2b49dea93fac145fcf90145b6eb031e6553b6105
Showing with 48 additions and 34 deletions.
  1. +4 −4 DESCRIPTION
  2. +42 −28 R/fitMk.R
  3. +2 −2 man/sim.history.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
Version: 0.6-54
Date: 2018-2-3
Version: 0.6-55
Date: 2018-2-7
Title: Phylogenetic Tools for Comparative Biology (and Other Things)
Author: Liam J. Revell
Maintainer: Liam J. Revell <liam.revell@umb.edu>
@@ -13,6 +13,6 @@ ZipData: no
Description: A wide range of functions for phylogenetic analysis. Functionality is concentrated in phylogenetic comparative biology, but also includes a diverse array of methods for visualizing, manipulating, reading or writing, and even inferring phylogenetic trees and data. Included among the functions in phylogenetic comparative biology are various for ancestral state reconstruction, model-fitting, simulation of phylogenies and data, and multivariate analysis. There are a broad range of plotting methods for phylogenies and comparative data which include, but are not restricted to, methods for mapping trait evolution on trees, for projecting trees into phenotypic space or a geographic map, and for visualizing correlated speciation between trees. Finally, there are a number of functions for reading, writing, analyzing, inferring, simulating, and manipulating phylogenetic trees and comparative data not covered by other packages. For instance, there are functions for randomly or non-randomly attaching species or clades to a phylogeny, for estimating supertrees or consensus phylogenies from a set, for simulating trees and phylogenetic data under a range of models, and for a wide variety of other manipulations and analyses that phylogenetic biologists might find useful in their research.
License: GPL (>= 2)
URL: http://github.com/liamrevell/phytools
Packaged: 2018-2-3 12:00:00 EST
Packaged: 2018-2-7 12:00:00 EST
Repository:
Date/Publication: 2018-2-3 12:00:00 EST
Date/Publication: 2018-2-7 12:00:00 EST
View
@@ -265,48 +265,62 @@ EXPM<-function(x,...){
## function to simulate multiple-rate Mk multiMk
## written by Liam J. Revell 2018
sim.multiMk<-function(tree,Q,anc=NULL,...){
sim.multiMk<-function(tree,Q,anc=NULL,nsim=1,...){
if(hasArg(as.list)) as.list<-list(...)$as.list
else as.list<-FALSE
ss<-rownames(Q[[1]])
tt<-map.to.singleton(reorder(tree))
P<-vector(mode="list",length=nrow(tt$edge))
for(i in 1:nrow(tt$edge))
P[[i]]<-expm(Q[[names(tt$edge.length)[i]]]*tt$edge.length[i])
if(is.null(anc)) anc<-sample(ss,1)
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-anc
for(i in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[i]][STATES[i,1],])[,1]==1)]
STATES[i,2]<-new
ii<-which(tt$edge[,1]==tt$edge[i,2])
if(length(ii)>0) STATES[ii,1]<-new
if(nsim>1) X<- if(as.list) vector(mode="list",length=nsim) else
data.frame(row.names=tt$tip.label)
for(i in 1:nsim){
a<-if(is.null(anc)) sample(ss,1) else anc
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-a
for(j in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[j]][STATES[j,1],])[,1]==1)]
STATES[j,2]<-new
ii<-which(tt$edge[,1]==tt$edge[j,2])
if(length(ii)>0) STATES[ii,1]<-new
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
if(nsim>1) X[,i]<-x else X<-x
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
x
X
}
## constant-rate Mk model simulator
## written by Liam J. Revell 2018
sim.Mk<-function(tree,Q,anc=NULL,...){
sim.Mk<-function(tree,Q,anc=NULL,nsim=1,...){
if(hasArg(as.list)) as.list<-list(...)$as.list
else as.list<-FALSE
ss<-rownames(Q)
tt<-reorder(tree)
P<-vector(mode="list",length=nrow(tt$edge))
for(i in 1:nrow(tt$edge))
P[[i]]<-expm(Q*tt$edge.length[i])
if(is.null(anc)) anc<-sample(ss,1)
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-anc
for(i in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[i]][STATES[i,1],])[,1]==1)]
STATES[i,2]<-new
ii<-which(tt$edge[,1]==tt$edge[i,2])
if(length(ii)>0) STATES[ii,1]<-new
if(nsim>1) X<- if(as.list) vector(mode="list",length=nsim) else
data.frame(row.names=tt$tip.label)
for(i in 1:nsim){
a<-if(is.null(anc)) sample(ss,1) else anc
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-a
for(j in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[j]][STATES[j,1],])[,1]==1)]
STATES[j,2]<-new
ii<-which(tt$edge[,1]==tt$edge[j,2])
if(length(ii)>0) STATES[ii,1]<-new
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
if(nsim>1) X[[i]]<-x else X<-x
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
x
X
}
View
@@ -5,8 +5,8 @@
\title{Simulate character history or a discrete character at the tips of the tree under some model}
\usage{
sim.history(tree, Q, anc=NULL, nsim=1, ...)
sim.Mk(tree, Q, anc=NULL, ...)
sim.multiMk(tree, Q, anc=NULL, ...)
sim.Mk(tree, Q, anc=NULL, nsim=1, ...)
sim.multiMk(tree, Q, anc=NULL, nsim=1, ...)
}
\arguments{
\item{tree}{a phylogenetic tree as an object of class \code{"phylo"}. For the case of \code{sim.multiMk} \code{tree} should be an object of class \code{"simmap"} in which the regimes for simulation have been mapped onto the tree.}

0 comments on commit 2b49dea

Please sign in to comment.