diff --git a/R/contrast_stats.R b/R/contrast_stats.R
new file mode 100644
index 0000000..823f67a
--- /dev/null
+++ b/R/contrast_stats.R
@@ -0,0 +1,27 @@
+
+#
+contrast_stats <- function(AC_stats,
+                           BC_stats,
+                           CI = 0.95) {
+  upper <- 0.5 + CI/2
+  ci_range <- c(1-upper, upper)
+  
+  contrasts <- list(
+    AB = AC_stats$mean - BC_stats$mean,
+    AC = AC_stats$mean,
+    BC = BC_stats$mean)
+  
+  contrast_variances <- list(
+    AB = AC_stats$var + BC_stats$var,
+    AC = AC_stats$var,
+    BC = BC_stats$var)
+  
+  contrast_ci <- list(
+    AB = contrasts$AB + qnorm(ci_range)*as.vector(sqrt(contrast_variances$AB)),
+    AC = contrasts$AC + qnorm(ci_range)*as.vector(sqrt(contrast_variances$AC)),
+    BC = contrasts$BC + qnorm(ci_range)*as.vector(sqrt(contrast_variances$BC)))
+  
+  list(contrasts = contrasts,
+       variances = contrast_variances,
+       CI = contrast_ci)
+}
diff --git a/R/outstandR.R b/R/outstandR.R
index cdff129..4c94d0e 100644
--- a/R/outstandR.R
+++ b/R/outstandR.R
@@ -50,50 +50,13 @@ outstandR <- function(AC.IPD, BC.ALD, strategy, CI = 0.95, ...) {
   if (!inherits(strategy, "strategy"))
     stop("strategy argument must be a class strategy.")
   
-  # select data according to formula
-  ipd <- model.frame(strategy$formula, data = AC.IPD)
-  
-  term.labels <- attr(terms(strategy$formula), "term.labels")
-  mean_names <- paste0("mean.", term.labels)
-  sd_names <- paste0("sd.", term.labels)
-  term_names <- c(mean_names, sd_names)
-  
-  # remove treatment labels
-  term_names <- sort(term_names[!grepl(pattern = "trt", term_names)])
-  
-  # replace outcome variable name
-  response_var <- all.vars(strategy$formula)[1]
-  response_names <- gsub(pattern = "y", replacement = response_var,
-                         x = c("y.B.sum", "y.B.bar", "N.B", "y.C.sum", "y.C.bar", "N.C")) 
-  
-  keep_names <- c(term_names, response_names)
-  
-  ald <- BC.ALD[keep_names]
+  ipd <- prep_ipd(strategy$formula, AC.IPD)
+  ald <- prep_ald(strategy$formula, BC.ALD)
 
-  AC_outstandR <- IPD_stats(strategy, ipd = ipd, ald = ald, ...) 
-  BC_outstandR <- ALD_stats(ald = ald) 
-  
-  upper <- 0.5 + CI/2
-  ci_range <- c(1-upper, upper)
-  
-  contrasts <- list(
-    AB = AC_outstandR$mean - BC_outstandR$mean,
-    AC = AC_outstandR$mean,
-    BC = BC_outstandR$mean)
-  
-  contrast_variances <- list(
-    AB = AC_outstandR$var + BC_outstandR$var,
-    AC = AC_outstandR$var,
-    BC = BC_outstandR$var)
-  
-  contrast_ci <- list(
-    AB = contrasts$AB + qnorm(ci_range)*as.vector(sqrt(contrast_variances$AB)),
-    AC = contrasts$AC + qnorm(ci_range)*as.vector(sqrt(contrast_variances$AC)),
-    BC = contrasts$BC + qnorm(ci_range)*as.vector(sqrt(contrast_variances$BC)))
+  AC_stats <- IPD_stats(strategy, ipd = ipd, ald = ald, ...) 
+  BC_stats <- ALD_stats(ald = ald) 
   
-  stats <- list(contrasts = contrasts,
-                variances = contrast_variances,
-                CI = contrast_ci)
+  stats <- contrast_stats(AC_stats, BC_stats, CI)
   
   structure(stats,
             CI = CI,
diff --git a/R/prep_data.R b/R/prep_data.R
new file mode 100644
index 0000000..9f11e53
--- /dev/null
+++ b/R/prep_data.R
@@ -0,0 +1,28 @@
+# prepare data functions
+
+#
+prep_ipd <- function(form, data) {
+  # select data according to formula
+  model.frame(form, data = data)
+}
+
+#
+prep_ald <- function(form, data) {
+  
+  term.labels <- attr(terms(form), "term.labels")
+  mean_names <- paste0("mean.", term.labels)
+  sd_names <- paste0("sd.", term.labels)
+  term_names <- c(mean_names, sd_names)
+  
+  # remove treatment labels
+  term_names <- sort(term_names[!grepl(pattern = "trt", term_names)])
+  
+  # replace outcome variable name
+  response_var <- all.vars(form)[1]
+  response_names <- gsub(pattern = "y", replacement = response_var,
+                         x = c("y.B.sum", "y.B.bar", "N.B", "y.C.sum", "y.C.bar", "N.C")) 
+  
+  keep_names <- c(term_names, response_names)
+  
+  data[keep_names]
+}