Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache=FALSE)
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=75),tidy=TRUE)
## ---- echo=FALSE, message=FALSE, warning=FALSE--------------------------------
require(OUwie)
## ---- eval=TRUE---------------------------------------------------------------
data(tworegime)
pp <- OUwie(tree, trait, model="OUM", root.station=TRUE, scaleHeight=TRUE, shift.point=0, algorithm="invert", quiet=TRUE)
pp
## ---- eval=TRUE---------------------------------------------------------------
data(tworegime)
pp <- OUwie(tree, trait, model="OUM", root.station=FALSE, scaleHeight=TRUE, shift.point=0, algorithm="invert", quiet=TRUE)
pp
## ---- eval=TRUE, echo=FALSE, fig.height=10, fig.width=6.5, fig.cap = "The edges are painted by regime, assuming an optimum $\\theta_i$ for each color. As shown in Ho and Ane (2014a) the left shows an case of unidentifiability case because every regime forms a connected component. The tree on the right shows a case of identifiability because the black regime covers two completely disconnected parts in the tree."----
par(mfcol=c(1,2), mar=c(4,4,0.5,0.5), oma=c(1,2,1,1))
load("simUnidentifiable.Rsave")
regimes <- OUwie:::GetParameterPainting(phy=phy, data=data, rates=c(3,1,2), k.pars=3)
nb.tip <- Ntip(phy)
nb.node <- Nnode(phy)
comp <- numeric(Nedge(phy))
comp[match(1:(Ntip(phy)), phy$edge[,2])] <- as.factor(regimes[1:nb.tip])
comp[match((2+Ntip(phy)):(Nedge(phy)+1), phy$edge[,2])] <- as.factor(regimes[(nb.tip+1):(nb.tip+nb.node)][-1])
co <- c("black", "blue", "orange")
plot.phylo(phy, edge.color=co[comp], show.tip.label=FALSE)
load("simIdentifiable.Rsave")
regimes <- OUwie:::GetParameterPainting(phy=phy, data=data, rates=c(3,1,2), k.pars=3)
nb.tip <- Ntip(phy)
nb.node <- Nnode(phy)
comp <- numeric(Nedge(phy))
comp[match(1:(Ntip(phy)), phy$edge[,2])] <- as.factor(regimes[1:nb.tip])
comp[match((2+Ntip(phy)):(Nedge(phy)+1), phy$edge[,2])] <- as.factor(regimes[(nb.tip+1):(nb.tip+nb.node)][-1])
co <- c("black", "blue", "orange")
plot.phylo(phy, edge.color=co[comp], show.tip.label=FALSE)
## ---- eval=TRUE---------------------------------------------------------------
load("simUnidentifiable.Rsave")
check.identify(phy, data)
## ---- eval=TRUE---------------------------------------------------------------
load("simIdentifiable.Rsave")
check.identify(phy, data)
## ---- eval=FALSE--------------------------------------------------------------
# load("simsOUidentify_8")
# surfaceDatThetaR_2 <- OUwie.contour(oum.root, focal.params=c("theta_Root", "theta_2"), focal.params.upper=c(10,10), nreps=1000, n.cores=4)
# surfaceDatTheta1_2 <- OUwie.contour(oum.root, focal.params=c("theta_1", "theta_2"), focal.params.upper=c(10,10), nreps=1000, n.cores=4)
## ---- eval=FALSE--------------------------------------------------------------
# plot(surfaceDatThetaR_2, mle.point="red", levels=c(0:20*0.1), xlab=expression(theta[root]), ylab=expression(theta[2]) , xlim=c(0,5,1), ylim=c(0,5,1), col=grey.colors(21, start=0, end=1))
## ---- eval=TRUE, echo=FALSE, fig.height=3.5, fig.width=6.5, fig.cap = "A contour plot of a OUM model, with the $\\theta_{root}$ included in the model. (A) The contour for when $\\theta_{root}$ and $\\theta_{2}$ are specified as the focal parameters, and (B) shows the likelihood surface for when $\\theta_{1}$ and $\\theta_{2}$ are specified. In both cases, the likelihood surface appears as a ridge, indicating that the regimes are not separately identifiable."----
par(mfcol=c(1,2), mar=c(4,4,0.5,0.5), oma=c(1,2,1,1))
load("surfaceDatRIndentR_2.Rsave")
plot(surfaceDatRIndent, mle.point="red", levels=c(0:20*0.1), xlab=expression(theta[root]), ylab=expression(theta[2]) , xlim=c(0,5,1), ylim=c(0,5,1), col=grey.colors(21, start=0, end=1))
load("surfaceDatRIndent1_2.Rsave")
plot(surfaceDatRIndent, mle.point="red", levels=c(0:20*0.1), xlab=expression(theta[1]), ylab=expression(theta[2]) , xlim=c(0,5,1), ylim=c(0,5,1), col=grey.colors(21, start=0, end=1))
## ---- eval=TRUE, echo=FALSE, fig.height=3.5, fig.width=6.5, fig.cap = "A contour plot of a OUM model, with the $\\theta_{root}$ removed from the model. (A) The contour for when $\\theta_{1}$ and $\\theta_{2}$ are specified as the focal parameters, and (B) for when $\\theta_{1}$ and $\\theta_{3}$ are specified, the likelihood surface is sufficiently peaked. In other words, the likelihood surface no longer appears as a ridge and regimes are separately identifiable."----
par(mfcol=c(1,2), mar=c(4,4,0.5,0.5), oma=c(1,2,1,1))
load("surfaceDatNRIndent1_2.Rsave")
plot(surfaceDatNRIndent, mle.point="red", levels=c(0:20*0.1), xlab=expression(theta[1]), ylab=expression(theta[2]) , xlim=c(0,5,1), ylim=c(0,5,1), col=grey.colors(21, start=0, end=1))
load("surfaceDatNRIndent1_3.Rsave")
plot(surfaceDatNRIndent, mle.point="red", levels=c(0:20*0.1), xlab=expression(theta[1]), ylab=expression(theta[3]) , xlim=c(0,5,1), ylim=c(0,5,1), col=grey.colors(21, start=0, end=1))
## ---- eval=TRUE, echo=TRUE----------------------------------------------------
data(tworegime)
set.seed(42)
ouwiefit <- OUwie(tree, trait, model="OUM", scaleHeight=TRUE, root.station=FALSE, shift.point=0.5, algorithm="invert", quiet=TRUE, check.identify=FALSE)
## ---- eval=FALSE, echo=TRUE---------------------------------------------------
# recon <- OUwie.anc(ouwiefit)
## ---- eval=TRUE, echo=TRUE----------------------------------------------------
recon <- OUwie.anc(ouwiefit, knowledge=TRUE)
## ---- eval=TRUE, echo=FALSE, fig.height=10, fig.width=6.5, fig.cap = "A plot of the ancestral state reconstruction under an OUwie model."----
plot(recon, fsize=0.5)
## ---- eval=TRUE, echo=TRUE----------------------------------------------------
data(tworegime)
three.point <- OUwie(tree, trait, model="OUMV", root.station=FALSE, scaleHeight=TRUE, shift.point=0.5, algorithm="three.point", quiet=TRUE, check.identify=FALSE)
three.point
## ---- eval=TRUE, echo=TRUE----------------------------------------------------
data(tworegime)
invert <- OUwie(tree, trait, model="OUMV", root.station=FALSE, scaleHeight=TRUE, shift.point=0.5, algorithm="invert", quiet=TRUE)
invert
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.