R Under development (unstable) (2025-08-19 r88650 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.59)
     x1      x2 
0.99604 0.00396 

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

Mixture log(OR) (delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -2.82444    0.91352 -4.61490   -1.034 -3.0918  0.001989
psi1         1.58685    0.38477  0.83271    2.341  4.1241 3.721e-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.621)
     x1      x2 
0.99567 0.00433 

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

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

Scaled effect size (negative direction, sum of negative coefficients = -3.69)
   x1    x2 
0.973 0.027 

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

Prob(Y ~ count):
            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -1.85235    0.40220 -2.64064  -1.0641 -4.6056 4.113e-06
psi1         0.62129    0.19529  0.23852   1.0041  3.1813  0.001466
Prob(Y ~ zero):
            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept)  -9.0730   157.0545  -316.89   298.75 -0.0578   0.9539
psi1         -3.6874        NaN      NaN      NaN     NaN      NaN
Warning messages:
1: In sqrt(diag(object$vcov)) : NaNs produced
2: In sqrt(var) : NaNs produced
> 
> 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.7114858  0.3388613 -2.3756417 -1.0473299 -5.050697 4.402001e-07
psi1         0.5429245  0.1416279  0.2653389  0.8205102  3.833457 1.263551e-04

> 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.71    0.42    3.12