update tests

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: qgcomp
 Title: Quantile G-Computation
 Version: 2.8.6
-Date: 2022-01-23
+Date: 2022-01-24
 Authors@R: 
     person("Alexander", "Keil", , "akeil@unc.edu", role = c("aut", "cre"))
 Author: Alexander Keil [aut, cre]

README.md

@@ -1,5 +1,7 @@
 `qgcomp` v2.8.6
 
+[![R-CMD-check](https://github.com/alexpkeil1/qgcomp/workflows/R-CMD-check/badge.svg)](https://github.com/alexpkeil1/qgcomp/actions) [![](https://www.r-pkg.org:443/badges/version/qgcomp)](https://www.r-pkg.org:443/pkg/qgcomp) [![CRAN RStudio mirror downloads](http://cranlogs.r-pkg.org/badges/qgcomp)](https://www.r-pkg.org:443/pkg/qgcomp) [![codecov](https://codecov.io/gh/alexpkeil1/qgcomp/branch/master/graph/badge.svg)](https://app.codecov.io/gh/alexpkeil1/qgcomp) [![CodeFactor](https://www.codefactor.io/repository/github/alexpkeil1/qgcomp/badge)](https://www.codefactor.io/repository/github/alexpkeil1/qgcomp)
+
 #### QGcomp (quantile g-computation): estimating the effects of exposure mixtures. Works for continuous, binary, and right-censored survival outcomes.
 
 #### Flexible, unconstrained, fast and guided by modern causal inference principles

tests/test_asis.R

@@ -35,12 +35,16 @@ stopifnot(all.equal(
 
 stopifnot(all.equal(
   obj1$fit$coefficients, 
-  obj2$fit$coefficients, 
   obj3$fit$coefficients, check.names=FALSE))
+stopifnot(all.equal(
+  obj1$fit$coefficients, 
+  obj2$fit$coefficients, check.names=FALSE))
 stopifnot(all.equal(
   obj1$msmfit$coefficients, 
-  obj2$msmfit$coefficients, 
   obj3$msmfit$coefficients, check.names=FALSE))
+stopifnot(all.equal(
+  obj1$msmfit$coefficients, 
+  obj2$msmfit$coefficients, check.names=FALSE))
 
 f0b = d ~ x1 + x2
 f4 = d ~ x1*x2

tests/test_l_boot_vs_noboot.R

@@ -0,0 +1,58 @@
+cat("# boot vs no boot test\n")
+library("qgcomp")
+
+dgm <- function(N){
+ dat <- data.frame(id=1:N) 
+ dat <- within(dat, {
+     u = 0
+     x1 = runif(N)*4 + u
+     x2 = runif(N)*4 + u
+     x3 = runif(N)*4 + u
+     x4 = runif(N)*4 + u
+     x5 = runif(N)*4 + u
+     x6 = runif(N)*4 + u
+     y = rnorm(N, x1+x2, 2)
+})
+ dat[,c('y', paste0("x", 1:6))]
+}
+
+Xnm = c(paste0("x", 1:6))
+
+repit <- function(i){
+  dat = dgm(50)
+  m1 = qgcomp.noboot(y~., expnms=Xnm, data = dat, family=gaussian(), q=4)
+  m2 = qgcomp.boot(  y~., expnms=Xnm, data = dat, family=gaussian(), q=4, B=2, parallel=FALSE)
+  res = c(m1$coef, m1$var.coef, 1*(m1$pval>0.05), with(m1, ci.coef[1]<2 & ci.coef[2]>2), m2$coef, m2$var.coef, 1*(m2$pval>0.05), with(m2, ci.coef[2,1]<2 & ci.coef[2,2]>2))
+  names(res) <- c("psiint", "psi", "varint", "var",  "powint", "pow",  "cover", "b.psiint", "b.psi", "b.varint", "b.var", "b.powint", "b.pow", "b.cover")
+  res
+}
+
+
+#res = mclapply(1:1000, repit)
+res = lapply(1:2, repit)
+res = simplify2array(res)
+
+# equality within toleraance
+stopifnot(all.equal(res["psiint",],res["b.psiint",]))
+stopifnot(all.equal(res["psi",],res["b.psi",]))
+
+#' \dontest{
+#' # bootstrap and regular variance good
+#' repit2 <- function(i){
+#'   dat = dgm(500)
+#'   m1 = qgcomp.noboot(y~., expnms=c("x1", "x2"), data = dat, family=gaussian(), q=4)
+#'   m2 = qgcomp.boot(  y~., expnms=c("x1", "x2"), data = dat, family=gaussian(), q=4, B=5, parallel=TRUE)
+#'   res = c(m1$coef, m1$var.coef, 1*(m1$pval>0.05), with(m1, ci.coef[1]<2 & ci.coef[2]>2), m2$coef, m2$var.coef, 1*(m2$pval>0.05), with(m2, ci.coef[2,1]<2 & ci.coef[2,2]>2))
+#'   names(res) <- c("psiint", "psi", "varint", "var",  "powint", "pow",  "cover", "b.psiint", "b.psi", "b.varint", "b.var", "b.powint", "b.pow", "b.cover")
+#'   res
+#' }
+#' 
+#' res = lapply(1:2, repit2)
+#' res = simplify2array(res)
+#' 
+#' 
+#' stopifnot(all.equal(res["var",],res["b.var",], tolerance=sqrt(0.01)))
+#' }
+
+
+cat("done")

tests/test_l_cox_msmtest.R

@@ -0,0 +1,70 @@
+cat("# cox msm test\n")
+library(qgcomp)
+# does the simulation version of a MS cox model yield the same result as
+# the sum-of-random-variables version when there are no non-exposure covariates
+# THIS TEST TAKES A LONG TIME!
+#set.seed(50)
+#N=500
+#dat <- data.frame(time=(tmg <- pmin(.1,rweibull(N, 10, 0.1))), 
+#                d=1.0*(tmg<0.1), x1=runif(N), x2=runif(N), z=runif(N))
+#expnms=paste0("x", 1:2)
+#f = survival::Surv(time, d)~x1 + x2
+#f2 = survival::Surv(time, d)~x1 + x2 + z
+#f3 = survival::Surv(time, d)~x1 + x2 + I(z-mean(z))
+#survival::coxph(f, data = dat)
+#survival::coxph(f2, data = dat)
+#(obj <- qgcomp.cox.noboot(f, expnms = expnms, data = dat))
+#(obj2 <- qgcomp.cox.boot(f, expnms = expnms, data = dat, B=1000, MCiter=20000))
+#(obj3 <- qgcomp.cox.boot(f2, expnms = expnms, data = dat, B=1000, MCiter=20000, parallel=TRUE))
+#
+#stopifnot(all.equal(obj$psi, obj2$psi, tolerance = .01))
+#stopifnot(all.equal(obj$var.psi, obj2$var.psi, tolerance = .1))
+#
+##(r00 <- qgcomp.cox.noboot(f2, expnms = expnms, data = dat))
+##system.time(r0 <- qgcomp.cox.boot(f2, expnms = expnms, data = dat, 
+##                            B=16, MCiter=20000, parallel=FALSE))
+###  Expected time to finish: 1.54 minutes 
+###  user  system elapsed 
+###  90.438  11.711 102.804 
+##system.time(r1 <- qgcomp.cox.boot(f2, expnms = expnms, data = dat, 
+##                            B=16, MCiter=20000, parallel=TRUE))
+### fairly high overhead vs. mclapply (but works on windows)
+### Expected time to finish: 3.75 minutes 
+### user  system elapsed 
+### 10.922   1.846  51.252 
+#
+#
+#
+##(obj4 <- qgcomp.cox.boot(survival::Surv(time, d)~factor(x1) + splines::bs(x2) + z, 
+##                         expnms = expnms, data = dat, 
+##                         B=10, MCiter=20000, parallel=FALSE, degree=2))
+##
+##lapply(1:5, mean)
+##parallel::mclapply(1:5, mean)
+### windows friendly version
+##future.apply::future_sapply(1:5, mean)
+#
+## LATE ENTRY
+
+set.seed(50)
+N=50
+dat <- data.frame(stop=(tmg <- pmin(.1,rweibull(N, 10, 0.1))), 
+                  start = pmax(0, tmg-runif(N)),
+                  d=1.0*(tmg<0.1), x1=runif(N), x2=runif(N), z=runif(N))
+expnms=paste0("x", 1:2)
+f = survival::Surv(start,stop, d)~x1 + x2
+f = survival::Surv(start,stop, d)~x1 + x2
+suppressWarnings(obj <- qgcomp.cox.boot(f, expnms = expnms, data = dat, B=2, MCsize=200, parallel=FALSE))
+suppressWarnings(obj <- qgcomp.cox.boot(f, expnms = expnms, data = dat, deg=2, B=2, MCsize=200, q=NULL, parallel=FALSE))
+#plot(obj)
+summary(obj)
+
+ymat = obj$fit$y
+tval = grep("stop|time",colnames(ymat) , value=TRUE)
+stop = as.numeric(ymat[,tval])
+times = sort(-sort(-unique(stop))[-1])
+datatimes = with(dat, sort(unique(stop*d))[-1])
+
+stopifnot(all.equal(times, datatimes))
+
+cat("done")

tests/test_l_utilities.R

@@ -0,0 +1,41 @@
+cat('testing utility functions')
+
+
+qgcomp:::construction("msg", "test")
+
+qgcomp:::cox()
+
+qgcomp:::zi()
+
+
+# zero inflated models
+set.seed(50)
+n=100
+dat <- data.frame(y=rbinom(n, 1, 0.5)*rpois(n, 1.2), x1=runif(n), x2=runif(n), z=runif(n))
+
+# poisson count model, mixture in both portions
+qgcomp::qgcomp.zi.noboot(f=y ~ z + x1 + x2 | x1 + x2, expnms = c('x1', 'x2'), 
+    data=dat, q=2, dist="poisson")
+
+qgcomp::qgcomp.zi.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), 
+   data=dat, q=2, dist="negbin") # equivalent
+
+
+ftz = qgcomp::qgcomp.zi.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), 
+                            data=dat, q=3, dist="negbin", B=2, parallel=FALSE, MCsize = 100) # equivalent
+fth = qgcomp::qgcomp.hurdle.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), 
+                            data=dat, q=3, dist="poisson", B=2, parallel=FALSE, MCsize = 100) # equivalent
+
+
+res = try(qgcomp:::glance.qgcompfit(ftz), silent=TRUE)
+stopifnot(inherits(res,"try-error"))
+res = try(qgcomp:::glance.qgcompfit(fth), silent=TRUE)
+stopifnot(inherits(res,"try-error"))
+
+ft = qgcomp::qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), 
+                           data=dat, q=3, family="gaussian") # equivalent
+qgcomp:::glance.qgcompfit(ft)
+ft = qgcomp::qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), 
+                           data=dat, q=3, family="poisson") # equivalent
+
+qgcomp:::tidy.qgcompfit(ft)

