# Power/type I error calculation for the proportional hazards model with piecewise constant hazard and random a0 Power/type I error calculation using the normalized power prior for the proportional hazards model with piecewise constant hazard and random $a_0$ ```r power.phm.random.a0( historical, n.subjects, n.events, n.intervals, change.points = NULL, samp.prior.beta, samp.prior.lambda, dist.enroll, param.enroll, rand.prob = 0.5, prob.drop = 0, param.drop = 0, dist.csr = "Constant", param.csr = 10000, min.follow.up = 0, max.follow.up = 10000, prior.beta.mvn = NULL, prior.beta.mean = rep(0, 50), prior.beta.sd = rep(1000, 50), prior.lambda0.hp1 = rep(10^(-5), 50), prior.lambda0.hp2 = rep(10^(-5), 50), prior.a0.shape1 = rep(1, 10), prior.a0.shape2 = rep(1, 10), prior.lambda.hp1 = rep(10^(-5), 50), prior.lambda.hp2 = rep(10^(-5), 50), lower.limits = NULL, upper.limits = rep(100, 50), slice.widths = rep(0.1, 50), nMC = 10000, nBI = 250, delta = 0, nullspace.ineq = ">", gamma = 0.95, N = 10000 ) ``` ## Arguments - `historical`: List of historical dataset(s). East historical dataset is stored in a list which contains four **named** elements: `time`, `event`, `X` and `S`. * `time` is a vector of follow up times. * `event` is a vector of status indicators. Normally 0=alive and 1=dead. * `X` is a matrix of covariates. The first column must be the treatment indicator. * `S` is a vector of integers, where each integer represents the stratum that the subject belongs to. For example, if there are three strata, S can take values 1, 2 or 3. - `n.subjects`: Number of subjects enrolled. - `n.events`: Number of events at which the trial will stop. - `n.intervals`: Vector of integers, indicating the number of intervals for the baseline hazards for each stratum. The length of the vector should be equal to the total number of strata. - `change.points`: List of vectors. Each vector in the list contains the change points for the baseline hazards for each stratum. The length of the list should be equal to the total number of strata. For a given stratum, if there is only one interval, then `change.points` should be `NULL` for that stratum. By default, we assign the change points so that the same number of events are observed in all the intervals in the historical data. These change points are used for data generation. The change points used during model fitting are assigned by default so that the same number of events are observed in all the intervals in the pooled historical and generated current data. - `samp.prior.beta`: Matrix of possible values of $\beta$ to sample (with replacement) from. Each row is a possible $\beta$ vector (a realization from the sampling prior for $\beta$). - `samp.prior.lambda`: List of matrices, where each matrix represents the sampling prior for the baseline hazards for each stratum. The number of columns of each matrix should be equal to the number of intervals for that stratum. - `dist.enroll`: Distribution for enrollment times. The choices are "Uniform" or "Exponential". - `param.enroll`: Parameter for the distribution of enrollment times. If `dist.enroll` is "Uniform", the enrollment times follow Unif(0, `param.enroll`). If `dist.enroll` is "Exponential", the enrollment times follow Exponential(rate=`param.enroll`). - `rand.prob`: Randomization probability for the treated group. The default value is 0.5. - `prob.drop`: Probability of subjects dropping out of the study (non-administrative censoring). The default value is zero. - `param.drop`: Parameter for dropout time simulations. The dropout times follow Unif(0, `param.drop`). The default value is zero. - `dist.csr`: Distribution for (administrative) censorship times. The choices are "Uniform", "Constant" and "Exponential". The default choice is "Constant". - `param.csr`: Parameter for the (administrative) censorship times. If `dist.csr` is "Uniform", the censorship times follow Unif(0, `param.csr`). If `dist.csr` is "Constant", the censorship times of all subjects are equal to `param.csr`. If `dist.csr` is "Exponential", the censorship times follow Exponential(rate=`param.csr`). The default value is 10^4. - `min.follow.up`: Minimum amount of time for which subjects are followed up. The default value is zero. - `max.follow.up`: Maximum amount of time for which subjects are followed up. The default value is 10^4. - `prior.beta.mvn`: List of multivariate normal approximations of the normalized power prior for $\beta$. Each element in the list is a list with three components, the mean vector, the covariance matrix and the weight of the multivariate normal distribution. The normalized power prior for $\beta$ is approximated by the weighted mixture of the multivariate normal distributions provided. By default (`prior.beta.mvn=NULL`), a single multivariate normal distribution is assumed. The user can use the `approximate.prior.beta` function to obtain samples of $\beta$ from the normalized power prior, and use any mixture of multivariate normals to approximate the normalized power prior for $\beta$. - `prior.beta.mean`: (Only applies if `prior.beta.mvn=NULL`) Vector of means of the normal initial prior on $\beta$. The default value is zero for all the elements of $\beta$. - `prior.beta.sd`: (Only applies if `prior.beta.mvn=NULL`) Vector of standard deviations of the normal initial prior on $\beta$. The default value is 10^3 for all the elements of $\beta$. - `prior.lambda0.hp1`: (Only applies if `prior.beta.mvn=NULL`) Vector of first hyperparameters of the Gamma prior on $\lambda_0$. The length of the vector should be equal to the dimension of $\lambda_0$, i.e., the total number of intervals for all strata. The default value is 10^(-5) for all the elements of $\lambda_0$. - `prior.lambda0.hp2`: (Only applies if `prior.beta.mvn=NULL`) Vector of second hyperparameters of the Gamma prior on $\lambda_0$. The length of the vector should be equal to the dimension of $\lambda_0$, i.e., the total number of intervals for all strata. The default value is 10^(-5) for all the elements of $\lambda_0$. - `prior.a0.shape1`: (Only applies if `prior.beta.mvn=NULL`) Vector of the first shape parameters of the independent beta priors for $a_0$. The length of the vector should be equal to the number of historical datasets. The default is a vector of one's. - `prior.a0.shape2`: (Only applies if `prior.beta.mvn=NULL`) Vector of the second shape parameters of the independent beta priors for $a_0$. The length of the vector should be equal to the number of historical datasets. The default is a vector of one's. - `prior.lambda.hp1`: Vector of first hyperparameters of the Gamma prior on $\lambda$. The length of the vector should be equal to the dimension of $\lambda$, i.e., the total number of intervals for all strata. The default value is 10^(-5) for all the elements of $\lambda$. - `prior.lambda.hp2`: Vector of second hyperparameters of the Gamma prior on $\lambda$. The length of the vector should be equal to the dimension of $\lambda$, i.e., the total number of intervals for all strata. The default value is 10^(-5) for all the elements of $\lambda$. - `lower.limits`: Vector of lower limits for parameters ($\beta$ and $\lambda$, in this order) to be used by the slice sampler. The length of the vector should be equal to the total number of parameters. The default is -100 for $\beta$ and 0 for $\lambda$ (may not be appropriate for all situations). - `upper.limits`: Vector of upper limits for parameters ($\beta$ and $\lambda$, in this order) to be used by the slice sampler. The length of the vector should be equal to the total number of parameters. The default is 100 for all parameters (may not be appropriate for all situations). - `slice.widths`: Vector of initial slice widths for parameters ($\beta$ and $\lambda$, in this order) to be used by the slice sampler. The length of the vector should be equal to the total number of parameters. The default is 0.1 for all parameters (may not be appropriate for all situations). - `nMC`: Number of iterations (excluding burn-in samples) for the slice sampler. The default is 10,000. - `nBI`: Number of burn-in samples for the slice sampler. The default is 250. - `delta`: Prespecified constant that defines the boundary of the null hypothesis. The default is zero. - `nullspace.ineq`: Character string specifying the inequality of the null hypothesis. The options are ">" and "<". If ">" is specified, the null hypothesis is $H_0$: $\beta_1$ $\ge$ $\delta$. If "<" is specified, the null hypothesis is $H_0$: $\beta_1$ $\le$ $\delta$. The default choice is ">". - `gamma`: Posterior probability threshold for rejecting the null. The null hypothesis is rejected if posterior probability is greater `gamma`. The default is 0.95. - `N`: Number of simulated datasets to generate. The default is 10,000. ## Returns Power or type I error is returned, depending on the sampling prior used. The posterior probabilities of the alternative hypothesis are returned. The average posterior means of $\beta$ and $\lambda$ are also returned. ## Details The proportional hazards model with piecewise constant hazard is implemented. We assume $\beta$ is the regression coefficients. We assume the first column of the covariate matrix is the treatment indicator, and the corresponding parameter is $\beta_1$. Here $a_0$ is modeled as random with a normalized power prior. The normalized power prior for $\beta$ is approximated by a weighted mixture of multivariate normal distributions provided in `prior.beta.mvn`. The user can use the `approximate.prior.beta` function to obtain samples of $\beta$ from the normalized power prior, and use any mixture of multivariate normals to approximate the normalized power prior for $\beta$. By default, a single multivariate normal distribution is assumed. Baseline hazard parameters for the current and historical data are NOT shared. The baseline hazards of the current data are denoted by $\lambda$. The baseline hazards of the historical data are denoted by $\lambda_0$. We assume Gamma priors for $\lambda$ and $\lambda_0$. To perform sample size determination, we test the hypotheses $$ H_0: \beta_1 \ge \delta $$ and $$ H_1: \beta_1 < \delta. $$ The sampling prior for the treatment parameter can be generated from a normal distribution (see examples). For example, suppose one wants to compute the power for the hypotheses $H_0: \beta_1 \ge 0$ and $H_1: \beta_1 < 0.$ To approximate the sampling prior for $\beta_1$, one can simply sample from a normal distribution with negative mean, so that the mass of the prior falls in the alternative space. Conversely, to compute the type I error rate, one can sample from a normal distribution with positive mean, so that the mass of the prior falls in the null space. The sampling prior for the other parameters ($\beta_2$, ..., $\beta_p$ and $\lambda$) can be generated from the posterior based on the historical data. This can be achieved by the function phm.fixed.a0 with `current.data` set to `FALSE` (see the vignette). Posterior samples are obtained through slice sampling. The default lower limits are -100 for $\beta$ and 0 for $\lambda$. The default upper limits for the parameters are 100. The default slice widths for the parameters are 0.1. The defaults may not be appropriate for all situations, and the user can specify the appropriate limits and slice width for each parameter. If a sampling prior with support in the null space is used, the value returned is a Bayesian type I error rate. If a sampling prior with support in the alternative space is used, the value returned is a Bayesian power. ## Examples ```r # Simulate two historical datasets set.seed(1) n <- 100 P <- 4 time1 <- round(rexp(n, rate=0.5),1) event1 <- rep(1,n) X1 <- matrix(rbinom(n*P,prob=0.5,size=1), ncol=P) S1 <- c(rep(1,n/2),rep(2,n/2)) time2 <- round(rexp(n, rate=0.7),1) event2 <- rep(1,n) X2 <- matrix(rbinom(n*P,prob=0.5,size=1), ncol=P) S2 <- c(rep(1,n/2),rep(2,n/2)) historical <- list(list(time=time1, event=event1, X=X1, S=S1), list(time=time2, event=event2, X=X2, S=S2)) n.subjects <- 100 n.events <- 30 # We choose three intervals for the first stratum and two intervals for the second stratum n.intervals <- c(3,2) change.points <- list(c(1,2),1) # Generate sampling priors # The null hypothesis here is H0: beta_1 >= 0. To calculate power, # we can provide samples of beta_1 such that the mass of beta_1 < 0. # To calculate type I error, we can provide samples of beta_1 such that # the mass of beta_1 >= 0. samp.prior.beta1 <- rnorm(100, mean=-1, sd=1) # Here, mass is put on the alternative region, so power is calculated. samp.prior.beta <- cbind(samp.prior.beta1, matrix(rnorm(100*(P-1)), 100, P-1)) # Point mass sampling priors are used for lambda lambda_strat1 <- matrix(c(0.5, 0.5, 0.5), nrow=1) lambda_strat2 <- matrix(c(0.7, 0.7), nrow=1) samp.prior.lambda <- list(lambda_strat1, lambda_strat2) nMC <- 50 # nMC should be larger in practice nBI <- 50 N <- 5 # N should be larger in practice result <- power.phm.random.a0(historical=historical, n.subjects=n.subjects, n.events=n.events, n.intervals=n.intervals, change.points=change.points, samp.prior.beta=samp.prior.beta, samp.prior.lambda=samp.prior.lambda, prior.a0.shape1 = c(1,1), prior.a0.shape2 = c(1,1), dist.enroll="Uniform", param.enroll=0.5, nMC=nMC, nBI=nBI, delta=0, nullspace.ineq=">", N=N) result$`power/type I error` result$`average posterior mean of beta` result$`average posterior mean of lambda` ``` ## References Ibrahim, J. G., Chen, M.-H. and Sinha, D. (2001). Bayesian Survival Analysis. New York: Springer Science & Business Media. Psioda, M. A. and Ibrahim, J. G. (2019). Bayesian clinical trial design using historical data that inform the treatment effect. Biostatistics 20, 400–415. Shen, Y., Psioda, M. A., and Joseph, J. G. (2023). BayesPPD: an R package for Bayesian sample size determination using the power and normalized power prior for generalized linear models. The R Journal, 14(4). ## See Also `phm.random.a0` and `approximate.prior.beta`

power.phm.random.a0 function