version 0.0.1.0

DESCRIPTION

@@ -0,0 +1,24 @@
+Package: CoxBcv
+Type: Package
+Title: Bias-Corrected Sandwich Variance Estimators for Marginal Cox
+        Analysis of Cluster Randomized Trials
+Version: 0.0.1.0
+Authors@R: c(person("Xueqi","Wang", role = c("aut", "cre"),
+            email = "xueqi.wang@duke.edu"),
+            person("Elizabeth", "Turner", role = "aut",
+            email = "liz.turner@duke.edu"),
+            person("Fan", "Li", role = "aut",
+            email = "fan.f.li@yale.edu"))
+Maintainer: Xueqi Wang <xueqi.wang@duke.edu>
+Description: The implementation of bias-corrected sandwich variance estimators for the analysis of cluster randomized trials with time-to-event outcomes using the marginal Cox model, proposed by Wang et al. (under review).
+License: GPL (>= 2)
+RoxygenNote: 7.1.2
+Encoding: UTF-8
+Imports: pracma, survival
+NeedsCompilation: no
+Packaged: 2022-03-05 21:37:54 UTC; xueqiwang
+Author: Xueqi Wang [aut, cre],
+  Elizabeth Turner [aut],
+  Fan Li [aut]
+Repository: CRAN
+Date/Publication: 2022-03-09 08:10:02 UTC

MD5

@@ -0,0 +1,22 @@
+fe9296859b8b838aa684fb7b2d10ef0d *DESCRIPTION
+e12da9e98fa677978105bbd6ebabd8ec *NAMESPACE
+dd1a218f132c5b4f89effdd66ee2e2e2 *R/CoxBcv.fg.R
+eb277fbc789724883f0b48792064f476 *R/CoxBcv.fgmr.R
+3191e478f853c72b5444c109eb8e20b0 *R/CoxBcv.kc.R
+21d61eec8a2018452cd94a4e55ee2285 *R/CoxBcv.kcmr.R
+f2bb7b13d44e20ad5651bf98f4300624 *R/CoxBcv.mbn.R
+686c82c4cd435d8e90f9951823c54aa2 *R/CoxBcv.mbnmr.R
+f144da0442157fb1102ea78c4e08f6e6 *R/CoxBcv.md.R
+f5d6cb0b0a51dc06f0b57d36e905081a *R/CoxBcv.mdmr.R
+86ac3fda22c3b1d6c0a7e7db8110208b *R/CoxBcv.mr.R
+00baede614be0fa02eb454f31d0c435d *R/CoxBcv.rob.R
+198b3ff4eb6edb59d110a4861a180455 *man/CoxBcv.fg.Rd
+373f0afd532c9437869a175e1aa95eed *man/CoxBcv.fgmr.Rd
+400f8681fceff000894fead6ae5b60e8 *man/CoxBcv.kc.Rd
+6e0c2c192b0e2af8413c314aa35361ca *man/CoxBcv.kcmr.Rd
+519fa5985b860f4bf3b83f2b2988377f *man/CoxBcv.mbn.Rd
+3380c06484b26674377285f21e7730c9 *man/CoxBcv.mbnmr.Rd
+97643d60342eac16ac870e1b32d0c8e8 *man/CoxBcv.md.Rd
+fd57c8219c867fa561a028fa557cd932 *man/CoxBcv.mdmr.Rd
+02e9982cb90c40ca81ad8eccf3b3ee24 *man/CoxBcv.mr.Rd
+47f57e4c0a4c77f975191d48b259c0dd *man/CoxBcv.rob.Rd

NAMESPACE

@@ -0,0 +1,15 @@
+# Generated by roxygen2: do not edit by hand
+
+export(CoxBcv.fg)
+export(CoxBcv.fgmr)
+export(CoxBcv.kc)
+export(CoxBcv.kcmr)
+export(CoxBcv.mbn)
+export(CoxBcv.mbnmr)
+export(CoxBcv.md)
+export(CoxBcv.mdmr)
+export(CoxBcv.mr)
+export(CoxBcv.rob)
+import(pracma)
+import(survival)
+importFrom(stats,coef)

R/CoxBcv.fg.R

