logPower_ttest: Power Calculation For t-tests on Log-Transformed Data

In bozenne/butils: Various useful functions

Power Calculation For t-tests on Log-Transformed Data

Description

Power calculation for t-tests on log-transformed data.

Usage

logPower_ttest(
  type,
  equivalence = FALSE,
  gamma,
  sigma2,
  mu = NULL,
  rho = NULL,
  n = NULL,
  ratio.n = 1,
  sig.level = 0.05,
  power = 0.8,
  method.meanvar = "lognorm",
  method.cor = "uniroot",
  n.large = 1e+05,
  trace = TRUE,
  ncpus = NULL
)

Arguments

type	[character] type of study (parameter of interest): type="paired one.sample" (mean change within the group), or type="two.sample" (difference in mean between two independent groups), or type="paired two.sample" (difference in mean change between two independent groups)
equivalence	If TRUE, the alternative hypothesis is no change/difference. If FALSE, the null hypothesis is no change/difference.
gamma	[numeric] When considering two groups: relative difference in expected values between the two groups on the original scale. When considering a single group: expected value of the individual relative differences on the original scale.
sigma2	[numeric] When the argument mu is unspecified: coefficient of variation (i.e. standard error divided by expectation) of the outcome on the original scale. When the argument mu is specified: variance of the outcome on the original scale.
mu	[numeric] expected value of the outcome on the original scale for the control group or at baseline.
rho	[numeric 0-1] correlation coefficient between the paired values on the original scale.
n	[integer >1] Sample size.
ratio.n	[numeric >1] ratio between the sample size in larger and the smaller group. By default 1. Note that the larger group is the reference group (with expectation mu) and the smaller group the active group (with expectation gamma*mu).
sig.level	[numeric 0-1] type 1 error
power	[numeric 0-1] statistical power (i.e. complement to 1 of the type 2 error)
method.meanvar	[character] method used to identify the expected mean and variance on the log scale. Either assumes that the outcome is log-normally distributed on the original scale (method.meanvar = "lognorm"), or that it is normally distributed on the log-scale (method.meanvar = "lognorm2").
method.cor	[character] method used to identify the correlation on the log scale: one of "uniroot", "optim", or "taylor".
n.large	[integer, >0] sample size used to indentify the correlation coefficient or assess the error made when identifying the parameters. Should be large.
trace	[logical] Should a progress bar be displayed when estimating the power for various combinaisons of parameters?
ncpus	[integer, >0] Number of cores to be used, i.e., how many processes can be run simultaneously.

Details

The power calculation is performed under the assumption that on the log-scale, the variance of the outcome is the same between the two groups or over time.

References

Shein-Chung Chow , Jun Shao & Hansheng Wang (2002). A note on sample size calculation for mean comparisons based on noncentral t-statistics, Journal of Biopharmaceutical Statistics, 12:4, 441-456, DOI: 10.1081/BIP-120016229 Gerald Van Belle and Donald C. Martin. Sample size as a function of coefficient of variation and ratio of means (1993). The American Statistician 47(3) 165-167.

Examples

