R/independenceTest.R

Defines functions print.independenceTest independenceTest

Documented in independenceTest

#' Non-parametric independence test for bivariate survival data
#'
#' @title Non-parametric independence test for bivariate survival data
#' @param formula a formula object, with the response on the left of a ~
#'          operator, and the terms on the right.  The response must be a
#'          survival object as returned by the \code{Surv} function. The RHS must contain a 'cluster' term
#' @param data a data.frame containing the variables in the model
#' @param weight weight function for test. Default is 'independence'
#' which is optimal for the Frank copula. Other options are 'dabrowska', 'atRisk', a weight matrix
#' or a number (in which case the weight is constant).
#' @return test statistic, SE and p-value for independence test.
#' @seealso dabrowska
#' @references Shih, Joanna H., and Thomas A. Louis. "Tests of independence for bivariate survival data." Biometrics (1996): 1440-1449.
#' @examples
#' library(survival)
#' data("diabetic")
#' independenceTest(Surv(time,status)~cluster(id), data = diabetic)
#' @import stats
#' @importFrom survival Surv
#' @importFrom Rfast sort_unique
#' @export
#' @author Jeppe E. H. Madsen <jeppe.ekstrand.halkjaer@gmail.com>
independenceTest <- function(formula, data = NULL, weight = "independence"){
    Call <- match.call()
    d <- uniTrans(formula, data)
    if(ncol(d) != 4) 
        stop("RHS needs a 'cluster(id)' element")
    x <- d$x; y <- d$y; xstatus <- d$xstatus; ystatus <- d$ystatus
    eyy <- eyyfunc(x,y,sort_unique(x),sort_unique(y))
    nu1 <- nrow(eyy)
    nu2 <- ncol(eyy)
    m1 <- survfit(survival::Surv(x,xstatus)~1)
    m2 <- survfit(survival::Surv(y,ystatus)~1)
    if(class(weight) == "matrix"){}
    else if(weight == "independence"){
        weight <- m1$surv %o% m2$surv
    }
    else if(weight == "dabrowska"){
        weight <- dabrowska(formula, data)$surv
    }
    else if(weight == "atRisk"){
        weight <- eyy
    }
    else if(length(weight)==1){
        weight <- matrix(weight, nrow = nu1, ncol = nu2)
    }
    else{
        stop("'weight' has to be either 'independence', 'dabrowska', 'atRisk', a number or a matrix")
    }
    cumhaz1 <- m1$cumhaz
    cumhaz2 <- m2$cumhaz
    n <- length(x)
    M1 <- martin(x,xstatus,m1$time,cumhaz1,n,nu1)
    M2 <- martin(y,ystatus,m2$time,cumhaz2,n,nu2)
    dM1 <- rbind(M1[1,],diff(M1))
    dM2 <- rbind(M2[1,],diff(M2))
    testStatistic <- 0
    for(i in 1:n){
       testStatistic <- testStatistic + sum((dM1[,i] %o% dM2[,i]) * weight)
    }
    ## Variance estimate
    help <- c(cumhaz1[1],diff(cumhaz1)) %o% c(cumhaz2[1],diff(cumhaz2))
    varianceEst <- sum(weight^2 * eyy * help) * n
    ## Present wrapping
    test <- testStatistic / sqrt(varianceEst)
    out <- list(call = Call, est = testStatistic, se = sqrt(varianceEst), zval = test,
                p = 2 * pnorm(-abs(test)))
    class(out) <- "independenceTest"
    out
}

#' @export
print.independenceTest <- function(x, digits = max(3L, getOption("digits") - 3L), symbolic.cor = x$symbolic.cor, 
    signif.stars = getOption("show.signif.stars"), ...) 
{
    cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"), 
        "\n\n", sep = "")
    coefs <- cbind(Estimate = x$est, `Std. Error` = x$se, 
                   `z value` = x$zval, `Pr(>|z|)` = x$p)
    rownames(coefs) <- "Test"
    printCoefmat(coefs, digits = digits, signif.stars = signif.stars, 
                 na.print = "NA", ...)
    cat("\n")
    invisible(x)
}

#' @export
summary.independenceTest <- function(object, digits = max(3L, getOption("digits") - 3L), symbolic.cor = object$symbolic.cor, 
    signif.stars = getOption("show.signif.stars"), ...) 
{
    cat("\nCall:\n", paste(deparse(object$call), sep = "\n", collapse = "\n"), 
        "\n\n", sep = "")
    coefs <- cbind(Estimate = object$est, `Std. Error` = object$se, 
                   `z value` = object$zval, `Pr(>|z|)` = object$p)
    rownames(coefs) <- "Test"
    printCoefmat(coefs, digits = digits, signif.stars = signif.stars, 
                 na.print = "NA", ...)
    cat("\n")
    invisible(object)
}
Jeepen/biSurv documentation built on Sept. 30, 2023, 3:55 a.m.