Movatterモバイル変換


[0]ホーム

URL:


Revticulate

An Rpackage for interacting with RevBayes from R and RStudio

Installation

Revticulate can be installed in two ways. The first is via CRAN,using the defaultinstall.packages function in R:

install.packages("Revticulate")

The second is via the remotes package, a lightweight package enablinginstallation from GitHub repositories.

remotes::install_github("revbayes/Revticulate")

The GitHub repository for Revticulate contains cutting-edge featuresand may contain bugfixes, but the CRAN is known to be stable foreveryday use.

Upon first installation, Revticulate will run a package check.
This check searches for and .Renviron file that contains a RevBayespath. If the package doesn’t find this file, or finds it without thepath, the package prompts the user to useusethis::edit_r_environ(). This opens the .Renviron file,and the user will enterrb={absolute path to revbayes}.This can be edited at any time if there are multiple installs on thesystem, or if you recompile RevBayes and want to use a new version.

Before using Revticulate in knitr, make sure the following is in yoursetup chunk:

library(Revticulate)knitRev()

Using RevBayes in a KnitRChunk

RevBayes can be used in a KnitR chunk by changing the header torb instead ofr. In the below chunk, we createan object calledexample and use the assignment operator togive it the value 1. Then we print it.

example <- 1.0example

This is not an overtly useful thing to do, however. Let’s erase theprevious chunk using theclearRev() function. This removesprior code from the RevBayes environment. Very handy if you make amistake!

clearRev()

We could, instead, choose to do something a little more useful. Howabout reading in a data matrix and making a quick starting tree?

morpho <- readDiscreteCharacterData("bears.nex")num_taxa <- morpho.size()num_branches <- 2 * num_taxa - 2taxa <- morpho.names()br_len_lambda ~ dnExp(0.2)phylogeny ~ dnUniformTopologyBranchLength(taxa, branchLengthDistribution=dnExponential(br_len_lambda))phylogeny

Anything entered in anrb block will be interpreted asRev code, and all the normal Rev syntax will apply. For a nice overviewof Rev language and syntax, please seethistutorial.

One thing researchers are often interested in doing is making anobject in Rev and then viewing it in R. The best way to accomplish thisis with thedoRev() function. When using this function, theRevCode you’d like to run goes in the parentheses of thedoRev function. These are then exportable to R. In thisexample, we load the dataset used in the published tutorial “Estimatinga time-calibrated phylogeny of fossil and extant taxa using RevBayes”{@barido2020estimating}.

doRev(input = 'morpho <- readDiscreteCharacterData("bears.nex")num_taxa <- morpho.size()num_branches <- 2 * num_taxa - 2taxa <- morpho.names()br_len_lambda ~ dnExp(0.2)phylogeny ~ dnUniformTopologyBranchLength(taxa, branchLengthDistribution=dnExponential(br_len_lambda))phylogeny')

ThedoRev function is then used to extract the object.Note that knitr chunks can only have one language type. Thus, to use aRev Object in another chunk, it must be exported. In this case, aphylogeny is not a simple numeric type, and Revticualte automates thecoercion from a string to a Newick tree that can be read by Phytools orsimilar.

phylogeny <- doRev("phylogeny")phytools::plotTree(phylogeny)

We may choose to clear RevBayes objects out of memory so that theyare not being consistently echoed to the screen.

clearRev()

One nice facet of having RevBayes running in an R notebook is theability to flip to visualizations of the different distributions we use.For example, here is the code for a common parameterization of thediscrete Gamma distribution on site rates.

alpha_morpho ~ dnUniform( 0, 1E6 );rates_morpho := fnDiscretizeGamma( alpha_morpho, alpha_morpho, 4 )alpha_morpho

If you aren’t a big stats person, this might not mean much to you, interms of what this distribution actually looks like. But it is importantto develop intuitions for what common distributions look like and whatthis says about our data. So, we can use R’s built-in graphicscapabilities to have a look at what 1000 draws from this gamma will looklike.

doRev('alpha_morpho ~ dnUniform( 0, 1E6 );rates_morpho := fnDiscretizeGamma( alpha_morpho, alpha_morpho, 4 )')library(ggplot2)alpha_value <- doRev("alpha_morpho")alpha_valuedraws <- rgamma(1000, shape = alpha_value, rate = alpha_value)hist(draws, xlab = "Value")

Simulate coin flips

It’s adviseable if you’re switching gears to a new activity to clearthe Rev environment of workspace objects from old activities:

clearRev()

Note thatclearRev is an R function, and must beexecuted in an R chunk.

# The number of coin flipsn <- 100# The number of headsx <- 50x

Initialize Markov Chain

We have to start MCMC off with some initial parameter values. One wayto do this is to randomly draw values of the parameter (\(p\)) from the prior distribution. We assumea flat beta prior distribution (\(\alpha =1\) and\(\beta = 1\)).

alpha <- 1beta <- 1p <- rbeta (n=1, alpha, beta) [1]p

Likelihood Function

We next specify the likelihood function. We use the binomialprobability for the likelihood function. Since the likelihood is definedonly for\(p\) between 0 and 1, wereturn 0 if\(p\) is outside thisrange.

{rb include=FALSE} function likelihood(p) { if(p < 0 || p > 1) return 0 l<-dbinomial(x,p,n,log=false) return l }

The function can then be executed in the next cell:

likelihood(p)

The RepRev() environment

The functionrepRev() can be called in the console (orin non-RStudio versions of R) to use RevBayes directly to programinteractively. TherepRev() environment is denoted withrb>>>. To exit, type Ctrl + C. It is notcompatible with KnitR, being a console tool.

# repRev()# rb>>> 1+2# [1] 3

Running Long Computations

Perhaps we would like to run a longer computation from R usingRevticulate. For example, maybe we have made an MCMC script, and wouldlike to run it using Revticulate, and then automatically process theoutput in R. In the below example, we will run an MCMC and automaticallyevaluate it for convergence using the package conveience.

Included with the package, we have a script called ‘mcmc_mk.Rev,’which runs a short phylogenetic estimation for a small morphologicaldataset from bears using the Mk model {@Lewis2001}. Once the computation iscomplete, convergence is diagnosed with the R package convenience {@fabreti2021}. Pleasenote that this will take about 5 minutes if executed.

library(convenience)callRevFromTerminal("mcmc_mk.Rev")checkConvergence(path = "vignettes/output/")

[8]ページ先頭

©2009-2025 Movatter.jp