Nothing
`bpec` <- function(seqsFile, coordsFile, dims = 2, iter = 100000, postSamples = 100, maxMig = 5, ds = 0, colorCode = c(7,5,6,3,2,8)) {
rawSeqs = read.nexus.data(seqsFile)
coordsLocs = bpec.loadCoords(coordsFile)
#run the Markov chain Monte Carlo sampler
bpecout = bpec.mcmc(rawSeqs, coordsLocs, maxMig, iter, ds, postSamples, dims)
return(bpecout)
}
`plot.bpec` <- function(x, GoogleEarth = 0, colorCode = c(7,5,6,3,2,8), ... ) {
# plot geographical cluster contour map
bpec.contourPlot(x, GoogleEarth = 0, colorCode = colorCode)
# plot tree network with cluster indicators
bpec.Tree = bpec.treePlot(x, colorCode = colorCode)
# now also plot the environmental covariates
par(mfrow = c(2, 3)) #split the plot window into 2x3 to fit all the covariates
if (x$input$coordsDims > 2) {
bpec.covariatesPlot(x, colorCode = colorCode)
}
if (GoogleEarth == 1) {
bpec.Geo = bpec.geoTree(x,file="GoogleEarthTree.kml")
}
}
`summary.bpec` <- function(object,...) {
cat(paste('bpec ran with the following settings:\n'))
cat(paste('Number of input sequences',object$input$seqCountOrig,'\n'))
cat(paste('Input sequence length',object$input$seqLengthOrig,'\n'))
cat(paste('Number of iterations',object$input$iter,'\n'))
cat(paste('Number of saved iterations',dim(object$clust$sampleMeans)[3]/2,'\n'))
cat(paste('Dimensions',object$input$coordsDims,'\n'))
cat(paste('Parsimony relaxation',object$input$ds,'\n'))
cat(paste('Maximum number of migrations',length(object$tree$migProbs) - 1,'\n'))
cat(paste('\nThe results of bpec are: \n'))
cat(paste('Number of haplotypes (including missing)',dim(object$preproc$seq)[1],'\n'))
cat(paste('of which',dim(object$preproc$seq)[1]-object$input$seqCountOrig,'are missing\n'))
cat(paste('Effective sequence length',object$preproc$seqLength,'\n'))
cat(paste('The most likely number of migrations is',which.max(object$tree$migProbs)))
seqLabels = object$preproc$seqsFile[object$preproc$seqLabels]
maxMig = length(object$tree$migProbs) - 1
rootProbMean = (object$tree$rootProbs[1, ] + object$tree$rootProbs[2, ]) / 2
if (which.max(rootProbMean) <= length(seqLabels)) {
root = seqLabels[which.max(rootProbMean)]
writeLines(paste("\nThe most likely root node is ", root, sep=""))
}
if (which.max(rootProbMean) > length(seqLabels)) {
root = which.max(rootProbMean)
writeLines(paste("\nThe most likely root node is extinct", sep=""))
}
maxVar = 0
chainMeans = array(0, dim = c(dim(object$clust$sampleMeans)[1], dim(object$clust$sampleMeans)[2], 2))
flag = 0
postSamples = length(object$clust$sampleMeans[1, 1, ]) / 2
for(i in 1:(maxMig+1)) {
for(j in 1:object$input$coordsDims) {
for(l in 1:2) {
chainMeans[j, i, l] = mean(object$clust$sampleMeans[j, i, (1 + (l-1) * postSamples):(l * postSamples)], na.rm = TRUE)
}
if (sum(is.na(chainMeans[j, i, ])) == 0) {
if (abs(chainMeans[j, i, 1] - chainMeans[j, i, 2]) > 0.05 * (max(object$input$coordsLocs[, j]) - min(object$input$coordsLocs[, j])) && mean(is.na(object$clust$sampleMeans[j, i, ])) < 0.5) {
writeLines("NO CLUSTER CONVERGENCE: You need to re-run the sampler with more iterations")
flag = 1
break
}
}
}
if (flag == 1) {
break
}
}
rootProbs1 = object$tree$rootProbs[1:object$preproc$count]
rootProbs2 = object$tree$rootProbs[(object$preproc$count + 1):(object$preproc$count * 2)]
if (max(abs(rootProbs1 / sum(rootProbs1) - rootProbs2 / sum(rootProbs2))) > 0.5 / object$preproc$count) {
writeLines("NO ROOT CONVERGENCE: You need to re-run the sampler with more iterations");
}
}
`print.bpec` <- function(x,...) {
cat(paste('\nThe results of bpec are: \n'))
cat(paste('Number of haplotypes (including missing)',dim(x$preproc$seq)[1],'\n'))
cat(paste('of which',dim(x$preproc$seq)[1]-x$input$seqCountOrig,'are missing\n'))
cat(paste('Effective sequence length',x$preproc$seqLength,'\n'))
cat(paste('The most likely number of migrations is',which.max(x$tree$migProbs)))
seqLabels = x$preproc$seqsFile[x$preproc$seqLabels]
maxMig = length(x$tree$migProbs) - 1
rootProbMean = (x$tree$rootProbs[1, ] + x$tree$rootProbs[2, ]) / 2
if (which.max(rootProbMean) <= length(seqLabels)) {
root = seqLabels[which.max(rootProbMean)]
writeLines(paste("\nThe most likely root node is ", root, sep=""))
}
if (which.max(rootProbMean) > length(seqLabels)) {
root = which.max(rootProbMean)
writeLines(paste("\nThe most likely root node is extinct", sep=""))
}
maxVar = 0
chainMeans = array(0, dim = c(dim(x$clust$sampleMeans)[1], dim(x$clust$sampleMeans)[2], 2))
flag = 0
postSamples = length(x$clust$sampleMeans[1, 1, ]) / 2
for(i in 1:(maxMig+1)) {
for(j in 1:x$input$coordsDims) {
for(l in 1:2) {
chainMeans[j, i, l] = mean(x$clust$sampleMeans[j, i, (1 + (l-1) * postSamples):(l * postSamples)], na.rm = TRUE)
}
if (sum(is.na(chainMeans[j, i, ])) == 0) {
if (abs(chainMeans[j, i, 1] - chainMeans[j, i, 2]) > 0.05 * (max(x$input$coordsLocs[, j]) - min(x$input$coordsLocs[, j])) && mean(is.na(x$clust$sampleMeans[j, i, ])) < 0.5) {
writeLines("NO CLUSTER CONVERGENCE: You need to re-run the sampler with more iterations")
flag = 1
break
}
}
}
if (flag == 1) {
break
}
}
rootProbs1 = x$tree$rootProbs[1:x$preproc$count]
rootProbs2 = x$tree$rootProbs[(x$preproc$count + 1):(x$preproc$count * 2)]
if (max(abs(rootProbs1 / sum(rootProbs1) - rootProbs2 / sum(rootProbs2))) > 0.5 / x$preproc$count) {
writeLines("NO ROOT CONVERGENCE: You need to re-run the sampler with more iterations");
}
}
`mean.bpec` <- function(x, ...) {
cat('The mean cluster centres are:\n')
print(apply(x$clust$sampleMeans, c(1,2), 'mean', na.rm = TRUE))
cat('The mean cluster covariances are:\n')
print(apply(x$clust$sampleCovs, c(1,2,3), 'mean', na.rm = TRUE))
cat('The mean root posterior probabilities are:\n')
print(apply(x$tree$rootProbs, 2, 'mean', na.rm = TRUE))
}
input <- function(bpecout) UseMethod("input")
`input.bpec` <- function(bpecout) {
output = list()
output$seqCountOrig = bpecout$input$seqCountOrig
output$seqLengthOrig = bpecout$input$seqLengthOrig
output$iter = bpecout$input$iter
output$ds = bpecout$input$ds
output$coordsLocs = bpecout$input$coordsLocs
output$coordsDims = bpecout$input$coordsDims
output$locNo = bpecout$input$locNo
output$locData = bpecout$input$locData
return(output)
}
preproc <- function(bpecout) UseMethod("preproc")
`preproc.bpec` <- function(bpecout){
output = list()
output$seq = bpecout$preproc$seq
output$seqsFile = bpecout$preproc$seqsFile
output$seqLabels = bpecout$preproc$seqLabels
output$seqIndices = bpecout$preproc$seqIndices
output$seqLength = bpecout$preproc$seqLength
output$noSamples = bpecout$preproc$noSamples
output$count = bpecout$preproc$count
return(output)
}
#setGeneric("tree", function(bpecout) standardGeneric("tree"))
output.tree <- function(bpecout) UseMethod("output.tree")
`output.tree.bpec` <- function(bpecout){
output = list()
output$clado = bpecout$tree$clado
output$levels = bpecout$tree$levels
output$edgeTotalProb = bpecout$tree$edgeTotalProb
output$rootProbs = bpecout$tree$rootProbs
output$treeEdges = bpecout$tree$treeEdges
output$rootLocProbs = bpecout$tree$rootLocProbs
output$migProbs = bpecout$tree$migProbs
return(output)
}
output.clust <- function(bpecout) UseMethod("output.clust")
`output.clust.bpec` <- function(bpecout){
output = list()
output$sampleMeans = bpecout$clust$sampleMeans
output$sampleCovs = bpecout$clust$sampleCovs
output$sampleIndices = bpecout$clust$sampleIndices
output$clusterProbs = bpecout$clust$clusterProbs
return(output)
}
#setGeneric("mcmc", function(x) standardGeneric("mcmc"))
#mcmc <- function(x) UseMethod("mcmc")
#setMethod("mcmc", "bpec",
#function(x, ...) {
# output = list()
# output$MCMCparams = bpecout$mcmc$MCMCparams
# output$MCChainMeans = bpecout$mcmc$MCChainMeans
# output$codaInput = bpecout$mcmc$codaInput
#
# return(output)
# }
#)
#
#setMethod("mcmc", "list",
#function(x, ...) {
# mcmc.list(x)
# }
#)
#setGeneric("mcmc", function(x) standardGeneric("mcmc"))
output.mcmc <- function(bpecout) UseMethod("output.mcmc")
`output.mcmc.bpec` <- function(bpecout){
output = list()
output$MCMCparams = bpecout$mcmc$MCMCparams
output$MCChainMeans = bpecout$mcmc$MCChainMeans
output$codaInput = bpecout$mcmc$codaInput
return(output)
}
#`mcmc.list` <- function(bpecout){
# return(coda::mcmc(bpecout))
#}
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.