Use quantile g-computation with quantiles that are derived from weighted data, rather than the analytic sample

weighted_quantiles_in_qgcomp.R

#!/usr/local/env R
# a way to base the quantiles off of weighted data

library(qgcomp)
library(DescTools) # need this to calculate weighted quantiles
data(metals)

# new function to calculate weighted quantiles
quantize.wtd <- function (data, expnms, q = 4, breaks = NULL, weights=NULL) 
{
  e <- new.env()
  e$retbr <- list()
  qt <- function(i) {
    datmat <- as.numeric(unlist(data[, expnms[i]]))
    if (!is.null(breaks)) {
      br <- breaks[[i]]
      e$retbr[[i]] <- breaks[[i]]
    }
    else {
      if(is.null(weights)){
        br <- unique(quantile(datmat, probs = seq(0, 1, by = 1/q), 
                                  na.rm = TRUE))
      }
      if(!is.null(weights)){
        br <- unique(Quantile(datmat, probs = seq(0, 1, by = 1/q),
                              weights = weights,
                                  na.rm = TRUE))
      }
      br[1] <- -1e+64
      br[length(br)] <- 1e+64
      e$retbr[[i]] <- br
    }
    cut(datmat, breaks = br, labels = FALSE, include.lowest = TRUE) - 
      1
  }
  if (length(expnms) == 1) {
    data[, expnms] <- qt(1)
  }
  else {
    data[, expnms] <- sapply(1:length(expnms), qt)
  }
  return(list(data = data, breaks = e$retbr))
}



submetals = metals[,c("disease_time", "disease_state","arsenic", "barium")]
expnms = c("arsenic", "barium")

#find quantiles of weighted dataset (I'm just using randomly assigned weights)
submetals$wt = runif(nrow(submetals))*2

# use the new function to calculate the cutpoints/breaks
wtd.breaks = quantize.wtd(submetals, expnms, q=4, weights = submetals$wt)$breaks
# compare with standard breaks from qgcomp
unwtd.breaks = quantize(submetals, expnms, q=4)$breaks


# use the new cutpoints in either qgcomp function
qgcomp.cox.noboot(Surv(disease_time, disease_state) ~ arsenic + barium, data = submetals, breaks=wtd.breaks)
qgcomp.cox.boot(Surv(disease_time, disease_state) ~ arsenic + barium, data = submetals, breaks=unwtd.breaks, B=2)

#default
qgcomp.cox.noboot(Surv(disease_time, disease_state) ~ arsenic + barium, data = submetals, q=4)