Merge pull request #157 from DeclareDesign/mh_next_round

add sd

R/multi_arm_designer.R

@@ -4,8 +4,9 @@
 #'
 #' @param N An integer. Sample size.
 #' @param m_arms An integer. Number of arms.
-#' @param means A numeric vector of length \code{m_arms}.  Average outcome in each arm.
-#' @param sds A nonnegative numeric vector of length \code{m_arms}. Standard deviations for each of the arms.
+#' @param outcome_means A numeric vector of length \code{m_arms}.  Average outcome in each arm.
+#' @param sd A nonnegative scalar. Standard deviations for shock for each unit (common across arms).
+#' @param outcome_sds A nonnegative numeric vector of length \code{m_arms}. Standard deviations for additional shock for each unit for each of the arms.
 #' @param conditions A vector of length \code{m_arms}. The names of each arm. It can be numeric or a character without blank spaces. 
 #' @param fixed A character vector. Names of arguments to be fixed in design. By default \code{m_arms} and \code{conditions} are always fixed.
 #' @return A function that returns a design.
@@ -21,48 +22,50 @@
 #'
 #'
 #' # A design with different mean and sd in each arm
-#' design <- multi_arm_designer(means = c(0, 0.5, 2), sd =  c(1, 0.1, 0.5))
+#' design <- multi_arm_designer(outcome_means = c(0, 0.5, 2), sd =  c(1, 0.1, 0.5))
 #'
 # A design with fixed sds and means. N is the sole modifiable argument.
