R Under development (unstable) (2023-08-08 r84908 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 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.14)
x1 
 1 

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

Mixture log(OR) (Delta method CI):

            Estimate Std. Error Lower CI Upper CI Z value  Pr(>|z|)
(Intercept) -2.14509    0.84939  -3.8099 -0.48032 -2.5254   0.01156
psi1         1.91462    0.43341   1.0651  2.76409  4.4175 9.983e-06
> 
> 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.626)
    x1     x2 
0.9611 0.0389 

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

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

Scaled effect size (negative direction, sum of negative coefficients = -19.6)
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.79184    0.47733 -2.72738  -0.8563 -3.7539 0.0001741
psi1         0.62596    0.21293  0.20862   1.0433  2.9397 0.0032849
Prob(Y ~ zero):
            Estimate Std. Error Lower CI Upper CI Z value Pr(>|z|)
(Intercept) -12.7787    66.7284  -143.56   118.01 -0.1915   0.8481
psi1         -9.1228    62.9406  -132.48   114.24 -0.1449   0.8848
> 
> 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.7576035 0.11748334 -1.9878666 -1.5273404 -14.96045        0
psi1         0.5678248 0.03923117  0.4909331  0.6447165  14.47382        0

> 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.68    0.50    4.17