Nothing
R/assemblePriorNode.R
In ips: Interfaces to Phylogenetic Software in R
Defines functions
assemblePriorNode
## This code is part of the ips package
## © C. Heibl 2020 (last update 2020-03-13)
assemblePriorNode <- function(x){
# ## Loop over partitions to add childen to <prior>
# ## ----------------------------------------------
# for (i in seq_along(id)){
# prior <- assemblePriorNode(id, i, link.trees, prior)
# }
i <- 1
# ## Define taxonsets
# ## ----------------
# if (!missing(taxonset)) {
#
# ## Loop over taxonsets
# ## -------------------
# for (j in seq_along(taxonset)){
#
# taxonSetPrior <- assembleTaxonSetPrior(taxonset = taxonset[[j]],
# setname = names(taxonset)[j],
# id = id, i = i, link.trees = link.trees)
# if (!is.null(taxonSetPrior)){
# prior <- addChildren(prior, kids = list(taxonSetPrior))
# }
# }
# }
prior <- xmlNode("distribution",
attrs = c(id = "prior",
spec = "util.CompoundDistribution"))
if (x$tree == "Yule"){
kids <- list(
xmlNode("distribution",
attrs = c(id = paste0("YuleModel.t:", x$partition[i]),
spec = "beast.evolution.speciation.YuleModel",
birthDiffRate = paste0("@birthRate.t:", x$partition[i]),
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("YuleBirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@birthRate.t:", x$partition[i])),
.children = uniformPrior(x))
)
}
if (x$tree == "Calibrated Yule"){
kids <- list(
xmlNode("distribution",
attrs = c(birthRate = paste0("@birthRateY.t:", x$partition[i]),
id = paste0("CalibratedYuleModel.t:", x$partition[i]),
spec = "beast.evolution.speciation.CalibratedYuleModel",
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("CalibratedYuleBirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@birthRateY.t:", x$partition[i])),
.children = uniformPrior(x, upper = "1000.0"))
)
}
if (x$tree == "Birth Death"){
kids <- list(
xmlNode("distribution",
attrs = c(birthDiffRate = paste0("@BDBirthRate.t:", x$partition[i]),
id = paste0("BirthDeath.t:", x$partition[i]),
relativeDeathRate = paste0("@BDDeathRate.t:", x$partition[i]),
spec = "beast.evolution.speciation.BirthDeathGernhard08Model",
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("BirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@BDBirthRate.t:", x$partition[i])),
.children = uniformPrior(x, upper = "1000.0")),
xmlNode("prior",
attrs = c(id = paste0("DeathRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@BDDeathRate.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
)
}
if (x$tree == "Fossilized Birth Death"){
kids <- list(
xmlNode("distribution",
attrs = c(id = paste0("FBD.t:", x$partition[i]),
spec = "beast.evolution.speciation.SABirthDeathModel",
conditionOnRhoSampling = "true",
diversificationRate = paste0("@diversificationRateFBD.t:", x$partition[i]),
origin = paste0("@originFBD.t:", x$partition[i]),
samplingProportion = paste0("@samplingProportionFBD.t:", x$partition[i]),
tree = paste0("@Tree.t:", x$partition[i]),
turnover = paste0("@turnoverFBD.t:", x$partition[i])),
.children = lapply(c("rFBD.t", "rhoFBD.t"), parameter, i = i, x = x)),
xmlNode("prior",
attrs = c(id = paste0("diversificationRatePriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@diversificationRateFBD.t:", x$partition[i])),
.children = uniformPrior(x)),
xmlNode("prior",
attrs = c(id = paste0("originPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@originFBD.t:", x$partition[i])),
.children = uniformPrior(x)),
xmlNode("prior",
attrs = c(id = paste0("samplingProportionPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@samplingProportionFBD.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL)),
xmlNode("prior",
attrs = c(id = paste0("turnoverPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@turnoverFBD.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL)))
if (x$clock == "Strict Clock"){
kids <- c(kids, list(xmlNode("prior",
attrs = c(id = paste0("ClockPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@clockRate.c:", x$partition[i])),
.children = uniformPrior(x))))
}
if (x$clock == "Relaxed Clock Exponential"){
UCMeanRatePrior.c <- xmlNode("prior",
attrs = c(id = paste0("UCMeanRatePrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucedMean.c:", x$partition[i])),
.children = uniformPrior(x))
kids <- c(kids, list(UCMeanRatePrior.c))
}
}
## RELAXED CLOCK LOG NORMAL
## ------------------------
if (x$clock == "Relaxed Clock Log Normal"){
## Default: 'ucldStdevPrior.c'
clock.kids <- list(xmlNode("prior",
attrs = c(id = paste0("ucldStdevPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucldStdev.c:", x$partition[i])),
.children = list(gammaPrior(x, "0.5396", "0.3819"))))
## Have to add 'MeanRatePrior.c', when tip dating .children = list(0.3819)))))))
if (!is.na(x$tip.dates)){
clock.kids <- c(list(xmlNode("prior", attrs = c(id = paste0("MeanRatePrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucldMean.c:", x$partition[i])),
.children = uniformPrior(x))),
clock.kids)
}
kids <- c(kids, clock.kids)
}
## RANDOM LOCAL CLOCK
## ------------------------
if (x$clock == "Random Local Clock"){
clock.kids <- list(xmlNode("prior",
attrs = c(id = paste0("RRateChangesPrior.c:", x$partition[i]),
name = "distribution"),
.children = list(
xmlNode("x", attrs = c(id = paste0("RRateChanges.c:", x$partition[i]),
spec = "util.Sum"),
value = xmlNode("arg", attrs = c(idref = paste0("Indicators.c:", x$partition[i])))),
xmlNode("distr", attrs = c(id = "Poisson.0",
spec = "beast.math.distributions.Poisson"),
value = xmlNode("parameter",
attrs = c(id = realParameter(x),
spec = "parameter.RealParameter",
estimate = "false",
name = "lambda"),
value = "0.6931471805599453"))
)),
xmlNode("prior",
attrs = c(id = paste0("RRatesPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@clockrates.c:", x$partition[i])),
.children = list(gammaPrior(x, "0.5396", "0.3819"))))
kids <- c(kids, clock.kids)
}
if (x$subst.model == "HKY"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
KappaPrior.s <- list(
xmlNode("parameter",
attrs = c(estimate = "false",
id = "RealParameter.17",
name = "M",
spec = "parameter.RealParameter"),
value = "1.0"),
xmlNode("parameter",
attrs = c(estimate = "false",
id = "RealParameter.18",
name = "S",
spec = "parameter.RealParameter"),
value = "1.25"))
KappaPrior.s <- xmlNode("LogNormal",
attrs = c(id = "LogNormalDistributionModel.2", name = "distr"),
.children = KappaPrior.s)
KappaPrior.s <- xmlNode("prior",
attrs = c(id = paste0("KappaPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa.s:", x$partition[i])),
value = KappaPrior.s)
kids <- c(kids, list(FrequenciesPrior.s, KappaPrior.s))
}
if (x$subst.model == "TN93"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
kappa1Prior.s <- xmlNode("prior",
attrs = c(id = paste0("kappa1Prior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa1.s:", x$partition[i])),
value = LogNormalPrior(x, 1.0, 1.25))
kappa2Prior.s <- xmlNode("prior",
attrs = c(id = paste0("kappa2Prior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa2.s:", x$partition[i])),
value = LogNormalPrior(x, 1.0, 1.25))
kids <- c(kids, list(FrequenciesPrior.s, kappa1Prior.s, kappa2Prior.s))
}
if (x$subst.model == "GTR"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
RateACPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateACPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAC.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateAGPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateAGPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAG.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "20.0")) ## RLC beta of AC and AG exchanged in BEAUTI
RateATPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateATPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAT.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateCGPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateCGPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateCG.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateGTPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateGTPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateGT.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
kids <- c(kids, list(FrequenciesPrior.s, RateACPrior.s, RateAGPrior.s,
RateATPrior.s, RateCGPrior.s, RateGTPrior.s))
}
prior <- addChildren(prior, kids = kids)
prior
}
Try the ips package in your browser
Any scripts or data that you put into this service are public.
ips documentation built on May 29, 2024, 4:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions assemblePriorNode
## This code is part of the ips package
## © C. Heibl 2020 (last update 2020-03-13)
assemblePriorNode <- function(x){
# ## Loop over partitions to add childen to <prior>
# ## ----------------------------------------------
# for (i in seq_along(id)){
# prior <- assemblePriorNode(id, i, link.trees, prior)
# }
i <- 1
# ## Define taxonsets
# ## ----------------
# if (!missing(taxonset)) {
#
# ## Loop over taxonsets
# ## -------------------
# for (j in seq_along(taxonset)){
#
# taxonSetPrior <- assembleTaxonSetPrior(taxonset = taxonset[[j]],
# setname = names(taxonset)[j],
# id = id, i = i, link.trees = link.trees)
# if (!is.null(taxonSetPrior)){
# prior <- addChildren(prior, kids = list(taxonSetPrior))
# }
# }
# }
prior <- xmlNode("distribution",
attrs = c(id = "prior",
spec = "util.CompoundDistribution"))
if (x$tree == "Yule"){
kids <- list(
xmlNode("distribution",
attrs = c(id = paste0("YuleModel.t:", x$partition[i]),
spec = "beast.evolution.speciation.YuleModel",
birthDiffRate = paste0("@birthRate.t:", x$partition[i]),
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("YuleBirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@birthRate.t:", x$partition[i])),
.children = uniformPrior(x))
)
}
if (x$tree == "Calibrated Yule"){
kids <- list(
xmlNode("distribution",
attrs = c(birthRate = paste0("@birthRateY.t:", x$partition[i]),
id = paste0("CalibratedYuleModel.t:", x$partition[i]),
spec = "beast.evolution.speciation.CalibratedYuleModel",
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("CalibratedYuleBirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@birthRateY.t:", x$partition[i])),
.children = uniformPrior(x, upper = "1000.0"))
)
}
if (x$tree == "Birth Death"){
kids <- list(
xmlNode("distribution",
attrs = c(birthDiffRate = paste0("@BDBirthRate.t:", x$partition[i]),
id = paste0("BirthDeath.t:", x$partition[i]),
relativeDeathRate = paste0("@BDDeathRate.t:", x$partition[i]),
spec = "beast.evolution.speciation.BirthDeathGernhard08Model",
tree = paste0("@Tree.t:", x$partition[i]))),
xmlNode("prior",
attrs = c(id = paste0("BirthRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@BDBirthRate.t:", x$partition[i])),
.children = uniformPrior(x, upper = "1000.0")),
xmlNode("prior",
attrs = c(id = paste0("DeathRatePrior.t:", x$partition[i]),
name = "distribution",
x = paste0("@BDDeathRate.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
)
}
if (x$tree == "Fossilized Birth Death"){
kids <- list(
xmlNode("distribution",
attrs = c(id = paste0("FBD.t:", x$partition[i]),
spec = "beast.evolution.speciation.SABirthDeathModel",
conditionOnRhoSampling = "true",
diversificationRate = paste0("@diversificationRateFBD.t:", x$partition[i]),
origin = paste0("@originFBD.t:", x$partition[i]),
samplingProportion = paste0("@samplingProportionFBD.t:", x$partition[i]),
tree = paste0("@Tree.t:", x$partition[i]),
turnover = paste0("@turnoverFBD.t:", x$partition[i])),
.children = lapply(c("rFBD.t", "rhoFBD.t"), parameter, i = i, x = x)),
xmlNode("prior",
attrs = c(id = paste0("diversificationRatePriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@diversificationRateFBD.t:", x$partition[i])),
.children = uniformPrior(x)),
xmlNode("prior",
attrs = c(id = paste0("originPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@originFBD.t:", x$partition[i])),
.children = uniformPrior(x)),
xmlNode("prior",
attrs = c(id = paste0("samplingProportionPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@samplingProportionFBD.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL)),
xmlNode("prior",
attrs = c(id = paste0("turnoverPriorFBD.t:", x$partition[i]),
name = "distribution",
x = paste0("@turnoverFBD.t:", x$partition[i])),
.children = uniformPrior(x, upper = NULL)))
if (x$clock == "Strict Clock"){
kids <- c(kids, list(xmlNode("prior",
attrs = c(id = paste0("ClockPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@clockRate.c:", x$partition[i])),
.children = uniformPrior(x))))
}
if (x$clock == "Relaxed Clock Exponential"){
UCMeanRatePrior.c <- xmlNode("prior",
attrs = c(id = paste0("UCMeanRatePrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucedMean.c:", x$partition[i])),
.children = uniformPrior(x))
kids <- c(kids, list(UCMeanRatePrior.c))
}
}
## RELAXED CLOCK LOG NORMAL
## ------------------------
if (x$clock == "Relaxed Clock Log Normal"){
## Default: 'ucldStdevPrior.c'
clock.kids <- list(xmlNode("prior",
attrs = c(id = paste0("ucldStdevPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucldStdev.c:", x$partition[i])),
.children = list(gammaPrior(x, "0.5396", "0.3819"))))
## Have to add 'MeanRatePrior.c', when tip dating .children = list(0.3819)))))))
if (!is.na(x$tip.dates)){
clock.kids <- c(list(xmlNode("prior", attrs = c(id = paste0("MeanRatePrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@ucldMean.c:", x$partition[i])),
.children = uniformPrior(x))),
clock.kids)
}
kids <- c(kids, clock.kids)
}
## RANDOM LOCAL CLOCK
## ------------------------
if (x$clock == "Random Local Clock"){
clock.kids <- list(xmlNode("prior",
attrs = c(id = paste0("RRateChangesPrior.c:", x$partition[i]),
name = "distribution"),
.children = list(
xmlNode("x", attrs = c(id = paste0("RRateChanges.c:", x$partition[i]),
spec = "util.Sum"),
value = xmlNode("arg", attrs = c(idref = paste0("Indicators.c:", x$partition[i])))),
xmlNode("distr", attrs = c(id = "Poisson.0",
spec = "beast.math.distributions.Poisson"),
value = xmlNode("parameter",
attrs = c(id = realParameter(x),
spec = "parameter.RealParameter",
estimate = "false",
name = "lambda"),
value = "0.6931471805599453"))
)),
xmlNode("prior",
attrs = c(id = paste0("RRatesPrior.c:", x$partition[i]),
name = "distribution",
x = paste0("@clockrates.c:", x$partition[i])),
.children = list(gammaPrior(x, "0.5396", "0.3819"))))
kids <- c(kids, clock.kids)
}
if (x$subst.model == "HKY"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
KappaPrior.s <- list(
xmlNode("parameter",
attrs = c(estimate = "false",
id = "RealParameter.17",
name = "M",
spec = "parameter.RealParameter"),
value = "1.0"),
xmlNode("parameter",
attrs = c(estimate = "false",
id = "RealParameter.18",
name = "S",
spec = "parameter.RealParameter"),
value = "1.25"))
KappaPrior.s <- xmlNode("LogNormal",
attrs = c(id = "LogNormalDistributionModel.2", name = "distr"),
.children = KappaPrior.s)
KappaPrior.s <- xmlNode("prior",
attrs = c(id = paste0("KappaPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa.s:", x$partition[i])),
value = KappaPrior.s)
kids <- c(kids, list(FrequenciesPrior.s, KappaPrior.s))
}
if (x$subst.model == "TN93"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
kappa1Prior.s <- xmlNode("prior",
attrs = c(id = paste0("kappa1Prior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa1.s:", x$partition[i])),
value = LogNormalPrior(x, 1.0, 1.25))
kappa2Prior.s <- xmlNode("prior",
attrs = c(id = paste0("kappa2Prior.s:", x$partition[i]),
name = "distribution",
x = paste0("@kappa2.s:", x$partition[i])),
value = LogNormalPrior(x, 1.0, 1.25))
kids <- c(kids, list(FrequenciesPrior.s, kappa1Prior.s, kappa2Prior.s))
}
if (x$subst.model == "GTR"){
FrequenciesPrior.s <- xmlNode("prior",
attrs = c(id = paste0("FrequenciesPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@freqParameter.s:", x$partition[i])),
.children = uniformPrior(x, upper = NULL))
RateACPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateACPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAC.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateAGPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateAGPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAG.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "20.0")) ## RLC beta of AC and AG exchanged in BEAUTI
RateATPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateATPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateAT.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateCGPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateCGPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateCG.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
RateGTPrior.s <- xmlNode("prior",
attrs = c(id = paste0("RateGTPrior.s:", x$partition[i]),
name = "distribution",
x = paste0("@rateGT.s:", x$partition[i])),
value = gammaPrior(x, "0.05", "10.0"))
kids <- c(kids, list(FrequenciesPrior.s, RateACPrior.s, RateAGPrior.s,
RateATPrior.s, RateCGPrior.s, RateGTPrior.s))
}
prior <- addChildren(prior, kids = kids)
prior
}
Try the ips package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.