generation of zero-inflated Poisson data in R

Question

I want to generate data from a zero-inflated Poisson distribution in R using the mpath package in the following way:

n <- 50
p <- 20
r <- 0.5    

# true parameter values 
beta <- numeric(p)  
beta[1] <- 1.8
beta[2] <- 0.5
beta[5] <- -0.3

gamma <- numeric(p)
gamma[1] <- 0.4
gamma[2] <- -0.5
gamma[3] <- 0.9

S <- matrix(0, nrow = p, ncol = p)
for(i in 1:p) {
  for(j in i:p) {
    S[i, j] <- r^(j - i)
    S[j, i] <- S[i, j]
  }
} 
X <- mvrnorm(n, rep(0, length(beta)), S)
Z <- mvrnorm(n, rep(0, length(gamma)), S)
bet0 <- matrix(beta, ncol = 1)
gama <- matrix(gamma, ncol = 1)
b <- (exp(Z %*% gama) / (1 + exp(Z %*% gama)))
a <- exp(X %*% bet0)
Y <- rzi(n, X, Z, a, b, family = "poisson", infl = TRUE)

but instead of getting the simulated data an error of the from Error in x %*% a[-1] : non-conformable arguments occur but I am unable to understand or solve it.

Answer 1

The question as, you asked it, is really a programming question and would have been more suitable for StackOverflow. But as it is worth to clarify some of the underlying ideas I'm answering it here with a little bit more context.

rzi() from mpath expects you to supply regressor matrices without intercept but expects that the first element of the coefficient vector is the intercept. (The documentation explains this incorrectly, though.) The intercept is then always included in the model and a logit link is always used for the zero-inflation probabilities.

Your code however had no intercepts at all. And it already computed the full vector of Poisson means $\lambda_i$ and zero-inflation probabilities $\pi_i$. From there it is easy to generate a zero-inflated Poisson variable from first principles rather than a supplied function where you struggled with the correct usage. Hence, I'll explain the latter here using somewhat simplified R code.

First, we fix the dimensions: n is the number of observations and p the number of "real" regressors excluding the intercept:

n <- 50
p <- 20

Then we set up the true coefficients including a zero intercept in both parts:

beta <- c(0, 1.8, 0.5, -0.3, rep(0, p - 3))
gamma <- c(0, 0.4, -0.5, 0.9, rep(0, p - 3))

The real regressors are multivariate normal with standard normal marginals and an AR(1) correlation structure $S_{i, j} = 0.5^{|i - j|}$.

corr <- matrix(0, nrow = 20, ncol = 20)
corr <- 0.5^abs(col(corr) - row(corr))
x <- MASS::mvrnorm(n, rep(0, p), corr)
z <- MASS::mvrnorm(n, rep(0, p), corr)
x <- cbind(1, x)
z <- cbind(1, z)

To compute the expectation $\log(\lambda_i) = x_i^\top \beta$ of the Poisson component a log-link is used (inverse: exp), and for the probabilitiy $\mathrm{logit}(\pi_i) = z_i^\top \gamma$ for zero-inflation a logit link is used (inverse: cumulative distribution function of the logistic distribution).

lambda <- exp(x %*% beta)
pi <- plogis(z %*% gamma)

Then with probability $\pi_i$ we observe a zero and with probability $1 - \pi_i$ we observe something from a Poisson($\lambda_i$) distribution. This can be compactly written using rpois() for the Poisson distribution multiplied with 0 vs. 1 from a Bernoulli distribution rbinom(size = 1).

y <- rpois(n, lambda = lambda) * rbinom(n, prob = 1 - pi, size = 1)