Continuous Marker Logistic Power

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.

• Number of Simulations (1000): The number of simulations to run. For better approximations, increase the number of simulations.
• Total Sample Size
• Alpha Level (0.05): \(\alpha\), the 1 or 2-sided significance level
• Response Probability (marker value = 0): The probabiltiy of the event at marker value equal to 0
• Effect Size (Odds Ratio per Unit): The odds ratio associated with a one unit change in the marker
• Unit (Interquartile Range): Unit of marker effect is either Interquartile Range or Unstandardized
• Sided: One-sided or two-sided test
• Marker distribution (Uniform): Either standard uniform, or normal (standardized to a mean of zero), or lognormal
• Log-Normal Distribution Parameter (1): The spread/shape parameter (sigma) for each of the distributions. It is fixed to 1 for uniform and normal and defaulted to 1 for log-normal. The log-normal density is specified as \[ f(x) = \frac{e^{-(log(x) - \mu)^2/(2\sigma^2)}}{\sqrt{2\pi}{\sigma}x} \] where μ and σ are the mean and standard deviation of the logarithm of the normal random variable.

Output Items

• Estimated power
• Estimated effect and empirical histogram simulated expected confidence interval

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))
}