Fbms

.Rbuildignore

@@ -1,5 +1,9 @@
 .idea
 test
+Make_data
+tests_current
+README.md
 .DS_Store
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^CRAN-SUBMISSION$

.gitignore

@@ -1,3 +1,7 @@
 inst/doc
 run_algorithm.R
 .Rproj.user
+.idea
+test
+*.o
+*.so

DESCRIPTION

@@ -1,20 +1,27 @@
-Package: GMJMCMC
+Package: FBMS
 Type: Package
-Title: Genetically Modified Mode Jumping MCMC
-Version: 1.0
-Date: 2023-06-30
+Title: Flexible Bayesian Model Selection and Model Averaging
+Version: 1.1
+Date: 2025-02-26
 Encoding: UTF-8
-Author: Jon Lachamnn
-Maintainer: Jon Lachmann <[email protected]>
-Description: Implements the GMJMCMC and MJMCMC algorithms used to fit Bayesian Generalized Non-Linear Models (BGNLM).
+Language: en-US
+Authors@R: c(
+    person("Jon", "Lachmann", email = "[email protected]", role = c("cre", "aut")),
+    person("Aliaksandr", "Hubin", email = "[email protected]", role = "aut")
+    )
+Description: Implements the Mode Jumping Markov Chain Monte Carlo algorithm described in <doi:10.1016/j.csda.2018.05.020> and its Genetically Modified counterpart described in <doi:10.1613/jair.1.13047> as well as the sub-sampling versions described in <doi:10.1016/j.ijar.2022.08.018> for flexible Bayesian model selection and model averaging.
 License: GPL-2
 Depends:
     R (>= 3.5.0),
     fastglm,
     GenSA,
     parallel,
+    methods,
     stats,
-    graphics
+    graphics,
+    r2r,
+    BAS,
+    tolerance
 Imports:
     Rcpp
 LinkingTo:
@@ -24,6 +31,7 @@ Suggests:
     knitr,
     rmarkdown,
     markdown
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)
 VignetteBuilder: knitr
+LazyData: true

Make_data/Read_sanger_data.R

@@ -0,0 +1,51 @@
+#Read Sanger data
+#Download first from https://ftp.sanger.ac.uk/pub/genevar/ the files
+ #https://ftp.sanger.ac.uk/pub/genevar/CEU_children_norm_march2007.zip
+ #https://ftp.sanger.ac.uk/pub/genevar/CHB_unrelated_norm_march2007.zip
+ #https://ftp.sanger.ac.uk/pub/genevar/JPT_unrelated_norm_march2007.zip
+ #https://ftp.sanger.ac.uk/pub/genevar/YRI_parents_norm_march2007.zip
+#and unzip these
+
+#Specify right path
+path = "/mn/sarpanitu/ansatte-u2/geirs/prj/FBMS/data/"
+x1 = t(read.table(paste(path,"CEU_parents_norm_march2007.txt",sep=""),header=T))
+x2 = t(read.table(paste(path,"CHB_unrelated_norm_march2007.txt",sep=""),header=T))
+x3 = t(read.table(paste(path,"JPT_unrelated_norm_march2007.txt",sep=""),header=T))
+x4 = t(read.table(paste(path,"YRI_parents_norm_march2007.txt",sep=""),header=T))
+
+rnames = c(rownames(x1)[-1],rownames(x2)[-1],rownames(x3)[-1],rownames(x4)[-1])
+nam = x1[1,]
+x1 = apply(x1[-1,],2,as.numeric)
+x2 = apply(x2[-1,],2,as.numeric)
+x3 = apply(x3[-1,],2,as.numeric)
+x4 = apply(x4[-1,],2,as.numeric)
+df = rbind(x1,x2,x3,x4)
+colnames(df) = nam
+rownames(df) = rnames
+
+#Make columnn 24266 the first column, corresponding to CCT8 
+#(from the 	illumina_Human_WG-6_array_content.csv file)
+df <- df[,c(24266,1:24265,24267:ncol(df))]
+
+#Choose response variable
+SangerData = df
+usethis::use_data(SangerData,overwrite=TRUE)
+
+#Rename columns
+#colnames(df) = c("y",paste0("x",1:47292))
+
+#Reduced dataset first by those having maximum expression levels below the 
+#25-th percentile of all measured expression levels
+q = quantile(df[,-1],0.25)
+foo = apply(df[,-1],2,max)
+nC = ncol(df)-1
+df = df[,c(1,1+c(1:nC)[foo>q])]
+
+#Reduced dataset by deleting rows with range<2
+drange = function(x){diff(range(x))}
+foo2 = apply(df,2,drange)
+nC = ncol(df)-1
+df = df[,c(1,1+c(1:nC)[foo2>2])]
+SangerData2 = as.data.frame(df)
+usethis::use_data(SangerData2,overwrite=TRUE)
+

NAMESPACE

