R Under development (unstable) (2025-02-22 r87796 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 = 1.64)
x1 
 1 

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

Mixture log(OR) (delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -2.33970    0.79047 -3.88899 -0.79041 -2.9599  0.003077
psi1         1.57827    0.39551  0.80308  2.35346  3.9905 6.594e-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.59)
x1 
 1 

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

Prob(Y ~ zero/count):
Scaled effect size (positive direction, sum of positive coefficients = 0)
None

Scaled effect size (negative direction, sum of negative coefficients = -30.8)
   x1    x2 
0.506 0.494 

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

Prob(Y ~ count):
            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -1.68325    0.41696 -2.50047 -0.86603 -4.0370 5.414e-05
psi1         0.53519    0.20121  0.14083  0.92956  2.6599  0.007817
Prob(Y ~ zero):
            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)   7.8299    49.7571  -89.692   105.35  0.1574   0.8750
psi1        -30.8235   229.6727 -480.974   419.33 -0.1342   0.8932
> 
> 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.7631522 0.17862055 -2.1132421 -1.4130624 -9.870937 0.000000e+00
psi1         0.5640181 0.07600852  0.4150441  0.7129921  7.420459 1.167955e-13

> 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 
   1.89    0.32    2.15