[f2e496]: / R / study.R

Download this file

245 lines (203 with data), 10.3 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
#The Study class definition (used in predict from parameters)
#together with the exported methods. Also see study_constructors.R
#' @include lag.R common.R ctrlSpec.R sfn.R eventPrediction_package.R
NULL
##' Class defining the Study
##' @slot HR Hazard ratio to be detected
##' @slot alpha Significance level [0,1] (see also two-sided indicator)
##' @slot power Power [0,1]
##' @slot two.sided If TRUE, two sided test will be used (i.e. alpha/2).
##' @slot r Control:Experimental subject balance (1:r), i.e. nE/nC=r. r=1 corresponds to equally
##' many subjects in both arms. 2 means we have twice the number of subjects in the experimental arm.
##' Specifically \code{floor(r*N/(r+1))} subjects are
##' allocated to the experimental arm and all other subjects are allocated to the control arm.
##' @slot N Number of subjects to be recruited (integer)
##' @slot study.duration Number of months the study will be going.
##' @slot ctrlSpec A CtrlSpec object which calculates the control group median. This object will be created automatically
##' when calling a constructor for the Study class.
##' @slot dropout A list of CtrlSpec object which calculates the median drop out rate for the control arm (index 1) and
##' active arm (index 2).
##' This object will be created automatically when calling a constructor for the study class
##' @slot dropout.shape The Weibull shape parameter of the dropout hazard function
##' @slot k non-uniformity of accrual (integer, 1=uniform). Non-uniform accrual is allowed for
##' using the following distribution for the probability of a patient entering the trial at time \eqn{b}
##' within the accrual period \eqn{[0,B]}: \eqn{F(b)=b_k/B_k}; \eqn{f(b)=k b_{k-1}/B_k} where \eqn{k} is the
##' measure of non-uniformity (\eqn{k>0}). \eqn{k=1} indicates uniform accrual. This implies that during
##' the first half of the accrual period, \eqn{1/2^k} of the patients will be recruited. Half of the patients
##' will be recruited by time \eqn{B/2^{1/k}}.
##' @slot acc.period Accrual time in months
##' @slot shape The Weibull shape parameter
##' @slot followup The time a subject is followed after randomization, if Inf then there is no fixed time period
##' @slot type Character: The study type, either "Oncology" or "CRGI"
##' @slot lag.settings The \code{LaggedEffect} object describing any lag effect for the study
##' @export
setClass( "Study",
slots= list( HR = "numeric", # Hazard ratio
alpha = "numeric",
power = "numeric",
two.sided = "logical",
r = "numeric",
N = "numeric", # Patients to be recruited
study.duration = "numeric", # Length of study (In months)
ctrlSpec = "CtrlSpec", # Median of control arm
dropout = "list",
dropout.shape ="numeric",
k = "numeric",
acc.period = "numeric",
shape="numeric",
followup="numeric",
type="character",
lag.settings = "LagEffect"),
validity = function(object){
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
ans <- ""
if(object@shape <= 0) ans <- paste(ans,"Invalid shape.")
if(object@acc.period <= 0) ans <- paste(ans,"Invalid acc.period")
if(object@k <= 0) ans <- paste(ans,"Invalid k.")
if(object@N <= 0 || !is.wholenumber(object@N) || length(object@N)>1 ) ans <- paste(ans,"Invalid N.")
if(object@study.duration <= 0 ) ans <- paste(ans,"Invalid study.duration.")
if(object@r < 0 ) ans <- paste(ans,"Invalid r.")
if(object@study.duration <= object@acc.period)
ans <- paste(ans,"acc.period must be < study.duration")
if(object@r == 0){
if(!is.na(object@HR)) ans <- paste(ans,"HR must be NA if r = 0")
if(!is.na(object@power)) ans <- paste(ans,"power must be NA if r = 0")
if(!is.na(object@alpha)) ans <- paste(ans,"alpha must be NA if r = 0")
}
else{
if(is.na(object@HR) || object@HR >= 1 || object@HR <= 0) ans <- paste(ans,"Invalid HR.")
if(object@alpha >= 1 || object@alpha <= 0 ) ans <- paste(ans,"Invalid alpha.")
if(object@power >= 1 || object@power <= 0 ) ans <- paste(ans,"Invalid power.")
}
if(!object@type %in% c("Oncology","CRGI")) ans <- paste(ans,"Invalid type.")
if(length(object@followup) > 1 || object@followup <= 0 ) ans <- paste(ans,"Invalid followup.")
if(class(object@lag.settings)!="LagEffect")ans <- paste(ans,"Invalid lag.settings")
if(!isNullLag(object@lag.settings)){
if(isSingleArm(object) && !is.na(object@lag.settings@L.HazardRatio)){
ans <- paste(ans,"lag.settings@L.HazardRatio must be as.numeric(NA) if study has one arm")
}
if(!isSingleArm(object) && is.na(object@lag.settings@L.HazardRatio)){
ans <- paste(ans,"lag.settings@L.HazardRatio cannot be NA if study has more than one arm")
}
}
if(!is.infinite(object@followup) && !isNullLag(object@lag.settings)){
ans <- paste(ans, "Cannot use lagged settings with a study which has a finite follow up time")
}
if(ans=="") return(TRUE)
ans
}
)
##' @name show
##' @rdname show-methods
##' @aliases show,Study-method
##' @export
setMethod("show", signature(object="Study"),
function(object) {
cat("Study definition:\n")
cat(paste("Number of Patients (N):",object@N,"\n"))
cat(paste("Study duration:",object@study.duration,"months\n"))
cat(paste("Accrual period:",object@acc.period,"months\n"))
cat(paste("Accrual uniformity (k):",object@k,"\n"))
if(!isSingleArm(object)){
cat(paste("Control arm survival:"),object@ctrlSpec@text,"\n")
cat(paste("Hazard Ratio:",object@HR,"\n"))
cat(paste("Ratio of control to experimental 1:",object@r,"\n",sep=""))
cat(paste("alpha:",object@alpha))
if(object@two.sided){
cat("(two sided)\n")
}
else{
cat("(one sided)\n")
}
cat(paste("Power:",object@power,"\n"))
}
else{
cat(paste("Survival:"),object@ctrlSpec@text,"\n")
cat("Single Arm trial\n")
}
if(object@shape==1){
cat("Exponential survival\n")
}
else{
cat(paste("Weibull survival with shape parameter",object@shape,"\n"))
}
if(!is.infinite(object@followup)){
cat("Subject follow up period:",object@followup,"months\n")
}
if(isSingleArm(object)){
outputdropouttext("Subject",object@dropout[[1]],object@dropout.shape)
}
else{
outputdropouttext("Control Arm",object@dropout[[1]])
outputdropouttext("Active Arm",object@dropout[[2]],object@dropout.shape)
}
show(object@lag.settings)
})
##' Is the \code{Study} a single arm study
##' @rdname isSingleArm-methods
##' @name isSingleArm
##' @param study A \code{Study} object
##' @return TRUE if study is single arm, FALSE otherwise
##' @export
setGeneric("isSingleArm",
def=function(study)
standardGeneric("isSingleArm"))
##' @rdname isSingleArm-methods
##' @aliases isSingleArm,Study-method
##' @name isSingleArm
##' @export
setMethod("isSingleArm",representation(study="Study"),
function(study){
study@r==0
})
##' @name predict
##' @rdname predict-methods
##' @aliases predict,Study-method
##' @param step.size The resolution of the grid to be used to calculate the time v
##' expected events curves
##' @export
setMethod( "predict", representation( object="Study" ),
function( object, time.pred=NULL,
event.pred=NULL,step.size=0.5 ) {
study <- object
lagged <- study@lag.settings
validatePredictArguments(time.pred,event.pred,step.size,study@study.duration)
#times for the plotting function
grid.times <- seq( 0, study@study.duration, step.size )
#Next calculate the values of rate parameters
lambda <- lambda.calc( study@ctrlSpec@median, study@HR, study@shape )
if(isNullLag(lagged)){
lambdaot <- as.numeric(NA)
}
else{
lambdaot <- lambda.calc( lagged@ctrlSpec@median, lagged@L.HazardRatio,study@shape )
}
dropout.lambda <- GetDropoutLambda(study)
#Create the survival functions
sfns <- GetSurvivalFunctions(lambda,lambdaot,lagged@Lag.T,isSingleArm(study),study@shape,
study@followup,study@dropout.shape,dropout.lambda)
#Calculate event details for plot
grid <- CalculateEventsGivenTimes(grid.times,study,sfns,calc.at.risk=FALSE)
#Calculate average HR
av.HR <-if(isNullLag(lagged) || isSingleArm(study)) study@HR
else calculateAverageHR(sfns,study,lagged@Lag.T,lagged@L.HazardRatio,
lambdaot,lambda,study@shape)
#Calculate the critical value
critical.data <- CalculateCriticalValue(study,sfns,grid,av.HR)
#Calculate events at given times
predict.data <- if(!is.null(time.pred)) CalculateEventsGivenTimes(time.pred,study,sfns)
else predict.data <- data.frame()
#Calculate times for given number of events
if(!is.null(event.pred)){
predict.data <- rbind(predict.data,CalculateTimesGivenEvents(event.pred,study,sfns,grid))
}
return(AnalysisResults(
critical.HR = critical.data[["chr"]],
critical.data = critical.data[["critical.data"]],
critical.events.req= critical.data[["critical.events.req"]],
av.hr=av.HR,
grid=grid,
predict.data=predict.data,
study=study,
sfns=sfns))
})