@@ -4,6 +4,7 @@ S3method(plot,gmjmcmc)
 S3method(plot,gmjmcmc_merged)
 S3method(plot,mjmcmc)
 S3method(plot,mjmcmc_parallel)
+S3method(predict,bgnlm_model)
 S3method(predict,gmjmcmc)
 S3method(predict,gmjmcmc_merged)
 S3method(predict,gmjmcmc_parallel)
@@ -14,24 +15,34 @@ S3method(summary,gmjmcmc)
 S3method(summary,gmjmcmc_merged)
 S3method(summary,mjmcmc)
 S3method(summary,mjmcmc_parallel)
-export(cos.rad)
+export(compute_effects)
+export(cos_deg)
+export(diagn_plot)
 export(erf)
-export(exp.dbl)
+export(exp_dbl)
+export(fbms)
+export(fbms.mlik.master)
 export(gauss)
 export(gaussian.loglik)
 export(gaussian.loglik.alpha)
-export(gaussian.loglik.bic.irlssgd)
+export(gaussian.loglik.g)
+export(gaussian_tcch_log_likelihood)
 export(gelu)
 export(gen.params.gmjmcmc)
 export(gen.params.mjmcmc)
 export(gen.probs.gmjmcmc)
 export(gen.probs.mjmcmc)
+export(get.best.model)
+export(get.mpm.model)
+export(glm.loglik)
+export(glm.logpost.bas)
 export(gmjmcmc)
-export(gmjmcmc.iact)
 export(gmjmcmc.parallel)
 export(hs)
-export(linear.g.prior.loglik)
+export(lm.logpost.bas)
+export(log_prior)
 export(logistic.loglik)
+export(logistic.loglik.ala)
 export(logistic.loglik.alpha)
 export(marginal.probs)
 export(merge_results)
@@ -60,27 +71,67 @@ export(pm2)
 export(relu)
 export(set.transforms)
 export(sigmoid)
-export(sin.rad)
+export(sin_deg)
 export(sqroot)
 export(string.population)
 export(string.population.models)
+export(to23)
 export(to25)
 export(to35)
+export(to72)
 export(troot)
+importFrom(BAS,CCH)
+importFrom(BAS,EB.local)
+importFrom(BAS,Jeffreys)
+importFrom(BAS,TG)
+importFrom(BAS,beta.prime)
+importFrom(BAS,bic.prior)
+importFrom(BAS,g.prior)
+importFrom(BAS,hyper.g)
+importFrom(BAS,hyper.g.n)
+importFrom(BAS,hypergeometric1F1)
+importFrom(BAS,hypergeometric2F1)
+importFrom(BAS,intrinsic)
+importFrom(BAS,phi1)
+importFrom(BAS,tCCH)
+importFrom(BAS,uniform)
 importFrom(GenSA,GenSA)
 importFrom(Rcpp,sourceCpp)
 importFrom(fastglm,fastglm)
 importFrom(graphics,barplot)
+importFrom(graphics,lines)
 importFrom(graphics,text)
+importFrom(methods,is)
+importFrom(parallel,clusterEvalQ)
+importFrom(parallel,clusterExport)
+importFrom(parallel,detectCores)
+importFrom(parallel,makeCluster)
 importFrom(parallel,mclapply)
+importFrom(parallel,parLapply)
+importFrom(parallel,stopCluster)
+importFrom(r2r,has_key)
+importFrom(r2r,hashmap)
+importFrom(stats,Gamma)
 importFrom(stats,acf)
+importFrom(stats,as.formula)
 importFrom(stats,binomial)
 importFrom(stats,cor)
+importFrom(stats,gaussian)
+importFrom(stats,glm.control)
 importFrom(stats,lm)
 importFrom(stats,lm.fit)
+importFrom(stats,median)
+importFrom(stats,model.matrix)
+importFrom(stats,model.response)
 importFrom(stats,pnorm)
+importFrom(stats,poisson)
+importFrom(stats,predict)
+importFrom(stats,qnorm)
 importFrom(stats,rbinom)
 importFrom(stats,rnorm)
 importFrom(stats,runif)
+importFrom(stats,sd)
 importFrom(stats,var)
-useDynLib(GMJMCMC)
+importFrom(tolerance,F1)
+importFrom(utils,sessionInfo)
+useDynLib(FBMS)

R/RcppExports.R

@@ -2,10 +2,10 @@
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 set_alphas <- function(formula) {
-    .Call('_GMJMCMC_set_alphas', PACKAGE = 'GMJMCMC', formula)
+    .Call('_FBMS_set_alphas', PACKAGE = 'FBMS', formula)
 }
 
 vec_in_mat <- function(mat, vec, firstCol = 0L, lastCol = 0L) {
-    .Call('_GMJMCMC_vec_in_mat', PACKAGE = 'GMJMCMC', mat, vec, firstCol, lastCol)
+    .Call('_FBMS_vec_in_mat', PACKAGE = 'FBMS', mat, vec, firstCol, lastCol)
 }
 

