add cluster option to assocCoxPH

smgogarten · Feb 8, 2020 · 8a3fd58 · 8a3fd58
1 parent 45418ce
commit 8a3fd58
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 2 deletions.
diff --git a/R/assocCoxPH.R b/R/assocCoxPH.R
@@ -109,6 +109,7 @@ assocCoxPH <- function(genoData,
                        covar = NULL,
                        ivar = NULL,
                        strata = NULL,
+                       cluster = NULL,
                        scan.exclude = NULL,
                        LRtest = FALSE,
                        effectAllele = c("minor", "alleleA"),
@@ -136,13 +137,14 @@ assocCoxPH <- function(genoData,
 
     ## read in outcome and covariate data
     if (verbose) message("Reading in Phenotype and Covariate Data...")
-    dat <- .modelData(genoData, chr, event, c(time.to.event, covar, strata), ivar)
+    dat <- .modelData(genoData, chr, event, c(time.to.event, covar, strata, cluster), ivar)
     ## identify samples with any missing data
     keep <- keep & complete.cases(dat)
     dat <- dat[keep,,drop=FALSE]
     if (!is.null(ivar)) ivar <- paste0(ivar, ":genotype")
     if (!is.null(strata)) strata <- paste0("strata(", paste(strata, collapse=","), ")")
-    model.string <- paste("surv ~", paste(c(covar, ivar, strata, "genotype"), collapse=" + "))
+    if (!is.null(cluster)) cluster <- paste0("cluster(", paste(cluster, collapse=","), ")")
+    model.string <- paste("surv ~", paste(c(covar, ivar, strata, cluster, "genotype"), collapse=" + "))
 
     ## sample size, assuming no missing genotypes
     n <- nrow(dat)

diff --git a/inst/unitTests/assocCoxPH_test.R b/inst/unitTests/assocCoxPH_test.R
@@ -79,6 +79,21 @@ test_strata <- function() {
     checkTrue(!all(is.na(assoc$Est)))
 }
 
+test_cluster <- function() {
+    event <- "event"
+    time.to.event <- "time.to.event"
+    covar <- "age"
+    cluster <- "site"
+
+    genoData <- .cphGenoData()
+    assoc <- assocCoxPH(genoData,
+                       event,
+                       time.to.event,
+                       covar = covar,
+                       cluster = cluster)
+    checkTrue(!all(is.na(assoc$Est)))
+}
+
 test_GxE <- function() {
     event <- "event"
     time.to.event <- "time.to.event"

diff --git a/man/assocCoxPH.Rd b/man/assocCoxPH.Rd
@@ -11,6 +11,7 @@ assocCoxPH(genoData,
            covar = NULL,
            ivar = NULL,
            strata = NULL,
+           cluster = NULL,
            scan.exclude = NULL,
            LRtest = FALSE,
            effectAllele = c("minor", "alleleA"),
@@ -39,6 +40,7 @@ assocCoxPH(genoData,
   \item{ivar}{the name of the variable in \code{covar} to include as an interaction with genotype}
   \item{strata}{a vector of names of variables to stratify on for a
     stratified analysis}
+  \item{cluster}{the name of a column (in the scan annotation of \code{genoData}) which clusters the observations for robust variance. See \code{\link[survival]{cluster}}}.
   \item{scan.exclude}{a vector of scanIDs for scans to exclude}
   \item{LRtest}{logical for whether to perform Likelihood Ratio Tests in addition to Wald tests (which are always performed).}
   \item{effectAllele}{whether the effects should be returned in terms of the minor allele for the tested sample (\code{effectAllele="minor"}) or the allele returned by \code{getAlleleA(genoData)}  (\code{effectAllele="alleleA"}). If the minor allele is alleleB for a given SNP, the difference between these two options will be a sign change for the beta estimate.}