Susan L. Tucker
The goal of this analysis is to perform power calculations needed to select a method for assigning patients to high-risk versus lower-risk residual disease (RD) groups in our validation study.
In our validation study, we measured FABP4 expression levels by qRT-PCR for 139 tumor samples.
To test our hypothesis that RD occurs more frequently among patients with high FABP4 levels, we will compare RD incidence in subgroups of patients predicted to be at high risk versus lower risk of RD.
Specifically, we will assign a pre-determined number of patients (nHigh) to the high risk group, and the remaining patients (139 - nHigh) to the low risk group. A 1-sided Fisher's exact test will be used to compare the incidence of RD in the two risk groups.
Therefore, the specific goal of this report is to select a value for nHigh, based on power calculations.
The power calculation proceeds as follows.
Using a binomial random-number generator, we randomly generate outcomes (RD versus no RD) based on assumed probabilities of RD in the high-risk and low-risk groups.
The observed (randomly generated) RD rates in the high and low risk groups are compared using the 1-sided Fisher exact test. The test is considered a success if P<0.05.
The two steps described above are repeated for a total of 10,000 simulations, and the proportion of successful tests (the estimated power of the test) is recorded.
In performing the power calculations, we assume that the overall probability of RD in the validation cohort (pRD) is 75%, based on the observed RD rates of 77% (378/491) and 73% (136/186) in the TCGA and Tothill data sets, respectively.
To determine an appropriate value for the probability of RD in the high-risk group (pRDhigh), we compute the positive predictive value (PPV = proportion of patients with RD in the high-risk group) in the TCGA and Tothill data sets as a function of the call rate (the proportion of patients assigned to the high-risk group). PPV will serve as an estimate of pRDhigh.
We note that the probability of RD in the low-risk group (pRDlow) can be computed from pRD and pRDhigh, as described later in this report, so it need not be specified separately in performing the power calculations.
To compute PPV as a function of call rate in the TCGA and Tothill data sets, we use the RData objects containing clinical and gene expression data that were created in previous reports (assembleTCGAClinical, assembleTothillClinical). Patients are filtered as described previously (filterTCGASamples, filterTothillSamples).
The power calculations show that assigning 20%, 25% or 30% of patients to the high-risk group would likely provide very high power (>90%) of achieving success in our validation study provided PPV values are at the high end of the range seen in the TCGA and Tothill data sets.
For PPV values on the low end of the range seen for TCGA and Tothill, call rates of 20% or 30% would likely provide relatively low power (50-60%).
For a call rate of 25%, the calculations suggest that the power will be reasonably high (>76%) for any value of PPV falling within the range observed for TCGA and Tothill.
Based on the power calculations, we elect to assign the 25% of patients with the highest FABP4 values to our predicted high-risk group.
The TCGA and Tothill data objects created previously are loaded.
load(file.path("RDataObjects", "tcgaExpression.RData"))
load(file.path("RDataObjects", "tcgaClinical.RData"))
load(file.path("RDataObjects", "tcgaFilteredSamples.RData"))
load(file.path("RDataObjects", "tothillExpression.RData"))
load(file.path("RDataObjects", "tothillClinical.RData"))
load(file.path("RDataObjects", "tothillFilteredSamples.RData"))
Filtrations are applied to the TCGA data.
tcgaSampleUseLong <- rownames(tcgaFilteredSamples[which(tcgaFilteredSamples[,
"sampleUse"] == "Used"), ])
tcgaSampleUse <- substr(tcgaSampleUseLong, 1, 12)
tcgaRDUse <- tcgaRD[tcgaSampleUse]
tcgaExprUse <- tcgaExpression[, tcgaSampleUseLong]
colnames(tcgaExprUse) <- tcgaSampleUse
Filtrations are applied to the Tothill data.
tothillSampleUse <- rownames(tothillFilteredSamples[which(tothillFilteredSamples[,
"sampleUse"] == "Used"), ])
tothillRDUse <- tothillRD[tothillSampleUse]
tothillExprUse <- tothillExpression[, tothillSampleUse]
Before performing power calculations, we calculate PPV in the TCGA and Tothill data sets as a function of the proportion of patients assigned to the high risk group based on FABP4 expression levels.
We first obtain the vectors of FABP4 expression values for the two data sets.
library(hthgu133a.db)
probeFABP4 <- names(which(unlist(mget(rownames(tcgaExprUse), hthgu133aSYMBOL)) ==
"FABP4"))
probeFABP4
## [1] "203980_at"
tcgaFABP4 <- tcgaExprUse[probeFABP4, ]
tothillFABP4 <- tothillExprUse[probeFABP4, ]
We next sort the vectors by FABP4 expression level.
tcgaFABP4sorted <- tcgaFABP4[order(tcgaFABP4)]
tcgaRDsorted <- tcgaRDUse[order(tcgaFABP4)]
tothillFABP4sorted <- tothillFABP4[order(tothillFABP4)]
tothillRDsorted <- tothillRDUse[order(tothillFABP4)]
PPV is computed and plotted as a function of call rate.
tcgaNPts <- length(tcgaFABP4)
tcgaPctCall <- seq(from = tcgaNPts, to = 1, by = -1)
tcgaPPV <- c()
for (i in 1:tcgaNPts) {
tcgaPPV <- c(tcgaPPV, length(which(tcgaRDsorted[i:tcgaNPts] == "RD")))
}
tcgaPPV <- 100 * tcgaPPV/tcgaPctCall
tcgaPctCall <- 100 * tcgaPctCall/tcgaNPts
tothillNPts <- length(tothillFABP4)
tothillPctCall <- seq(from = tothillNPts, to = 1, by = -1)
tothillPPV <- c()
for (i in 1:tothillNPts) {
tothillPPV <- c(tothillPPV, length(which(tothillRDsorted[i:tothillNPts] ==
"RD")))
}
tothillPPV <- 100 * tothillPPV/tothillPctCall
tothillPCall <- 100 * tothillPctCall/tothillNPts
plot(tcgaPctCall, tcgaPPV, type = "l", xlab = "Patients Assigned to High-risk Group (%)",
ylab = "Positive Predictive Value (%)", ylim = c(70, 100), lwd = 2)
lines(tothillPctCall, tothillPPV, col = c("red"), lwd = 2)
legend(x = 50, y = 75, legend = c("TCGA", "Tothill et al."), lty = c(1, 1),
lwd = 2, col = c("black", "red"))
The following notation is used here:
We observe that pRD = pRDhigh * pCall + pRDlow * (1 - pCall). Therefore pRDlow can be determined once pRD, pCall and pRDhigh are specified.
We begin by defining the function that performs the power calculation.
powerCalc <- function(nSam, nCall, pRD, pRDhigh, sig, nits) {
# nSam = number of samples (to be fixed here at 139) nCall = number called
# with RD sig = significance level defining success of the 1-sided Fisher's
# exact test (to be fixed here at P=0.05) nits = number of iterations in the
# numerical simulation pCall and pRD as defined above
pCall <- nCall/nSam
pRDlow = (pRD - pRDhigh * pCall)/(1 - pCall)
callRD <- c(rep(1, nCall), rep(0, nSam - nCall))
fisherP <- c()
for (i in 1:nits) {
observedRD <- c()
for (j in 1:nCall) {
observedRD <- c(observedRD, rbinom(1, 1, pRDhigh))
}
for (j in 1:(nSam - nCall)) {
observedRD <- c(observedRD, rbinom(1, 1, pRDlow))
}
fisherP <- c(fisherP, fisher.test(callRD, observedRD, or = 1, alternative = "greater",
conf.int = FALSE)$p.value)
}
power <- sum(fisherP < sig)/nits
return(power)
}
We now run the power calculations. For each calculation, we fix the number of samples at nSam = 139, the overall probability of RD at 75%, and the number of iterations at 10,000.
We assume that the number of patients predicted to be in the high-risk group (nCall) is 28, 35 or 42, corresponding to call rates of about 20%, 25% or 30%, respectively.
Based on the PPV plots obtained from the TCGA and Tothill data, we assume that pRDhigh = 90, 92.5 or 95% when the call rate is 20-25%, and we assume that pRDhigh = 85, 90 or 95% when the call rate is 30%. These represent approximately upper and lower limits for the PPV, as well as an intermediate value.
powerEst <- matrix(0, nrow = 9, ncol = 3)
colnames(powerEst) <- c("nCall", "pRDhigh", "Power")
rownames(powerEst) <- rep("", nrow(powerEst))
set.seed(5913)
powerEst[1, ] <- c(28, 0.95, powerCalc(139, 28, 0.75, 0.95, 0.05, 10000))
powerEst[2, ] <- c(28, 0.925, powerCalc(139, 28, 0.75, 0.925, 0.05, 10000))
powerEst[3, ] <- c(28, 0.9, powerCalc(139, 28, 0.75, 0.9, 0.05, 10000))
powerEst[4, ] <- c(35, 0.95, powerCalc(139, 35, 0.75, 0.95, 0.05, 10000))
powerEst[5, ] <- c(35, 0.925, powerCalc(139, 35, 0.75, 0.925, 0.05, 10000))
powerEst[6, ] <- c(35, 0.9, powerCalc(139, 35, 0.75, 0.9, 0.05, 10000))
powerEst[7, ] <- c(42, 0.95, powerCalc(139, 42, 0.75, 0.95, 0.05, 10000))
powerEst[8, ] <- c(42, 0.9, powerCalc(139, 42, 0.75, 0.9, 0.05, 10000))
powerEst[9, ] <- c(42, 0.85, powerCalc(139, 42, 0.75, 0.85, 0.05, 10000))
We list the results of the power calculations.
powerEst
## nCall pRDhigh Power
## 28 0.950 0.9098
## 28 0.925 0.7925
## 28 0.900 0.6224
## 35 0.950 0.9750
## 35 0.925 0.8981
## 35 0.900 0.7622
## 42 0.950 0.9916
## 42 0.900 0.8606
## 42 0.850 0.5015
getwd()
## [1] "/Users/slt/SLT WORKSPACE/EXEMPT/OVARIAN/Ovarian residual disease study 2012/RD manuscript/Web page for paper/Webpage"
sessionInfo()
## R version 3.0.2 (2013-09-25)
## Platform: x86_64-apple-darwin10.8.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] hthgu133a.db_2.9.0 org.Hs.eg.db_2.9.0 RSQLite_0.11.4
## [4] DBI_0.2-7 AnnotationDbi_1.22.6 Biobase_2.20.1
## [7] BiocGenerics_0.6.0 knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.1 formatR_0.9 IRanges_1.18.4 stats4_3.0.2
## [5] stringr_0.6.2 tools_3.0.2