tests/test_l_zi.R

@@ -0,0 +1,36 @@
+cat("# zi test\n")
+library(qgcomp)
+set.seed(1231)
+n=300
+dat <- data.frame(y=rbinom(n, 1, 0.5)*rpois(n, 1.2), x1=runif(n), x2=runif(n), z=runif(n))
+
+# should not cause error
+qgcomp.zi.noboot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=2, dist="negbin")
+qgcomp.zi.noboot(f=y ~ z + x1 + x2 | x1 + x2 + z, expnms = c('x1', 'x2'), data=dat, q=2, dist="negbin")
+qgcomp.hurdle.noboot(f=y ~ z + x1 + x2 | x1 + x2 + z, expnms = c('x1', 'x2'), data=dat, q=2, dist="negbin")
+
+qgcomp.zi.boot(f=y ~ z + x1 + x2 | x1 + x2 + z, B=1, expnms = c('x1', 'x2'), msmcontrol =zimsm.fit.control(predmethod="catprobs") , data=dat, q=2, dist="negbin")
+qgcomp.hurdle.boot(f=y ~ z + x1 + x2 | x1 + x2 + z, B=1, expnms = c('x1', 'x2'), msmcontrol =zimsm.fit.control(predmethod="components") , data=dat, q=2, dist="negbin")
+
+
+
+qgcomp.zi.boot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=4, B = 10, MCsize = 100, 
+               dist="poisson")
+
+# qgcomp.hurdle.boot is iffy when mixture is not in zero model
+#rr = qgcomp.hurdle.boot(f=y ~ z + x1 + x2 | x1 + x2 + z, expnms = c('x1', 'x2'), data=dat, q=2, B = 10, MCsize = 1000)
+#summary(rr)
+#'\dontrun{
+#' qgcomp.zi.boot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=4, B = 10, MCsize = 1000, 
+#'                dist="negbin", msmcontrol = zimsm.fit.control(predmethod="catprobs"))
+#' ff = qgcomp.zi.boot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=4, B = 10, MCsize = 100, 
+#'                dist="geometric")
+#' plot(ff)
+#'}              
+
+
+# should cause error
+res <- try(qgcomp.zi.noboot(f=y ~ z + x1 + x2 | x1 + z, expnms = c('x1', 'x2'), data=dat, q=2, dist="negbin"), silent = TRUE)
+stopifnot(class(res)=="try-error")
+
+cat("done")