-#' design <- multi_arm_designer(N = 80, m_arms = 4, means = 1:4,
-#'                              fixed = c("means", "sds"))
+#' design <- multi_arm_designer(N = 80, m_arms = 4, outcome_means = 1:4,
+#'                              fixed = c("outcome_means", "outcome_sds"))
 #'
 
 multi_arm_designer <- function(N = 30,
                                m_arms = 3,
-                               means = rep(0, m_arms),
-                               sds = rep(1, m_arms),
+                               outcome_means = rep(0, m_arms),
+                               sd = 1,
+                               outcome_sds = rep(0, m_arms),
                                conditions = 1:m_arms,
                                fixed = NULL) {
   # Housekeeping
   Y_Z_1 <- Z <- NULL
-  sds_ <- sds 
-  means_ <- means
+  outcome_sds_ <- outcome_sds 
+  outcome_means_ <- outcome_means
   N_ <- N
   if (m_arms <= 1 || round(m_arms) != m_arms)
     stop("m_arms should be an integer greater than one.")
-  if (length(means) != m_arms ||
-      length(sds) != m_arms ||
+  if (length(outcome_means) != m_arms ||
+      length(outcome_sds) != m_arms ||
       length(conditions) != m_arms)
-    stop("means, sds and conditions arguments must be of length m_arms.")
-  if (any(sds <= 0)) stop("sds should be nonnegative")
-  if (!"sds" %in% fixed) sds_ <-  sapply(1:m_arms, function(i) expr(sds[!!i]))
-  if (!"means" %in% fixed) means_ <-  sapply(1:m_arms, function(i) expr(means[!!i]))
+    stop("outcome_means, outcome_sds and conditions arguments must be of length m_arms.")
+  if (sd < 0) stop("sd should be nonnegative")
+  if (any(outcome_sds < 0)) stop("outcome_sds should be nonnegative")
+  if (!"outcome_sds" %in% fixed) outcome_sds_ <-  sapply(1:m_arms, function(i) expr(outcome_sds[!!i]))
+  if (!"outcome_means" %in% fixed) outcome_means_ <-  sapply(1:m_arms, function(i) expr(outcome_means[!!i]))
   if (!"N" %in% fixed) N_ <- expr(N)
 
   # Create helper vars to be used in design
-  errors <- sapply(1:m_arms, function(x) quos(rnorm(!!N_, 0, !!!sds_[x])))
+  errors <- sapply(1:m_arms, function(x) quos(rnorm(!!N_, 0, !!!outcome_sds_[x])))
   error_names <- paste0("u_", 1:m_arms)
   names(errors) <- error_names
-  population_expr <- expr(declare_population(N = !!N_, !!!errors))
+  population_expr <- expr(declare_population(N = !!N_, !!!errors, u = rnorm(!!N_)*sd))
 
   conditions <- as.character(conditions)
 
   f_Y <- formula(
     paste0("Y ~ ",paste0(
-      "(", means_, " + ", error_names,
+      "(", outcome_means_, " + ", error_names,
       ")*( Z == '", conditions, "')",
-      collapse = " + ")))
+      collapse = " + "), "+ u"))
 
   potential_outcomes_expr <-
     expr(

R/simple_factorial_designer.R

@@ -22,7 +22,8 @@
 #' @param mean_A0B1 A number. Mean outcome in A=0, B=1 condition.
 #' @param mean_A1B0 A number. Mean outcome in A=1, B=0 condition.
 #' @param mean_A1B1 A number. Mean outcome in A=1, B=1 condition.
-#' @param outcome_sds A non-negative 4-vector.  Standard deviation in each condition, in order AB = 00, 01, 10, 11.
+#' @param sd A nonnegative scalar. Standard deviations for shock for each unit (common across arms).
+#' @param outcome_sds A non-negative 4-vector.Standard deviation of (additional) unit level shock in each condition, in order AB = 00, 01, 10, 11.
 #' @return A two-by-two factorial design.
 #' @author \href{https://declaredesign.org/}{DeclareDesign Team}
 #' @concept experiment factorial
@@ -63,23 +64,24 @@ simple_factorial_designer <- function(N = 100,
                                       mean_A0B1 = outcome_means[2],
                                       mean_A1B0 = outcome_means[3],
                                       mean_A1B1 = outcome_means[4],
-                                      outcome_sds = rep(1,4)
+                                      sd = 1,
+                                      outcome_sds = rep(0,4)
 ){
   Y_A_0_B_0 <- Y_A_0_B_1 <- Y_A_1_B_0 <- Y_A_1_B_1 <- A <- B <- Y <- NULL
   if((w_A < 0) || (w_B < 0) || (w_A > 1) || (w_B > 1)) stop("w_A and w_B must be in 0,1")
-  if(max(outcome_sds < 0) )      stop("sd must be non-negative")
+  if(max(c(sd, outcome_sds) < 0) )      stop("sd must be non-negative")
   if(max(c(prob_A, prob_B) < 0)) stop("prob_ arguments must be non-negative")
   if(max(c(prob_A, prob_B) > 1))  stop("prob_ arguments must not exceed 1")
   {{{
 
     # M: Model
-    population <- declare_population(N)
+    population <- declare_population(N, u = rnorm(N, sd=sd))
 
     potentials <- declare_potential_outcomes(
-      Y_A_0_B_0 = mean_A0B0 + rnorm(N, sd = outcome_sds[1]),  
-      Y_A_0_B_1 = mean_A0B1 + rnorm(N, sd = outcome_sds[2]),  
-      Y_A_1_B_0 = mean_A1B0 + rnorm(N, sd = outcome_sds[3]),
-      Y_A_1_B_1 = mean_A1B1 + rnorm(N, sd = outcome_sds[4]))
+      Y_A_0_B_0 = mean_A0B0 + u + rnorm(N, sd = outcome_sds[1]),  
+      Y_A_0_B_1 = mean_A0B1 + u + rnorm(N, sd = outcome_sds[2]),  
+      Y_A_1_B_0 = mean_A1B0 + u + rnorm(N, sd = outcome_sds[3]),
+      Y_A_1_B_1 = mean_A1B1 + u + rnorm(N, sd = outcome_sds[4]))
 
 
     # I: Inquiry
@@ -97,7 +99,7 @@ simple_factorial_designer <- function(N = 100,
     # Factorial assignments
     assign_A <- declare_assignment(prob = prob_A, assignment_variable = A)
     assign_B <- declare_assignment(prob = prob_B, assignment_variable = B, blocks = A)
-    reveal_Y   <- declare_reveal(Y_variables = Y, assignment_variables = c(A,B))
+    reveal_Y <- declare_reveal(Y_variables = Y, assignment_variables = c(A,B))
 
     # A: Answer Strategy
     estimator_1 <- declare_estimator(Y ~ A + B,

tests/testthat/test_designers.R

@@ -163,9 +163,10 @@ test_that(desc = "cluster_sampling_designer errors when it should",
 
 test_that(desc = "multi_arm_designer errors when it should",
           code = {
-            expect_error(multi_arm_designer(means = rep(1,2),m_arms = 10))
-            expect_error(multi_arm_designer(m_arms = .5,means = 2))
-            expect_error(multi_arm_designer(sds = c(-10,-10),means = c(2,2), m_arms = 2))
+            expect_error(multi_arm_designer(outcome_means = rep(1,2), m_arms = 10))
+            expect_error(multi_arm_designer(m_arms = .5,outcome_means = 2))
+            expect_error(multi_arm_designer(outcome_sds = c(-10,-10),outcome_means = c(2,2), m_arms = 2))
+            expect_error(multi_arm_designer(sd = -1))
           })
 
 

vignettes/multi_arm.Rmd

@@ -32,7 +32,7 @@ In settings of multiple treatment arms, we could do a number of pairwise compari
 - **A**nswer strategy: We fit a linear regression model with individual indicators for each of the treatments as covariates. The average treatment effect is equal to the regression coefficients, which are computed by subtracting the mean of the comparison group from the mean of each treatment group.
 
 
-```{r, code = get_design_code(multi_arm_designer( means = c(0.5, 1, 2))), eval=TRUE}
+```{r, code = get_design_code(multi_arm_designer(outcome_means = c(0.5, 1, 2))), eval=TRUE}
 ```
 
 ```{r}
@@ -53,10 +53,10 @@ In R, you can generate a multi-arm design using the template function `multi_arm
 library(DesignLibrary)
 ```
 
-We can then create specific designs by defining values for each argument. For example, we can create a design called `my_multi_arm_design` where `N`, `m_arms`, and `means` set to 80, 4, and `c(-0.2, 0.2, 0.1, 0)`, respectively, and other parameters use default values. To do so, we run the lines below.
+We can then create specific designs by defining values for each argument. For example, we can create a design called `my_multi_arm_design` where `N`, `m_arms`, and `outcome_means` set to 80, 4, and `c(-0.2, 0.2, 0.1)`, respectively, and other parameters use default values. To do so, we run the lines below.
 
 ```{r, eval=FALSE}
-my_multi_arm_design <- multi_arm_designer(N = 80, m_arms = 4, means = c(-0.2, 0.2, 0.1, 0))
+my_multi_arm_design <- multi_arm_designer(N = 80, m_arms = 4, outcome_means = c(-0.2, 0.2, 0.1))
 ```
 
 You can see more details on the `multi_arm_designer()` function, its arguments, and default values, by running the following line of code: