Power Calculation for Interaction Term Binary and Continuous Variable Cox Regression

Appendix 1

A1.1. Exponential distribution of event times (covariate effect is a linear effect of time)

If event times follow an exponential distribution, the baseline hazard function is equal to ${\lambda }_0\left(t\right)=\lambda$ . Then the cumulative hazard function is equal to $\begin{array}{rl}H\left(t\right|\mathbf{\text{x}},z)& ={\int }_0^t{\lambda }_0\left(s\right)e^{\gamma \mathbf{\text{x}}+\beta \left(s\right)z}\text{d}s={\int }_0^t\lambda e^{\gamma \mathbf{\text{x}}}{e}^{\beta \left(s\right)z}\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{\int }_0^t{e}^{\beta \left(s\right)z}ds\\ & =\lambda e^{\gamma \mathbf{\text{x}}}{\int }_0^t{e}^{({\beta }_0+{\beta }_{1}s)z}\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{\int }_0^t{e}^{{\beta }_{0}z}{e}^{{\beta }_{1}sz}ds=\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t{e}^{{\beta }_{1}sz}ds\\ & =\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{{\beta }_{1}z}[e^{{\beta }_{1}sz}{]}_0^t=\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{{\beta }_{1}z}[e^{{\beta }_{1}tz}-1\left]\right(\text{provided}{\beta }_{1}z\ne 0).\end{array}$ The distribution function of the event times is equal to $F\left(t\right|\mathbf{\text{x}},z)=1-\text{exp}(-H\left(t\right|\mathbf{\text{x}},z)=1-\text{exp}\left(-\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{{\beta }_{1}z}[e^{{\beta }_{1}tz}-1]\right).$ Therefore, $u=\text{exp}\left(-\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{{\beta }_{1}z}[e^{{\beta }_{1}tz}-1]\right)\sim U[0,1].$ Solving for t gives: $\text{log}\left(u\right)=-\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{{\beta }_{1}z}(e^{{\beta }_{1}tz}-1),$ thus $\frac{-{\beta }_{1}z\text{log}\left(u\right)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}=e^{{\beta }_{1}tz}-1,$ thus $1-\frac{{\beta }_{1}z\text{log}\left(u\right)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}=e^{{\beta }_{1}tz},$ thus $\text{log}\left(1-\frac{{\beta }_{1}z\text{log}\left(u\right)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\right)={\beta }_{1}zt,$ resulting in $t=\frac{1}{{\beta }_{1}z}\text{log}\left(1-\frac{{\beta }_{1}z\text{log}\left(u\right)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\right),\text{provided that}{\beta }_{1}z\ne 0.$

If β ₁ z = 0, then the effect of the covariate is time-invariant and methods based on those described by Bender et al. can be used to generate an event time.

A1.2. Weibull distribution of event times (covariate effect is a linear effect of time)

The baseline hazard function is equal to ${\lambda }_0\left(t\right)=\lambda \upsilon t^{\upsilon -1}$ . Then the cumulative hazard function is $\begin{array}{rl}H\left(t\right|\mathbf{\text{x}},z)& ={\int }_0^t{\lambda }_0\left(s\right)e^{\gamma \mathbf{\text{x}}+\beta \left(s\right)z}\text{d}s={\int }_0^t\lambda \upsilon s^{\upsilon -1}{e}^{\gamma \mathbf{\text{x}}}{e}^{\beta \left(s\right)z}\text{d}s=\lambda \upsilon e^{\gamma \mathbf{\text{x}}}{\int }_0^t{s}^{\upsilon -1}{e}^{({\beta }_0+{\beta }_{1}s)z}\text{d}s\\ & =\lambda \upsilon e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t{s}^{\upsilon -1}{e}^{{\beta }_{1}sz}\text{d}s=\lambda \upsilon e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}\left[t^{\upsilon }(-{\beta }_{1}tz{)}^{-\upsilon }\right(\mathrm{\Gamma }\left(\upsilon \right)-\mathrm{\Gamma }(\upsilon ,-{\beta }_{1}tz)\left)\right].\end{array}$ The cumulative hazard function involves a difference of a Gamma function and an Incomplete Gamma function. Then, the cumulative distribution of event times is $F\left(t\right|\mathbf{\text{x}},z)=1-\text{exp}[-H\left(t\right|\mathbf{\text{x}},z),$ and $u=\text{exp}[-\lambda \upsilon e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}(t^{\upsilon }(-{\beta }_{1}tz{)}^{-\upsilon }\left(\mathrm{\Gamma }\right(\upsilon )-\mathrm{\Gamma }(\upsilon ,-{\beta }_{1}tz\left)\right)\left)\right]\sim U(0,1).$

An event time can be generated by solving the above expression for t and evaluating the expression at $u\sim U(0,1)$ . However, a closed-form expression for the inverse of the above expression cannot be obtained. As an alternative approach, one could use numerical methods to find an approximate inverse to the above expression for a given value of u.

A1.3. Gompertz distribution of event times (covariate effect is a linear effect of time)

The baseline hazard function is equal to ${\lambda }_0\left(t\right)=\lambda e^{\alpha t}$ . Then the cumulative hazard function is $\begin{array}{rl}H\left(t\right|\mathbf{\text{x}}\text{, z)}& ={\int }_0^t{\lambda }_0\left(s\right)e^{\gamma \mathbf{\text{x}}+\beta \left(s\right)z}\text{d}s={\int }_0^t\lambda e^{\alpha s}{e}^{\gamma \mathbf{\text{x}}}{e}^{({\beta }_0+{\beta }_{1}s)z}\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{\int }_0^t{e}^{\alpha s}{e}^{{\beta }_{0}z}{e}^{{\beta }_{1}zs}\text{d}s\\ & =\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t{e}^{\alpha s}{e}^{{\beta }_{1}zs}\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t{e}^{(\alpha +{\beta }_{1}z)s}\text{d}s\\ & =\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{\alpha +{\beta }_{1}z}[e^{(\alpha +{\beta }_{1}z)s}{]}_0^t=\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{\alpha +{\beta }_{1}z}[e^{(\alpha +{\beta }_{1}z)t}-1]\text{provided that}\alpha +{\beta }_{1}z\ne 0.\end{array}$ Then the cumulative distribution of event times is $F\left(t\right|\mathbf{\text{x}},z)=1-\text{exp}(-H\left(t\right|\mathbf{\text{x}},z\left)\right)$ and $u=\text{exp}(-H(t|\mathbf{\text{x}},z))\sim U(0,1).$ To invert this function, we have that $u=\text{exp}\left(-\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{\alpha +{\beta }_{1}z}[e^{(\alpha +{\beta }_{1}z)t}-1]\right),$ thus $\text{log}\left(u\right)=-\frac{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}{\alpha +{\beta }_{1}z}[e^{(\alpha +{\beta }_{1}z)t}-1],$ $\frac{-(\alpha +{\beta }_{1}z)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\text{log}\left(u\right)=e^{(\alpha +{\beta }_{1}z)t}-1,$ $1-\frac{(\alpha +{\beta }_{1}z)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\text{log}\left(u\right)=e^{(\alpha +{\beta }_{1}z)t},$ so that $\text{log}\left(1-\frac{(\alpha +{\beta }_{1}z)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\text{log}\left(u\right)\right)=(\alpha +{\beta }_{1}z)t,$ so that $t=\frac{1}{\alpha +{\beta }_{1}z}\text{log}\left(1-\frac{(\alpha +{\beta }_{1}z)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}\text{log}\left(u\right)\right),\text{provided that}\alpha +{\beta }_{1}z\ne 0.$ If $\alpha +{\beta }_{1}z=0$ then the calculation of the cumulative hazard function must be modified as follows: $H\left(t\right|\mathbf{\text{x}},z)=\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t{e}^0\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}{\int }_0^t\text{d}s=\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}t.$ Then, we have that $u=\text{exp}(-H(t|\mathbf{\text{x}},z))=\text{exp}(-\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}t)$ . To invert this expression, we have that $\text{log}\left(u\right)=-\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}t$ and hence that $t=\frac{-\text{log}\left(u\right)}{\lambda e^{\gamma \mathbf{\text{x}}}{e}^{{\beta }_{0}z}}.$

Appendix 2

A2.1. R and SAS code for simulating event time data with an exponential distribution when the log-hazard ratio of the binary variable changes linear with time

# The following software code is provided for illustrative purposes only and comes with absolutely no warranty.

# Simulate data from a Cox-exponential model when the effect of

# treatment

# varies as a linear function of time.

N <- 10000

# Number of subjects in each simulated dataset.

# Parameter of the exponential distribution.

lambda <- 0.01

g1 <- log(1.5)

# Log-hazard ratio for the continuous covariate.

b0 <- log(1.1)

# Log-hazard ratio for the binary covariate at t = 0.

hr <- 1.005

b1 <- log(hr)

# Effect of time on the log-hazard ratio for the binary covariate.

prev <- 0.25

# Prevalence of the binary covariate.

set.seed(1)

# Set random number seed for reproducibility.

x <- rnorm(N)

z <- rbinom(N,size=1,prob=prev)

u <- runif(N,0,1)

event.time <- ifelse(b1*z==0,

 -log(u)/(lambda*exp(g1*x)),

 (1/b1*z) * log(1-b1*z*log(u)/(lambda*exp(g1*x + b0*z))))

simdata <- cbind(x,z,event.time)

write(t(simdata),ncol=3,file="exp.dat",append=T)

%macro dgp_exp(N=,lambda=,hr_cont=,b0=,b1=,prev=0,ranseed=);

*** N is the size of the simulated dataset;

*** lambda is the parameter for the exponential distribution of

event times;

*** hr_cont is the hazard ratio for the continuous covariate;

*** b0 is the hazard ratio for the binary variable at time t = 0;

*** b1 is the relative change in the effect of binary covariate

as a function;

**   of time;

*** prev is the prevalence of the binary covariate;

*** ranseed is a seed for random number generation to ensure ;

***   reproducibility of the results.;

data randata;

 call streaminit(&ranseed);

 do i = 1 to &N;

  /* Random uniform variable for generating event times */

  u = rand("Uniform");

  /* The binary covariate with a time-varying covariate effect */

  z = rand("binomial",&prev,1);

  /* The continuous covariate for which we are adjusting */

  x = rand("normal",0,1);

  /* Parameter for the exponential distribution */

  lambda = &lambda;

  gamma = log(&hr_cont);

  /* A one-standard deviation increase in x increases the hazard

    of the outcome by this amount */

  b0 = log(&b0);

  b1 = log(&b1);

  if (b1*z=0) then

   event_time = -log(u)/(lambda*exp(gamma*x));

  else

   event_time = (1/b1*z) * log(1 - b1*z*log(u)/(lambda*

    exp(gamma*x + b0*z)));

  event = 1;

  output;

 end;

run;

%mend dgp_exp;

%dgp_exp(N=10000,lambda=0.01,hr_cont=1.5,b0=1.1,b1=1.005,

 ranseed=1,prev=0.25);

A2.2. R code for simulating event time data with a Weibull distribution when the log-hazard ratio of the binary variable changes linear with time

# The following software code is provided for illustrative purposes only and comes with absolutely no warranty.

# Simulate data from a Cox-Weibull model when the effect of

# treatment varies as a linear function of time.

# There is not a closed-form expression for the inverse of the

# survival function (S).

# Numerical methods are used to invert the survival function.

# We use Brent's method for root finding to invert the function.

library(pracma)

# pracma package has a function for Brent's method for finding

# roots and a function for the incomplete Gamma function.

N <- 10000

# Number of subjects in each simulated dataset.

# Parameters of the Weibull distribution.

lambda <- 0.001

nu <- 1

g1 <- log(1.5)

# Log-hazard ratio for the continuous covariate.

b0 <- log(1.1)

# Log-hazard ratio for the binary covariate at t=0.

hr <- 1.001

b1 <- log(hr)

# Effect of time on the log-hazard ratio for the binary covariate.

set.seed(1)

# Set random number seed for reproducibility.

simdata.matrix <- NULL

# Define the survival function

S2 <- function(t){

 exp(-lambda*nu*exp(g1*x + b0*z) * ((t^nu) * (-b1*t*z)^(-nu) *

  (gamma(nu) - gammainc(-b1*t*z,nu)[2]) ) ) - u

}

for (i in 1:N){

 u <- runif(1,0,1)

 x <- rnorm(1)

 z <- rbinom(1,size=1,prob=0.25)

 if (z==0){

  event.time <- (-log(u)/(lambda*exp(g1*x)))^(1/nu)

 } else {

  brent.sol <- brent(S2,0.001,250000)

  event.time <- brent.sol$root

 }

 simdata <- c(x,z,event.time)

 simdata.matrix <- rbind(simdata.matrix,simdata)

}

 write(t(simdata.matrix),ncol=3,file="weibull.dat",append=T)

A2.3. R and SAS code for simulating event time data with a Gompertz distribution when the log-hazard ratio of the binary variable changes linear with time

# The following software code is provided for illustrative purposes only and comes with absolutely no warranty.

# Simulate data from a Cox-Gompertz model when the effect of

# treatment varies as a linear function of time.

N <- 10000

# Number of subjects in each simulated dataset.

# Parameters of the Gompertz distribution.

lambda <- 0.0001

alpha <- 0.025

g1 <- log(1.5)

# Log-hazard ratio for the continuous covariate.

b0 <- log(1.1)

# Log-hazard ratio for the binary covariate at t=0.

hr <- 1.005

b1 <- log(hr)

# Effect of time on the log-hazard ratio for the binary covariate.

prev <- 0.25

# Prevalence of the binary covariate.

set.seed(1)

# Set random number seed for reproducibility.

x <- rnorm(N)

z <- rbinom(N,size=1,prob=prev)

u <- runif(N,0,1)

event.time <- ifelse(alpha+b1*z==0,

 -log(u)/(lambda*exp(g1*x + b0*z)),

 (1/(alpha + b1*z)) * log(1 - (alpha + b1*z)*log(u)

/(lambda*exp(g1*x + b0*z))))

simdata <- cbind(x,z,event.time)

write(t(simdata),ncol=3,file="gompertz.dat",append=T)

%macro dgp_gompertz(N=,lambda=,alpha=,hr_cont=,b0=,b1=,

prev=0,ranseed=);

*** N is the size of the simulated dataset;

*** lambda and alpha are the parameters of the Gompertz

distribution;

*** hr_cont is the hazard ratio for the continuous covariate;

*** b0 is the hazard ratio for the binary variable at time t = 0;

*** b1 is the relative change in the effect of binary covariate

as a function;

**   of time;

*** prev is the prevalence of the binary covariate;

*** ranseed is a seed for random number generation to ensure;

***   reproducibility of the results.;

data randata;

 call streaminit(&ranseed);

 do i = 1 to &N;

  /* Random uniform variable for generating event times */

  u = rand("Uniform");

  /* The binary covariate with a time-varying covariate effect */

  z = rand("binomial",&prev,1);

  /* The continuous covariate for which we are adjusting */

  x = rand("normal",0,1);

  /* Parameters for the Gompertz distribution */

  lambda = &lambda;

  alpha = &alpha;

  gamma = log(&hr_cont);

  b0 = log(&b0);

  b1 = log(&b1);

  if alpha + b1*z = 0 then

   event_time = -log(u)/(lambda * exp(gamma*x + b0*z));

  else

   event_time = (1/(alpha + b1*z)) * log(1 - (alpha + b1*z)

   *log(u)/(lambda*exp(gamma*x + b0*z)));

  event = 1;

  output;

 end;

run;

%mend dgp_gompertz;

%dgp_gompertz(N=10000,lambda=0.0001,alpha=0.025,hr_cont=1.5,

 b0=1.1,b1=1.005,prev=0.25,ranseed=1);

savageprommented.blogspot.com

Source: https://www.tandfonline.com/doi/full/10.1080/00949655.2017.1397151

Share this post