Running 'test.R' [168s]
  Running 'test_catmod.R' [158s]
Running the tests in 'tests/test.R' failed.
Complete output:
  > 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.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 = 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 
  
  
  ## 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.0099575 (-0.34389, 0.32398)
  > (qfitemmboot <- qgcomp.emm.boot(f=y ~ z + x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=2, family=gaussian()))
  ## Mixture slope parameters (bootstrap CI):
  
              Estimate Std. Error Lower CI Upper CI t value Pr(>|t|)
  (Intercept)  0.58062    0.13609  0.31389  0.84735  4.2665 0.000107
  psi1        -0.27807    0.20677 -0.68334  0.12720 -1.3448 0.185739
  z           -0.10410    0.16748 -0.43234  0.22415 -0.6216 0.537516
  z:mixture    0.26811    0.25828 -0.23812  0.77434  1.0380 0.305050
  > (qfitemmb <- qgcomp.emm.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 = 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 
  
  
  ## 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.014683 (-0.34524, 0.31587)
  > # 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
  > 
  > dat$z=as.factor(dat$z)
  > (qfitemmf <- qgcomp.emm.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 = 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 
  
  
  ## 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.0099575 (-0.0099575, -0.0099575)
  > 
  > 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.00185)
  x2 
   1 
  
  Scaled effect size (negative direction, sum of negative coefficients = -0.0229)
  x1 
   1 
  
  Mixture slope parameters (Delta method CI):
  
               Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
  (Intercept)  0.548202   0.059166  0.43224  0.66417  9.2654 5.329e-15
  psi1        -0.021011   0.076533 -0.17101  0.12899 -0.2745    0.7843
  > (qfitemm <- qgcomp.emm.noboot(f=y ~ x1 + x2, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=4, 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.0592)
      x1     x2 
  0.9237 0.0763 
  
  
  ## Mixture slope parameters (Delta method CI):
  
               Estimate Std. Error  Lower CI Upper CI t value  Pr(>|t|)
  (Intercept)  0.616013   0.096419  0.427035 0.804990  6.3889 6.754e-09
  psi1        -0.059161   0.061513 -0.179724 0.061402 -0.9618    0.3387
  zb          -0.189183   0.159770 -0.502326 0.123960 -1.1841    0.2394
  zb:mixture   0.123677   0.093073 -0.058744 0.306098  1.3288    0.1872
  zm          -0.142767   0.136765 -0.410820 0.125287 -1.0439    0.2993
  zm:mixture   0.080101   0.083141 -0.082851 0.243053  0.9634    0.3379
  > (qfitemmboot <- qgcomp.emm.boot(f=y ~ x1 + x2, degree=1, emmvar="z", expnms = c('x1', 'x2'), data=dat, q=4, family=gaussian()))
  ## Mixture slope parameters (bootstrap CI):
  
               Estimate Std. Error  Lower CI Upper CI t value  Pr(>|t|)
  (Intercept)  0.616013   0.089368  0.440855 0.791170  6.8930 7.265e-10
  psi1        -0.059161   0.059253 -0.175295 0.056973 -0.9984    0.3207
  zb          -0.189183   0.130213 -0.444396 0.066030 -1.4529    0.1497
  zm          -0.142767   0.136974 -0.411231 0.125697 -1.0423    0.3001
  zb:mixture   0.123677   0.085641 -0.044177 0.291531  1.4441    0.1522
  zm:mixture   0.080101   0.082308 -0.081221 0.241423  0.9732    0.3331
  > (qfitemmb <- qgcomp.emm.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 = 0
  Scaled effect size (positive direction, sum of positive effects = 0)
  None
  Scaled effect size (negative direction, sum of negative effects = -0.129)
     x1    x2 
  0.885 0.115 
  
  
  ## Mixture slope parameters (Delta method CI):
  
               Estimate Std. Error Lower CI Upper CI t value  Pr(>|t|)
  (Intercept)  0.595054   0.073765  0.45048  0.73963  8.0669 2.614e-12
  psi1        -0.129387   0.122968 -0.37040  0.11163 -1.0522    0.2955
  zb          -0.136649   0.125718 -0.38305  0.10975 -1.0869    0.2799
  zb:mixture   0.256635   0.189988 -0.11573  0.62900  1.3508    0.1801
  zm          -0.082950   0.105563 -0.28985  0.12395 -0.7858    0.4340
  zm:mixture   0.110889   0.178245 -0.23847  0.46024  0.6221    0.5354
  > # 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.0592)
      x1     x2 
  0.9237 0.0763 
  
  > getstrateffects(qfitemm, emmval=0.0)
  Joint effect at z=0
           Estimate Std. Error Lower CI Upper CI z value Pr(>|z|)
  Mixture -0.059161   0.061513 -0.17972 0.061402 -0.9618   0.3362
  > getstratweights(qfitemm, emmval=1.0)
  ## Qgcomp weights/partial effects at z = 1
  Scaled effect size (positive direction, sum of positive effects = 0)
  None
  Scaled effect size (negative direction, sum of negative effects = -0.0592)
      x1     x2 
  0.9237 0.0763 
  
  > getstrateffects(qfitemm, emmval=1.0)
  Joint effect at z=1
           Estimate Std. Error Lower CI Upper CI z value Pr(>|z|)
  Mixture -0.059161   0.061513 -0.17972 0.061402 -0.9618   0.3362
  > getstratweights(qfitemmb, emmval=1.0)
  ## Qgcomp weights/partial effects at z = 1
  Scaled effect size (positive direction, sum of positive effects = 0)
  None
  Scaled effect size (negative direction, sum of negative effects = -0.129)
     x1    x2 
  0.885 0.115 
  
  > 
  > 
  > 
  > ##### 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.0527)
     x5    x3    x1 
  0.487 0.363 0.150 
  
  Scaled effect size (negative direction, sum of negative coefficients = -0.0171)
     x4    x2 
  0.524 0.476 
  
  Mixture slope parameters (Delta method CI):
  
              Estimate Std. Error  Lower CI Upper CI t value  Pr(>|t|)
  (Intercept) 0.404928   0.082956  0.242336  0.56752  4.8812 2.196e-06
  psi1        0.035567   0.040379 -0.043574  0.11471  0.8808    0.3795
  > (qfitemm <- qgcomp.emm.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.291)
      x5     x1     x4     x3 
  0.4814 0.3382 0.1208 0.0596 
  Scaled effect size (negative direction, sum of negative effects = -0.0606)
  x2 
   1 
  
  
  ## Mixture slope parameters (Delta method CI):
  
               Estimate Std. Error Lower CI Upper CI t value Pr(>|t|)
  (Intercept)  0.390147   0.143550  0.10879  0.67150  2.7178 0.007185
  psi1         0.230021   0.267845 -0.29494  0.75499  0.8588 0.391552
  z            0.016956   0.124277 -0.22662  0.26053  0.1364 0.891622
  z:mixture   -0.059779   0.231724 -0.51395  0.39439 -0.2580 0.796710
  > 
  > 
  > qfitemm$fit
  
  Call:  glm(formula = newform, family = ..1, data = qdata, weights = weights)
  
  Coefficients:
  (Intercept)           x1           x2           x3           x4           x5  
      0.39015      0.09829     -0.06065      0.01732      0.03513      0.13993  
            z         x1:z         x2:z         x3:z         x4:z         x5:z  
      0.01696     -0.08144      0.05026      0.06546     -0.02742     -0.06665  
  
  Degrees of Freedom: 199 Total (i.e. Null);  188 Residual
  Null Deviance:	    15.11 
  Residual Deviance: 14.4 	AIC: 67.4
  > getstratweights(qfitemm, emmval=0.0)
  ## Qgcomp weights/partial effects at z = 0
  Scaled effect size (positive direction, sum of positive effects = 0.291)
      x5     x1     x4     x3 
  0.4814 0.3382 0.1208 0.0596 
  Scaled effect size (negative direction, sum of negative effects = -0.0606)
  x2 
   1 
  
  > getstratweights(qfitemm, emmval=1.0)
  ## Qgcomp weights/partial effects at z = 1
  Scaled effect size (positive direction, sum of positive effects = 0.181)
      x3     x5     x1     x4 
  0.4583 0.4057 0.0933 0.0427 
  Scaled effect size (negative direction, sum of negative effects = -0.0104)
  x2 
   1 
  
  > getstratweights(qfitemm, emmval=2.0)
  ## Qgcomp weights/partial effects at z = 2
  Scaled effect size (positive direction, sum of positive effects = 0.195)
     x3    x2    x5 
  0.761 0.205 0.034 
  Scaled effect size (negative direction, sum of negative effects = -0.0843)
     x1    x4 
  0.766 0.234 
  
  > getstrateffects(qfitemm, emmval=0.0)
  Joint effect at z=0
          Estimate Std. Error Lower CI Upper CI z value Pr(>|z|)
  Mixture  0.23002    0.26784 -0.29494  0.75499  0.8588   0.3905
  > getstrateffects(qfitemmb, emmval=1.0)
  Joint effect at z=1
          Estimate Std. Error Lower CI Upper CI z value Pr(>|z|)
  Mixture -0.12939    0.12297  -0.3704  0.11163 -1.0522   0.2927
  > getstratweights(qfitemm, emmval=-2.0)
  ## Qgcomp weights/partial effects at z = -2
  Scaled effect size (positive direction, sum of positive effects = 0.624)
     x5    x1    x4 
  0.438 0.418 0.144 
  Scaled effect size (negative direction, sum of negative effects = -0.275)
     x2    x3 
  0.587 0.413 
  
  > getstrateffects(qfitemmb, emmval=-2.0)
  Joint effect at z=-2
          Estimate Std. Error Lower CI Upper CI z value Pr(>|z|)
  Mixture -0.12939    0.12297  -0.3704  0.11163 -1.0522   0.2927
  > 
  > 
  > # -----------------------------------
  > # 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.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.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.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.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.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.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.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.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 
     6.39    0.51    6.90 
Running the tests in 'tests/test_catmod.R' failed.
Complete output:
  > library(qgcomp)
  > library(qgcompint)
  > set.seed(23)
  > dat2 <- simdata_quantized_emm(
  +   outcometype="logistic",
  +   # sample size
  +   n = 100,
  +   # correlation between x1 and x2,x3,...
  +   corr=c(0.8,0.6,0.3,-0.3,-0.3,-0.3),
  +   # model intercept
  +   b0=-2,
  +   # linear model coefficients for x1,x2,... at referent level of interacting variable
  +   mainterms=c(0.3,-0.1,0.1,0.0,0.3,0.1,0.1),
  +   # linear model coefficients for product terms between x1,x2,... and interacting variable
  +   prodterms = c(1.0,0.0,0.0,0.0,0.2,0.2,0.2),
  +   # type of interacting variable
  +   ztype = "categorical",
  +   # number of levels of exposure
  +   q = 4,
  +   # residual variance of y
  +   yscale = 2.0
  + )
  Warning message:
  In .dgm_quantized_logistic_emm(N = n, b0 = b0, mainterms = mainterms,  :
    model implies > 10% of observations with very high/low (<0.5%) outcome probability, which may distort estimates
  > 
  > head(dat2)
    z x1 x2 x3 x4 x5 x6 x7 y
  1 0  1  1  0  3  0  3  2 0
  2 2  2  3  2  1  3  1  1 1
  3 0  3  0  0  3  0  0  2 1
  4 2  1  1  1  0  3  2  2 1
  5 2  2  2  2  2  0  0  1 1
  6 0  0  0  0  0  1  3  3 0
  > table(dat2$z)
  
   0  1  2 
  33 33 34 
  > table(dat2$y)
  
   0  1 
  35 65 
  > dat2$z = as.factor(dat2$z)
  > qfit2 <- qgcomp.emm.noboot(y~x1,
  +                            data = dat2,
  +                            expnms = paste0("x",1:1),
  +                            emmvar = "z",
  +                            q = 4)
  > qfit2
  ## Qgcomp weights/partial effects at z = 0
  Scaled effect size (positive direction, sum of positive effects = 0.112)
  x1 
   1 
  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.038714   0.096934 -0.151274  0.22870  0.3994   0.69051
  psi1        0.112205   0.049806  0.014587  0.20982  2.2529   0.02657
  z1          0.573912   0.138179  0.303086  0.84474  4.1534 7.155e-05
  z1:mixture  0.043462   0.072059 -0.097771  0.18470  0.6031   0.54785
  z2          0.615288   0.138944  0.342962  0.88761  4.4283 2.539e-05
  z2:mixture  0.046243   0.075029 -0.100811  0.19330  0.6163   0.53915
  > qfit2$fit
  
  Call:  glm(formula = newform, data = qdata, weights = weights)
  
  Coefficients:
  (Intercept)           x1           z1           z2        x1:z1        x1:z2  
      0.03871      0.11220      0.57391      0.61529      0.04346      0.04624  
  
  Degrees of Freedom: 99 Total (i.e. Null);  94 Residual
  Null Deviance:	    22.75 
  Residual Deviance: 10.77 	AIC: 74.93
  > 
  > proc.time()
     user  system elapsed 
     2.10    0.35    2.45