EasyWeibullParms()
An easy way for mere humans to specify the Weibull(a,b) distribution as parameterized in base R/stats by dweibull(), pweibull(), etc. Comments within this Rfunc's code summarize the theory. The final example (#5) deals with a Monte Carlo study to perform a sample-size analysis of a Cox-model survival analysis with right-censoring.
EasyWeibullParms() is included within these Rfuncs: none yet released.
An easy way for mere humans to specify the Weibull(a,b) distribution as parameterized in base R/stats by dweibull(), pweibull(), etc. Comments within this Rfunc's code summarize the theory. The final example (#5) deals with a Monte Carlo study to perform a sample-size analysis of a Cox-model survival analysis with right-censoring.
- Current version released 09 January 2017.
- Download EasyWeibullParms.Rfunc170109a.R
EasyWeibullParms() is included within these Rfuncs: none yet released.
Arguments
Median=NA
Median=NA
- median of X, paired with RS.P=, the P*100% relative spread.
- The P*100% Relative Spread is the ratio of the upper versus lower limits of the middle P*100% of the distribution of X (defined by P=). Used only with Median=. This "Median+RS.P" method is isomorphic to using the mean (central tendency) and SD (spread) to set the Normal distribution ("Mean+SD").
- Q=Q.P or Q=c(Q.P1, Q.P2), one or two quantiles of X (defined by P=).
- P= defines RS.P or Q.P (depending on specifications for other arguments). A single value (P=0.xx) is required if using RS.P; P = 0.xx defines Q=Q.P. P=c(P1, P2) defines Q=c(Q.P1, Q.P2).
- = c(HR.1, HR.2, ..., HR.r). If supplied, determines r Weibull(a,b) distributions that differ according to the proportional hazards model of Cox regression. See Examples 3, 4, and 5.
- The exponential distribution is an oft-used special case of the Weibull in which the hazard is constant. Exponential=TRUE sets a to 1.0. b is set using Median=, Q=, Mean=, or Rate=.
- Used only with Exponential=TRUE
- Used only with Exponential=TRUE
Object Returned
- A data.frame (unnamed) with one row per hazard ratio and three columns: HazRatio, a, and b. See examples.
Examples
These examples are extracted from the "RunExamples" block of code in EasyWeibullParms(). Each example was verified using large N Monte Carlo simulations. See the Rfunc.
Example 1: Exponential(Median=30).
These examples are extracted from the "RunExamples" block of code in EasyWeibullParms(). Each example was verified using large N Monte Carlo simulations. See the Rfunc.
Example 1: Exponential(Median=30).
> (Ex1 <- EasyWeibullParms(Exponential=TRUE, Median=30))
HazRatio a b
1 NA 1 43.280855
Note: The Rfuncs file shows how the same problem is handled using the Q=, Mean=, and Rate= methods.
Example 2: Weibull(Median=30, RS.90=8)
> (Ex2 <- EasyWeibullParms(Median=30, RS.P=8,P=0.90))
HazRatio a b
1 NA 1.955998 36.18253
> # The following gives the same a and b.
> (Q.90 <- qweibull(0.90,Ex2$a,Ex2$b))
[1] 55.42184
> EasyWeibullParms(Q=c(30,Q.90),P=c(0.50,0.90))
HazRatio a b
1 NA 1.955998 36.18253
Example 3: Base Weibull(Q=c(15,55), P=c(0.05,0.95)) & HazRatio=c(0.50,1,2,3)
> (Ex3 <- EasyWeibullParms(Q=c(15,55), P=c(0.05,0.95),
+ HazRatio=c(0.50,1,2,3)))
HazRatio a b
1 0.5 3.130484 48.34034
2 1.0 3.130484 38.73904
3 2.0 3.130484 31.04474
4 3.0 3.130484 27.27329
> # Verifying through a large-N Monte Carlo simulation.
> # The true hazard ratio of group 3 vs, group 2, is HR.4/HR.3 = 3/2 = 1.5
> group4vs3 <- rep(c(0,1),each=100000)
> set.seed(123455)
> x <- c(rweibull(100000,Ex3[3,"a"],Ex3[3,"b"]),
+ rweibull(100000,Ex3[4,"a"],Ex3[4,"b"]))
> tapply(x,group4vs3+3,median)
3 4
27.63146 24.25857
> # Verifying proportional hazards theory and coding by using cph() in rms package
> # to perform Cox proportional hazard regression.
> rmsloaded <- require(rms)
> if (!rmsloaded) {
+ stop("The rms package needs to be installed on this computer.")
+ }
> cat("\nThe rms package has been loaded.")
The rms package has been loaded.
> options(datadist="dd")
> dd <- datadist(group4vs3)
> x.surv <- Surv(x)
> CPHfit = cph(x.surv ~ group4vs3)
> (HR.est <- exp(CPHfit$coefficients)) # compare to true HR = 3/2 = 1.5
group4vs3
1.506613
Example 4. Two-group problem with right censoring. Random censoring is relatively minor. All observations are deterministically censored at 90.
- time.to.event ~ Base Weibull(Q.50=30, RS.90=8), HazRatio = c(1,0.60).
- time.to.censor = min(X~exponential(median=600), 90).
- time = min(Time-to-event, Time-to-censor)
- censored = TRUE, if time.to.event > time.to.censor
= FALSE, otherwise - n = 250 per group
> (Ex4.event <- EasyWeibullParms(Median=30, P=0.95, RS.P=8,
+ HazRatio=c(1.00, 0.60)))
HazRatio a b
1 1.0 2.395629 34.95948
2 0.6 2.395629 43.26838
> (Ex4.censor <- EasyWeibullParms(Exponential=TRUE,Median=600))
HazRatio a b
1 NA 1 865.617
> n <- 250
> group2 <- rep(c(0,1),each=n)
> set.seed(17563)
> time.to.event <- round(c(rweibull(n,Ex4.event[1,"a"],Ex4.event[1,"b"]),
+ rweibull(n,Ex4.event[2,"a"],Ex4.event[2,"b"])),2)
> time.to.censor <- pmin(round(c(rweibull(2*n,1,Ex4.censor[1,"b"]))),90)
> time <- pmin(time.to.event, time.to.censor)
> censored <- (time.to.censor < time.to.event)
> data.frame(group2,time.to.event,time.to.censor,time,censored)[
+ c(1:5,n:(n+5),2*n),]
group2 time.to.event time.to.censor time censored
1 0 31.89 90 31.89 FALSE
2 0 18.62 90 18.62 FALSE
3 0 23.80 36 23.80 FALSE
4 0 65.12 90 65.12 FALSE
5 0 9.92 90 9.92 FALSE
250 0 29.08 90 29.08 FALSE
251 1 43.00 90 43.00 FALSE
252 1 25.47 90 25.47 FALSE
253 1 37.53 90 37.53 FALSE
254 1 28.79 17 17.00 TRUE
255 1 29.53 90 29.53 FALSE
500 1 23.72 90 23.72 FALSE
> sum(censored)/n # proportion of cases censored
[1] 0.108
> rmsloaded <- require(rms)
> if (!rmsloaded) {
+ stop("The rms package needs to be installed on this computer.")
+ }
> cat("\nThe rms package has been loaded.")
The rms package has been loaded.
> options(datadist="dd")
> dd <- datadist(group2)
> TimeSurv <- Surv(time, censored==0)
> CPHfit = cph(TimeSurv ~ group2)
> (HR.est <- exp(CPHfit$coefficients)) # compare to true (population) value: 0.60
group2
0.5548066
> (HR.CI95 <- exp(confint(CPHfit)))
2.5 % 97.5 %
group2 0.4596883 0.6696067
Example 5. Continuing Example 4, conduct a sample-size analysis to address the "crucial" power question, What is the chance that the upper 95% confidence limit for the hazard ratio is less than 0.80? In other words, what is the chance that a study with N = 250+250 (and given the scenarios) will support concluding that the hazard rate for group 2 is at least 20% less than for group 1?
> n <- 250
> Ntrials=5000
> set.seed(1631)
> group2 <- rep(c(0,1),each=n)
> HR.est <- HR.CI95.upper <- numeric(Ntrials)
>
> for (i in 1:Ntrials) {
+ time.to.event <- round(c(rweibull(n,Ex4.event[1,"a"],Ex4.event[1,"b"]),
+ rweibull(n,Ex4.event[2,"a"],Ex4.event[2,"b"])),2)
+ time.to.censor <- pmin(round(c(rweibull(2*n,1,Ex4.censor[1,"b"]))),90)
+ time <- pmin(time.to.event, time.to.censor)
+ censored <- (time.to.censor < time.to.event)
+ TimeSurv <- Surv(time, censored==0)
+ CPHfit = cph(TimeSurv ~ group2)
+ HR.est[i] <- exp(CPHfit$coefficients) # True (population) value: 0.60
+ HR.CI95.upper[i] <- exp(confint(CPHfit))[2]
+ {
+ if ((i==1) & (Ntrials>500)) {
+ cat(paste("\n\nRunning ",Ntrials," Monte Carlo trials: ",sep=""))
+ }
+ if (trunc(i/100) == i/100) { cat(".") }
+ if (trunc(i/500) == i/500) { cat(i) }
+ }
+ }
Running 5000 Monte Carlo trials: .....500.....1000.....1500.....2000.....2500
.....3000.....3500.....4000.....4500.....5000
> mean(HR.est)
[1] 0.6020945
> # estimated traditional power: % with upper 95% conf. limit < 1.00
> sum(HR.CI95.upper < 1.00)/Ntrials
[1] 1
> # estimated "crucial" power: % with upper 95% conf. limit < 0.80
> sum(HR.CI95.upper < 0.80)/Ntrials
[1] 0.87