Design Analysis - Frequentist Code

Anna Vesely

2019-12-04

Cohen’s d

Given the vectors \(x_1\) and \(x_2\) with the observed values of two groups, the function cohen_d computes Cohen’s \(d\) effect size statistic with the corresponding t-statistic and p-value:

\[ d=\frac{\bar{\mathbf{x}_1}-\bar{\mathbf{x}_2}}{s_{\text{pool}}},\quad\quad t=d\,\sqrt{\frac{n_1\, n_2}{n_1+n_2}},\quad\quad p=2\,P(T_{r}\geq |t|) \]

where \(T\) follows a Student’s t-distribution with \(r=n_1+n_2-2\) degrees of freedom.

cohen_d <- function(x1, x2) {
  
  (m1 <- mean(x1, na.rm=TRUE))
  (m2 <- mean(x2, na.rm=TRUE))
  S1 <- sum((x1 - m1)^2, na.rm=TRUE)
  S2 <- sum((x2 - m2)^2, na.rm=TRUE)
  (n1 <- length(na.omit(x1)))
  (n2 <- length(na.omit(x2)))
  (df <- n1 + n2 - 2)
  
  (pool.s <- sqrt((S1 + S2) / df))
  (d <- (m2 - m1) / pool.s)
  
  s1 <- sd(x1, na.rm=TRUE)
  s2 <- sd(x2, na.rm=TRUE)
  (tstat <- (m2 - m1) / (pool.s * sqrt((n1 + n2) / (n1 * n2))))
  (pval <- pt(abs(tstat), df, lower.tail=FALSE) * 2)
  
  return(list(d=d, tstat=tstat, pval=pval))
}

Design Analysis Given a Plausible Value of d

Simulations of Expected Power and Errors

In the comparison between two groups, given a plausible effect size \(d_{\text{sim}}\), the sample sizes \(n_1\) and \(n_2\) of the groups and the significance level \(\alpha\), the function sim\_estimate computes the expected power \(\beta\), the type S error \(e_S\) (the risk that a significant effect is estimated in the wrong direction) and the type M error \(e_M\) (the average overestimation of an effect that emerges as significant).

In particular, the function simulates \(\mathbf{x}_1\) as \(n_1\) observations from a \(\mathcal{N}(0,1)\) distribution, and \(\mathbf{x}_2\) as \(n_2\) observations from a \(\mathcal{N}(d_{\text{sim}},1)\) distribution. Then the values of Cohen’s \(d\) and the corresponding test statistic \(t\) and p-value \(p\) are estimated. Such procedure is repeated \(B\) times (default \(10^4\)), and leads to:

\[\begin{align*} \hat{\beta}=\frac{\#\{i\,:\,p_i<\alpha\}}{B},\quad\quad \hat{e}_S=\frac{\#\{i\,:\,p_i<\alpha,\;t_i<0\}}{\#\{i\,:\,p_i<\alpha\}},\quad\quad \hat{e}_M=\frac{\text{mean}\{|d_i|\,:\,p_i<\alpha\}}{d_{\text{sim}}}. \end{align*}\]

Prospective Analysis: from Power to Sample Size

In the comparison between two groups, given a plausible value of effect size \(d_{\text{sim}}\), the significance level \(\alpha\) and the required power \(\beta\), the function estimates the sample size per group \(n\) needed to reach such power, as well as the type S error \(e_S\) and the type M error \(e_M\). The function also requires the minimum and the maximum size values (\(n_{\min}\) and \(n_{\max}\)), the number \(B\) of iterations and a tolerance \(\epsilon\).

Firstly, the function evaluates \(\hat{\beta}_{\max}\), the expected power given the effect size \(d_{\text{sim}}\) and the sample sizes \(n_1=n_2=n_{\max}\). If \(\hat{\beta}_{\max}<\beta\), the function stops, because a bigger value \(n_{\max}\) is needed.

Otherwise, let \((l^{(1)},\, u^{(1)})=(n_{\min},\, n_{\max})\). For each step \(j=1,2,\ldots\), the function defines \[\hat{n}^{(j)}=\lfloor\text{median}\{l^{(j)},\,l^{(j)}+1,\,\ldots,\, u^{(j)}\}\rfloor ,\] and \(\hat{\beta}^{(j)}\) as the expected power given the effect size \(d_{\text{sim}}\) and the sample sizes \(n_1=n_2=\hat{n}^{(j)}\).

Retrospective Analysis: from Sample Size to Power

In the comparison between two groups, given a plausible value of effect size \(d_{\text{sim}}\), the significance level \(\alpha\) and the sample size per group \(n\), the function estimates the power \(\beta\), the type S error \(e_S\) and the type M error \(e_M\). It makes use of the function .

Prospective or Retrospective Design Analysis

In the comparison between two groups, given a plausible value of effect size \(d_{\text{sim}}\) and the significance level \(\alpha\), the function performs a design analysis.

Design Analysis Given a Plausible Distribution for d

Sampling of Cohen’s d

The function requires the package . It samples a vector \(\mathbf{d}\) of \(B_0\) values from a given distribution. Such distribution is bounded on an interval \((d_{\min},d_{\max})\), and may be

Extended Retrospective Design Analysis

The function requires the package . In the comparison between two groups, it performs a retrospective design analysis. In particular, given the sizes \(n_1\) and \(n_2\) of the two groups, a plausible effect size (either a fixed value or a distribution) and the significance level \(\alpha\), it estimates the power \(\beta\), the type S error \(e_S\) and the type M error \(e_M\).

design_est<- function(n1=5, n2=n1, target_d=NULL,
                       target_d_limits=NULL,
                       distribution=c("uniform", "normal"), k=1/6,
                       sig.level=0.05, B=500, B0=500,
                       return_data=FALSE) {
  
  distribution <- match.arg(distribution)
  
  if (sum(sapply(list(target_d, target_d_limits), is.null)) != 1) {
    stop("Exactly one of 'target_d' and 'target_d_limits' must be NULL")
  }
  
  ########### fixed d
  if (!is.null(target_d)) {
    if (target_d <= 0) stop("Target_d must be greater than 0")
    temp <- sim_estimate(sim_d=target_d, target_d=target_d, n1=n1,
                         n2=n2, B=B)
    results <- list(power=temp[1,"power"], typeS=temp[1,"typeS"],
                    typeM=temp[1,"typeM"])
    outList <- list(call=list(n1=n1, n2=n2, target_d=target_d, B=B),
                    results=results)
    return(outList)
  }
  
  ############## plausible interval and distribution for d
  if (!is.null(target_d_limits)) {
    if (min(target_d_limits)<0)
      stop("Minimum of target_d_limits must be greater than 0")
    (call_est <- list(n1=n1, n2=n2, target_d_limits=target_d_limits,
                      distribution=distribution, k=k, B=B, B0=B0))
    (my <- (sum(target_d_limits) / 2))
    (sy <- diff(target_d_limits) * k)
    y <- sampling_d(target_d_limits, distribution=distribution,
                    B0=B0, k=k) ###
    sim_cohen <- y$y

    
    data <- t(sapply(1:B0, function(b){
      temp <- sim_estimate(sim_cohen[b], target_d=my, n1=n1, n2=n2,
                           B=B)
      c(temp[1,"power"], temp[1,"typeS"], temp[1,"typeM"])
    }))
    data <- data.frame(data)
    colnames(data) <- c("power", "typeS", "typeM")
    
    (power <- mean(data$power))
    (typeS <- mean(data$typeS, na.rm=TRUE))
    (typeM <- mean(data$typeM, na.rm=TRUE))
    
    results <- list(power=power, typeS=typeS, typeM=typeM)
    if (return_data) {
      results$data <- data
      
      outList <-  list(call=call_est, results=results)
    } else {
      outList <-  list(call=call_est, results=results)
    }
    return(outList)
  }
}