R/abalone-data.R

@@ -0,0 +1,31 @@
+#' Physical measurements of 4177 abalones, a species of sea snail.
+#'
+#' %% ~~ A concise (1-5 lines) description of the dataset. ~~
+#'
+#' See the web page \url{https://archive.ics.uci.edu/ml/datasets/Abalone} for
+#' more information about the data set.
+#'
+#' @name abalone
+#' @docType data
+#' @format A data frame with 4177 observations on the following 9 variables.
+#' \describe{
+#' \item{Diameter}{Diameter Perpendicular to length, continuous}
+#' \item{Height}{Height with with meat in shell, continuous.}
+#' \item{Length}{Longest shell measurement, continuous}
+#' \item{Rings}{+1.5 gives the age in years, integer}
+#' \item{Sex}{Sex of the abalone, \code{F} is female, \code{M} male, and \code{I} infant, categorical.}
+#' \item{Weight_S}{Grams after being dried, continuous.}
+#' \item{Weight_Sh}{Grams weight of meat, continuous.}
+#' \item{Weight_V}{Grams gut weight (after bleeding), continuous.}
+#' \item{Weight_W}{Grams whole abalone, continuous.} }
+#' @source Dua, D. and Graff, C. (2019). UCI Machine Learning Repository
+#' \url{https://archive.ics.uci.edu/ml/}. Irvine, CA: University of California,
+#' School of Information and Computer Science.
+#'
+#' @keywords datasets
+#' @examples
+#'
+#' data(abalone)
+#' ## maybe str(abalone) ; plot(abalone) ...
+#'
+NULL

R/alpha_generation.R

@@ -3,24 +3,29 @@
 # Created by: jonlachmann
 # Created on: 2021-03-16
 
-gen.alphas <- function (strategy, feature, data, loglik) {
-  if (strategy == 1) stop("Not implemented.")
-  else if (strategy == 2) stop("Not implemented.")
-  else if (strategy == 3) feature <- alpha_3(feature, data, loglik)
+gen.alphas <- function (strategy, feature, data, loglik, verbose) {
+  if (strategy == "deep") feature <- alpha_3(feature, data, loglik, verbose)
+  else if (strategy == "random") feature <- alpha_4(feature) else stop("Not implemented.")
   return(feature)
 }
 
 #' Alpha generator using strategy 1 as per Hubin et. al.
 #' TODO: This is just a placeholder.
+#'
+#' @noRd
+#'
 #' @param feature The feature to generate alphas for
 alpha_1 <- function (feature) {
   return(feature)
 }
 
 #' Alpha generator using strategy 2 as per Hubin et. al.
 #' TODO: This is just a placeholder.
+#'
+#' @noRd
+#'
 #' @param feature The feature to generate alphas for
 alpha_2 <- function (feature) {
   return(feature)
 }
 
@@ -29,31 +34,48 @@
 #' @param feature The feature to generate alphas for
 #' @param data The dataset used
 #' @param loglik log likelihood function to use
-alpha_3 <- function (feature, data, loglik) {
+#'
+#' @noRd
+#'
+alpha_3 <- function (feature, data, loglik, verbose) {
   # Create the string representation of the feature with variable alphas
   featfun <- print.feature(feature, dataset = TRUE, alphas = TRUE)
   featfun <- set_alphas(featfun)
   # Return if there are no alphas to set
   if (featfun$count == 0) return(feature)
 
   # Set initial range for Simulated Annealing
-  cat("Generating alphas\n")
   range <- 10
   done <- FALSE
   while (!done) {
     # Run simulated annealing on current range
     sares <- GenSA::GenSA(rnorm(featfun$count), loglik,
                       rep(-range / 2, featfun$count), rep(range / 2, featfun$count),
-                      control = list(max.call = 5e3), data, featfun$formula)
+                      control = list(max.call = 5e3, maxit = 20,nb.stop.improvement = 5), data, featfun$formula)
     # Check if any estimate is on the edge of the range, if so, extend the range and run again
     if (sum((sares$par == (-range / 2)) + (sares$par == (range / 2))) != 0) range <- range*2
     else done <- TRUE
   }
   if (sum(sares$par == 0) == featfun$count) {
-    cat("All zero feature occured.\n")
+    if (verbose) cat("All zero feature occured.\n")
     return(NULL)
   }
   # Inject the new alphas into the feature
   feature <- update.alphas(feature, sares$par)
   return(feature)
 }
+
+#' Alpha generator using strategy 4 as per Hubin et. al.
+#'
+#' @param feature The feature to generate alphas for
+#'
+#' @noRd
+#'
+alpha_4 <- function (feature) {
+  # Create the string representation of the feature with variable alphas
+  featfun <- print.feature(feature, dataset = TRUE, alphas = TRUE)
+  featfun <- set_alphas(featfun)
+
+  feature <- update.alphas(feature, rnorm(featfun$count,0,1))
+  return(feature)
+}