#!/bin/sh
# This is a shell archive (produced by shar 3.49)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#
# made 02/17/1995 15:06 UTC by feh@biostat
# Source directory /usr/local/slibrary/local
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 983 -rw-r--r-- cpower.README
# 2909 -rw-r--r-- cpower.s
# 2566 -rw-r--r-- ciapower.s
# 3655 -rw-r--r-- .Data/.Help/cpower
# 2039 -rw-r--r-- .Data/.Help/ciapower
#
# ============= cpower.README ==============
if test -f 'cpower.README' -a X"$1" != X"-c"; then
echo 'x - skipping cpower.README (File already exists)'
else
echo 'x - extracting cpower.README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cpower.README' &&
cpower - Compute power of two-sample Cox or logrank survival test
ciapower- Compute power of interaction test for exponential survival
X
X Created by Frank Harrell (feh@biostat.mc.duke.edu)
X
To install the programs and help files, place this file in your local library
directory $SHOME/library/local. (Type !echo $SHOME while under S if you
need to find out what $SHOME is). Then type sh fff, where fff is the name
you stored this shar file under. Use "sh fff -c" to overwrite existing
files. The sh command will create these files:
cpower.README, cpower.s, ciapower.s, and .Data/.Help/cpower,ciapower. Under S
or S-Plus while in the /local area, source the cpower.s and ciapower.s
files.
X
COPYRIGHT NOTICE
X
You may distribute these functions freely as long as you do so without
financial gain, that you acknowledge the source, and that you
communicate improvements to the author.
X
The author is willing to help with problems. Send E-mail to
feh@biostat.mc.duke.edu.
X
SHAR_EOF
chmod 0644 cpower.README ||
echo 'restore of cpower.README failed'
Wc_c="`wc -c < 'cpower.README'`"
test 983 -eq "$Wc_c" ||
echo 'cpower.README: original size 983, current size' "$Wc_c"
fi
# ============= cpower.s ==============
if test -f 'cpower.s' -a X"$1" != X"-c"; then
echo 'x - skipping cpower.s (File already exists)'
else
echo 'x - extracting cpower.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cpower.s' &&
#tref time at which mortalities estimated
#n total sample size
#mc tref-year mortality, control
#r % reduction in m1c by intervention
#accrual duration of accrual period
#tmin minimum follow-up time
#noncomp.c % non-compliant in control group (drop-ins)
#noncomp.i % non-compliant in intervention group (non-adherers)
#alpha type I error
#nc Sample size for control (if not n/2)
#ni Sample size for intervention (if not n/2)
#pr set to T to print intermediate results
#
#non-compliance handled by an approximation of Eq. 5.4 of
#Lachin JM, Foulkes MA (1986): Evaluation of sample size and power for
#analyses of survival with allowance for nonuniform patient entry,
#losses to follow-up, noncompliance, and stratification.
#Here we're using log hazard ratio instead of their hazard difference
X
cpower <- function(tref,
X n,
X mc,
X r,
X accrual,
X tmin,
X noncomp.c=0,
X noncomp.i=0,
X alpha=.05,
X nc, ni,
X pr=T) {
X
if(mc>1) stop("mc should be a fraction")
#Find mortality in intervention group
mi <- (1-r/100)*mc
X
if(missing(nc) | missing(ni)) {nc <- n/2; ni <- n/2} else n <- nc+ni
X
if(pr) {
X cat("\nAccrual duration:",accrual,"y Minimum follow-up:",tmin,"y\n")
X cat("\nTotal sample size:",n,"\n")
X cat("\nAlpha=",alpha,"\n")
X d <- c("Control","Intervention")
X m <- c(mc,mi)
X names(m) <- d
X cat("\n",tref,"-year Mortalities\n",sep=""); print(m)
}
X
#Find exponential hazards for all groups
lamc <- -log(1-mc)/tref
lami <- -log(1-mi)/tref
X
if(pr) {
X lam <- c(lamc,lami)
X names(lam) <- d
X cat("\nHazard Rates\n"); print(lam)
}
X
#Find probability that a subject will have her event observed during
#the study, for all groups
tmax <- tmin+accrual
pc <- 1-1/accrual/lamc*(exp(-tmin*lamc)-exp(-tmax*lamc))
pi <- 1-1/accrual/lami*(exp(-tmin*lami)-exp(-tmax*lami))
X
if(pr) {
X p <- c(pc,pi)
X names(p) <- d
X cat("\nProbabilities of an Event During Study\n")
X print(p)
}
X
#Find expected number of events, all groups
mc <- pc*nc
mi <- pi*ni
X
if(pr) {
X m <- c(mc,mi)
X names(m) <- d
X cat("\nExpected Number of Events\n")
X print(round(m,1))
}
X
#Find expected value of observed log hazard ratio
delta <- log(lami/lamc)
if(pr) cat("\nHazard ratio:",format(exp(delta)),"\n")
X
if(noncomp.c+noncomp.i>0) {
X if(pr) cat("\nDrop-in rate (controls):",noncomp.c,
X "%\nNon-adherence rate (intervention):",noncomp.i,"%\n",sep="")
X delta <- delta * (1 - (noncomp.c+noncomp.i)/100)
X if(pr) cat("Effective hazard ratio with non-compliance:",
X format(exp(delta)),"\n")
}
X
#Find its variance
v <- 1/mc + 1/mi
#Get same as /sasmacro/samsizc.sas if use 4/(mc+mi)
X
sd <- sqrt(v)
if(pr) cat("Standard deviation of log hazard ratio:",format(sd),"\n")
X
z <- -qnorm(alpha/2)
X
1 - ( pnorm(z - abs(delta)/sd) - pnorm(-z - abs(delta)/sd) )
}
X
X
SHAR_EOF
chmod 0644 cpower.s ||
echo 'restore of cpower.s failed'
Wc_c="`wc -c < 'cpower.s'`"
test 2909 -eq "$Wc_c" ||
echo 'cpower.s: original size 2909, current size' "$Wc_c"
fi
# ============= ciapower.s ==============
if test -f 'ciapower.s' -a X"$1" != X"-c"; then
echo 'x - skipping ciapower.s (File already exists)'
else
echo 'x - extracting ciapower.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ciapower.s' &&
#tref time at which mortalities estimated
#n1 total sample size, stratum 1
#n2 total sample size, stratum 2
#m1c tref-year mortality, stratum 1 control
#m2c " " 2 "
#r1 % reduction in m1c by intervention, stratum 1
#r2 % reduction in m2c by intervention, stratum 2
#accrual duration of accrual period
#tmin minimum follow-up time
#alpha type I error
#pr set to T to print intermediate results
X
ciapower <- function(tref,
X n1,
X n2,
X m1c,
X m2c,
X r1,
X r2,
X accrual,
X tmin,
X alpha=.05,
X pr=T) {
X
#Find mortality in intervention groups
if(m1c>1 | m2c>1) stop("m1c and m2c must be fractions")
m1i <- (1-r1/100)*m1c
m2i <- (1-r2/100)*m2c
X
if(pr) {
X cat("\nAccrual duration:",accrual,"y Minimum follow-up:",tmin,"y\n")
X cat("\nSample size Stratum 1:",n1," Stratum 2:",n2,"\n")
X cat("\nAlpha=",alpha,"\n")
X d <- list(c("Stratum 1","Stratum 2"), c("Control","Intervention"))
X m <- cbind(c(m1c,m2c),c(m1i,m2i))
X dimnames(m) <- d
X cat("\n",tref,"-year Mortalities\n",sep=""); print(m)
}
X
#Find exponential hazards for all groups
lam1c <- -log(1-m1c)/tref
lam2c <- -log(1-m2c)/tref
lam1i <- -log(1-m1i)/tref
lam2i <- -log(1-m2i)/tref
X
if(pr) {
X lam <- cbind(c(lam1c,lam2c),c(lam1i,lam2i))
X dimnames(lam) <- d
X cat("\nHazard Rates\n"); print(lam)
}
X
#Find probability that a subject will have her event observed during
#the study, for all groups
tmax <- tmin+accrual
p1c <- 1-1/accrual/lam1c*(exp(-tmin*lam1c)-exp(-tmax*lam1c))
p2c <- 1-1/accrual/lam2c*(exp(-tmin*lam2c)-exp(-tmax*lam2c))
p1i <- 1-1/accrual/lam1i*(exp(-tmin*lam1i)-exp(-tmax*lam1i))
p2i <- 1-1/accrual/lam2i*(exp(-tmin*lam2i)-exp(-tmax*lam2i))
X
if(pr) {
X p <- cbind(c(p1c,p2c), c(p1i,p2i))
X dimnames(p) <- d
X cat("\nProbabilities of an Event During Study\n")
X print(p)
}
X
#Find expected number of events, all groups
m1c <- p1c*n1/2
m2c <- p2c*n2/2
m1i <- p1i*n1/2
m2i <- p2i*n2/2
X
if(pr) {
X m <- cbind(c(m1c,m2c), c(m1i,m2i))
X dimnames(m) <- d
X cat("\nExpected Number of Events\n")
X print(round(m,1))
}
X
#Find expected value of observed log hazard ratio
delta <- log((lam1i/lam1c)/(lam2i/lam2c))
if(pr) cat("\nRatio of hazard ratios:",format(exp(delta)),"\n")
X
#Find its variance
v <- 1/m1c + 1/m2c + 1/m1i + 1/m2i
sd <- sqrt(v)
if(pr) cat("Standard deviation of log ratio of ratios:",format(sd),"\n")
X
z <- -qnorm(alpha/2)
#if(pr) cat("\nCritical value:",format(z),"\n")
X
1 - ( pnorm(z - abs(delta)/sd) - pnorm(-z - abs(delta)/sd) )
}
X
X
SHAR_EOF
chmod 0644 ciapower.s ||
echo 'restore of ciapower.s failed'
Wc_c="`wc -c < 'ciapower.s'`"
test 2566 -eq "$Wc_c" ||
echo 'ciapower.s: original size 2566, current size' "$Wc_c"
fi
# ============= .Data/.Help/cpower ==============
if test ! -d '.Data'; then
echo 'x - creating directory .Data'
mkdir '.Data'
fi
if test ! -d '.Data/.Help'; then
echo 'x - creating directory .Data/.Help'
mkdir '.Data/.Help'
fi
if test -f '.Data/.Help/cpower' -a X"$1" != X"-c"; then
echo 'x - skipping .Data/.Help/cpower (File already exists)'
else
echo 'x - extracting .Data/.Help/cpower (Text)'
sed 's/^X//' << 'SHAR_EOF' > '.Data/.Help/cpower' &&
.tr @@
.BG
.FN cpower
.TL
Power of Cox/log-rank Two-Sample Test
.DN
Assumes exponential distributions for both treatment groups.
Uses the George-Desu method along with
formulas of Schoenfeld that allow estimation of the expected number of
events in the two groups.
To allow for drop-ins (noncompliance to control therapy, crossover to
intervention) and noncompliance of the intervention, the method of
Lachin and Foulkes is used.
.CS
cpower(tref, n, mc, r, accrual, tmin, noncomp.c=0, noncomp.i=0,
X alpha=0.05, nc, ni, pr=T)
.RA
.AG tref
time at which mortalities estimated
.AG n
total sample size (both groups combined). If allocation is unequal
so that there are not `n/2' observations in each group, you may specify
the sample sizes in `nc' and `ni'.
.AG mc
tref-year mortality, control
.AG r
% reduction in `mc' by intervention
.AG accrual
duration of accrual period
.AG tmin
minimum follow-up time
.OA
.AG noncomp.c
% non-compliant in control group (drop-ins)
.AG noncomp.i
% non-compliant in intervention group (non-adherers)
.AG alpha
type I error probability. A 2-tailed test is assumed.
.AG nc
number of subjects in control group
.AG ni
number of subjects in intervention group. `nc' and `ni' are specified
exclusive of `n'.
.AG pr
set to `F' to suppress printing of details
.RT
power
.SE
prints
.DT
For handling noncompliance, uses a modification of formula (5.4) of
Lachin and Foulkes. Their method is based on a test for the difference
in two hazard rates, whereas `cpower' is based on testing the difference
in two log hazards. It is assumed here that the same correction factor
can be approximately applied to the log hazard ratio as Lachin and Foulkes applied to
the hazard difference.
.sp \\
Note that Schoenfeld approximates the variance
of the log hazard ratio by `4/m', where `m' is the total number of events,
whereas the George-Desu method uses the slightly better `1/m1 + 1/m2'.
Power from this function will thus differ slightly from that obtained with
the SAS `samsizc' program.
.SH AUTHOR
Frank Harrell
.sp 0
Division of Biometry
.sp 0
Duke University Medical Center
.sp 0
feh@biostat.mc.duke.edu
.SH REFERENCES
Peterson B, George SL: Controlled Clinical Trials 14:511-522; 1993.
.sp \\
Lachin JM, Foulkes MA: Biometrics 42:507-519; 1986.
.sp \\
Schoenfeld D: Biometrics 39:499-503; 1983.
.SA
`ciapower', `bpower'
.EX
#In this example, 4 plots are drawn on one page, one plot for each
#combination of noncompliance percentage. Within a plot, the
#5-year mortality % in the control group is on the x-axis, and
#separate curves are drawn for several % reductions in mortality
#with the intervention. The accrual period is 1.5y, with all
#patients followed at least 5y and some 6.5y.
X
par(mfrow=c(2,2),oma=c(3,0,3,0))
X
morts <- seq(10,25,length=50)
red <- c(10,15,20,25)
X
for(noncomp in c(0,10,15,-1)) {
X if(noncomp>=0) nc.i <- nc.c <- noncomp else {nc.i <- 25; nc.c <- 15}
X z <- paste("Drop-in ",nc.c,"%, Non-adherence ",nc.i,"%",sep="")
X plot(0,0,xlim=range(morts),ylim=c(0,1),
X xlab="5-year Mortality in Control Patients (%)",
X ylab="Power",type="n")
X title(z)
X cat(z,"\n")
X lty <- 0
X for(r in red) {
X lty <- lty+1
X power <- morts
X i <- 0
X for(m in morts) {
X i <- i+1
X power[i] <- cpower(5, 14000, m/100, r, 1.5, 5, nc.c, nc.i, pr=F)
X }
X lines(morts, power, lty=lty)
X }
X if(noncomp==0)legend(18,.55,rev(paste(red,"% reduction",sep="")),
X lty=4:1,bty="n")
}
mtitle("Power vs Non-Adherence for Main Comparison",
X ll="alpha=.05, 2-tailed, Total N=14000",cex.l=.8)
.KW sample size
.KW design
.KW survival
.WR
SHAR_EOF
chmod 0644 .Data/.Help/cpower ||
echo 'restore of .Data/.Help/cpower failed'
Wc_c="`wc -c < '.Data/.Help/cpower'`"
test 3655 -eq "$Wc_c" ||
echo '.Data/.Help/cpower: original size 3655, current size' "$Wc_c"
fi
# ============= .Data/.Help/ciapower ==============
if test -f '.Data/.Help/ciapower' -a X"$1" != X"-c"; then
echo 'x - skipping .Data/.Help/ciapower (File already exists)'
else
echo 'x - extracting .Data/.Help/ciapower (Text)'
sed 's/^X//' << 'SHAR_EOF' > '.Data/.Help/ciapower' &&
.tr @@
.BG
.FN ciapower
.TL
Power of Interaction Test for Exponential Survival
.DN
Uses the method of Peterson and George to compute the power of an
interaction test in a 2 x 2 setup in which all 4 distributions are
exponential. This will be the same as the power of the Cox model
test if assumptions hold. The test is 2-tailed.
The duration of accrual is specifed
(constant accrual is assumed), as is the minimum follow-up time.
The maximum follow-up time is then `accrual + tmin'. Treatment
allocation is assumed to be 1:1.
.CS
ciapower(tref, n1, n2, m1c, m2c, r1, r2, accrual, tmin,
X alpha=0.05, pr=T)
.RA
.AG tref
time at which mortalities estimated
.AG n1
total sample size, stratum 1
.AG n2
total sample size, stratum 2
.AG m1c
tref-year mortality, stratum 1 control
.AG m2c
tref-year mortality, stratum 2 control
.AG r1
% reduction in `m1c' by intervention, stratum 1
.AG r2
% reduction in `m2c' by intervention, stratum 2
.AG accrual
duration of accrual period
.AG tmin
minimum follow-up time
.OA
.AG alpha
type I error probability
.AG pr
set to `F' to suppress printing of details
.RT
power
.SE
prints
.SH AUTHOR
Frank Harrell
.sp 0
Division of Biometry
.sp 0
Duke University Medical Center
.sp 0
feh@biostat.mc.duke.edu
.SH REFERENCES
Peterson B, George SL: Controlled Clinical Trials 14:511-522; 1993.
.SA
`cpower'
.EX
# Find the power of a race x treatment test. 25% of patients will
# be non-white and the total sample size is 14000.
# Accrual is for 1.5 years and minimum follow-up is 5y.
# Reduction in 5-year mortality is 15% for whites, 0% or -5% for
# non-whites. 5-year mortality for control subjects if assumed to
# be 0.18 for whites, 0.23 for non-whites.
n <- 14000
for(nonwhite.reduction in c(0,-5)) {
X cat("\n\n\n% Reduction in 5-year mortality for non-whites:",
X nonwhite.reduction, "\n\n")
X pow <- ciapower(5, .75*n, .25*n, .18, .23, 15, nonwhite.reduction,
X 1.5, 5)
X cat("\n\nPower:",format(pow),"\n")
}
.KW power
.KW sample size
.KW survival
.KW study design
.WR
SHAR_EOF
chmod 0644 .Data/.Help/ciapower ||
echo 'restore of .Data/.Help/ciapower failed'
Wc_c="`wc -c < '.Data/.Help/ciapower'`"
test 2039 -eq "$Wc_c" ||
echo '.Data/.Help/ciapower: original size 2039, current size' "$Wc_c"
fi
exit 0