Merge pull request #5 from m-freitag/dev_mf

Dev mf

Author: m-freitag

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: cjpowR
 Type: Package
 Title: A Priori Power Analyses for Conjoint Experiments
-Version: 1.0.1
+Version: 1.0.2
 Author: Markus Freitag, Julian Schuessler
 Maintainer: Markus Freitag <[email protected]>
 Description: Implements the formulae provided by Schuessler/Freitag (2020) and provides some plotting functions.

R/amce.R

@@ -24,54 +24,53 @@
 #' @param levels the number of levels of the treatment.
 #' @param n effective sample size, default is NULL.
 #' @param treat.prob default is 0.5.
-#' @param sims number of simulation runs to compute the expected type M error ("exaggeration ratio"), default 100k. If NULL, E(type M) is not computed.
+#' @param sims number of simulation runs to compute the expected type M error ("exaggeration ratio"), default 100k. If NULL or, E(type M) is not computed (much faster).
 #' @keywords conjoint, power analysis, AMCE
 #' @export
 #' @export
 #' @export
 #' @export
 #' @examples
 #' # This gives the minimum required effective sample size (type S, E(type M)):
-#' df = cjpowr_amce(amce = 0.05, power = 0.8, levels = 5)
+#' df <- cjpowr_amce(amce = 0.05, power = 0.8, levels = 5)
 #'
 #' # For example, for a conjoint with 2 profiles and 4 tasks, n becomes:
 #'
-#' df$n/(2*4)
+#' df$n / (2 * 4)
 #'
-#' #This gives the power (type S, E(type M)):
+#' # This gives the power (type S, E(type M)):
 #' cjpowr_amce(amce = 0.05, n = 7829.258, levels = 5)
 #'
-#' #Generating an interactive plot for type M error:
+#' # Generating an interactive plot for type M error:
 #'
-#' d <- expand.grid(amce = c(0.01, 0.02, 0.03, 0.05), n = seq(from = 100, to = 50000, length.out = 1000))
+#' d <- expand.grid(
+#'     amce = c(0.01, 0.02, 0.03, 0.05), 
+#'     n = seq(from = 100, to = 50000, length.out = 1000), 
+#'     alpha = 0.05, 
+#'     levels = 2,
+#'     treat.prob = 0.5,
+#'     sims = 10000) #set to 0 if you want to make an interactive plot for something other than Type M error
 #'
-#' # Purrr Style:
-#' library(purrr)
-#' set.seed(123)
-#' df <- pmap_df(d, function(amce, n) cjpowr_amce(amce = amce, n = n, sims = 1000, levels = 5, alpha = 0.05, treat.prob = 0.5))
-#'
-#' # Base R:
-#' set.seed(123)
-#' cjpowr_amce_vec <- Vectorize(cjpowr_amce)
-#' df2 <- t(cjpowr_amce_vec(amce = d$amce, n = d$n, sims = 1000, levels = 5, alpha = 0.05, treat.prob = 0.5))
-#'
-#' df2 <- data.frame(df2)
-#' df2[] <- lapply(df2, unlist)
+#' df <- list2DF(do.call(cjpowr_amce, d))
 #'
 #' library(plotly)
-#' 
-#' plot_ly(df, x = ~n, y = ~exp_typeM, type = 'scatter', mode = 'lines', linetype = ~amce) %>%
-#'  layout(
-#'    xaxis = list(title = "Effective Sample Size",
-#'                 zeroline = F,
-#'                 hoverformat = '.0f'),
-#'    yaxis = list(title = "Exaggeration Ratio",
-#'                 range = c(0,10),
-#'                 zeroline = F,
-#'                 hoverformat = '.2f'),
-#'    legend=list(title=list(text='<b> AMCE </b>')),
-#'    hovermode = "x unified"
-#'  )
+#'
+#' plot_ly(df, x = ~n, y = ~exp_typeM, type = "scatter", mode = "lines", linetype = ~amce) %>%
+#'   layout(
+#'     xaxis = list(
+#'       title = "Effective Sample Size",
+#'       zeroline = F,
+#'       hoverformat = ".0f"
+#'     ),
+#'     yaxis = list(
+#'       title = "Exaggeration Ratio",
+#'       range = c(0, 10),
+#'       zeroline = F,
+#'       hoverformat = ".2f"
+#'     ),
+#'     legend = list(title = list(text = "<b> AMCE </b>")),
+#'     hovermode = "x unified"
+#'   )
 #' @section Literature:
 #'
 #' 1. Schuessler, J. and M. Freitag (2020). Power Analysis for Conjoint Experiments.
