Updated R package to pass CRAN checks (had to manually remove some pragmas from json header)

Author: andrewherren

NAMESPACE

@@ -71,7 +71,7 @@ export(saveBARTModelToJsonString)
 export(saveBCFModelToJsonFile)
 export(saveBCFModelToJsonString)
 importFrom(R6,R6Class)
-importFrom(stats, coef)
+importFrom(stats,coef)
 importFrom(stats,lm)
 importFrom(stats,model.matrix)
 importFrom(stats,qgamma)

R/bart.R

@@ -575,7 +575,7 @@ bart <- function(X_train, y_train, W_train = NULL, group_ids_train = NULL,
                 }
                 if (has_rfx) {
                     resetRandomEffectsModel(rfx_model, rfx_samples, forest_ind, sigma_alpha_init)
-                    resetRandomEffectsTracker(rfx_tracker_train, rfx_model, rfx_dataset_train, outcome_train, rfx_samples, forest_ind)
+                    resetRandomEffectsTracker(rfx_tracker_train, rfx_model, rfx_dataset_train, outcome_train, rfx_samples)
                 }
                 if (sample_sigma_global) current_sigma2 <- global_var_samples[forest_ind + 1]
             } else if (has_prev_model) {

R/bcf.R

@@ -158,9 +158,11 @@
 #' tau_train <- tau_x[train_inds]
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, pi_train = pi_train, 
 #'                  X_test = X_test, Z_test = Z_test, pi_test = pi_test)
-#' # plot(rowMeans(bcf_model$mu_hat_test), mu_test, xlab = "predicted", ylab = "actual", main = "Prognostic function")
+#' # plot(rowMeans(bcf_model$mu_hat_test), mu_test, xlab = "predicted", 
+#' #      ylab = "actual", main = "Prognostic function")
 #' # abline(0,1,col="red",lty=3,lwd=3)
-#' # plot(rowMeans(bcf_model$tau_hat_test), tau_test, xlab = "predicted", ylab = "actual", main = "Treatment effect")
+#' # plot(rowMeans(bcf_model$tau_hat_test), tau_test, xlab = "predicted", 
+#' #      ylab = "actual", main = "Treatment effect")
 #' # abline(0,1,col="red",lty=3,lwd=3)
 bcf <- function(X_train, Z_train, y_train, pi_train = NULL, group_ids_train = NULL, 
                 rfx_basis_train = NULL, X_test = NULL, Z_test = NULL, pi_test = NULL, 
@@ -872,7 +874,7 @@ bcf <- function(X_train, Z_train, y_train, pi_train = NULL, group_ids_train = NU
                 }
                 if (has_rfx) {
                     resetRandomEffectsModel(rfx_model, rfx_samples, forest_ind, sigma_alpha_init)
-                    resetRandomEffectsTracker(rfx_tracker_train, rfx_model, rfx_dataset_train, outcome_train, rfx_samples, forest_ind)
+                    resetRandomEffectsTracker(rfx_tracker_train, rfx_model, rfx_dataset_train, outcome_train, rfx_samples)
                 }
                 if (adaptive_coding) {
                     current_b_1 <- b_1_samples[forest_ind + 1]
@@ -1190,6 +1192,8 @@ bcf <- function(X_train, Z_train, y_train, pi_train = NULL, group_ids_train = NU
         "num_gfr" = num_gfr, 
         "num_burnin" = num_burnin, 
         "num_mcmc" = num_mcmc, 
+        "keep_every" = keep_every,
+        "num_chains" = num_chains,
         "has_rfx" = has_rfx, 
         "has_rfx_basis" = has_basis_rfx, 
         "num_rfx_basis" = num_basis_rfx, 
@@ -1290,9 +1294,11 @@ bcf <- function(X_train, Z_train, y_train, pi_train = NULL, group_ids_train = NU
 #' tau_train <- tau_x[train_inds]
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, pi_train = pi_train)
 #' preds <- predict(bcf_model, X_test, Z_test, pi_test)
-#' # plot(rowMeans(preds$mu_hat), mu_test, xlab = "predicted", ylab = "actual", main = "Prognostic function")
+#' # plot(rowMeans(preds$mu_hat), mu_test, xlab = "predicted",
+#' #      ylab = "actual", main = "Prognostic function")
 #' # abline(0,1,col="red",lty=3,lwd=3)
-#' # plot(rowMeans(preds$tau_hat), tau_test, xlab = "predicted", ylab = "actual", main = "Treatment effect")
+#' # plot(rowMeans(preds$tau_hat), tau_test, xlab = "predicted", 
+#' #      ylab = "actual", main = "Treatment effect")
 #' # abline(0,1,col="red",lty=3,lwd=3)
 predict.bcf <- function(bcf, X_test, Z_test, pi_test = NULL, group_ids_test = NULL, rfx_basis_test = NULL){
     # Preprocess covariates
@@ -1475,13 +1481,14 @@ predict.bcf <- function(bcf, X_test, Z_test, pi_test = NULL, group_ids_test = NU
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' rfx_samples <- getRandomEffectSamples(bcf_model)
 getRandomEffectSamples.bcf <- function(object, ...){
     result = list()
@@ -1561,13 +1568,14 @@ getRandomEffectSamples.bcf <- function(object, ...){
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' # bcf_json <- convertBCFModelToJson(bcf_model)
 convertBCFModelToJson <- function(object){
     jsonobj <- createCppJson()
@@ -1617,6 +1625,8 @@ convertBCFModelToJson <- function(object){
     jsonobj$add_scalar("num_burnin", object$model_params$num_burnin)
     jsonobj$add_scalar("num_mcmc", object$model_params$num_mcmc)
     jsonobj$add_scalar("num_samples", object$model_params$num_samples)
+    jsonobj$add_scalar("keep_every", object$model_params$keep_every)
+    jsonobj$add_scalar("num_chains", object$model_params$num_chains)
     jsonobj$add_scalar("num_covariates", object$model_params$num_covariates)
     if (object$model_params$sample_sigma_global) {
         jsonobj$add_vector("sigma2_samples", object$sigma2_samples, "parameters")
@@ -1700,13 +1710,14 @@ convertBCFModelToJson <- function(object){
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' # saveBCFModelToJsonFile(bcf_model, "test.json")
 saveBCFModelToJsonFile <- function(object, filename){
     # Convert to Json
@@ -1773,13 +1784,14 @@ saveBCFModelToJsonFile <- function(object, filename){
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' # saveBCFModelToJsonString(bcf_model)
 saveBCFModelToJsonString <- function(object){
     # Convert to Json
@@ -1848,13 +1860,14 @@ saveBCFModelToJsonString <- function(object){
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' # bcf_json <- convertBCFModelToJson(bcf_model)
 #' # bcf_model_roundtrip <- createBCFModelFromJson(bcf_json)
 createBCFModelFromJson <- function(json_object){
@@ -1993,13 +2006,14 @@ createBCFModelFromJson <- function(json_object){
 #' rfx_basis_train <- rfx_basis[train_inds,]
 #' rfx_term_test <- rfx_term[test_inds]
 #' rfx_term_train <- rfx_term[train_inds]
+#' bcf_params <- list(sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
 #'                  Z_test = Z_test, pi_test = pi_test, group_ids_test = group_ids_test,
 #'                  rfx_basis_test = rfx_basis_test, 
 #'                  num_gfr = 100, num_burnin = 0, num_mcmc = 100, 
-#'                  sample_sigma_leaf_mu = TRUE, sample_sigma_leaf_tau = FALSE)
+#'                  params = bcf_params)
 #' # saveBCFModelToJsonFile(bcf_model, "test.json")
 #' # bcf_model_roundtrip <- createBCFModelFromJsonFile("test.json")
 createBCFModelFromJsonFile <- function(json_filename){
@@ -2100,24 +2114,55 @@ createBCFModelFromJsonString <- function(json_string){
 #' @examples
 #' n <- 100
 #' p <- 5
-#' X <- matrix(runif(n*p), ncol = p)
-#' f_XW <- (
-#'     ((0 <= X[,1]) & (0.25 > X[,1])) * (-7.5) + 
-#'     ((0.25 <= X[,1]) & (0.5 > X[,1])) * (-2.5) + 
-#'     ((0.5 <= X[,1]) & (0.75 > X[,1])) * (2.5) + 
-#'     ((0.75 <= X[,1]) & (1 > X[,1])) * (7.5)
-#' )
-#' noise_sd <- 1
-#' y <- f_XW + rnorm(n, 0, noise_sd)
+#' x1 <- rnorm(n)
+#' x2 <- rnorm(n)
+#' x3 <- rnorm(n)
+#' x4 <- rnorm(n)
+#' x5 <- rnorm(n)
+#' X <- cbind(x1,x2,x3,x4,x5)
+#' p <- ncol(X)
+#' g <- function(x) {ifelse(x[,5] < -0.44,2,ifelse(x[,5] < 0.44,-1,4))}
+#' mu1 <- function(x) {1+g(x)+x[,1]*x[,3]}
+#' mu2 <- function(x) {1+g(x)+6*abs(x[,3]-1)}
+#' tau1 <- function(x) {rep(3,nrow(x))}
+#' tau2 <- function(x) {1+2*x[,2]*(x[,4] > 0)}
+#' mu_x <- mu1(X)
+#' tau_x <- tau2(X)
+#' pi_x <- 0.8*pnorm((3*mu_x/sd(mu_x)) - 0.5*X[,1]) + 0.05 + runif(n)/10
+#' Z <- rbinom(n,1,pi_x)
+#' E_XZ <- mu_x + Z*tau_x
+#' snr <- 3
+#' group_ids <- rep(c(1,2), n %/% 2)
+#' rfx_coefs <- matrix(c(-1, -1, 1, 1), nrow=2, byrow=TRUE)
+#' rfx_basis <- cbind(1, runif(n, -1, 1))
+#' rfx_term <- rowSums(rfx_coefs[group_ids,] * rfx_basis)
+#' y <- E_XZ + rfx_term + rnorm(n, 0, 1)*(sd(E_XZ)/snr)
+#' X <- as.data.frame(X)
+#' X$x4 <- factor(X$x4, ordered = TRUE)
+#' X$x5 <- factor(X$x5, ordered = TRUE)
 #' test_set_pct <- 0.2
 #' n_test <- round(test_set_pct*n)
 #' n_train <- n - n_test
 #' test_inds <- sort(sample(1:n, n_test, replace = FALSE))
 #' train_inds <- (1:n)[!((1:n) %in% test_inds)]
 #' X_test <- X[test_inds,]
 #' X_train <- X[train_inds,]
+#' pi_test <- pi_x[test_inds]
+#' pi_train <- pi_x[train_inds]
+#' Z_test <- Z[test_inds]
+#' Z_train <- Z[train_inds]
 #' y_test <- y[test_inds]
 #' y_train <- y[train_inds]
+#' mu_test <- mu_x[test_inds]
+#' mu_train <- mu_x[train_inds]
+#' tau_test <- tau_x[test_inds]
+#' tau_train <- tau_x[train_inds]
+#' group_ids_test <- group_ids[test_inds]
+#' group_ids_train <- group_ids[train_inds]
+#' rfx_basis_test <- rfx_basis[test_inds,]
+#' rfx_basis_train <- rfx_basis[train_inds,]
+#' rfx_term_test <- rfx_term[test_inds]
+#' rfx_term_train <- rfx_term[train_inds]
 #' bcf_model <- bcf(X_train = X_train, Z_train = Z_train, y_train = y_train, 
 #'                  pi_train = pi_train, group_ids_train = group_ids_train, 
 #'                  rfx_basis_train = rfx_basis_train, X_test = X_test, 
@@ -2177,6 +2222,7 @@ createBCFModelFromCombinedJsonString <- function(json_string_list){
     model_params[["sample_sigma_leaf_mu"]] <- json_object_default$get_boolean("sample_sigma_leaf_mu")
     model_params[["sample_sigma_leaf_tau"]] <- json_object_default$get_boolean("sample_sigma_leaf_tau")
     model_params[["include_variance_forest"]] <- include_variance_forest
+    model_params[["propensity_covariate"]] <- json_object_default$get_string("propensity_covariate")
     model_params[["has_rfx"]] <- json_object_default$get_boolean("has_rfx")
     model_params[["has_rfx_basis"]] <- json_object_default$get_boolean("has_rfx_basis")
     model_params[["num_rfx_basis"]] <- json_object_default$get_scalar("num_rfx_basis")
@@ -2263,7 +2309,7 @@ createBCFModelFromCombinedJsonString <- function(json_string_list){
         output[["rfx_samples"]] <- loadRandomEffectSamplesCombinedJson(json_object_list, 0)
     }
 
-    class(output) <- "bartmodel"
+    class(output) <- "bcf"
     return(output)
 }
 

R/calibration.R

@@ -1,13 +1,13 @@
-#' Calibrate the scale parameter on an inverse gamma prior for the global error variance as in Chipman et al (2022) [1]
+#' Calibrate the scale parameter on an inverse gamma prior for the global error variance as in Chipman et al (2022)
 #' 
-#' [1] Chipman, H., George, E., Hahn, R., McCulloch, R., Pratola, M. and Sparapani, R. (2022). Bayesian Additive Regression Trees, Computational Approaches. In Wiley StatsRef: Statistics Reference Online (eds N. Balakrishnan, T. Colton, B. Everitt, W. Piegorsch, F. Ruggeri and J.L. Teugels). https://doi.org/10.1002/9781118445112.stat08288
+#' Chipman, H., George, E., Hahn, R., McCulloch, R., Pratola, M. and Sparapani, R. (2022). Bayesian Additive Regression Trees, Computational Approaches. In Wiley StatsRef: Statistics Reference Online (eds N. Balakrishnan, T. Colton, B. Everitt, W. Piegorsch, F. Ruggeri and J.L. Teugels). https://doi.org/10.1002/9781118445112.stat08288
 #'
 #' @param y Outcome to be modeled using BART, BCF or another nonparametric ensemble method.
 #' @param X Covariates to be used to partition trees in an ensemble or series of ensemble.
-#' @param W [Optional] Basis used to define a "leaf regression" model for each decision tree. The "classic" BART model assumes a constant leaf parameter, which is equivalent to a "leaf regression" on a basis of all ones, though it is not necessary to pass a vector of ones, here or to the BART function. Default: `NULL`.
+#' @param W (Optional) Basis used to define a "leaf regression" model for each decision tree. The "classic" BART model assumes a constant leaf parameter, which is equivalent to a "leaf regression" on a basis of all ones, though it is not necessary to pass a vector of ones, here or to the BART function. Default: `NULL`.
 #' @param nu The shape parameter for the global error variance's IG prior. The scale parameter in the Sparapani et al (2021) parameterization is defined as `nu*lambda` where `lambda` is the output of this function. Default: `3`.
-#' @param quant [Optional] Quantile of the inverse gamma prior distribution represented by a linear-regression-based overestimate of `sigma^2`. Default: `0.9`.
-#' @param standardize [Optional] Whether or not outcome should be standardized (`(y-mean(y))/sd(y)`) before calibration of `lambda`. Default: `TRUE`.
+#' @param quant (Optional) Quantile of the inverse gamma prior distribution represented by a linear-regression-based overestimate of `sigma^2`. Default: `0.9`.
+#' @param standardize (Optional) Whether or not outcome should be standardized (`(y-mean(y))/sd(y)`) before calibration of `lambda`. Default: `TRUE`.
 #'
 #' @return Value of `lambda` which determines the scale parameter of the global error variance prior (`sigma^2 ~ IG(nu,nu*lambda)`)
 #' @export 

R/forest.R

@@ -256,7 +256,7 @@ ForestSamples <- R6::R6Class(
         }, 
 
         #' @description
-        #' Add a numeric (i.e. X[,i] <= c) split to a given tree in the ensemble
+        #' Add a numeric (i.e. `X[,i] <= c`) split to a given tree in the ensemble
         #' @param forest_num Index of the forest which contains the tree to be split
         #' @param tree_num Index of the tree to be split
         #' @param leaf_num Leaf to be split
@@ -695,7 +695,7 @@ Forest <- R6::R6Class(
         }, 
 
         #' @description
-        #' Add a numeric (i.e. X[,i] <= c) split to a given tree in the ensemble
+        #' Add a numeric (i.e. `X[,i] <= c`) split to a given tree in the ensemble
         #' @param tree_num Index of the tree to be split
         #' @param leaf_num Leaf to be split
         #' @param feature_num Feature that defines the new split

R/model.R

@@ -153,6 +153,7 @@ createRNG <- function(random_seed = -1){
 #' @param alpha Root node split probability in tree prior
 #' @param beta Depth prior penalty in tree prior
 #' @param min_samples_leaf Minimum number of samples in a tree leaf
+#' @param max_depth Maximum depth of any tree in the ensemble in the mean model. Setting to ``-1`` does not enforce any depth limits on trees.
 #'
 #' @return `ForestModel` object
 #' @export