# R/rmean.R In mbrueckner/AHR: Estimation and Testing of Average Hazard Ratios

#### Documented in print.rmdrmeanDiffrmeanDiff.ahr

#' Estimate difference of restricted mean survival based on (weighted) Kaplan-Meier estimates
#' of the survival functions in each group.
#'
#' @title rmeanDiff
#' @param L time-limit specifying up to which time restricted mean will be calculated
#' @param formula an object of class '"formula"' specifying the conditional survival model
#' @param data data frame containing the variables in formula
#' @param rr.subset logical vector defining subset of observations to use for response rate estimation (default: use all observations)
#' @return An object of class '"rmd"', i.e. a list containing:
#'  \item{L}{time limit, i.e. restricted mean up to time L is calculated}
#'  \item{rmean1}{restricted mean in group 1}
#'  \item{rmean2}{restricted mean in group 2}
#'  \item{rmean.diff}{estimated restricted mean difference}
#'  \item{var.rmean1}{an estimate of the asymptotic variance of the restricted mean in group 1}
#'  \item{var.rmean2}{an estimate of the asymptotic variance of the restricted mean in group 2}
#'  \item{var.rmean.diff}{an estimate of the asymptotic variance of the restricted mean difference}
#'  \item{Z.rmean}{the standardized test statistic for testing rmean.diff=0}
#'  \item{p.value}{p-value corresponding to Z.rmean}
#' @export
#' @examples
#' T <- c(rexp(100, 1), rexp(100, 2))
#' C <- c(rexp(100, 1), rexp(100, 2))
#' Y <- pmin(T, C)
#' D <- T <= C
#' Z <- rep(c(0,1), c(100, 100))
#' fit <- rmeanDiff(2, formula=Surv(Y, D) ~ Z, data.frame(Y=Y, D=D, Z=Z))
rmeanDiff <- function(L, formula, data, rr.subset=rep(TRUE, nrow(data))) {
if(!is.null(formula)) data <- parseFormula(formula, data)

grps <- levels(data$Trt) if(nlevels(data$Trt) != 2) stop("Difference in restricted means can only be estimated in two-sample setting")

times <- getTimes(L, data)

data1 <- data[data$Trt == grps[1], ] data2 <- data[data$Trt == grps[2], ]

trt.sub <- data$Trt[rr.subset] rr.subset1 <- rr.subset[trt.sub == grps[1]] rr.subset2 <- rr.subset[trt.sub == grps[2]] n <- c(nrow(data1), nrow(data2)) param1 <- list(alpha=1, var=TRUE, cov=TRUE, left.limit=FALSE, rr.subset=rr.subset1) param2 <- list(alpha=1, var=TRUE, cov=TRUE, left.limit=FALSE, rr.subset=rr.subset2) fit1 <- wkm(times, data1, param1) fit2 <- wkm(times, data2, param2) rmeanDiffFit(L, times, fit1, fit2, n, n/sum(n), FALSE) } #' Estimate difference of restricted mean survival (based on ahr object as returned by ahr) #' #' This function is usefull if the function 'ahr' has already been called, since the survival estimates #' in the object returned by 'ahr' can be reused. #' #' @title rmeanDiff.ahr #' @param ahr.obj object of class '"ahr"' #' @return An object of class '"rmd"', i.e. a list containing: #' \item{L}{time limit, i.e. restricted mean up to time L is calculated} #' \item{rmean1}{restricted mean in group 1} #' \item{rmean2}{restricted mean in group 2} #' \item{rmean.diff}{estimated restricted mean difference} #' \item{var.rmean1}{an estimate of the asymptotic variance of the restricted mean in group 1} #' \item{var.rmean2}{an estimate of the asymptotic variance of the restricted mean in group 2} #' \item{var.rmean.diff}{an estimate of the asymptotic variance of the restricted mean difference} #' \item{Z.rmean}{the standardized test statistic for testing rmean.diff=0} #' \item{p.value}{p-value corresponding to Z.rmean} #' @export #' @examples #' T <- c(rexp(100, 1), rexp(100, 2)) #' C <- c(rexp(100, 1), rexp(100, 2)) #' Y <- pmin(T, C) #' D <- T <= C #' Z <- rep(c(0,1), c(100, 100)) #' fit <- avgHR(2, data.frame(Y=Y, D=D, Z=Z), formula=Surv(Y, D) ~ Z) #' rmd <- rmeanDiff.ahr(fit) rmeanDiff.ahr <- function(ahr.obj) { if(class(ahr.obj) != "ahr") stop("Input must be of class 'ahr'") if(nlevels(ahr.obj$groups) != 2) stop("Difference in restricted means can only be estimated in two-sample setting")

