---
title: "COMPoissonReg Demo"
date: "`r format(Sys.time(), '%B %-e, %Y')`"
output: pdf_document
---

## COM Poisson Regression

Generate some data.

```{r}
library(COMPoissonReg, quietly = TRUE)
set.seed(1235)

n = 400
x = runif(n, 1, 4)
X = model.matrix(~ x)
beta.true = c(1, 0.5)
lambda.true = exp(X %*% beta.true)
nu.true = 0.75
y = rcmp(n, lambda = lambda.true, nu = nu.true)
dat = data.frame(y = y, x = x)
```

Fit the model. The package notices that we specify only an intercept for $\nu$, and also displays an estimate for $\nu$ directly by using the appropriate transformation.

```{r}
cmp.out = glm.cmp(y ~ x, data = dat)
print(cmp.out)
```

Plot the randomized quantile residuals.

```{r, fig.width=4, fig.height=4}
y.hat = predict(cmp.out)
res = resid(cmp.out, type = "quantile")
qqnorm(res); qqline(res, lty = 2, col = "red", lwd = 2)
plot(y.hat, res)
```

## Zero-Inflated COM Poisson Regression

```{r, echo=FALSE}
rm(list = ls())
```


Generate some data.

```{r}
library(COMPoissonReg, quietly = TRUE)
set.seed(1235)

n = 400
x = runif(n, 1, 4)
X = model.matrix(~ x)
S = matrix(1, n, 1)
W = model.matrix(~ x)
beta.true = c(1, 2)
gamma.true = 1
zeta.true = c(0.05, -1)
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 the model. The package notices that we specify only an intercept for $\nu$, and also displays an estimate for $\nu$ directly by using the appropriate transformation.

```{r}
options(COMPoissonReg.optim.method = 'BFGS')
zicmp.out = glm.cmp(y ~ x, formula.nu = ~ 1, formula.p = ~ x, data = dat)
print(zicmp.out)
```

Plot the randomized quantile residuals.

```{r, fig.width=4, fig.height=4}
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)
```