version 1.1.0

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: optrefine
 Title: Optimally Refine Strata
-Version: 1.0.0
+Version: 1.1.0
 Authors@R: 
     person("Katherine", "Brumberg", , "kbrum@wharton.upenn.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-5193-6250"))
@@ -12,13 +12,13 @@ BugReports: https://github.com/kkbrum/optrefine/issues
 License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Depends: R (>= 2.10), MASS, Rglpk, sampling, ggplot2
 Suggests: covr, gurobi, testthat (>= 3.0.0)
 Config/testthat/edition: 3
 NeedsCompilation: no
-Packaged: 2022-11-07 21:31:34 UTC; katherine
+Packaged: 2023-04-18 16:12:11 UTC; katherine
 Author: Katherine Brumberg [aut, cre] (<https://orcid.org/0000-0002-5193-6250>)
 Maintainer: Katherine Brumberg <kbrum@wharton.upenn.edu>
 Repository: CRAN
-Date/Publication: 2022-11-08 14:20:08 UTC
+Date/Publication: 2023-04-18 19:20:09 UTC

MD5

@@ -1,17 +1,17 @@
-1008c86f87f1cec8489d3621ffb84a0b *DESCRIPTION
+5f01bd23ab6c637bb4c4df7272effddb *DESCRIPTION
 8fb9529835f8f530d0475b23c06fde16 *NAMESPACE
-745d8da0e01e6eb832f37a6401ffa41c *NEWS.md
+e475b40f61dfe76e7404c67cf2025553 *NEWS.md
 d36c1f69a84e2dc2b8396dceb074a056 *R/best_split.R
 ac87b09049dcc3a95849952f67492d50 *R/calc_smds.R
 e435e6b633e5dc5bb8822d7927b9b518 *R/data.R
 916a3b2236caf42e29db5981958e2d5f *R/plot.strat.R
 0897c5867eb877de28dc55d73ea883a2 *R/print.strat.R
 a402effda4b4550380e95d57a0e6f10f *R/prop_strat.R
 f187a75e17d28b2417852c02a0c9e198 *R/rand_pvals.R
-d54e57cbc65f12c2df1c8a32c4568c0f *R/refine.R
-033c3d830a9ba8467a215eff540e898d *R/split_stratum.R
+26198bc38be3c05c8bc76c126ae3efe8 *R/refine.R
+ddb8ab1279a199b9dbc45765ae4839cf *R/split_stratum.R
 5deb71c328e61bab62bf2c6268ece6b9 *R/strat.R
-4a8c621c5d1869b340fcc229c6157ee1 *README.md
+58fb912aabf8d94bbd96a0e1952972e8 *README.md
 c055510b95ea37125a153dd198e94f10 *data/rhc_X.rda
 490380f001dd2b8f071631d722c121bb *man/best_split.Rd
 c6789311b1ce886a17750f3ed1b94c1a *man/calc_smds.Rd
@@ -20,7 +20,7 @@ bccacc1bedbb368d34b2c3bb7730e19a *man/plot.strat.Rd
 86c2349e601409a289abbe144a38c0d9 *man/print.strat.Rd
 20b714e0ee6f267b55387dbcc7babb3b *man/prop_strat.Rd
 4ae39dc84751eafdafe804c1920e221c *man/rand_pvals.Rd
-a42526b0bcf091e4e720fe08c0477e43 *man/refine.Rd
+dab86b6158bb8a3109be80bb38d9b6c6 *man/refine.Rd
 d06871afcd45d60be2be4b609000a12c *man/rhc_X.Rd
 6943b1cca35a0a32deeb832336a7e104 *man/split_stratum.Rd
 aeb5396c85e22aea95e68e48b08e0ce1 *man/strat.Rd
@@ -33,6 +33,6 @@ dc2e1ff0e3a24aefa3df6e1fcfdf8d78 *tests/testthat/test-best_split.R
 5c4836e1537ba854c33f06fc7d23947c *tests/testthat/test-print.R
 ca92e980837ce86f954bd3ec8ee308ef *tests/testthat/test-prop_strat.R
 209d6646b949f662d5ff9dc16ce26d20 *tests/testthat/test-rand_pvals.R
-2b58ce5caf2e99e031a834bac9541e26 *tests/testthat/test-refine.R
+fc63b5bd72280c6cff98624ed441cccf *tests/testthat/test-refine.R
 5df7e229f1e70f6631903c98183dbf3d *tests/testthat/test-split_stratum.R
 43ef1c7a8282079f99513045cb670436 *tests/testthat/test-strat.R

NEWS.md

@@ -1,3 +1,8 @@
 # optrefine 1.0.0
 
+* Added a couple more possible elements to the `options` argument to the `refine` function. These are `minsplit` for the minimum number of treated or control units to allow in a split stratum and `threads` for the number of threads to use in the optimization.
+* Changed the value of `details$valueLP` to `NA` if the integer program was solved instead of giving the same value as `details$valueIP`.
+
+# optrefine 1.0.0
+
 * Added a `NEWS.md` file to track changes to the package.

R/refine.R

@@ -37,6 +37,9 @@
 #' Only used if criterion is set to "combo". Default is 5}
 #'     \item{ist: }{which strata to split. Should be a level from the specified `strata` or
 #' a vector of multiple levels. Default is to split all strata}
+#'     \item{minsplit: }{The minimum number of treated and control units to allow in a refined stratum.
+#' Default is 10}
+#'     \item{threads: }{How many threads you'd like the optimization to use if using the "gurobi" solver. Uses all available threads by default}
 #' }
 #'
 #' Note that setting a seed before using this function will ensure that the results are reproducible
@@ -62,7 +65,6 @@
 #'         \item{wMax: }{weight placed on the maximum standardized mean difference in the optimization
 #'         (see the `details` about the `options` for the optimization)}
 #'         \item{X_std: }{standardized version of `X`}
-#'         \item{threads: }{how many threads you'd like the optimization to use if using the "gurobi" solver. Uses all available threads by default}
 #'         }
 #'     }
 #' }
@@ -93,7 +95,8 @@ refine <- function(object = NULL, z = NULL, X = NULL, strata = NULL,
                    options = list()) {
 
   stopifnot(is.list(options))
-  stopifnot(all(names(options) %in% c("solver", "standardize", "criterion", "integer", "wMax", "ist", "threads")))
+  stopifnot(all(names(options) %in% c("solver", "standardize", "criterion",
+                                      "integer", "wMax", "ist", "threads", "minsplit")))
   if (is.null(options$solver)) {
     if(requireNamespace("gurobi", quietly = TRUE)) {
       solver <- "gurobi"
@@ -109,12 +112,12 @@ refine <- function(object = NULL, z = NULL, X = NULL, strata = NULL,
   if (is.null(options$wMax)) {wMax <- 5} else {wMax <- options$wMax}
   if (is.null(options$ist)) {ist <- NULL} else {ist <- options$ist}
   if (is.null(options$threads)) {threads <- NULL} else {threads <- options$threads}
+  if (is.null(options$minsplit)) {min_split <- 10} else {min_split <- options$minsplit}
 
   stopifnot(solver %in% c("Rglpk", "gurobi"))
   stopifnot(is.logical(standardize))
   stopifnot(is.logical(integer))
   stopifnot(criterion %in% c("combo", "max", "sum"))
-  stopifnot(is.null(ist) || all(ist %in% levels(object$base_strata)))
   stopifnot(!is.null(object) || (!is.null(z) && !is.null(X)))
 
   if (is.null(object)) {
@@ -128,6 +131,7 @@ refine <- function(object = NULL, z = NULL, X = NULL, strata = NULL,
   } else {
     X <- object$X
   }
+  stopifnot(is.null(ist) || all(ist %in% levels(object$base_strata)))
 
   # Only split the strata indicated
   if (is.null(ist)) {
@@ -150,7 +154,6 @@ refine <- function(object = NULL, z = NULL, X = NULL, strata = NULL,
 
   s <- rep(NA,length(object$z))
 
-  min_split <- 10
   if (criterion == "sum") {
     wMax <- 0
     wEach <- 1

R/split_stratum.R

@@ -179,10 +179,10 @@ split_stratum <- function(z, X, strata, ist, nc, nt,
   rand_t_prop <- sapply(1:S, function(x) {(sum(st_mat[x, z == 1]) - sum(st_mat[x, z == 1] == 1)) / sum(st_mat[x, z == 1])})
 
   # Randomized rounding
-  n_rr <- 10
-  for (rr_i in 1:n_rr) {
-    st_mat <- round(pr, 5)
-    if(!integer) {
+  if(!integer) {
+    n_rr <- 10
+    for (rr_i in 1:n_rr) {
+      st_mat <- round(pr, 5)
       if (S > 2) {
         # (sample sizes only correct in expectation since we draw independently for each unit)
         for (i in 1:I) {
@@ -220,45 +220,51 @@ split_stratum <- function(z, X, strata, ist, nc, nt,
           st_mat[2, ] <- 1 - st_mat[1, ]
         }
       }
-    }
-    st <- which(st_mat == 1, arr.ind = TRUE)[, 1]
-    Sadj <- length(unique(st))
+      st <- which(st_mat == 1, arr.ind = TRUE)[, 1]
+      Sadj <- length(unique(st))
 
-    # Evaluate covariate balance of best integer solution
-    eval <- NULL
-    for (s in 1:S) {
-      if (sum(z == 1 & st == s) > 0) {
-        tmean <- colMeans(X[z == 1 & st == s, , drop = FALSE])
-      } else {
-        tmean <- rep(NA, J)
-      }
-      if (sum(z == 0 & st == s) > 0) {
-        cmean <- colMeans(X[z == 0 & st == s, , drop = FALSE])
-      } else {
-        cmean <- rep(NA, J)
+      # Evaluate covariate balance of best integer solution
+      eval <- NULL
+      for (s in 1:S) {
+        if (sum(z == 1 & st == s) > 0) {
+          tmean <- colMeans(X[z == 1 & st == s, , drop = FALSE])
+        } else {
+          tmean <- rep(NA, J)
+        }
+        if (sum(z == 0 & st == s) > 0) {
+          cmean <- colMeans(X[z == 0 & st == s, , drop = FALSE])
+        } else {
+          cmean <- rep(NA, J)
+        }
+        eval <- rbind(eval, tmean, cmean)
       }
-      eval <- rbind(eval, tmean, cmean)
-    }
-    rownames(eval)<-c(paste0(c("T", "C"), rep(1:S, each = 2)))
+      rownames(eval)<-c(paste0(c("T", "C"), rep(1:S, each = 2)))
 
-    epsIP <- matrix(NA,2*S+1,J)
-    for (s in 1:S) {
-      epsIP[(2*(s-1)+1), ] <- pmax(0, eval[2*s, ] - eval[2*s-1, ])
-      epsIP[2*s, ] <- pmax(0, eval[2*s-1, ] - eval[2*s, ])
-    }
-    epsIP[(2*S + 1),]<-apply(epsIP[1:(2*S),, drop = FALSE], 2, max, na.rm = TRUE)
-    rownames(epsIP)<-c(paste(c("Pos eps","Neg eps"), rep(1:S, each = 2)), "max")
-    if (!is.null(colnames(X))) colnames(epsIP)<-colnames(X)
-    epsIP_nona <- epsIP
-    epsIP_nona[is.na(epsIP)] <- 0
+      epsIP <- matrix(NA,2*S+1,J)
+      for (s in 1:S) {
+        epsIP[(2*(s-1)+1), ] <- pmax(0, eval[2*s, ] - eval[2*s-1, ])
+        epsIP[2*s, ] <- pmax(0, eval[2*s-1, ] - eval[2*s, ])
+      }
+      epsIP[(2*S + 1),]<-apply(epsIP[1:(2*S),, drop = FALSE], 2, max, na.rm = TRUE)
+      rownames(epsIP)<-c(paste(c("Pos eps","Neg eps"), rep(1:S, each = 2)), "max")
+      if (!is.null(colnames(X))) colnames(epsIP)<-colnames(X)
+      epsIP_nona <- epsIP
+      epsIP_nona[is.na(epsIP)] <- 0
 
-    valueIP <- (wEach * sum(epsIP_nona[1:(2*S),]) + wMax * max(epsIP_nona[2*S + 1, ], na.rm = TRUE)) / (wEach * J * Sadj + wMax)
-    if (valueIP < best_valueIP) {
-      best_valueIP <- valueIP
-      best_selection <- st
+      valueIP <- (wEach * sum(epsIP_nona[1:(2*S),]) + wMax * max(epsIP_nona[2*S + 1, ], na.rm = TRUE)) / (wEach * J * Sadj + wMax)
+      if (valueIP < best_valueIP) {
+        best_valueIP <- valueIP
+        best_selection <- st
+      }
     }
+  } else {
+    best_valueIP <- v
+    v <- NA
+    st <- which(st_mat == 1, arr.ind = TRUE)[, 1]
+    best_selection <- st
   }
 
+
   return(list(valueIP = best_valueIP, valueLP = v, n_smds = wEach * J * Sadj + wMax, n_fracs = n_fracs,
               rand_c_prop = rand_c_prop, rand_t_prop = rand_t_prop, pr = pr, selection = best_selection))
 }

README.md

@@ -16,6 +16,13 @@ covariate balance.
 
 ## Installation
 
+You can install the release version of optrefine from
+[CRAN](https://cran.r-project.org/) with:
+
+``` r
+install.packages("optrefine")
+```
+
 You can install the development version of optrefine from
 [GitHub](https://github.com/) with:
 

tests/testthat/test-refine.R

@@ -74,7 +74,7 @@ test_that("Arguments integer, wMax, and ist work as expected", {
                                wMax = 10,
                                ist = 2))
   expect_false(is.null(ref$details$valueLP))
-  expect_equal(ref$details$valueLP, ref$details$valueIP, tolerance = 10e-6)
+  expect_true(is.na(ref$details$valueLP))
   expect_equal(ref$details$valueIP, 0.249396391, tolerance = 10e-6)
 })