Nothing
## -----------------------------------------------------------------------------
## Load the package
library( ratematrix )
## Load the data
data( "centrarchidae" )
## A list with data, phy.map, and pred
names( centrarchidae )
## -----------------------------------------------------------------------------
resp.data <- centrarchidae$data
class( resp.data )
head( resp.data )
## -----------------------------------------------------------------------------
pred.data <- centrarchidae$pred
names( pred.data )
## Table show the distribution of the data for the predictor regime.
table( pred.data )
## -----------------------------------------------------------------------------
phy <- centrarchidae$phy.map
## Drop the regimes of the stochastic map.
phy <- mergeSimmap(phy, drop.regimes = TRUE)
## ---- echo=FALSE--------------------------------------------------------------
load( system.file("extdata", "intro_vignette.RData", package = "ratematrix") )
## ---- eval=FALSE--------------------------------------------------------------
# handle <- ratematrixJointMCMC(data_BM = resp.data, data_Mk = pred.data, phy = phy
# , gen = 100000, dir = tempdir())
# mcmc <- readMCMC(handle)
## -----------------------------------------------------------------------------
plotRatematrix(mcmc)
## ---- eval=FALSE--------------------------------------------------------------
# handle2 <- ratematrixJointMCMC(data_BM = resp.data, data_Mk = pred.data, phy = phy
# , gen = 100000, dir = tempdir())
# mcmc2 <- readMCMC(handle2)
## -----------------------------------------------------------------------------
Rfactor <- checkConvergence(mcmc, mcmc2)
Rfactor
## -----------------------------------------------------------------------------
merged.mcmc <- mergePosterior(mcmc, mcmc2)
## -----------------------------------------------------------------------------
corr <- extractCorrelation(merged.mcmc)
dim( corr ) ## 1500 samples and 2 regimes.
head( corr )
## -----------------------------------------------------------------------------
hist(x = corr[,1], xlim = c(-1,1), main = "Generalist", col = "grey"
, border = "white", breaks = 20, freq = FALSE)
hist(x = corr[,2], xlim = c(-1,1), main = "Specialist", col = "grey"
, border = "white", breaks = 20, freq = FALSE)
## -----------------------------------------------------------------------------
boxplot(x = corr)
## -----------------------------------------------------------------------------
testRatematrix(chain = merged.mcmc, par = "rates")
## -----------------------------------------------------------------------------
testRatematrix(chain = merged.mcmc, par = "correlation")
## -----------------------------------------------------------------------------
names( merged.mcmc )
head( merged.mcmc$root )
## -----------------------------------------------------------------------------
names( merged.mcmc$matrix )
class( merged.mcmc$matrix$generalist )
length( merged.mcmc$matrix$generalist )
gen.rate <- merged.mcmc$matrix$generalist
spec.rate <- merged.mcmc$matrix$specialist
## The first element of the list:
gen.rate[[1]]
spec.rate[[1]]
## -----------------------------------------------------------------------------
## Extracts the rates for the first trait from the posterior distribution of generalists.
rate.tr1.gen <- sapply(gen.rate, function(x) x[1,1])
summary( rate.tr1.gen )
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.