added in function to return the mll of the root node of a hierarchy

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: fastbaps
 Title: A fast genetic clustering algorithm that approximates a Dirichlet Process Mixture model
 Version: 1.0.0
-Authors@R: person("Gerry", "Tonkin-Hill", email = "[email protected]", role = c("aut", "cre"))
+Authors@R: person("Gerry", "Tonkin-Hill", email = "[email protected]", role = c("aut", "cre"))
 Description: Takes a multiple sequence alignment as input and clusters according to the 'no-admixture' model. 
   It combines ideas from the Bayesian Hierarchical Clustering algorithm of Heller et al. <doi:10.1145/1102351.1102389>
   and hierBAPS <doi:10.1093/molbev/mst028> to produce a rapid and accurate clustering algorithm.

NAMESPACE

@@ -9,6 +9,7 @@ export(import_fasta_sparse_nt)
 export(multi_level_best_baps_partition)
 export(multi_res_baps)
 export(optimise_prior)
+export(root_mll)
 export(snp_dist)
 export(snp_similarity)
 import(Matrix)

R/boot_fast_baps.R

@@ -16,12 +16,13 @@
 #'
 #' @examples
 #'
+#' \dontrun{
 #' fasta.file.name <- system.file("extdata", "seqs.fa", package = "fastbaps")
 #' sparse.data <- import_fasta_sparse_nt(fasta.file.name)
 #' boot.result <- boot_fast_baps(sparse.data)
 #' dendro <- as.dendrogram(fast_baps(sparse.data))
 #' heatmap(boot.result, dendro, dendro)
-#'
+#' }
 #'
 #' @export
 boot_fast_baps <- function(sparse.data, k.init=NULL, n.replicates=100, hc.method='ward',

R/fix_clusters.R

@@ -13,13 +13,15 @@
 #'
 #' @examples
 #'
+#' \dontrun{
 #' fasta.file.name <- system.file("extdata", "seqs.fa", package = "fastbaps")
 #' sparse.data <- import_fasta_sparse_nt(fasta.file.name)
 #' sim.matrix <- snp_similarity(sparse.data)
 #' x <- as.dist(1-sim.matrix/max(sim.matrix))
 #' phylo <- ape::as.phylo(hclust(x, method="average"))
 #' clusters <- best_baps_partition(sparse.data, phylo)
 #' clusters <- fix_clusters(sparse.data, clusters)
+#' }
 #'
 #' @export
 fix_clusters <- function(sparse.data, clusters, n.iterations=1, quiet=FALSE){

R/root_mll.R

@@ -0,0 +1,35 @@
+#' root_mll
+#'
+#' Function to calculate the marginal log likelihood of the root of a hierarchy for use in model comparison.
+#'
+#' @import Matrix
+#'
+#'
+#' @param sparse.data a sparse SNP data object returned from import_fasta_sparse_nt
+#' @param h a hclust object representing the hierarchical clustering
+#' @param quiet suppress the printing of extra information (default=FALSE)
+#'
+#' @return the marginal log likelihood of the root node
+#'
+#' @examples
+#' fasta.file.name <- system.file("extdata", "seqs.fa", package = "fastbaps")
+#' fasta.file.name <- system.file("extdata", "seqs.fa", package = "fastbaps")
+#' sparse.data <- import_fasta_sparse_nt(fasta.file.name)
+#' baps.hc <- fast_baps(sparse.data)
+#' root_mll(sparse.data, baps.hc)
+#'
+#' @export
+root_mll <- function(sparse.data, h, quiet=FALSE){
+
+  # Check inputs
+  if(!is.list(sparse.data)) stop("Invalid value for sparse.data! Did you use the import_fasta_sparse_nt function?")
+  if(!(class(sparse.data$snp.matrix)=="dgCMatrix")) stop("Invalid value for sparse.data! Did you use the import_fasta_sparse_nt function?")
+  if(!is.numeric(sparse.data$consensus)) stop("Invalid value for sparse.data! Did you use the import_fasta_sparse_nt function?")
+  if(!is.matrix(sparse.data$prior)) stop("Invalid value for sparse.data! Did you use the import_fasta_sparse_nt function?")
+  if(!(class(h) %in% c("hclust", "phylo"))) stop("Invalid value for h! Should be a hclust or phylo object!")
+
+  llks <- tree_llk(sparse.data, h$merge)
+  root.mll <- llks$ptree[length(llks$ptree)]
+
+  return(root.mll)
+}

fastbaps.Rproj

@@ -18,4 +18,5 @@ StripTrailingWhitespace: Yes
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageCheckArgs: --as-cran
 PackageRoxygenize: rd,collate,namespace