# R/internal.R In MRTSampleSize: A Sample Size Calculator for Micro-Randomized Trials

#### Defines functions PowerCalculationSampleSizePowerCalparavalidateEffectAvailParametersvalidateProximalEffectParametersvalidateParametersgenerateTaugenerateBetagetTrueBeta

```#' @title  Internal MRT functions
#' @description   These internal functions are not usually called
#' directly by the user. Their names and capabilities may change
#' without warning from one version of \pkg{MRT} to the next.
#' @importFrom stats qf pf
#' @noRd
#'

getTrueBeta <- function(q1 = 0, q2 = 0.3, q3 = 28, days){

### output (coefficients of) standardized treatment effect or availability (quadratic
### form) given q1, q2, q3 ###
### only used when q3 (in days) is larger than half of the study

beta = matrix(0, 3)
beta = q1
a = sum(c(1:(days-1)))
b = sum(c(1:(days-1))^2)
mat = t(matrix(c(1, 2*(q3-1),a,b), 2, 2))
beta[2:3] = solve(mat) %*% matrix(c(0, q2 * days))

return(beta)
}

### Quadratic class of proximal treatment effect ###
N <- days*occ_per_day
M <- beta_mean
I <- beta_initial

beta1 <- getTrueBeta(q1 = 0, q2 = M-I, q3 = beta_quadratic_max , days = days)
b <- beta1
c <- beta1
sequence <- c(0:(days-1))
beta_input <- I + b*sequence + c * sequence^2
} else if(beta_shape == "constant"){

### Constant class of proximal treatment effect ###
beta_input <- replicate(days, beta_mean)

} else if(beta_shape == "linear"){

### Linear class of proximal treatment effect ###
initial <- beta_initial
mean <- beta_mean
days <- days

range <- (initial - mean) * 2

if(days == 1) {
beta_input <- mean
} else if (range == 0) {
beta_input <- replicate(days, mean)
} else {
num <- range / (days - 1)
beta_input <- seq(from = mean + range / 2,
to = mean - range / 2,
by = -num)
}
}
return(beta_input)
}

### Quadratic class of proximal treatment effect ###
N <- days*occ_per_day
M <- tau_mean
I <- 2*M - tau_initial

beta1 <- getTrueBeta(q1 = 0, q2 = M-I, q3 = tau_quadratic_max , days = days)
b <- beta1
c <- beta1
sequence = c(0:(days-1))
tau_input <- I + b*sequence + c * sequence^2
tau_input <- 2*M-tau_input
} else if(tau_shape == "constant"){

### Constant class of expected availability ###
tau_input <- replicate(days, tau_mean)

} else if(tau_shape == "linear"){

### Linear class of expected availability ###
initial <- tau_initial
mean <- tau_mean
days <- days

range <- (initial - mean) * 2

if(days == 1) {
tau_input <- mean
} else if (range == 0) {
tau_input <- replicate(days, mean)
} else {
num <- range / (days - 1)
tau_input <- seq(from = mean + range / 2,
to = mean - range / 2,
by = -num)
}
}
return(tau_input)
}

validateParameters <- function(days,occ_per_day,prob,
sigLev){

if(days != round(days)){
stop("Error: Please enter integer values for the number of days")
}
if(days <= 0){
stop("Error: Please specify the number of days greater than 0")
}
if(occ_per_day != round(occ_per_day)){
stop("Error:Please enter integer for the number of occasions per day")
}
if(occ_per_day <= 0){
stop("Error: Please specify the number of occasions per day greater than 0")
}

if(length(prob) == 1){
if(prob <= 0){
stop("Error: Please specify the randomization probability greater than 0")
}
if(prob >= 1){
stop("Error: Please specify the randomization probability less than 1")
}
} else {
if(length(prob) != days * occ_per_day)
{
stop("Error: The number of prob does not match the number of decision times provided in Study Setup.")
}
if(max(prob) > 1)
{
stop("Error: The provided randomization probability is greater than 1 for one or more decision times.")
}
if(min(prob) < 0)
{
stop("Error: The provided randomization probability is less than 0 for one or more decision times.")
}
}

if(!((dimB == as.integer(dimB)) & (dimB > 0))){

stop("Error: Please specify the parameters in the main effect to be a positive integer")

}

if(sigLev < 0){
stop("Error: Please specify the significance level greater than or equal to 0")
}
if(sigLev > 1){
stop("Error: Please specify the significance level less than or equal to 1")
}

validateProximalEffectParameters(beta_shape,beta_mean,beta_initial)

validateEffectAvailParameters(tau_shape,tau_mean)

}

if(beta_shape != "constant"
&& beta_shape !="linear"
stop("Please select a pattern for proximal treatment effect")
}

if(beta_mean <= 0){
stop("Error: Please specify an average standardized effect greater than 0")
}

if(beta_mean >= 1){
stop("Error: Please specify an average standardized effect less than 1")
}

if(beta_initial < 0){
stop("Error: Please specify a Standardized Initial Effect greater than or equal to 0")
}

if(beta_initial > beta_mean){
stop("Error: Please specify the standardized Initial Effect less than or equal to average standardized effect")
}
}
}

validateEffectAvailParameters <- function(tau_shape,tau_mean){
if(tau_shape != "constant"
&& tau_shape !="linear"
stop("Please select a pattern for expected availability")
}

if(tau_mean <= 0){
stop("Error: Please specify the mean of availability greatert than 0")
}

}

para <- function(order, Total, days, occ_per_day){

if(order == 0){

return(matrix(1, Total, 1))

}else{

Z <- matrix(0, Total, order + 1)
Z[,1] <- 1;
for(i in 1:order){

Z[, i+1] <- (rep(c(0:(days-1)), each =occ_per_day))^i

}

}

return(Z)

}

PowerCal <- function(p, q, N, d, beta, alpha){

C <- N * t(d) %*% beta %*% d ## noncentral para.
df2 <- N-(p+q) ## degrees of freedom
adj_c <- qf(1-alpha, df1 = p, df2 = df2) ### critical value

# output the power
power <- 1 - pf(q = adj_c, df1 = p, df2 = df2 , ncp = C)

return(power)
}

SampleSize <- function(days, occ_per_day, beta_t, tau, delta, alpha0, beta0, p, q, Nmax=1000){

### MRT Sample Size Calculation (Liao et. al 2015) ###

### INPUT:
### standardized treatment effect (beta_t), expected availability (tau),
### randomization probability (delta), significance level (alpha0)
### desired power (beta0), study setup (days, occ.per.days)
### p, q = number of parameters in Z and B (assume Z are quadratic, linear or constant)
### Nmax is the maximum number of participants to acheive the power

### OUTPUT: Sample size

Total <- days * occ_per_day;
Z <- para(order = p-1, Total, days, occ_per_day);

stopifnot(length(beta_t) == Total)
stopifnot(length(tau) == Total);

if(length(delta) == 1){

delta <- rep(delta, Total);

}else if(length(delta) == days){

delta <- rep(delta, each = occ_per_day);

}else{

stopifnot(length(delta) == Total)
}

beta <- t(Z) %*% diag(tau * delta * (1-delta)) %*% Z

d <- solve(t(Z) %*% diag(tau) %*% Z) %*% t(Z)  %*% diag(tau) %*% beta_t

N.all <- c(10:Nmax);

N <- Nmax
for( k in N.all){
power <- PowerCal(p,q,k,d,beta,alpha0)
if(power >= beta0){
N <- k
break
}
}

if(N == Nmax){

stop(paste("Cannot attain",beta0,"power when sample size is below", Nmax))

}

return(N)
}

PowerCalculation <- function(days, occ_per_day, N, beta_t, tau, delta, alpha0, p, q){

### MRT Power Calculation (Liao et. al 2015) ###

### INPUT:
### sample size (N)
### standardized treatment effect (beta_t), expected availability (tau),
### randomization probability (delta), significance level (alpha0)
### p, q = number of parameters in Z and B (assume Z are quadratic, linear or constant)

### OUTPUT: (Estimated) power

Total <- days * occ_per_day ;
Z <- para(order = p-1, Total, days, occ_per_day);

stopifnot(length(beta_t) == Total)
stopifnot(length(tau) == Total);

if(length(delta) == 1){

delta <- rep(delta, Total);

}else if(length(delta) == days){

delta <- rep(delta, each = occ_per_day);

}else{

stopifnot(length(delta) == Total)
}

beta <- t(Z) %*% diag(tau * delta * (1-delta)) %*% Z

d <- solve(t(Z) %*% diag(tau) %*% Z) %*% t(Z)  %*% diag(tau) %*% beta_t

power <- PowerCal(p,q,N,d,beta,alpha0)

return(power)
}
```

## Try the MRTSampleSize package in your browser

Any scripts or data that you put into this service are public.

MRTSampleSize documentation built on Sept. 13, 2020, 5:16 p.m.