# NOT RUN {
set.seed(50)
n=100
dat <- data.frame(y=rbinom(n, 1, 0.5)*rpois(n, 1.2), x1=runif(n), x2=runif(n), z=runif(n))
# poisson count model, mixture in both portions
# }
# NOT RUN {
# warning: the examples below can take a long time to run
res = qgcomp.zi.boot(f=y ~ x1 + x2 | x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="poisson", B=1000, MCsize=10000, parallel=TRUE)
qgcomp.zi.noboot(f=y ~ x1 + x2 | x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="poisson")
res
# accuracy for small MCsize is suspect (compare coefficients between boot/noboot versions),
# so re-check with MCsize set to larger value (this takes a long time to run)
res2 = qgcomp.zi.boot(f=y ~ x1 + x2 | x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="poisson", B=1000, MCsize=50000, parallel=TRUE)
res2
plot(density(res2$bootsamps[4,]))
# negative binomial count model, mixture and covariate in both portions
qgcomp.zi.boot(f=y ~ z + x1 + x2 | z + x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="negbin", B=10, MCsize=10000)
# weighted analysis (NOTE THIS DOES NOT WORK WITH parallel=TRUE!)
dat$w = runif(n)*5
qgcomp.zi.noboot(f=y ~ z + x1 + x2 | x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="poisson", weights=w)
# Expect this:
# Warning message:
# In eval(family$initialize) : non-integer #successes in a binomial glm!
qgcomp.zi.boot(f=y ~ x1 + x2 | x1 + x2, expnms = c('x1', 'x2'),
data=dat, q=4, dist="poisson", B=5, MCsize=50000, parallel=FALSE, weights=w)
# Log rr per one IQR change in all exposures (not on quantile basis)
dat$x1iqr <- dat$x1/with(dat, diff(quantile(x1, c(.25, .75))))
dat$x2iqr <- dat$x2/with(dat, diff(quantile(x2, c(.25, .75))))
# note that I(x>...) now operates on the untransformed value of x,
# rather than the quantized value
res2 = qgcomp.zi.boot(y ~ z + x1iqr + x2iqr + I(x2iqr>0.1) + I(x2>0.4) + I(x2>0.9) | x1iqr + x2iqr,
family="binomial", expnms = c('x1iqr', 'x2iqr'), data=dat, q=NULL, B=2,
degree=2, MCsize=200, dist="poisson")
res2
# }
Run the code above in your browser using DataLab