This is an implementation of the power calculation formula
derived by Schmoor et al. (2000) for
the following Cox proportional hazards regression in the epidemiological
studies:
$$h(t|x_1, x_2)=h_0(t)\exp(\beta_1 x_1+\beta_2 x_2 + \gamma (x_1 x_2)),$$
where both covariates $X_1$ and $X_2$ are binary variables.Suppose we want to check if
the hazard ratio of the interaction effect $X_1 X_2=1$ to $X_1 X_2=0$ is equal to $1$
or is equal to $\exp(\gamma)=\theta$.
Given the type I error rate $\alpha$ for a two-sided test, the power
required to detect a hazard ratio as small as $\exp(\gamma)=\theta$ is
$$power=\Phi\left(-z_{1-\alpha/2}+\sqrt{\frac{n}{G}[\log(\theta)]^2 p (1-p) \psi (1-\rho^2)}\right),$$
where $\psi$ is the proportion of subjects died of
the disease of interest, and
$$\rho=corr(X_1, X_2)=(p_1-p_0)\times\sqrt{\frac{q(1-q)}{p(1-p)}},$$
and
$p=Pr(X_1=1)$, $q=Pr(X_2=1)$, $p_0=Pr(X_1=1|X_2=0)$,
and $p_1=Pr(X_1=1 | X_2=1)$, and
$$G=\frac{[(1-q)(1-p_0)p_0+q(1-p_1)p_1]^2}{(1-q)q (1-p_0)p_0 (1-p_1) p_1},$$
and
$p0=Pr(X_1=1 | X_2=0)=myc/(mya+myc)$,
$p1=Pr(X_1=1 | X_2=1)=myd/(myb+myd)$,
$p=Pr(X_1=1)=(myc+myd)/n_{obs}$,
$q=Pr(X_2=1)=(myb+myd)/n_{obs}$,
$n_{obs}=mya+myb+myc+myd$.
$p_{00}=Pr(X_1=0,\mbox{and}, X_2=0)$,
$p_{01}=Pr(X_1=0,\mbox{and}, X_2=1)$,
$p_{10}=Pr(X_1=1,\mbox{and}, X_2=0)$,
$p_{11}=Pr(X_1=1,\mbox{and}, X_2=1)$.