tests/test_l_zzzplots.R

@@ -0,0 +1,77 @@
+cat("# plot test\n")
+library(qgcomp)
+set.seed(112312)
+n=100
+
+#base
+# binomial
+   dat <- data.frame(y=rbinom(n, 1, 0.5), x1=runif(n), x2=runif(n), z=runif(n))
+   ee = qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, family=binomial())
+   plot(ee)
+   plot(ee)
+   ff = qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, B=5, family=binomial(), rr=FALSE)
+   plot(ff)
+   ff$msmfit$family$link
+   pointwisebound.boot(ff)
+   qgcomp:::pointwisebound.boot_old(ff)
+   modelbound.boot(ff)
+   qgcomp:::modelbound.boot_old(ff)
+
+   gg = qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, B=5, family=binomial(), rr=TRUE)
+   plot(gg)
+   pointwisebound.boot(gg)
+   modelbound.boot(gg)
+
+
+# gaussian
+   dat <- data.frame(y=rnorm(n), x1=runif(n), x2=runif(n), z=runif(n))
+   ee = qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, family=gaussian())
+   plot(ee)
+   pointwisebound.noboot(ee) 
+   qgcomp:::pointwisebound.noboot_old(ee) 
+   ff = qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, B=3, family=gaussian())
+   plot(ff, flexfit = FALSE)
+   modelbound.boot(ff)
+   qgcomp:::modelbound.boot_old(ff)
+
+# poisson
+   dat <- data.frame(y=rpois(n, 1.2), x1=runif(n), x2=runif(n), z=runif(n))
+   ee = qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, family=poisson())
+   plot(ee)
+   pointwisebound.noboot(ee) 
+   qgcomp:::pointwisebound.noboot_old(ee) 
+   ff = qgcomp.boot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=7, B=5, family=poisson())
+   modelbound.boot(ff)
+   qgcomp:::modelbound.boot_old(ff)
+   plot(ff)
+
+#cox
+   dat <- data.frame(stop=(tmg <- pmin(.1,rweibull(n, 10, 0.1))), 
+                     start = pmax(0, tmg-runif(n)),
+                     d=1.0*(tmg<0.1), x1=runif(n), x2=runif(n), z=runif(n))
+   expnms=paste0("x", 1:2)
+   f = survival::Surv(start,stop, d)~x1 + x2
+   suppressWarnings(ee <- qgcomp.cox.noboot(f, expnms = expnms, data = dat))
+   res = try(pointwisebound.noboot(ee), silent = TRUE) # not working
+   stopifnot(class(res)=="try-error")   
+   plot(ee)
+   suppressWarnings(ff <- qgcomp.cox.boot(f, expnms = expnms, data = dat, B=2, MCsize=500))
+   plot(ff)
+   res = try(modelbound.boot(ff, pwonly=TRUE), silent = TRUE) # not working
+   stopifnot(class(res)=="try-error")   
+   curvelist = qgcomp.survcurve.boot(ff)
+   stopifnot(inherits(curvelist, "list"))
+
+# zi
+  dat <- data.frame(y=rbinom(n, 1, 0.5)*rpois(n, 1.2), x1=runif(n), x2=runif(n), z=runif(n))
+  ee = qgcomp.zi.noboot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=7, dist="negbin")
+  plot(ee)
+  res = try(pointwisebound.noboot(ee), silent = TRUE) # not working
+  stopifnot(class(res)=="try-error")   
+  ffz = qgcomp.zi.boot(f=y ~ z + x1 + x2 | z, expnms = c('x1', 'x2'), data=dat, q=7, B=2, MCsize=500, dist="negbin")
+  pointwisebound.boot(ffz)
+  modelbound.boot(ffz, pwonly=TRUE)
+  plot(ffz)
+
+
+cat("done")