##' @include semiParametric.R survivalModels.R
NULL
##' Method to calculate model restricted means
##' @name calcModelRmst
##' @rdname calcModelRmst-methods
##' @param object (SurvivalModel) A survival model - note there cannot be
##' any covariates and armAsFactor must be FALSE
##' @param ... additional arguments for this generic
##' @return (data.frame or FlexTable)
##' @export
setGeneric("calcModelRmst", function(object, ...){
standardGeneric("calcModelRmst")
})
##' @name calcModelRmst
##' @aliases calcModelRmst,SurvivalModel-method
##' @rdname calcModelRmst-methods
##' @param model (character) The name of the model for which to calculate the restricted mean
##' @param times (nuermic vector) times to calculate the restricted mean
##' @param class ('data.frame' or "FlexTable' (default)) type of output required
##' @param digits (numeric default=3) if outputting a FlexTable then the number of digits
##' to round the entries to
##' @export
setMethod("calcModelRmst", "SurvivalModel",
function(object, model, times, class=c("data.frame","FlexTable")[2], digits=3, ...){
#validation
if(object@armAsFactor){
stop("Cannot calculate restricted means if armAsFactor is TRUE")
}
if(length(object@covariates)!=0){
stop("Cannot calculate restricted means if covariates fitted in model")
}
if(any(!is.numeric(times) | times < 0)){
stop("Times must be numeric and non-negative")
}
if(length(class) != 1 || !class %in% c("data.frame","FlexTable")){
stop("Invalid class argument, should be 'data.frame' or 'FlexTable")
}
if(length(digits)!=1 || !is.numeric(digits) || !digits > 0 || is.infinite(digits) ||
is.na(digits)){
stop("Invalid digits argument")
}
if(length(model)!=1 || !model %in% names(object@models)){
stop("Invalid model argument must be one of ",
paste(names(object@models),collapse=", "))
}
#for each arm
rmsts <- lapply(object@models[[model]],function(oneModel){
#get the cdf function
tempF <- oneModel$dfns$p
#if spline need to add the knots argument for the dfns$p function to work
args <- list()
if(!is.null(oneModel$knots)){
args$knots <- oneModel$knots
}
#survival function
survFn <- function(x){
1 - do.call("tempF", c(args, list(q=x), oneModel$res[,"est"]))
}
#calculate restricted means (an optimization would be to
#not handle times indpendently but calculate [0,t1], [t1, t2], ... and
#sum them up as needed)
vapply(times, function(time){
tryCatch(
integrate(survFn, 0, time)$value,
error=function(cond) NA
)
},
FUN.VALUE = numeric(1))
})
rmsts <- as.data.frame(do.call("rbind",rmsts))
colnames(rmsts) <- NULL
#if two arms add a difference row
if(nrow(rmsts)==2){
rmsts <- rbind(rmsts,difference=as.numeric(rmsts[2,])-as.numeric(rmsts[1,]) )
}
#Add row of times
rmsts <- rbind(time=times, rmsts)
if(class=="data.frame") return(rmsts)
#create FlexTable
numRows <- nrow(rmsts)
numCols <- 1+ncol(rmsts)
MyFTable <- MyFTable <- FlexTable(numrow=numRows,numcol=numCols,
body.par.props=parProperties(text.align="right"),
header.text.props = textProperties(font.weight = "bold"),
body.cell.props = cellProperties(padding.right=1))
#Set borders
MyFTable[1:numRows,1:numCols,side='bottom'] <- borderProperties(width=0)
MyFTable[1:numRows,1:numCols,side='left'] <- borderProperties(width=0)
MyFTable[1:numRows,1:numCols,side='top'] <- borderProperties(width=0)
MyFTable[1:numRows,1:numCols,side='right'] <- borderProperties(width=0)
MyFTable[numRows,1:numCols,side='bottom'] <- borderProperties(width=3)
MyFTable[1,1:numCols,side='top'] <- borderProperties(width=3)
MyFTable[2,1:numCols,side='top'] <- borderProperties(width=3)
#Add in data to table
MyFTable[2:numRows,2:numCols] <- round(rmsts[2:numRows,], digits=digits)
MyFTable[1,2:numCols] <- times
MyFTable[1:numRows, 1] <- rownames(rmsts)
#Add header denoting which distribution
hR <- FlexRow(paste(getDistributionDisplayNames(model),"\nrestricted means"),
colspan = numCols,
par.properties=parProperties(text.align="center",padding=1),
cell.properties = cellProperties(border.width = 0),
text.properties = textProperties(font.weight = "bold"))
MyFTable <- addHeaderRow(MyFTable,hR)
MyFTable
}
)
##' Method to calculate RMST on subset of data contained in a
##' SemiParametricModel object
##' @name calcRmst
##' @rdname calcRmst-methods
##' @param object (SemiParametricModel) The object which was created when
##' fitting the Cox model
##' @param ... additional parameters needed for specific instances of this
##' generic
##' @return (rmst object) contains list of RMST values, differences and call or
##' a FlexTable for output into a word document (depending on the class variable)
##' @export
setGeneric("calcRmst", function(object, ...){
standardGeneric("calcRmst")
})
##' @name calcRmst
##' @aliases calcRmst,SemiParametricModel-method
##' @rdname calcRmst-methods
##' @param class ('rmst' or "FlexTable' (default)) type of output required
##' @param digits (numeric default=3) if outputting a FlexTable then the number of digits
##' to round the entries to
##' @param pval_digits (numeric default = same value as digits argument) decimal place rounding of p value
##' @export
setMethod("calcRmst", "SemiParametricModel", function(object,class=c("rmst","FlexTable")[2], digits=3, pval_digits=digits, ...){
if(length(class) != 1 || !class %in% c("rmst","FlexTable")){
stop("Invalid class argument, should be 'rmst' or 'FlexTable")
}
if(isSingleArm(object)){
stop("Cannot calculate rmst for a single arm trial")
}
# Create formula for Kaplan-Meier estimator
formulaToFit <- survivalFormula(armAsFactor=!isSingleArm(object),
covariates=character(0),
timeCol = object@endPointDef[["timeCol"]],
censorCol = object@endPointDef[["censorCol"]])
# Call RMST
result <- rmst(formula = formulaToFit,
data = object@survData@subject.data, ...)
if(class=="rmst"){
return(result)
}
#create FlexTable (note rmst function only works with two arms so size of table is fixed)
numRows <- 3
numCols <- 6
MyFTable <- FlexTable(numrow=numRows,numcol=numCols,
body.par.props=parProperties(text.align="right"),
header.text.props = textProperties(font.weight = "bold"),
body.cell.props = cellProperties(padding.right = 1))
#Check if the estimate is less than 0 and if so reverse the signs and switch the lower and upper CI's to show the
#absolute difference
if(result$diff[,1] < 0){
storedDiff <- result$diff
result$diff[,1] = -storedDiff[,1]
result$diff[,3] = -storedDiff[,4]
result$diff[,4] = -storedDiff[,3]
}
#Add data
MyFTable[3,2:5] <- round(result$diff[1,1:4], digits)
#P value rounding
pval <- result$diff[1,5]
p_digits <- as.numeric(pval_digits)
MyFTable[3,6] <- if(pval >= 10^(-p_digits)) round(pval, p_digits) else paste(" < 0." , paste(rep(0, p_digits -1),collapse=""), sep="","1")
MyFTable[1:2,2:5] <- round(result$RMST[,1:4],digits)
#Add 1st column (the arm names)
MyFTable[1:3,1] <- c(as.character(getArmNames(object@survData)),"Difference")
MyFTable[1:numRows,1] <- parProperties(text.align="left")
MyFTable[1:numRows,1] <- textProperties(font.weight = "bold")
#Add header
hR <- FlexRow(c("Arm","RMST", "SE","Lower CI", "Upper CI", "p-value"),
par.properties=parProperties(text.align="left"),
cell.properties =cellProperties(padding.right = 1),
text.properties = textProperties(font.weight = "bold"))
MyFTable <- addHeaderRow(MyFTable,hR)
MyFTable
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.