Zero Inflated Poisson model

Cédric NOEL - Jang SHILTZ

2026-08-03

The ZIP (Zero Inflated Poisson) employ two different process : a binary distribution that generate structural zero, i.e. the excess zeros, and and Poisson distribution that generate the counts.

Under this model, zero can occur for two reasons. The first is the counts is equal to zero with a probability \(P(Y_{it}=0)\) for \(Y_{it}\sim \mathcal{P}(\lambda_{ikt})\) and the second become of the fact the binary produce zero values with a probability \(\rho_{ikt}\).

We have \[\begin{align} P\left(Y_{it}=y_{it}|W_i=wi, C_i=ci\right) = \left\lbrace \begin{array}{l} \rho_{ikt}+(1-\rho_{ikt})e^{-\lambda_{ikt}},\ y_{it}=0 \\ (1-\rho_{ikt})\frac{\lambda_{ikt}^{y_{it}}e^{-\lambda_{ikt}}}{y_{it}!},\ y_{it}> 0 \end{array} \right. \end{align}\]

Parameters

We use artificial data contained in the library. There are 500 longitudinal data generated by the parameters above : + 2 groups ; + probability is \(\pi_1=0.3\) and \(\pi_2=0.7\) ; + period is 5 ; + we use 2 polynomial shape to calculate the parameters \(\lambda\) of the Poisson state - degree 2 and \(\beta_1=(1.2, 0.5, -0.06)\) ; - degree 2 and \(\beta_2=(0.89, 0.01, 0.01)\) ; + we use 2 polynomial shape to calculate the parameters \(\nu\) of the zero state - degree 1 and \(\nu_1=(-0.2,-0.1)\) ; - degree 1 and \(\nu_2=(-1,0)\) ;

The data contained the time variable dependent \(Y_i\) in data_ZIP[,2:6] , the time variable in data_ZIP[,7:11], the time covariate in data_ZIP[,13:17] and a covariate that influence the belonging probability in data_ZIP[,12].

library(trajeR)

data("data_ZIP")

matplot(
  t(data_ZIP[, 7:11]),
  t(data_ZIP[, 2:6]),
  pch = 1,
  type = 'b',
  col = "black",
  lty = 1,
  xlab = "Times",
  ylab = "Values",
  main = "Plot of the individual's trajectories"
)

We use each method to fit the model. For all methods, we specify the number of group of our model, ng=2 and the degree of the polynomial shape of the trajectories for the two state. For the Poisson state we choice a 3 degree polynomial and for the zero state a 2 degree polynomial. Thus for \(\beta\) degree is a vector \((3,3)\) and for \(\nu\) degree is a vector \((2,2)\). We specify hessian=TRUE to ask the calculus of the hessian matrix.

