R version 4.5.0 beta (2025-03-27 r88065 ucrt) -- "How About a Twenty-Six"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> cat("# bayesqgcomp test\n")
# bayesqgcomp test
> #devtools::install_github("alexpkeil1/qgcomp@dev")
> library("qgcomp")
> data("metals", package="qgcomp")
> 
> Xnm <- c(
+   'arsenic','barium','cadmium','calcium','chromium','copper',
+   'iron','lead','magnesium','manganese','mercury','selenium','silver',
+   'sodium','zinc'
+ )
> 
> # continuous outcome, example with perfect collinearity
> metals$cadmium2 = metals$cadmium
> Xnm2 = c(Xnm, "cadmium2")
> res = try(qc.fit <- qgcomp.noboot(y~.,dat=metals[,c(Xnm2, 'y')], family=gaussian()), silent=TRUE)
Including all model terms as exposures of interest

> stopifnot(class(res)=="try-error")
> # error
> 
> res = try(qc.fit <- qgcomp.noboot(y~.,dat=metals[,c(Xnm2, 'y')], family=gaussian(), bayes=TRUE))
Including all model terms as exposures of interest

> stopifnot(inherits(res, "qgcompfit"))
> 
> # compare to results with colinear exposure removed
> res = try(qc.fit2 <- qgcomp.noboot(y~.,dat=metals[,c(Xnm, 'y')], family=gaussian()))
Including all model terms as exposures of interest

> stopifnot(inherits(res, "qgcompfit"))
> 
> #hitting code coverage just to check
> dat = qgcomp:::.dgm_quantized(N=100)
> dat$y = as.numeric(dat$y>median(dat$y))
> 
> qgcomp(y~.,dat=dat, expnms=c("x1", "x2"), family="binomial", rr=FALSE)
Scaled effect size (positive direction, sum of positive coefficients = 2.04)
x1 
 1 

Scaled effect size (negative direction, sum of negative coefficients = -0.224)
x2 
 1 

Mixture log(OR) (delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -3.33386    0.89064  -5.0795  -1.5882 -3.7432 0.0001817
psi1         1.81475    0.41166   1.0079   2.6216  4.4084 1.041e-05
> 
> res = try(qgcomp(y~x1+x2,dat=dat, expnms=c("x1", "x2"), family=NULL, rr=FALSE), silent=TRUE)
NULL
> stopifnot(inherits(res, "try-error"))
> 
> qgcomp(y~x1+x2 | x1+x2,dat=dat, expnms=c("x1", "x2"), family="poisson")
Prob(Y ~ count):
Scaled effect size (positive direction, sum of positive coefficients = 0.319)
    x1     x2 
0.9419 0.0581 

Scaled effect size (negative direction, sum of negative coefficients = 0)
None

Prob(Y ~ zero/count):
Scaled effect size (positive direction, sum of positive coefficients = 8.76)
x2 
 1 

Scaled effect size (negative direction, sum of negative coefficients = -35.7)
x1 
 1 

Mixture log(OR/RR) (Delta method CI):

Prob(Y ~ count):
            Estimate Std. Error Lower CI  Upper CI Z value Pr(>|z|)
(Intercept) -0.99917    0.47862 -1.93725 -0.061096 -2.0876  0.03683
psi1         0.31907    0.21738 -0.10699  0.745131  1.4678  0.14216
Prob(Y ~ zero):
            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)   10.372    107.379  -200.09   220.83  0.0966   0.9230
psi1         -26.914    161.838  -344.11   290.28 -0.1663   0.8679
> 
> res = try(qgcomp(y~x1+x2 | x1+x2 +I(x2^2), B=2,dat=dat, expnms=c("x1", "x2"), bayes=TRUE, family="poisson"), silent=TRUE)
> stopifnot(inherits(res, "try-error"))
> 
> 
> ft = qgcomp(y~x1+x2 + I(x2^2), B=2,dat=dat, expnms=c("x1", "x2"), bayes=TRUE, family="poisson")
> pp = msm.predict(ft)
> ft = qgcomp(y~x1+x2,dat=dat, expnms=c("x1", "x2"), bayes=TRUE, family="poisson")
> res = try(msm.predict(ft), silent = TRUE)
> stopifnot(inherits(res, "try-error"))
> 
> 
> ft = qgcomp(y~x1+x2 + I(x2^2), B=2, q=NULL,dat=dat, expnms=c("x1", "x2"), bayes=TRUE, family=binomial())

Note: using all possible values of exposure as the
              intervention values

> summary(ft)
Mixture log(RR) (bootstrap CI):

$coefficients
              Estimate Std. Error   Lower CI   Upper CI   Z value     Pr(>|z|)
(Intercept) -1.6884449  0.3136663 -2.3032196 -1.0736702 -5.382933 7.328173e-08
psi1         0.5389419  0.1037504  0.3355949  0.7422889  5.194602 2.051578e-07

> qc.fit2 <- qgcomp.boot(y~.,expnms = Xnm, dat=metals[,c(Xnm, 'y')],B=2, q=NULL, family=gaussian())

Note: using quantiles of all exposures combined in order to set
          proposed intervention values for overall effect (25th, 50th, 75th %ile)
        You can ensure this is valid by scaling all variables in expnms to have similar ranges.
> 
> 
> ## slight shrinkage with heavily non-linear models that don't converge
> ## non-Bayesian
> #results5 = qgcomp(y~. + .^2 + .^3 + arsenic*cadmium,
> #                  expnms=Xnm,
> #                  metals[,c(Xnm, 'y')], family=gaussian(), q=10, B=10, 
> #                  seed=125, degree=3)
> #
> #print(results5)
> #plot(results5)
> #results5$bootsamps # samples are wack
> #
> ## Bayesian (takes a few mins)
> #results5bayes = qgcomp.boot(y~. + .^2 + .^3 + arsenic*cadmium,
> #                  expnms=Xnm,
> #                  metals[,c(Xnm, 'y')], family=gaussian(), q=10, B=10, 
> #                  seed=125, degree=3, 
> #                  bayes=TRUE)
> #
> ## note p-values are not computed due to negative
> ## degrees of freedom (can fix by using a z-statistic instead of a t-statistic)
> #print(results5bayes)
> #plot(results5bayes)
> #results5bayes$bootsamps # samples are reasonable
> 
> 
> cat("done")
done> 
> proc.time()
   user  system elapsed 
   3.26    0.34    3.57