rmeanDiffFit(ahr.obj$L, ahr.obj$times, ahr.obj$surv.fit[[1]], ahr.obj$surv.fit[[2]], ahr.obj$n, ahr.obj$p, ahr.obj$log.iis) } rmeanDiffFit <- function(L, times, fit1, fit2, n, p, log.iis) { snt <- 1:length(times) ## restricted means rmean1 <- stepIntegrate(fit1$S, times)
rmean2 <- stepIntegrate(fit2$S, times) rmean.diff <- rmean1 - rmean2 if((!is.null(fit1$COV) && !is.null(fit2$COV)) || log.iis) { if(log.iis) { v1 <- fit1$V / p[1]
v2 <- fit2$V / p[2] ##tmp <- simplify2array(lapply(snt, function(l) stepIntegrate(v1[pmin(snt, l)] + v2[pmin(snt, l)], times))) tmp1 <- simplify2array(lapply(snt, function(l) stepIntegrate(v1[pmin(snt, l)], times))) tmp2 <- simplify2array(lapply(snt, function(l) stepIntegrate(v2[pmin(snt, l)], times))) } else { rho1 <- fit1$COV / p[1]
rho2 <- fit2$COV / p[2] ##tmp <- simplify2array(lapply(snt, function(l) stepIntegrate(rho1[,l] + rho2[,l], times))) tmp1 <- simplify2array(lapply(snt, function(l) stepIntegrate(rho1[,l], times))) tmp2 <- simplify2array(lapply(snt, function(l) stepIntegrate(rho2[,l], times))) } var.rmean1 <- stepIntegrate(tmp1, times) / sum(n) ## / n1 var.rmean2 <- stepIntegrate(tmp2, times) / sum(n) ## / n2 var.rmean.diff <- var.rmean1 + var.rmean2 ## standardized difference Z <- rmean.diff / sqrt(var.rmean.diff) res <- list(rmean1=rmean1, rmean2=rmean2, rmean.diff=rmean.diff, var.rmean.diff=var.rmean.diff, var.rmean1=var.rmean1, var.rmean2=var.rmean2, Z.rmean=Z, p.value=1-pnorm(Z), L=L) } else res <- list(rmean.diff=rmean.diff, L=L) class(res) <- "rmd" res } #' Print rmd object #' #' @title print.rmd #' @param x an object of class '"rmd"'. #' @param digits minimal number of significant digits. #' @param ... further arguments passed to or from other methods. #' @method print rmd #' @export print.rmd <- function(x, digits=3, ...) { cat("\n") cat("\t Restricted mean difference", paste0("(L=", x$L, ")"), "\n")
cat("\n\n")

cat(paste0("Restricted mean in group 1: ", round(x$rmean1, digits), " (Std.dev: ", round(sqrt(x$var.rmean1), digits), ")\n"))
cat(paste0("Restricted mean in group 2: ", round(x$rmean2, digits), " (Std.dev: ", round(sqrt(x$var.rmean2), digits), ")"))

cat("\n\n")

cat(paste0("Restricted mean difference: ", round(x$rmean.diff, digits), " (Std.dev: ", round(sqrt(x$var.rmean.diff), digits), ", Z: ", round(x$Z.rmean, digits), ", p: ", round(x$p.value, digits), ")"))

cat("\n")

}

mbrueckner/AHR documentation built on May 22, 2017, 12:16 a.m.