R Under development (unstable) (2025-03-11 r87944 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. > library(qgcomp) > library(qgcompint) > set.seed(40) > # linear model > dat <- data.frame(y=runif (50), + x1=runif (50), + x2=runif (50), + z=rbinom(50, 1, 0.5), + r=rbinom(50, 1, 0.5)) > (qfit <- qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) Scaled effect size (positive direction, sum of positive coefficients = 0) None Scaled effect size (negative direction, sum of negative coefficients = -0.117) x2 x1 0.692 0.308 Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.516037 0.086514 0.34647 0.68560 5.9648 3.043e-07 psi1 -0.117418 0.130342 -0.37288 0.13805 -0.9008 0.3723 > (qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = -0.278) x2 x1 0.662 0.338 ## Qgcomp weights/partial effects at z = Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = 0) None ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.58062 0.11142 0.36224 0.79900 5.2112 4.787e-06 psi1 -0.27807 0.20757 -0.68490 0.12876 -1.3397 0.1872 z -0.10410 0.15683 -0.41148 0.20329 -0.6637 0.5103 z:mixture 0.26811 0.26854 -0.25822 0.79444 0.9984 0.3235 Estimate (CI), z=1: -0.27807 (-0.6849, 0.12876) > (qfitemmboot <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian(), B=10)) ## Mixture slope parameters (bootstrap CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.58062 0.12482 0.33598 0.82526 4.6517 3.139e-05 psi1 -0.27807 0.15194 -0.57587 0.01973 -1.8301 0.07417 z -0.10410 0.14908 -0.39628 0.18809 -0.6983 0.48877 z:mixture 0.26811 0.21010 -0.14367 0.67989 1.2761 0.20876 > (qfitemmee <- qgcomp.emm.glm.ee(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) ## Mixture slope parameters (robust CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.58062 0.12470 0.33621 0.825039 4.6560 3.096e-05 psi1 -0.27807 0.18653 -0.64366 0.087517 -1.4908 0.1433 z -0.10410 0.16457 -0.42665 0.218464 -0.6325 0.5304 mixture:z 0.26811 0.24757 -0.21712 0.753348 1.0830 0.2849 Estimate (CI), z=1: -0.27807 (-0.64366, 0.087517) > (qfitemmb <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2, bayes=TRUE, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = -0.266) x2 x1 0.664 0.336 ## Qgcomp weights/partial effects at z = Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = 0) None ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.575366 0.109514 0.36072 0.79001 5.2538 4.153e-06 psi1 -0.266227 0.203004 -0.66411 0.13165 -1.3114 0.1965 z -0.095933 0.153072 -0.39595 0.20408 -0.6267 0.5341 z:mixture 0.251544 0.260955 -0.25992 0.76301 0.9639 0.3403 Estimate (CI), z=1: -0.26623 (-0.66411, 0.13165) > (qfitemmbootb <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z", bayes=TRUE, expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian(), B=10)) ## Mixture slope parameters (bootstrap CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.575366 0.125877 0.328651 0.822080 4.5709 4.071e-05 psi1 -0.266227 0.144178 -0.548811 0.016358 -1.8465 0.07171 z -0.095933 0.107493 -0.306615 0.114749 -0.8925 0.37711 z:mixture 0.251544 0.153631 -0.049566 0.552654 1.6373 0.10886 > # check > getstratweights(qfitemm, emmval=0.0) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = -0.278) x2 x1 0.662 0.338 > getstrateffects(qfitemm, emmval=0.0) Joint effect at z=0 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.27807 0.20757 -0.6849 0.12876 -1.3397 0.1804 > getstratweights(qfitemm, emmval=1.0) ## Qgcomp weights/partial effects at z = 1 Scaled effect size (positive direction, sum of positive effects = 0.0028) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.0128) x2 1 > getstrateffects(qfitemm, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.0099575 0.1703789 -0.34389 0.32398 -0.0584 0.9534 > getstratweights(qfitemmb, emmval=1.0) ## Qgcomp weights/partial effects at z = 1 Scaled effect size (positive direction, sum of positive effects = 0.000944) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.0156) x2 1 > getstrateffects(qfitemmb, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.014683 0.168655 -0.34524 0.31587 -0.0871 0.9306 > getstrateffects(qfitemmee, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.0099575 0.1627897 -0.32902 0.3091 -0.0612 0.9512 > > #null > (qfitemmn <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2, bayes=TRUE, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=NULL, family=gaussian())) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0.143) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.255) x2 1 ## Qgcomp weights/partial effects at z = Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = 0) None ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.50495 0.15530 0.20057 0.80932 3.2515 0.002206 psi1 -0.11245 0.28832 -0.67755 0.45265 -0.3900 0.698405 z -0.17684 0.23980 -0.64685 0.29316 -0.7375 0.464759 z:mixture 0.37529 0.43261 -0.47261 1.22319 0.8675 0.390375 Estimate (CI), z=1: -0.11245 (-0.67755, 0.45265) > (qfitemmbootn <- qgcomp.emm.glm.boot(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=NULL, family=gaussian(), B=10)) 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. ## Mixture slope parameters (bootstrap CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.51583 0.15767 0.20681 0.82485 3.2716 0.002112 psi1 -0.13586 0.23968 -0.60563 0.33391 -0.5668 0.573764 z -0.20593 0.33595 -0.86438 0.45252 -0.6130 0.543123 z:mixture 0.43153 0.54698 -0.64053 1.50360 0.7889 0.434481 > (qfitemmeen <- qgcomp.emm.glm.ee(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, 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. ## Mixture slope parameters (robust CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.51583 0.15830 0.20558 0.82608 3.2587 0.002191 psi1 -0.13586 0.25573 -0.63708 0.36535 -0.5313 0.597957 z -0.20593 0.23730 -0.67102 0.25917 -0.8678 0.390316 mixture:z 0.43153 0.40062 -0.35368 1.21674 1.0771 0.287421 Estimate (CI), z=1: -0.13586 (-0.63708, 0.36535) > getstratweights(qfitemmn, emmval=1.0) ## Qgcomp weights/partial effects at z = 1 Scaled effect size (positive direction, sum of positive effects = 0.263) x2 x1 0.613 0.387 Scaled effect size (negative direction, sum of negative effects = 0) None > getstrateffects(qfitemmn, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture 0.26284 0.34164 -0.40676 0.93244 0.7694 0.4417 > getstrateffects(qfitemmbootn, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture 0.29567 0.35531 -0.40073 0.99206 0.8321 0.4053 > getstrateffects(qfitemmeen, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture 0.29567 0.30839 -0.30876 0.9001 0.9588 0.3377 > > > dat$z=as.factor(dat$z) > (qfitemmf <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = -0.278) x2 x1 0.662 0.338 ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0) None Scaled effect size (negative direction, sum of negative effects = -0.278) x2 x1 0.662 0.338 ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.58062 0.11142 0.36224 0.79900 5.2112 4.787e-06 psi1 -0.27807 0.20757 -0.68490 0.12876 -1.3397 0.1872 z1 -0.10410 0.15683 -0.41148 0.20329 -0.6637 0.5103 z1:mixture 0.26811 0.26854 -0.25822 0.79444 0.9984 0.3235 Estimate (CI), z=1: -0.27807 (-0.6849, 0.12876) > > checkres <- function(qfitemm) { + q0 <- q1 <- qfitemm$fit$data + #q0$z <- q1$z <- 1 + q0$x1 <- q0$x2 <- 0 + q1$x1 <- q1$x2 <- 1 + #qfitemm$fit$coefficients + qfitemm + mean(predict(qfitemm, newdata = q1[q1$z==0, ]))- + mean(predict(qfitemm, newdata = q0[q0$z==0, ])) + mean(predict(qfitemm, newdata = q1[q1$z==1, ]))- + mean(predict(qfitemm, newdata = q0[q0$z==1, ])) + } > > > > # categorical modifier > n <- 100 > dat <- data.frame(y=runif (n), + x1=runif (n), + x2=runif (n), + z=as.factor(sample(c("a", "b", "m"), n, replace=TRUE)), + r=rbinom(n, 1, 0.5)) > (qfit <- qgcomp.noboot(f=y ~ z + x1 + x2, expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) Scaled effect size (positive direction, sum of positive coefficients = 0.017) x1 1 Scaled effect size (negative direction, sum of negative coefficients = -0.0142) x2 1 Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.4299178 0.0704197 0.2919 0.56794 6.1051 2.119e-08 psi1 0.0027983 0.0807145 -0.1554 0.16100 0.0347 0.9724 > (qfitemm <- qgcomp.emm.glm.noboot(f=y ~ x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian())) ## Qgcomp weights/partial effects at z = a Scaled effect size (positive direction, sum of positive effects = 0.0215) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.0692) x2 1 ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.5163707 0.1192040 0.282735 0.750006 4.3318 3.766e-05 psi1 -0.0476841 0.0622652 -0.169722 0.074353 -0.7658 0.4457 zb -0.0082289 0.1601608 -0.322138 0.305681 -0.0514 0.9591 zb:mixture 0.0923514 0.0911490 -0.086297 0.271000 1.0132 0.3136 zm -0.0205183 0.1572054 -0.328635 0.287599 -0.1305 0.8964 zm:mixture 0.0560174 0.0898595 -0.120104 0.232139 0.6234 0.5346 > (qfitemmboot <- qgcomp.emm.glm.boot(f=y ~ x1 + x2, degree=1, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian(), B=10)) ## Mixture slope parameters (bootstrap CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.5163707 0.0867415 0.34636 0.686381 5.9530 4.994e-08 psi1 -0.0476841 0.0441180 -0.13415 0.038786 -1.0808 0.2827 zb -0.0082289 0.1437691 -0.29001 0.273553 -0.0572 0.9545 zm -0.0205183 0.1379753 -0.29094 0.249908 -0.1487 0.8821 zb:mixture 0.0923514 0.0743080 -0.05329 0.237992 1.2428 0.2172 zm:mixture 0.0560174 0.0887852 -0.11800 0.230033 0.6309 0.5297 > (qfitemmee <- qgcomp.emm.glm.ee(f=y ~ x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian())) ## Mixture slope parameters (robust CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.5163707 0.1113176 0.298192 0.734549 4.6387 1.184e-05 psi1 -0.0476841 0.0595845 -0.164468 0.069099 -0.8003 0.4257 zb -0.0082289 0.1370382 -0.276819 0.260361 -0.0600 0.9523 zm -0.0205183 0.1544373 -0.323210 0.282173 -0.1329 0.8946 mixture:zb 0.0923514 0.0770389 -0.058642 0.243345 1.1988 0.2338 mixture:zm 0.0560174 0.0858192 -0.112185 0.224220 0.6527 0.5156 > (qfitemmb <- qgcomp.emm.glm.noboot(f=y ~ x1 + x2, bayes=TRUE, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) ## Qgcomp weights/partial effects at z = a Scaled effect size (positive direction, sum of positive effects = 0.122) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.129) x2 1 ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.4363649 0.0909072 0.25819 0.61454 4.8001 6.112e-06 psi1 -0.0073751 0.1264609 -0.25523 0.24048 -0.0583 0.9536 zb 0.0728693 0.1301237 -0.18217 0.32791 0.5600 0.5768 zb:mixture 0.1387017 0.2093916 -0.27170 0.54910 0.6624 0.5094 zm 0.0995805 0.1180172 -0.13173 0.33089 0.8438 0.4010 zm:mixture -0.0573486 0.1817763 -0.41362 0.29893 -0.3155 0.7531 > # check > (levs = levels(dat$z)) [1] "a" "b" "m" > getstratweights(qfitemm, emmval="b") ## Qgcomp weights/partial effects at z = b Scaled effect size (positive direction, sum of positive effects = 0.0447) x1 x2 0.517 0.483 Scaled effect size (negative direction, sum of negative effects = 0) None > getstratweights(qfitemm, emmval="a") ## Qgcomp weights/partial effects at z = a Scaled effect size (positive direction, sum of positive effects = 0.0215) x1 1 Scaled effect size (negative direction, sum of negative effects = -0.0692) x2 1 > getstrateffects(qfitemm, emmval="m") Joint effect at z=m Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture 0.0083333 0.0647902 -0.11865 0.13532 0.1286 0.8977 > getstratweights(qfitemm, emmval="b") ## Qgcomp weights/partial effects at z = b Scaled effect size (positive direction, sum of positive effects = 0.0447) x1 x2 0.517 0.483 Scaled effect size (negative direction, sum of negative effects = 0) None > getstrateffects(qfitemm, emmval="b") Joint effect at z=b Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture 0.044667 0.066567 -0.085802 0.17514 0.671 0.5022 > getstratweights(qfitemmb, emmval="b") ## Qgcomp weights/partial effects at z = b Scaled effect size (positive direction, sum of positive effects = 0.131) x2 x1 0.641 0.359 Scaled effect size (negative direction, sum of negative effects = 0) None > > > > ##### continuous emm, assessing type 1 error only > n = 200 > dat <- data.frame(y=runif (n), + x1=runif (n), + x2=runif (n), + x3=runif (n), + x4=runif (n), + x5=runif (n), + z=rnorm(n, 1, 0.5), + r=rbinom(n, 1, 0.5)) > (qfit <- qgcomp.noboot(f=y ~ ., expnms = paste0("x", 1:5), data=dat, q=4, family=gaussian())) Scaled effect size (positive direction, sum of positive coefficients = 0.0206) x2 x5 0.78 0.22 Scaled effect size (negative direction, sum of negative coefficients = -0.0986) x1 x3 x4 0.361 0.353 0.286 Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.573594 0.082174 0.41254 0.7346518 6.9803 4.556e-11 psi1 -0.078014 0.041476 -0.15931 0.0032783 -1.8809 0.06148 > (qfitemm <- qgcomp.emm.glm.noboot(f=y ~x1+x2+x3+x4+x5, emmvar="z", expnms = paste0("x", 1:5), data=dat, q=2, family=gaussian())) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0.0474) x4 1 Scaled effect size (negative direction, sum of negative effects = -0.345) x2 x1 x5 x3 0.304 0.262 0.236 0.197 ## Mixture slope parameters (delta method CI): Estimate Std. Error Lower CI Upper CI t value Pr(>|t|) (Intercept) 0.598460 0.120495 0.36229 0.83463 4.9667 1.524e-06 psi1 -0.297682 0.211454 -0.71212 0.11676 -1.4078 0.1608 z -0.018558 0.105917 -0.22615 0.18904 -0.1752 0.8611 z:mixture 0.089297 0.182371 -0.26814 0.44674 0.4896 0.6250 > > > qfitemm$fit Call: glm(formula = newform, family = ..1, data = qdata, weights = weights) Coefficients: (Intercept) x1 x2 x3 x4 x5 0.598460 -0.090397 -0.105022 -0.068099 0.047417 -0.081582 z x1:z x2:z x3:z x4:z x5:z -0.018558 0.020035 0.150466 -0.007306 -0.120276 0.046378 Degrees of Freedom: 199 Total (i.e. Null); 188 Residual Null Deviance: 18.15 Residual Deviance: 17 AIC: 100.5 > getstratweights(qfitemm, emmval=0.0) ## Qgcomp weights/partial effects at z = 0 Scaled effect size (positive direction, sum of positive effects = 0.0474) x4 1 Scaled effect size (negative direction, sum of negative effects = -0.345) x2 x1 x5 x3 0.304 0.262 0.236 0.197 > getstratweights(qfitemm, emmval=1.0) ## Qgcomp weights/partial effects at z = 1 Scaled effect size (positive direction, sum of positive effects = 0.0454) x2 1 Scaled effect size (negative direction, sum of negative effects = -0.254) x3 x4 x1 x5 0.297 0.287 0.277 0.139 > getstratweights(qfitemm, emmval=2.0) ## Qgcomp weights/partial effects at z = 2 Scaled effect size (positive direction, sum of positive effects = 0.207) x2 x5 0.946 0.054 Scaled effect size (negative direction, sum of negative effects = -0.326) x4 x3 x1 0.592 0.254 0.154 > getstrateffects(qfitemm, emmval=0.0) Joint effect at z=0 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.29768 0.21145 -0.71212 0.11676 -1.4078 0.1592 > getstrateffects(qfitemmb, emmval=1.0) Joint effect at z=1 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.0073751 0.1264609 -0.25523 0.24048 -0.0583 0.9535 > getstratweights(qfitemm, emmval=-2.0) ## Qgcomp weights/partial effects at z = -2 Scaled effect size (positive direction, sum of positive effects = 0.288) x4 1 Scaled effect size (negative direction, sum of negative effects = -0.764) x2 x5 x1 x3 0.531 0.228 0.171 0.070 > getstrateffects(qfitemmb, emmval=-2.0) Joint effect at z=-2 Estimate Std. Error Lower CI Upper CI z value Pr(>|z|) Mixture -0.0073751 0.1264609 -0.25523 0.24048 -0.0583 0.9535 > > > # ----------------------------------- > # long run simulations (none of these are actually run in testing) > # ----------------------------------- > > ## type 1 error true linear model, binary modifier ## > rft <- function(i, n=50) { + dat <- data.frame(y=runif (n), + x1=runif (n), + x2=runif (n), + z=rbinom(n, 1, 0.5), + r=rbinom(n, 1, 0.5)) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + r, expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) + (qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + r, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) + c(qfit$pval[2], qfitemm$pval[c(2, 4)], qfit$psi, qfitemm$psi, qfitemm$psiint, qfit$covmat.psi, qfitemm$covmat.psi, qfitemm$covmat.psiint) + } > > runsims <- function() { + resl = lapply(1:1000, rft) + res = do.call(rbind, resl) + + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + apply(res[, 4:6], 2, mean) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + } > #runsims() > > ## type 1 error for logistic model, binary modifier ## > rftb <- function(i, n=50) { + dat <- data.frame(y=rbinom(n, 1, 0.5), + x1=runif (n), + x2=runif (n), + z=rbinom(n, 1, 0.5), + r=rbinom(n, 1, 0.5)) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + r, expnms = c('x1', 'x2'), data=dat, q=2, family=binomial())) + (qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + r, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=binomial())) + c(qfit$pval[2], qfitemm$pval[c(2, 4)], qfit$psi, qfitemm$psi, qfitemm$psiint, qfit$covmat.psi, qfitemm$covmat.psi, qfitemm$covmat.psiint) + } > > runsimsb <- function() { + resl = lapply(1:1000, rftb, n=100) + res = do.call(rbind, resl) + + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + apply(res[, 4:6], 2, mean) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + } > > > ## type 1 error for true linear model, continuous modifier ## > rft2 <- function(i, n=50) { + dat <- data.frame(y=runif (n), + x1=runif (n), + x2=runif (n), + z=rnorm(n, 1, 0.5), + r=rbinom(n, 1, 0.5)) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + r, expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + r, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian())) + c(qfit$pval[2], qfitemm$pval[c(2, 4)], qfit$psi, qfitemm$psi, qfitemm$psiint, qfit$covmat.psi, qfitemm$covmat.psi, qfitemm$covmat.psiint) + } > > #not run > runsims2 <- function() { + resl = lapply(1:10000, rft2, n=100) + res = do.call(rbind, resl) + + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + apply(res[, 4:6], 2, mean) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + } > #runsims2() > > ## bias, coverage, power for true linear model, binary modifier ## > rft3 <- function(i, n=100) { + dat <- qgcompint:::.dgm_quantized_linear_emm( + N = n, + b0=0, + mainterms= c(0.3, 0.0, 0.3, 0.0), + prodterms = c(0.3, 0.1, 0.1, 0.0), + ztype = "binary", + ncor=0, + corr=0.75 + ) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + x3 + x4, expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + x3 + x4, emmvar="z", expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[2], + pz = qfitemm$pval[c(2, 4)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint, + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint, + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > > runsims3 <- function() { + resl = lapply(1:1000, rft3, n=100) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + (bias <- est - truth) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > #runsims3() > > > ## bias, coverage, power for true linear model, continuous modifier ## > rft4 <- function(i, n=100) { + dat <- qgcompint:::.dgm_quantized_linear_emm( + N = n, + b0=0, + mainterms= c(0.3, 0.0, 0.3, 0.0), + prodterms = c(0.3, 0.1, 0.1, 0.0), + ztype = "cont", + ncor=0, + corr=0.75 + ) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + x3 + x4, expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + x3 + x4, emmvar="z", expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[2], + pz = qfitemm$pval[c(2, 4)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint, + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint, + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > > runsims4 <- function() { + resl = lapply(1:1000, rft4, n=100) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + (bias <- est - truth) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > > #runsims4() > > > > ## bias, coverage, power for true logistic model, binary modifier ## > rft5 <- function(i, n=100) { + dat <- qgcompint:::.dgm_quantized_logistic_emm( + N = n, + b0=0, + mainterms= c(0.3, 0.0, 0.3, 0.0), + prodterms = c(0.3, 0.1, 0.1, 0.0), + ztype = "binary", + ncor=0, + corr=0.75 + ) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + x3 + x4, expnms = paste0("x", 1:4), data=dat, q=NULL, family=binomial())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + x3 + x4, emmvar="z", expnms = paste0("x", 1:4), data=dat, q=NULL, family=binomial())) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[2], + pz = qfitemm$pval[c(2, 4)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint, + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint, + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > > runsims5 <- function() { + resl = lapply(1:1000, rft5, n=200) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + apply(res[, 4:6], 2, median) + (bias <- est - truth) + plot(density(res[, 4])) + plot(density(res[, 5])) + plot(density(res[, 6])) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > #runsims5() > > > ## bias, coverage, power for true logistic model, binary modifier, balancing modification ## > rft6 <- function(i, n=100) { + dat <- qgcompint:::.dgm_quantized_logistic_emm( + N = n, + b0=0, + mainterms= c(0.3, 0.0, 0.3, 0.0), + prodterms = c(-0.3, 0.15, 0.15, 0.0), + ztype = "binary", + ncor=0, + corr=0.75 + ) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + x3 + x4, expnms = paste0("x", 1:4), data=dat, q=NULL, family=binomial())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + x3 + x4, emmvar="z", expnms = paste0("x", 1:4), data=dat, q=NULL, family=binomial())) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[2], + pz = qfitemm$pval[c(2, 4)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint, + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint, + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > > runsims6 <- function() { + resl = lapply(1:1000, rft6, n=200) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + apply(res[, 4:6], 2, median) + (bias <- est - truth) + #plot(density(res[, 4])) + #plot(density(res[, 5])) + #plot(density(res[, 6])) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > #runsims6() > > > > ## bias, coverage, power for true Cox model, binary modifier, balancing modification ## > rft7 <- function(i, n=1000) { + dat = qgcompint:::.dgm_quantized_survival_emm( + N = n, + b0=0, + mainterms=c(1, 0, 0, 0), + prodterms = c(1, 0, 0, 0), + ztype = "bin", + ncor=0, + corr=0.75 + ) + f = survival::Surv(time, d)~x1 + x2 + x3 + x4+z + expnms = c("x1", "x2", "x3", "x4") + #(fit1 <- survival::coxph(f, data = dat)) + (qfit <- qgcomp.cox.noboot(f, expnms = expnms, data = dat, q=4)) + suppressMessages(qfitemm <- qgcomp.emm.cox.noboot(f, expnms = expnms, emmvar="z", data = dat, q=4)) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[1], + pz = qfitemm$pval[c(1, 3)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint, + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint, + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > runsims7 <- function() { + resl = lapply(1:1000, rft7, n=200) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + apply(res[, 4:6], 2, median) + (bias <- est - truth) + #plot(density(res[, 4])) + #plot(density(res[, 5])) + #plot(density(res[, 6])) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > #runsims7() > > > > > ## bias, coverage, power for true linear model, categorical modifier ## > rft8 <- function(i, n=100) { + dat <- qgcompint:::.dgm_quantized_linear_emm( + N = n, + b0=0, + mainterms= c(0.3, 0.0, 0.3, 0.0), + prodterms = c(0.3, 0.1, 0.1, 0.0), + ztype = "cat", + ncor=0, + corr=0.75 + ) + dat$z = as.factor(dat$z) + (qfit <- qgcomp.noboot(f=y ~ x1 + x2 + x3 + x4, expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + suppressMessages(qfitemm <- qgcomp.emm.glm.noboot(f=y ~ z + x1 + x2 + x3 + x4, emmvar="z", expnms = paste0("x", 1:4), data=dat, q=NULL, family=gaussian())) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[2], + pz = qfitemm$pval[c(2, 4)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint[1], + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint[1], + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > > runsims8 <- function() { + resl = lapply(1:1000, rft8, n=200) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + (bias <- est - truth) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > > #runsims8() > > > ## bias, coverage, power for true Cox model, categorical modifier, balancing modification ## > rft9 <- function(i, n=1000) { + dat = qgcompint:::.dgm_quantized_survival_emm( + N = n, + b0=0, + mainterms=c(1, 0, 0, 0), + prodterms = c(1, 0, 0, 0), + ztype = "cat", + ncor=0, + corr=0.75 + ) + dat$z = as.factor(dat$z) + f = survival::Surv(time, d)~x1 + x2 + x3 + x4 + expnms = c("x1", "x2", "x3", "x4") + #(fit1 <- survival::coxph(f, data = dat)) + (qfit <- qgcomp.cox.noboot(f, expnms = expnms, data = dat, q=4)) + suppressMessages(qfitemm <- qgcomp.emm.cox.noboot(f, expnms = expnms, emmvar="z", data = dat, q=4)) + truth = attr(dat, "truecoef") + c(pmarg = qfit$pval[1], + pz = qfitemm$pval[c(1, 3)], + psixmarg = qfit$psi, + psix0 = qfitemm$psi, + psix1 = qfitemm$psiint[1], + psivarxmarg = qfit$covmat.psi, + psivarx0 = qfitemm$covmat.psi, + psivarx0 = qfitemm$covmat.psiint[1], + truth0 = sum(truth$mainterms), + truthint = sum(truth$prodterms) + ) + } > runsims9 <- function() { + resl = lapply(1:1000, rft9, n=200) + res = as.data.frame(do.call(rbind, resl)) + (truth <- c(NA, mean(res$truth0), mean(res$truthint))) + pow = (res[, 1:3]<0.05)*1.0 + apply(pow, 2, mean) + (est <- apply(res[, 4:6], 2, mean)) + apply(res[, 4:6], 2, median) + (bias <- est - truth) + #plot(density(res[, 4])) + #plot(density(res[, 5])) + #plot(density(res[, 6])) + apply(res[, 4:6], 2, var) + apply(res[, 7:9], 2, mean) + # + varest = res[, 7:9] + ll = res[, 4:6] + sqrt(varest)*qnorm(.05/2) + ul = res[, 4:6] + sqrt(varest)*qnorm(1-.05/2) + cover = (t(ll) < truth) && (t(ul) > truth) + apply(cover, 1, mean) + } > #runsims9() > > proc.time() user system elapsed 4.45 0.67 5.12