For the Likelihood method we call ̀trajeR with option `Method =“L”̀ .

# Likelihood
solL <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "L",
  hessian = TRUE
)
solL
## Call TrajeR with 2 groups and a 2,2 degrees of polynomial shape of trajectory.
## Model : Zero Inflated Poisson
## Method : Likelihood 
## Convergence : Successful (Code 4)
## 
##   group   Parameter   Estimate   Std. Error   T for H0:   Prob>|T|  
##                                                param.=0             
## --------------------------------------------------------------------
##       1   Intercept    0.90782      0.09666     9.39219          0  
##              Linear    0.01872      0.07171     0.26111    0.79403  
##           Quadratic    0.00618      0.01156      0.5348    0.59283  
##                Nu11   -1.08146      0.16253    -6.65409          0  
##                Nu12    0.00454      0.04744     0.09568    0.92378  
##                                                                     
##       2   Intercept    1.04873      0.09549    10.98323          0  
##              Linear    0.61507      0.06494     9.47072          0  
##           Quadratic   -0.07773      0.01006    -7.72673          0  
##                Nu21    0.03162      0.16837      0.1878    0.85105  
##                Nu22   -0.16776      0.05219    -3.21428    0.00132  
## --------------------------------------------------------------------
##       1         pi1    0.65132       0.0227           0          0  
##       2         pi2    0.34868       0.0227   -27.53258          0  
## --------------------------------------------------------------------
## Likelihood : -5162.009

For use EM method we write the same syntax but change the variable Method to EM or ̀EMIRLS. In the first case, we find the parameters with quasi newton method and in the second one with Iterative Reweighted Least Square.

# EM
solEM <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EM",
  hessian = TRUE
)
solEMIRLS <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EMIRLS",
  hessian = TRUE
)
SolL
SolEM
SolEMIRLS
parameters sd parameters sd parameters sd
Beta 1
0.90782 0.09666 0.90736 0.07865 0.90782 0.07864
0.01872 0.07171 0.01897 0.05853 0.01872 0.05853
0.00618 0.01156 0.00615 0.00947 0.00618 0.00947
Beta 2
1.04873 0.09549 1.04876 0.07417 1.04873 0.07417
0.61507 0.06494 0.61499 0.04979 0.61507 0.04979
-0.07773 0.01006 -0.07772 0.00769 -0.07773 0.00769
Nu 1
-1.08146 0.16253 -1.08176 0.16536 -1.08146 0.16534
0.00454 0.04744 0.00461 0.04874 0.00454 0.04874
Nu 2
0.03162 0.16837 0.03159 0.18424 0.03162 0.18424
-0.16776 0.05219 -0.16775 0.05663 -0.16776 0.05663
Pi
0.65132 0.02270 0.65132 0.02151 0.65132 0.02151
0.34868 0.02270 0.34868 0.02151 0.34868 0.02151

Plot

We plot only the trajectories on the graph. We have just to use the native plot function. It use the class of the object to draw the polynomial shape. The values in abscissa are the first row of the time variable.

plotrajeR(solL)

We can add longitudinal data to this graph. If we want this on the plot we have to specify Y and A in the function plot(). By default colors are gray scale, but we can specify colors we want.

# colour's defintion
trans <- "70"
col1 <- "#034569"
col1.1 <- paste0("#64AAD0", trans)
col2 <- "#750062"
col2.1 <- paste0("#D962C7", trans)
cols1 <- c(col1.1, col2.1)
cols2 <- c(col1, col2)
vcol <- c(cols1, cols2)

plotrajeR(solEM, Y = data_ZIP[, 2:6], A = data_ZIP[, 7:11], dec = 5, col = vcol)

Adding covariate

In the same way that in the section below we can add covariate in the calculus of the parameters.

If we add risk variable we use option Risk =.

solLRisk <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  Risk = data_ZIP[, 12],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "L",
  hessian = TRUE
)
solLRisk
## Call TrajeR with 2 groups and a 2,2 degrees of polynomial shape of trajectory.
## Model : Zero Inflated Poisson
## Method : Likelihood 
## Convergence : Successful (Code 1)
## 
##   group   Parameter   Estimate   Std. Error   T for H0:   Prob>|T|  
##                                                param.=0             
## --------------------------------------------------------------------
##       1   Intercept    0.90742      0.09678     9.37572          0  
##              Linear    0.01908      0.07143     0.26703    0.78946  
##           Quadratic     0.0061      0.01152      0.5293    0.59664  
##                Nu11   -1.08237      0.16313    -6.63501          0  
##                Nu12    0.00462      0.04723     0.09787    0.92204  
##                                                                     
##       2   Intercept    1.04842      0.09471     11.0698          0  
##              Linear    0.61501      0.06492     9.47334          0  
##           Quadratic   -0.07772      0.01012    -7.67753          0  
##                Nu21    0.03215      0.16975     0.18942    0.84978  
##                Nu22   -0.16784      0.05211    -3.22118    0.00129  
## --------------------------------------------------------------------
##       1    Baseline          0           NA          NA         NA  
##                                                                     
##       2   Intercept   -0.81597      0.19942    -4.09172      4e-05  
##                   X    0.37652      0.33782     1.11455    0.26515  
## --------------------------------------------------------------------
## Likelihood : -5161.373

We can use all method too.

solEMRisk <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  Risk = data_ZIP[, 12],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EM",
  hessian = TRUE
)
solEMIRLSRisk <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  Risk = data_ZIP[, 12],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EMIRLS",
  hessian = TRUE
)

We have

SolLRisk
SolEMRisk
SolEMIRLSRisk
parameters sd parameters sd parameters sd
Beta 1
0.90742 0.09678 1.04813 0.0741151 1.04842 0.07411
0.01908 0.07143 0.61515 0.0497785 0.61501 0.04978
0.00610 0.01152 -0.07773 0.0076921 -0.07772 0.00769
Beta 2
1.04842 0.09471 0.90700 0.0786137 0.90742 0.07861
0.61501 0.06492 0.01937 0.0585259 0.01908 0.05852
-0.07772 0.01012 0.00605 0.0094679 0.00610 0.00947
Nu 1
-1.08237 0.16313 0.03212 0.1842311 0.03215 0.18423
0.00462 0.04723 -0.16783 0.0565987 -0.16784 0.05660
Nu 2
0.03215 0.16975 -1.08249 0.1655697 -1.08237 0.16556
-0.16784 0.05211 0.00465 0.0487714 0.00462 0.04877
Theta
0.00000 NA 0.00000 0.0000000 0.00000 0.00000
0.00000 NA 0.00000 0.0000000 0.00000 0.00000
0.81597 0.19942 0.81595 0.1978894 -0.81597 0.19789
-0.37652 0.33782 -0.37653 0.3341130 0.37652 0.33411

Adding time depdendent covariate

If we add a time-dependent covariate, we use the TCOV option.

solLTCOV <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  TCOV = data_ZIP[, 13:17],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "L",
  hessian = TRUE
)
solLTCOV 
## Call TrajeR with 2 groups and a 2,2 degrees of polynomial shape of trajectory.
## Model : Zero Inflated Poisson
## Method : Likelihood 
## Convergence : Successful (Code 4)
## 
##   group   Parameter   Estimate   Std. Error   T for H0:   Prob>|T|  
##                                                param.=0             
## --------------------------------------------------------------------
##       1   Intercept    0.88074      0.10032     8.77911          0  
##              Linear    0.01538      0.07234     0.21258    0.83167  
##           Quadratic    0.00677      0.01171     0.57809    0.56326  
##                Nu11   -1.07079      0.15832    -6.76342          0  
##                Nu12    0.00234      0.04563     0.05122    0.95915  
##               TCOV1    0.06375      0.04011     1.58937     0.1121  
##                                                                     
##       2   Intercept    1.04295      0.09688    10.76531          0  
##              Linear    0.61469      0.06573     9.35249          0  
##           Quadratic   -0.07767      0.01021    -7.60871          0  
##                Nu21    0.02789      0.16772     0.16629    0.86794  
##                Nu22   -0.16707      0.05278    -3.16564    0.00157  
##               TCOV1    0.01297       0.0334     0.38836    0.69778  
## --------------------------------------------------------------------
##       1         pi1    0.65198      0.02298           0          0  
##       2         pi2    0.34802      0.02298   -27.31879          0  
## --------------------------------------------------------------------
## Likelihood : -5160.531

We can use all method too.

solEMTCOV <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  TCOV = data_ZIP[, 13:17],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EM",
  hessian = TRUE
)
solEMIRLSTCOV <- trajeR(
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  TCOV = data_ZIP[, 13:17],
  ng = 2,
  degre = c(2, 2),
  degre.nu = c(1, 1),
  Model = "ZIP",
  Method = "EMIRLS",
  hessian = TRUE
)

We have

SolLTCOV SolEMTCOV SoEMIRLSTCOV
Beta 1
0.88074 0.88058 0.88074
0.01538 0.01546 0.01538
0.00677 0.00676 0.00677
Beta 2
1.04295 1.04209 1.04295
0.61469 0.61497 0.61469
-0.07767 -0.07769 -0.07767
Nu 1
-1.07079 -1.07096 -1.07079
0.00234 0.00237 0.00234
Nu 2
0.02789 0.02777 0.02789
-0.16707 -0.16703 -0.16707
TCOV
0.06375 0.06379 0.06375
0.01297 0.01319 0.01297
Pi
0.65198 0.65195 0.65198
0.34802 0.34805 0.34802

We can add the effect of other time covariate.

trans <- "70"
col1 <- "#034569"
col1.1 <- paste0("#64AAD0", trans)
col2 <- "#750062"
col2.1 <- paste0("#D962C7", trans)
col3 <- "#A68900"
col3.1 <- paste0("#FFE773", trans)
cols1 <- c(col1.1, col2.1, col3.1)
cols2 <- c(col1, col2, col3)
cols1 <- c(col1.1, col2.1)
cols2 <- c(col1, col2)
vcol <- c(cols1, cols2)
plotrajeR(
  solLTCOV,
  Y = data_ZIP[, 2:6],
  A = data_ZIP[, 7:11],
  TCOV = data_ZIP[, 13:17],
  col = vcol,
  plotcov = c(1, 1, 0, 0, 1, 0, 0, 1, 1, 1),
  mean = TRUE,
  alpha = 0.75
)