@@ -81,126 +80,121 @@
 #'    British Journal of Mathematical and Statistical Psychology 72(1), 1–17.
 
 cjpowr_amce <- function(amce, alpha = 0.05, power = NULL, levels = 2,
-                        treat.prob = 0.5, n = NULL, sims = 100000){
-
-  if (sum(sapply(list(power,n), is.null)) != 1)
+                        treat.prob = 0.5, n = NULL, sims = 100000) {
+  if (sum(sapply(list(power, n), is.null)) != 1) {
     stop("either 'n' or 'power' must be provided")
-  if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha > 1))
+  }
+  if (!is.null(alpha) && !is.numeric(alpha) || any(0 > alpha | alpha > 1)) {
     stop("'sig.level' must be numeric in [0, 1]")
+  }
 
-  if (is.null(sims)) {
+  lengths <- c(length(amce), length(alpha), length(power), length(levels), length(treat.prob), length(n), length(sims))
+  non_null_lengths <- lengths[lengths > 0]
 
-    if (is.null(n)) {
+  if (length(unique(non_null_lengths)) > 1) {
+    stop("All non-NULL input parameters must be of the same length.")
+  }
 
-      if (all(sapply(list(length(amce), length(alpha), length(power), length(levels)), function(x) x == 1)) == FALSE) {
-        stop("Please provide scalar values or apply a functional.")
-      }
+  length_output <- unique(non_null_lengths)
 
-      delta0 = 0.5 - (amce*treat.prob)
+  if (is.null(sims) | any(sims == 0) == TRUE) {
+    if (is.null(n)) {
+      delta0 <- 0.5 - (amce * treat.prob)
 
-      n = (levels/2)/(amce^2) * (qnorm(p = 1- alpha/2) + qnorm(p = power))^2*
+      n <- (levels / 2) / (amce^2) * (qnorm(p = 1 - alpha / 2) + qnorm(p = power))^2 *
         (
-          ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-            ((delta0)*(1 - delta0))/0.5
+          ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+            ((delta0) * (1 - delta0)) / 0.5
         )
 
-      se = (sqrt( (
-        ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-          ((delta0)*(1 - delta0))/0.5
-      ))/sqrt(2*(n/levels)))
+      se <- (sqrt((
+        ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+          ((delta0) * (1 - delta0)) / 0.5
+      )) / sqrt(2 * (n / levels)))
 
-      z=amce/se
+      z <- amce / se
 
-      pow=pnorm(z-qnorm(1-alpha/2)) + pnorm(-z-qnorm(1-alpha/2))
+      pow <- pnorm(z - qnorm(1 - alpha / 2)) + pnorm(-z - qnorm(1 - alpha / 2))
 
-      type_s = pnorm(-z-qnorm(1-alpha/2))/pow
+      type_s <- pnorm(-z - qnorm(1 - alpha / 2)) / pow
 
-      return(data.frame(amce = amce, n=n, type_s=type_s, power = power, alpha = alpha, levels = levels, delta0 = delta0))
+      return(data.frame(amce = amce, n = n, type_s = type_s, power = power, alpha = alpha, levels = levels, delta0 = delta0))
     }
 
     else if (is.null(power)) {
+      delta0 <- 0.5 - (amce * treat.prob)
 
-      if (all(sapply(list(length(amce), length(alpha), length(n), length(levels)), function(x) x == 1)) == FALSE) {
-        stop("Please provide scalar values or apply a functional.")
-      }
+      se <- (sqrt((
+        ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+          ((delta0) * (1 - delta0)) / 0.5
+      )) / sqrt(2 * (n / levels)))
 
-      delta0 = 0.5 - (amce * treat.prob)
+      z <- amce / se
 
-      se = (sqrt( (
-        ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-          ((delta0)*(1 - delta0))/0.5
-      ))/sqrt(2*(n/levels)))
+      power <- pnorm(z - qnorm(1 - alpha / 2)) + pnorm(-z - qnorm(1 - alpha / 2))
 
-      z=amce/se
+      type_s <- pnorm(-z - qnorm(1 - alpha / 2)) / power
 
-      power=pnorm(z-qnorm(1-alpha/2)) + pnorm(-z-qnorm(1-alpha/2))
-
-      type_s = pnorm(-z-qnorm(1-alpha/2))/power
-
-      return(data.frame(power=power, type_s=type_s, amce=amce, n=n, alpha = alpha, levels = levels, delta0 = delta0))
+      return(data.frame(power = power, type_s = type_s, amce = amce, n = n, alpha = alpha, levels = levels, delta0 = delta0))
     }
   }
   else if (is.null(n)) {
+    delta0 <- 0.5 - (amce * treat.prob)
 
-    if (all(sapply(list(length(amce), length(alpha), length(power), length(levels)), function(x) x == 1)) == FALSE) {
-      stop("Please provide scalar values or apply a functional.")
-    }
-
-    delta0 = 0.5 - (amce * treat.prob)
-
-    n = (levels/2)/(amce^2) * (qnorm(p = 1- alpha/2) + qnorm(p = power))^2*
+    n <- (levels / 2) / (amce^2) * (qnorm(p = 1 - alpha / 2) + qnorm(p = power))^2 *
       (
-        ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-          ((delta0)*(1 - delta0))/0.5
+        ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+          ((delta0) * (1 - delta0)) / 0.5
       )
 
-    se = (sqrt( (
-      ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-        ((delta0)*(1 - delta0))/0.5
-    ))/sqrt(2*(n/levels)))
+    se <- (sqrt((
+      ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+        ((delta0) * (1 - delta0)) / 0.5
+    )) / sqrt(2 * (n / levels)))
 
-    z=amce/se
+    z <- amce / se
 
-    pow=pnorm(z-qnorm(1-alpha/2)) + pnorm(-z-qnorm(1-alpha/2))
+    pow <- pnorm(z - qnorm(1 - alpha / 2)) + pnorm(-z - qnorm(1 - alpha / 2))
 
-    type_s = pnorm(-z-qnorm(1-alpha/2))/pow
+    type_s <- pnorm(-z - qnorm(1 - alpha / 2)) / pow
 
-    est <-  amce + se*rnorm(sims)
+    exp_typeM <- vector(mode = "numeric", length = length_output)
 
-    sig <- abs(est) > se*qnorm(1-alpha/2)
+    for (i in seq(1, length_output)) {
+      est <- amce[i] + se[i] * rnorm(sims[i])
 
-    exp_typeM <- abs(mean(abs(est)[sig])/amce)
+      sig <- abs(est) > se[i] * qnorm(1 - alpha[i] / 2)
 
-    return(data.frame(n=n, type_s=type_s, exp_typeM=exp_typeM, amce=amce, power=power, alpha = alpha, levels = levels, delta0 = delta0))
+      exp_typeM[i] <- abs(mean(abs(est)[sig]) / amce[i])
+    }
 
+    return(data.frame(n = n, type_s = type_s, exp_typeM = exp_typeM, amce = amce, power = power, alpha = alpha, levels = levels, delta0 = delta0))
   }
 
-  else if (is.null(power)){
+  else if (is.null(power)) {
+    delta0 <- 0.5 - (amce * treat.prob)
 
-    if (all(sapply(list(length(amce), length(alpha), length(n), length(levels)), function(x) x == 1)) == FALSE) {
-      stop("Please provide scalar values or apply a functional.")
-    }
-    
-    delta0 = 0.5 - (amce * treat.prob)
+    se <- (sqrt((
+      ((delta0 + amce) * (1 - (delta0 + amce))) / 0.5 +
+        ((delta0) * (1 - delta0)) / 0.5
+    )) / sqrt(2 * (n / levels)))
 
-    se = (sqrt( (
-      ((delta0 + amce)*(1 - (delta0 + amce)))/0.5 +
-        ((delta0)*(1 - delta0))/0.5
-    ))/sqrt(2*(n/levels)))
+    z <- amce / se
 
-    z=amce/se
+    power <- pnorm(z - qnorm(1 - alpha / 2)) + pnorm(-z - qnorm(1 - alpha / 2))
 
-    power=pnorm(z-qnorm(1-alpha/2)) + pnorm(-z-qnorm(1-alpha/2))
+    type_s <- pnorm(-z - qnorm(1 - alpha / 2)) / power
 
-    type_s = pnorm(-z-qnorm(1-alpha/2))/power
+    exp_typeM <- vector(mode = "numeric", length = length_output)
 
-    est <-  amce + se*rnorm(sims)
+    for (i in seq(1, length_output)) {
+      est <- amce[i] + se[i] * rnorm(sims[i])
 
-    sig <- abs(est) > se*qnorm(1-alpha/2)
+      sig <- abs(est) > se[i] * qnorm(1 - alpha[i] / 2)
 
-    exp_typeM <- abs(mean(abs(est)[sig])/amce)
+      exp_typeM[i] <- abs(mean(abs(est)[sig]) / amce[i])
+    }
 
-    return(data.frame(power=power, type_s=type_s, exp_typeM=exp_typeM,  amce=amce, n=n, alpha = alpha, levels = levels, delta0 = delta0 ))
+    return(data.frame(power = power, type_s = type_s, exp_typeM = exp_typeM, amce = amce, n = n, alpha = alpha, levels = levels, delta0 = delta0))
   }
-
-}
+}