library(COMPoissonReg)
set.seed(1234)
# ----- Generate the data -----
n = 400
x = runif(n, 2, 4)
X = model.matrix(~ x)
S = matrix(1, n, 1)
W = model.matrix(~ x)
beta.true = c(0.05, 0.50)
gamma.true = -0.3
zeta.true = c(-2.5, 0.05)
lambda.true = exp(X %*% beta.true)
nu.true = exp(S %*% gamma.true)
p.true = plogis(W %*% zeta.true)
y = rzicmp(n, lambda = lambda.true, nu = nu.true, p = p.true)
dat = data.frame(y = y, x = x)
# ----- Fit ZICMP model -----
options(COMPoissonReg.optim.method = "BFGS")
zicmp.out = glm.cmp(y ~ x, formula.nu = ~ 1, formula.p = ~ x, data = dat)
print(zicmp.out)
# ----- Residuals -----
y.hat = predict(zicmp.out)
res = resid(zicmp.out, type = "quantile")
qqnorm(res); qqline(res, lty = 2, col = "red", lwd = 2)
plot(y.hat, res)
# ----- Test for equidispersion -----
equitest(zicmp.out)
# ----- Deviance -----
tryCatch({
deviance(zicmp.out)
}, error = function(e) {
print(e)
})
# ----- Bootstrap -----
boot.out = parametric.bootstrap(zicmp.out, reps = 50, report.period = 10)
plot(density(boot.out[,1])); abline(v = beta.true[1], lty = 2, lwd = 2, col = "red")
plot(density(boot.out[,2])); abline(v = beta.true[2], lty = 2, lwd = 2, col = "red")
plot(density(boot.out[,3])); abline(v = gamma.true[1], lty = 2, lwd = 2, col = "red")
plot(density(boot.out[,4])); abline(v = zeta.true[1], lty = 2, lwd = 2, col = "red")
plot(density(boot.out[,5])); abline(v = zeta.true[2], lty = 2, lwd = 2, col = "red")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.