Nothing
# utility functions that are not part of the base functions but serve miscellaneous purposes
# allow dependencies that don't require installation at initial installation of qgcomp
.qgc.require <- function (package, message = paste("loading required package (",
package, ") failed", sep = "")){
if (!requireNamespace(package, quietly = FALSE)) {
stop(message, call. = FALSE)
}
invisible(TRUE)
}
.rmvnorm <- function(n,mu,Sigma){
# draw from multivariate normal distribution - this is a thin
# wrapper for either mvtnorm::rmvnorm or MASS::mvrnorm,
# depending on the moods of those package developers
#.qgc.require("mvtnorm")
#draw = mvtnorm::rmvnorm(n, mu, Sigma)
.qgc.require("MASS")
draw = MASS::mvrnorm(n, mu, Sigma, tol = 1e-6, empirical = FALSE, EISPACK = FALSE)
draw
}
construction <- function(i="", j=""){
msg ="This function is not yet fully implemented"
if(j != "") msg = paste0(msg, ": ", j)
if(i %in% c("msg", "message", "warning", "wrn", "warn")) warning(msg)
else stop(msg)
}
# fake cox family function
cox <- function(){
obj = binomial(link="log")
obj$family="cox"
obj
}
# fake zi family function
zi <- function(){
obj = binomial(link="log")
obj$family="zi"
obj
}
# fake multinomial family function
multinom_family <- function(){
obj = binomial(link="log")
obj$family="multinomial"
obj
}
se_comb <- function(expnms, covmat, grad=NULL){
#' @title Calculate standard error of weighted linear combination of random variables
#' @description This function uses the Delta method to calculate standard errors of linear
#' functions of variables (similar to `lincom` in Stata). Generally, users will not need to
#' call this function directly.
#' @details This function takes inputs of a set of exposure names (character vector)
#' and a covariance matrix (with colnames/rownames that contain the full set
#' of exposure names), as well as a possible `grad` parameter to calculate
#' the variance of a weighted combination of the exposures in `expnms`, where the
#' weights are based off of `grad` (which defaults to 1, so that this function
#' yields the variance of a sum of all variables in `expnms`)
#'
#' Here is simple version of the delta method for a linear combination of
#' three model coefficients:
#'
#' \eqn{f(\beta) = \beta_1 + \beta_2 + \beta_3}
#' given gradient vector
#' \deqn{G = [d(f(\beta))/d\beta_1 = 1,
#' d(f(\beta))/d\beta_2 = 1,
#' d(f(\beta))/d\beta_3 = 1]}
#' \eqn{t(G) Cov(\beta) G} = delta method variance, where t() is the transpose operator
#'
#' @param expnms a character vector with the names of the columns to be
#' of interest in the covariance matrix for a which a standard error will be
#' calculated (e.g. same as expnms in qgcomp fit)
#' @param covmat covariance matrix for parameters, e.g. from a model or
#' bootstrap procedure
#' @param grad the "weight" vector for calculating the contribution of each variable
#' in expnms to the final standard error. For a linear combination, this is equal
#' to a vector of ones (and is set automatically). Or can be calculated via the
#' grad.poly procedure, in the case of coming up with proper weights when the combination
#' of expnms derives from a polynomial function (as in qgcomp.glm.boot with degree>1).
#'
#' @export
#' @examples
#' vcov = rbind(c(1.2, .9),c(.9, 2.0))
#' colnames(vcov) <- rownames(vcov) <- expnms <- c("x1", "x2")
#' se_comb(expnms, vcov, c(1, 0))^2 # returns the given variance
#' se_comb(expnms, vcov, c(1, 1)) # default linear MSM fit: all exposures
#' # have equal weight
#' se_comb(expnms, vcov, c(.3, .1)) # used when one exposure contributes
#' # to the overall fit more than others = d(msmeffect)/dx
if(!is.matrix(covmat)) {
nm <- names(covmat)
covmat = matrix(covmat)
colnames(covmat) <- nm
}
weightvec <- rep(0, dim(covmat)[1])
#if(is.null(grad)) weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- 1
#if(!is.null(grad)) weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- grad
if (is.null(grad)){
weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- 1
} else if (!is.null(grad) && length(grad)==1){
weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- grad
} else if (!is.null(grad) && length(grad)==length(weightvec)){
weightvec <- grad
}
var <- weightvec %*% covmat %*% weightvec
#var <- sum(wcovmat)
sqrt(var)[1,,drop=TRUE] # should be a scalar
}
vc_comb <- function(aname="(Intercept)", expnms, covmat, grad=NULL){
#' @title Calculate covariance matrix between one random variable and a linear combination of
#' random variables
#' @description This function uses the Delta method to calculate a covariance matrix of linear
#' functions of variables and is used internally in qgcomp. Generally, users will not need to
#' call this function directly.
#' @details This function takes inputs of a name of random variable (character), as
#' set of exposure names (character vector) and a covariance matrix (with colnames/rownames
#' that contain the indepdendent variable and the full set
#' of exposure names). See \code{\link[qgcomp]{se_comb}} for details on variances of sums
#' of random variables. Briefly, for variables A, B and C with covariance matrix Cov(A,B,C),
#' we can calculate the covariance Cov(A,B+C) with the formula Cov(A,B) + Cov(A,C), and
#' Cov(A,B+C+D) = Cov(A,B) + Cov(A,C) + Cov(A,D), and so on.
#'
#' @param aname character scalar with the name of the first column of interest (e.g. variable
#' A in the examples given in the details section)
#' @param expnms a character vector with the names of the columns to be
#' of interest in the covariance matrix for a which a standard error will be
#' calculated (e.g. same as expnms in qgcomp fit)
#' @param covmat covariance matrix for parameters, e.g. from a model or
#' bootstrap procedure
#' @param grad not yet used
#'
#' @return A covariance matrix
#
#' @export
#' @examples
#' vcov = rbind(c(0.010051348, -0.0039332248, -0.0036965571),
#' c(-0.003933225, 0.0051807876, 0.0007706792),
#' c(-0.003696557, 0.0007706792, 0.0050996587))
#' colnames(vcov) <- rownames(vcov) <- c("(Intercept)", "x1", "x2")
#' expnms <- rownames(vcov)[2:3]
#' aname = rownames(vcov)[1]
#' vc_comb(aname, expnms, vcov) # returns the given covariance matrix
if(!is.matrix(covmat)) {
nm <- names(covmat)
covmat = matrix(covmat)
colnames(covmat) <- nm
}
weightvec <- rep(0, dim(covmat)[1])
# eventual extension: allow non-unity 'weights' such that the intervention
# could correspond to 1 unit increases in some variables, and < 1 unit increases
# in others
#if(!is.null(grad)) weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- grad
#if(!is.null(grad)) grad = NULL # not yet used
#if(is.null(grad)) weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- 1
if (is.null(grad)){
weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- 1
} else if (!is.null(grad) && length(grad)==1){
weightvec[which(colnames(as.matrix(covmat)) %in% expnms)] <- grad # need more testing with non-unity weights
} else if (!is.null(grad) && length(grad)==length(weightvec)){
weightvec <- grad
}
outcov = matrix(NA, nrow=2, ncol=2)
acol = which(colnames(as.matrix(covmat)) %in% aname)
bcol = which(colnames(as.matrix(covmat)) %in% expnms)
outcov[1,1] <- covmat[acol,acol]
outcov[1,2] <- outcov[2,1] <- sum(covmat[acol, bcol])
outcov[2,2] <- weightvec %*% covmat %*% weightvec # delta method
outcov
}
grad.poly <- function(intvals, degree){
# returns matrix with each column referring
if(degree==1){
mat <- matrix(1, nrow=length(intvals), ncol=1)
}else{
mat <- matrix(1, nrow=length(intvals), ncol=degree)
for(d in 2:degree){
mat[,d] <- d*poly(intvals, degree = d-1, simple = TRUE, raw = TRUE)[,d-1]
}
}
mat
}
quantize <- function (data, expnms, q=4, breaks=NULL) {
#' @title Quantizing exposure data
#' @description Create variables representing indicator functions with cutpoints defined
#' by quantiles. Output a list that includes: 1) a dataset that is a copy of data,
#' except that the variables whose names are included in the `expnms` variable are
#' transformed to their quantized version and 2) an unnamed list of the quantile cutpoints
#' that are used for each of the variables that were quantized
#'
#' @details This function creates categorical variables in place of the
#' exposure variables named in 'expnms.' For example, a continuous exposure
#' 'x1' will be replaced in the output data by another 'x1' that takes on values
#' 0:(q-1), where, for example, the value 1 indicates that the original x1 value
#' falls between the first and the second quantile.
#' @return A list containing the following fields
#' \describe{
#' \item{data}{a quantized version of the original dataframe}
#' \item{breaks}{a list of the quantile cutpoints used to create the quantized variables which
#' includes a very small number for the minimum and a very large number for the maximum to avoid
#' causing issues when using these breaks to quantize new data.}
#' }
#' @param data a data frame
#' @param expnms a character vector with the names of the columns to be
#' quantized
#' @param q integer, number of quantiles used in creating quantized variables
#' @param breaks (optional) list of (equal length) numeric vectors that
#' characterize the minimum value of each category for which to
#' break up the variables named in expnms. This is an alternative to using 'q'
#' to define cutpoints.
#' @concept variance mixtures
#' @import stats
#' @export
#' @examples
#' set.seed(1232)
#' dat = data.frame(y=runif(100), x1=runif(100), x2=runif(100), z=runif(100))
#' qdata = quantize(data=dat, expnms=c("x1", "x2"), q=4)
#' table(qdata$data$x1)
#' table(qdata$data$x2)
#' summary(dat[c("y", "z")]);summary(qdata$data[c("y", "z")]) # not touched
#' dat = data.frame(y=runif(100), x1=runif(100), x2=runif(100), z=runif(100))
#' # using 'breaks' requires specifying min and max (the qth quantile)
#' # example with theoretical quartiles (could be other relevant values)
#' qdata2 = quantize(data=dat, expnms=c("x1", "x2"),
#' breaks=list(c(-1e64, .25, .5, .75, 1e64),
#' c(-1e64, .25, .5, .75, 1e64)
#' ))
#' table(qdata2$data$x1)
#' table(qdata2$data$x2)
e <- new.env()
e$retbr <- list()
qt <- function(i){
# not exported
datmat <- as.numeric(unlist(data[, expnms[i]]))
if(!is.null(breaks)){
# prioritize breaks if given by user
br <- breaks[[i]]
e$retbr[[i]] <- breaks[[i]]
}else{
br <- unique(quantile(datmat, probs = seq(0, 1, by = 1 / q), na.rm = TRUE))
br[1] <- -1e64
br[length(br)] <- 1e64
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(seq_len(length(expnms)), qt)
data[, expnms] <- vapply(seq_len(length(expnms)), qt, rep(0.0, nrow(data)))
}
return(list(data=data, breaks=e$retbr))
}
checknames <- function(terms){
#' @title Check for valid model terms in a qgcomp fit
#' @description This is an internal function called by \code{\link[qgcomp]{qgcomp}},
#' \code{\link[qgcomp]{qgcomp.glm.boot}}, and \code{\link[qgcomp]{qgcomp.glm.noboot}},
#' but is documented here for clarity. Generally, users will not need to call
#' this function directly. This function tries to determine whether there are
#' non-linear terms in the underlying model, which helps infer whether the
#' appropriate function is called, and whether more explicit function calls
#' are needed.
#' @param terms model terms from attr(terms(modelfunction, data), "term.labels")
nonlin <- ifelse(sum(grep("\\(|\\:|\\^", terms))>0, TRUE, FALSE)
if(nonlin){
return(FALSE)
}else{
return(TRUE)
}
}
.qgcomp_object <- function(...){
res = list(...)
nms = names(res)
if(is.na(match("hasintercept", nms))) res$hasintercept = TRUE
if(is.na(match("bootstrap", nms))) res$bootstrap=FALSE
if(is.na(match("cov.yhat", nms))) res$cov.yhat=NULL
if(is.na(match("degree", nms))) res$degree=1
if(is.na(match("pos.psi", nms))) res$pos.psi = NULL
if(is.na(match("neg.psi", nms))) res$neg.psi = NULL
if(is.na(match("pos.weights", nms))) res$pos.weights = NULL
if(is.na(match("neg.weights", nms))) res$neg.weights = NULL
if(is.na(match("pos.size", nms))) res$pos.size = NULL
if(is.na(match("neg.size", nms))) res$neg.size = NULL
if(is.na(match("df", nms))) res$df = NULL
if(is.na(match("qx", nms))) res$qx = NULL
attr(res, "class") <- c("qgcompfit", "list")
res
}
.qgcompmult_object <- function(...){
res = .qgcomp_object(...)
attr(res, "class") <- c("qgcompmultfit", attr(res, "class"))
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.