A New Look at p-Values for Randomized Clinical Trials 23

# https://t.co/iIu46aY3mZ

# http://www.stat.columbia.edu/~gelman/research/published/pval_RCTs_rev3.pdf

# Evaluating a shrinkage estimator for the treatment effect in clinical trials 23

# https://onlinelibrary.wiley.com/doi/10.1002/sim.9992

# b= β+N(0, s²), SNR= β/σ, z=SNR+N(0, 1), β true effect

# We compute the shrinkage estimator β̂ and its SD σ from the unbiased estimator b and its SD s with the following R code:

shrinkage <- function (b, s) {

z <- b/s

p<- c(0.32, 0.31, 0.3, 0.07) # from Table 3

sigma <- c(0.61, 1.42, 2.16, 5.64)

sigma2 <- sigma^2

q <-p*dnorm (z, 0, sqrt (sigma2+1))

q <-q/sum (q) # conditional mixing probs

pm <-b*sigma2/(sigma2+1) # conditional means

pv <- s^2*sigma2/(sigma2+1) # conditional variances

data.frame(q, pm, pv)

}

# For example, if we observe b = 0.4 and s = 0.3 then the usual 95% confidence interval is -0.2 to 1.0. We can compute the shrinkage estimator as follows:

shrink <- shrinkage (b=0.4, s=0.3)

betahat <- sum (shrink$q*shrink$pm)

# We find β hat = 0.23. We can compute SD σ of the mixture distribution as follows:

pm2 <- sum(shrink$q*shrink$pm^2)

ps2 <- sum (shrink$q* shrink$pv)

sigma <- sqrt (ps2+ pm2 - betahat^2)

# We find σ =0.25, so the 95% interval is -0.25 to 0.72

calculate_power <- function(SNR) {

  power = pnorm(-1.96 - SNR) + 1 - pnorm(1.96 - SNR)

  return(power)

}

set.seed(123)  # Monte carlo simulation

# Define the parameters of the mixture distribution

proportions <- c(0.32, 0.31, 0.30, 0.07)

means <- c(0.00, 0.00, 0.00, 0.00)

std_devs <- c(0.61, 1.42, 2.16, 5.64)

# Function to generate data from a mixture distribution

generate_mixture <- function(n) {

# Determine the number of samples from each component

  samples_per_component <- sample(1:4, size = n, replace = TRUE, prob = proportions)

# Generate the samples

  rnorm(n, mean = means[samples_per_component], sd = std_devs[samples_per_component])

}

# Generate a sample from the mixture distribution

sample_size <- 10000

snr_sample <- generate_mixture(sample_size)

# Step 2: Add standard normal noise

snr_sample_with_noise <- snr_sample + rnorm(sample_size)

# Step 3: Compute two-sided p-value

# Wrong: p_values <- 2 * (1 - abs(snr_sample_with_noise))

p_values <- 2 * (1 - pnorm(abs(snr_sample_with_noise)))

# Step 4: Add another set of noise for replication study

replication_z <- snr_sample_with_noise + rnorm(sample_size)

# Step 5: Calculate statistical power (Transformation of SNR)

# Example transformation, replace with the specific transformation you need

power <- pnorm(-1.96 - replication_z) + 1 - pnorm(1.96 - replication_z)

# Result is a set of SNR, p-values, z-statistics, and power estimates

results <- data.frame(snr_sample, p_values, replication_z, power)

# View the first few rows of the results

head(results)

if(FALSE){

alpha <- 0.05

power <- 0.8

# For a two-tailed test, the critical value for alpha is the 97.5th percentile

z_alpha <- qnorm(1 - alpha / 2)

# The z-score for power is the value such that the area to its left under the standard normal curve is 0.9

z_beta <- qnorm(power)

# Effect size in standard errors

(effect_size_se <- z_alpha + z_beta)

}