Description

This program uses simulation to approximate the power of a test of a marker association with binary outcome data assuming an logistic regression model. The calculations include a test for the log odds ratio using the model where the logit of the probability of the event is linear in a continuous marker distribution. The choice of the marker distributions include uniform, normal, and lognormal. All marker distributions are mean centered for the analysis. The magnitude of the effect size can be specified based on a unit of the untransformed marker or the interquartile range of the marker.

It is useful to note that if the underlying distribution is assumed to be uniform, and the logistic regression model is linear in the marker, then effect size associated with the interquartile range will correspond to approximately same effect as a marker dichotomized at the median (or other quantile).

In addition to power calculations, the approximate confidence intervals of parameter effects are also presented.

Input Items

The user is prompted to enter values for the following items. For items that have initial default values set, the default values are given in parentheses.

Output Items

Note

This is a simulation-based calculation and therefore may take several seconds to calculate depending both on the number of simulations and sample size.

Statistical Code

The program is written in R.

View Code


function(n = 200, alpha = .05, sided = 2, p0 = .1, OR1, unit, dist, dist_param,
         nsim = 1000) {
    
    # biomarkers are mean centered in calculations 
    # n = total sample size, with n/2 in each arm
    # p0 = response probability at biomarker=0
    # OR1 = odds ratio associated with unit change in biomarker
    # unit = if "IQR" the unit is the interquartile range
    # split = split at median (or other quantile)
    # nsim = number of simulations used in the power calculaton
    
    beta0 = log(p0 / (1 - p0))
    beta1 = log(OR1)
    
    # This is where the biomarker distibution is chosen
    if (dist == "unif") { dist.generate = runif }
    else if (dist == "norm") { dist.generate = rnorm }
    else if (dist == "lnorm") { dist.generate = rlnorm }
    
    testscorea = coefa = nevent = rep(NA, nsim)
    
    # Loop over simulations 
    for (j in 1:nsim) {
        
        # Generate Data #
        #---------------#
        
        biomarker = dist.generate(n, 0, dist_param)
        
        biomarker0 = biomarker - mean(biomarker)
        if (unit == "IQR") {
            biomarker0 = biomarker0 / IQR(biomarker0)
        }
        
        # generate response probabilities                
        eta = beta0 + beta1 * biomarker0
        resp.prob = exp(eta) / (1 + exp(eta))
        resp4 = 1 * (runif(n) < resp.prob)  # create response outcomes
        
        # Do Analysis #
        #-------------#
        
        glm.conta = glm(resp4 ~ biomarker0, family = binomial(link = "logit"))      
        testscorea[j] = summary(glm.conta)$coef[2, 3]**2 # score test
        coefa[j] = glm.conta$coefficients[2]
        
    }
    
    powercont = mean(testscorea > qchisq(1 - alpha * sided, df = 1))
    coefcont = quantile(coefa, c(alpha / sided, .5, 1 - alpha / sided))
    
    if (dist == "lnorm") {
        plot_name = "log_normal_distribution.png"
        png(filename = plot_name)
            x1 = qlnorm((1:50) / 50, meanlog = 0, sdlog = dist_param)
            y1 = dlnorm(x1, meanlog = 0, sdlog = dist_param, log = FALSE)
            plot(x1, y1, xlab = "x", ylab = "log normal density", type = "n")
            lines(x1, y1, lwd = 2)
            title("Log Normal Distribution") 
        dev.off()
    } else {
        plot_name = ""
    }

    result =
        list(power_continuous = round(powercont, digits = 3),
             effect = round(exp(coefcont[2]), digits = 3),
             lower_ci = round(exp(coefcont[1]), digits = 3),
             upper_ci = round(exp(coefcont[3]), digits = 3),
             session_key = strsplit(getwd(), "/")[[1]][4],
             log_normal_plot_name = plot_name)

    return(jsonlite::toJSON(result, pretty = TRUE))
}