R Under development (unstable) (2025-01-21 r87610 ucrt) -- "Unsuffered Consequences"
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.31)
x1 
 1 

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

Mixture log(OR) (delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -3.95360    1.16647  -6.2398  -1.6674 -3.3894 0.0007005
psi1         2.11854    0.50142   1.1358   3.1013  4.2251 2.388e-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.425)
    x1     x2 
0.9549 0.0451 

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

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

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

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

Prob(Y ~ count):
            Estimate Std. Error   Lower CI Upper CI Z value Pr(>|z|)
(Intercept) -1.25586    0.48307 -2.2026540 -0.30907 -2.5998 0.009329
psi1         0.42526    0.22004 -0.0060109  0.85654  1.9326 0.053280
Prob(Y ~ zero):
            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)   9.1889    88.8890  -165.03   183.41  0.1034   0.9177
psi1        -31.0897   122.0096  -270.22   208.04 -0.2548   0.7989
> 
> 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.6986481  0.3160082 -2.3180128 -1.0792835 -5.375330 7.644274e-08
psi1         0.5497356  0.1032018  0.3474639  0.7520074  5.326805 9.995541e-08

> 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 
   2.95    0.42    3.32