@@ -0,0 +1,180 @@
+#' Fay and Graubard (FG) bias-corrected sandwich variance estimator
+#' 
+#' Calculate the Fay and Graubard (FG; 2001) bias-corrected sandwich variance estimator, for marginal Cox analysis of cluster randomized trials, 
+#' proposed by Wang et al. (under review).
+#' 
+#' @param Y vector of observed time-to-event data.
+#' @param Delta vector of censoring indicators.
+#' @param X matrix of marginal mean covariates with one column for one covariate (design matrix excluding intercept).
+#' @param ID vector of cluster identifiers.
+#' 
+#' @return
+#' \itemize{
+#'   \item coef - estimate of coefficients.
+#'   \item exp(coef) - estimate of hazard ratio.
+#'   \item FG-var - FG bias-corrected sandwich variance estimate of coef.
+#' }
+#' 
+#' @export
+#' 
+#' @references 
+#' Fay, M. P., & Graubard, B. I. (2001). 
+#' Small‐sample adjustments for Wald‐type tests using sandwich estimators. 
+#' Biometrics, 57(4), 1198-1206.
+#' 
+#' Wang, X., Turner, E. L., & Li, F. 
+#' Improving sandwich variance estimation for marginal Cox analysis of cluster randomized trials. 
+#' Under Review.
+#' 
+#' @examples 
+#' Y <- c(11,19,43,100,7,100,100,62,52,1,7,6)
+#' Delta <- c(1,1,1,0,1,0,0,1,1,1,1,1)
+#' X1 <- c(0,0,0,0,0,0,1,1,1,1,1,1)
+#' X2 <- c(-19,6,-25,48,10,-25,15,22,17,-9,45,12)
+#' ID <- c(1,1,2,2,3,3,4,4,5,5,6,6)
+#' 
+#' X <- X1
+#' CoxBcv.fg(Y,Delta,X,ID)
+#' 
+#' X <- cbind(X1,X2)
+#' CoxBcv.fg(Y,Delta,X,ID)
+#' 
+#' @importFrom stats coef
+#' @import pracma
+#' @import survival
+
+CoxBcv.fg <- function(Y,Delta,X,ID){
+
+  ######################################
+  # Step 1: prepare data elements
+  ######################################
+  # point estimate
+  test.cox_cluster <- coxph(Surv(Y,Delta)~X+cluster(ID))
+  beta <- as.matrix(coef(test.cox_cluster))
+
+  # sort observations by time
+  b <- order(Y)
+  Y <- sort(Y)
+  X <- as.matrix(X)[b,,drop=FALSE]
+  ID <- ID[b]
+  Delta <- Delta[b]
+  ny <- length(Y)
+  nbeta <- dim(as.matrix(X))[2]
+  UID <- sort(unique(ID))
+  n <- length(UID)
+
+  IDind <- zeros(length(UID), ny)
+  for (i in 1:length(UID)){
+    IDind[i, ID==UID[i]] <- 1
+  }
+
+  ######################################################
+  # Step 2: Model-based variance estimates
+  ######################################################
+  # the rate of counting process of event, or dN(t)
+  NN <- diag(Delta)
+  # use the following trick to obtain the at-risk process Y(t)
+  # each row is an individual
+  # each column is a specific time point (recall the counting process notation)
+  IndYY <- (t(repmat(t(Y),ny,1))>=repmat(t(Y),ny,1))
+  Xbeta <- c(X%*%beta)
+
+  # Three S matrices for variance calculation
+  S0beta <- colSums(IndYY*exp(Xbeta))
+  S1beta <- zeros(nbeta,ny)
+  for(k in 1:nbeta){
+    S1beta[k,] <- colSums(IndYY*exp(Xbeta)*X[,k])
+  }
+  S2beta <- array(0, c(ny,1,nbeta,nbeta))
+  for(k in 1:nbeta){
+    for(s in 1:nbeta){
+      S2beta[,,k,s] <- colSums(IndYY*exp(Xbeta)*X[,k]*X[,s])
+    }
+  }
+
+  Omega <- array(0,c(nbeta,nbeta,n))
+
+  # obtain cluster-specific matrices
+  for(i in UID){
+    # subset observations from each cluster
+    S0beta_c <- S0beta[IDind[i,]==1]
+    S1beta_c <- S1beta[,IDind[i,]==1,drop=FALSE]
+    Delta_c <- Delta[IDind[i,]==1]
+
+    # components for A matrix
+    for (k in 1:nbeta){
+      for (s in 1:nbeta){
+        Omega[k,s,i] <- sum(Delta_c*(S2beta[IDind[i,]==1,,k,s]/S0beta_c-
+                                       S1beta_c[k,]*S1beta_c[s,]/S0beta_c^2))
+      }
+    }
+  }
+
+  Ustar <- apply(Omega,c(1,2),sum)
+  naive <- solve(Ustar)
+
+  ######################################################
+  # Step 3: FG bias correction
+  ######################################################
+  # Breslow estimator of baseline hazard
+  dHY <- colSums(NN)/c(S0beta)
+  HY <- cumsum(dHY)
+
+  # obtain martingale increment: 
+  # recall that the martingale is the "residual" in survival context
+  epsilon <- NN-IndYY*repmat(t(dHY),ny,1)*exp(Xbeta)
+  nom <- nomFG <- zeros(nbeta,n)
+  Omega_m <- array(0,c(nbeta,nbeta,n))
+
+  # obtain cluster-specific matrices
+  for(i in UID){
+    # subset observations from each cluster
+    X_c <- X[IDind[i,]==1,,drop=FALSE]
+    epsilon_c <- epsilon[IDind[i,]==1,,drop=FALSE]
+    epsilon_c_all <- colSums(epsilon_c)
+    S0beta_c <- S0beta[IDind[i,]==1]
+    S1beta_c <- S1beta[,IDind[i,]==1,drop=FALSE]
+    ny_c <- sum(IDind[i,]==1)
+    Delta_c <- Delta[IDind[i,]==1]
+    IndYY_c <- IndYY[IDind[i,]==1,,drop=FALSE]
+    Xbeta_c <- Xbeta[IDind[i,]==1]
+    ylxb_c <- IndYY_c*exp(Xbeta_c)*repmat(dHY,ny_c,1)
+
+    # the trick is to loop through the dimension of the coefficients
+    # otherwise need to deal with multi-dimensional array, very complex
+    for (k in 1:nbeta){
+      # components for B matrix
+      tempk <- repmat(X_c[,k,drop=FALSE],1,ny)-repmat(S1beta[k,,drop=FALSE]/t(S0beta),ny_c,1)
+      nom[k,i] <- sum(as.matrix(tempk*epsilon_c)%*%repmat(1,ny,1))
+
+      for (s in 1:nbeta){
+        # true Omega
+        Omega_m[k,s,i] <- sum(Delta_c*(S2beta[IDind[i,]==1,,k,s]/S0beta_c-S1beta_c[k,]*S1beta_c[s,]/S0beta_c^2)) -
+          sum((repmat(S2beta[,,k,s]/S0beta-S1beta[k,]*S1beta[s,]/S0beta^2,ny_c,1)*ylxb_c)%*%repmat(1,ny,1)) +
+          sum((tempk*repmat(X_c[,s,drop=FALSE],1,ny)*ylxb_c)%*%repmat(1,ny,1))
+      }
+    }
+
+    # components for FG type correction
+    Hi <- zeros(nbeta,nbeta)
+    tempov <- Omega_m[,,i]%*%naive
+    diag(Hi) <- 1/sqrt(1-pmin(0.75,c(diag(tempov))))
+    nomFG[,i] <- Hi%*%nom[,i]
+  }
+
+  # variance estimator of FG type
+  UUFG <- tcrossprod(nomFG)
+  varFG <- naive%*%UUFG%*%naive
+
+  #############################################
+  # Output
+  # varFG: FG bias-corrected sandwich var
+  #############################################
+  bvarFG <- diag(varFG)
+  outbeta <- cbind(matrix(summary(test.cox_cluster)$coefficients[,1:2],ncol=2),bvarFG)
+  colnames(outbeta) <- c("coef","exp(coef)","FG-var")
+
+  return(list(outbeta=outbeta))
+}
+
+