if(require(MESS) && require(mvtnorm)){

rho <- 0.6
mu <- c(0,0.275)
Sigma <- matrix(c(0.5^2,0.5^2*rho,0.5^2*rho,0.5^2),2,2)

#### 1- one sample with paired data ####
X <- exp(rmvnorm(1e5, mean = mu, sigma = Sigma))

## 1a) H0: no difference
# proposed solution
logPower_ttest(mu = mean(X[,1]),
               sigma2 = var(X[,1]),
               gamma = 0.3, rho = cor(X[,1],X[,2]),
               type = "paired one.sample")

logPower_ttest(sigma2 = sd(X[,1]) / mean(X[,1]),
               gamma = 0.3, rho = cor(X[,1],X[,2]),
               type = "paired one.sample")

# no log-transform
beta <- mean(X[,1])*0.3
sigma <- sd(X[,2]-X[,1])
power_t_test(delta = beta/sigma, power = 0.80, type = "one.sample")

# using log-transform
beta <- mean(log(X[,2])-log(X[,1]))
sigma <- sd(log(X[,2])-log(X[,1]))
power_t_test(delta = beta/sigma, power = 0.80, type = "one.sample")

# using simulation
warper <- function(i, n){
X <- rmvnorm(n, mean = mu, sigma = Sigma)
t.test(X[,2]-X[,1])$p.value
}
mean(unlist(lapply(1:1e3,warper, n = 39))<=0.05)

## 1b) H0: >10% difference
# proposed solution
logPower_ttest(mu = mean(X[,1]),
               sigma2 = var(X[,1]),
               gamma = 0.1, rho = cor(X[,1],X[,2]),
               type = "paired one.sample", equivalence = TRUE)

## using log-transform
power_t_test(delta = log(1 + 0.1)/sd(log(X[,2]/X[,1])), power = 1 - (1-0.80)/2,
             type = "one.sample", alternative = "one.sided")
power_t_test(delta = -log(1 - 0.1)/sd(log(X[,2]/X[,1])), power = 1 - (1-0.80)/2,
             type = "one.sample", alternative = "one.sided")

## using simulation
warper <- function(i, n){
X <- exp(rmvnorm(n, mean = c(0,0), sigma = Sigma))
Z <- log(X[,2]/X[,1])
tt1 <- t.test(Z, mu = log(0.9), alternative = "greater")$p.value
tt2 <- t.test(Z, mu = log(1.1), alternative = "less")$p.value
return(c(colMeans(X),mean(Z),max(tt1,tt2)))
}
ls.res <- do.call(rbind,lapply(1:1e3,warper, n = 305))
mean(ls.res[,4]<=0.05)

#### 2- two independent samples ####
X <- rlnorm(1e5, meanlog = mu[1], sdlog = 0.5)
Y <- rlnorm(1e5, meanlog = mu[2], sdlog = 0.5)

## 2a) H0: no difference
# proposed solution
logPower_ttest(mu = mean(X),
               sigma2 = var(X),
               gamma = 0.3, type = "two.sample")

logPower_ttest(sigma2 = sd(X)/mean(X),
               gamma = 0.3, type = "two.sample")

## no log-transform
beta <- mean(X)*0.3
sigma <- sd(X)
power_t_test(delta = beta/sigma, power = 0.80, type = "two.sample")

## using log-transform
beta <- mean(log(Y))-mean(log(X))
sigma <- sqrt(var(log(Y))/2+var(log(X))/2)
power_t_test(delta = beta/sigma, power = 0.80, type = "two.sample")

## using simulation
warper <- function(i, n){
X <- rlnorm(n, meanlog = mu[1], sdlog = 0.5)
Y <- rlnorm(n, meanlog = mu[2], sdlog = 0.5)
t.test(log(X),log(Y))$p.value
}
mean(unlist(lapply(1:1e3,warper, n = 53))<=0.05)

## graphical display
df.power <- logPower_ttest(mu = mean(X),
               sigma2 = var(X),
               gamma = seq(0.1,0.4, length.out = 30),
               n = seq(20, 75, length.out = 30),
               type = "two.sample")

if(require(ggplot2) && require(metR)){
gg <- ggplot(df.power, aes(x = n1, y = gamma))
gg <- gg + geom_tile(aes(fill = power))
gg <- gg + geom_contour(aes(z = power))
gg <- gg + scale_fill_gradientn(colours = terrain.colors(10), limits = c(0,1))
gg <- gg + geom_label_contour(aes(z = power), skip = 0)
gg <- gg +  scale_y_continuous(labels = scales::percent)
gg <- gg + xlab("sample size in each group") + ylab("mean difference between the groups")
gg
}

#### 3- two independent samples with paired data ####
X <- exp(rmvnorm(1e5, mean = rep(mu[1],2), sigma = Sigma))
Y <- exp(rmvnorm(1e5, mean = mu, sigma = Sigma))

## 3a) H0: no difference
# proposed solution
logPower_ttest(mu = mean(X[,1]),
               sigma2 = var(X[,1]),
               gamma = 0.3, rho = cor(X[,1],X[,2]),
               type = "paired two.sample")

logPower_ttest(sigma2 = sd(X[,1])/mean(X[,1]),
               gamma = 0.3, rho = cor(X[,1],X[,2]),
               type = "paired two.sample")

logPower_ttest(sigma2 = sd(X[,1])/mean(X[,1]),
               gamma = 0.3, rho = cor(X[,1],X[,2]),
               type = "paired two.sample", ratio.n = 2)

## no log-transform
beta <- mean(X[,1])*0.3
sigma <- sd(X[,1])
power_t_test(delta = beta/sigma, power = 0.80, type = "two.sample")

## using log-transform
beta <- mean(log(Y[,2])-log(Y[,1]))-mean(log(X[,2])-log(X[,1]))
sigma <- sqrt(var(log(Y[,2])-log(Y[,1]))/2+var(log(X[,2])-log(X[,1]))/2)
power_t_test(delta = beta/sigma, power = 0.80, type = "two.sample")

## using simulation
warper <- function(i, n){
X <- exp(rmvnorm(n, mean = rep(mu[1],2), sigma = Sigma))
Y <- exp(rmvnorm(n, mean = mu, sigma = Sigma))
t.test(log(X[,2])-log(X[,1]),log(Y[,2])-log(Y[,1]))$p.value
}
mean(unlist(lapply(1:1e3,warper, n = 73))<=0.05)
}

bozenne/butils documentation built on Oct. 14, 2023, 6:19 a.m.