Zero Inflated Poisson model
In trajeR: Group Based Modeling Trajectory

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].

data_ZIP[,2:6] is a matrix with real.
data_ZIP[,7:11] is a matrix with time 1 to 10.
data_ZIP[,13:17] is a time covariate that influence the shape of the trajectories. It is a matrix with 0 and 1 value, like the presence or not of a characteristic on the individual.
data_ZIP[,12] is a matrix with real.

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

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
)

library(kableExtra)
t1 <- solL$tab[, 1:2]
t1[11:12, 1] <- exp(t1[11:12, 1]) / sum(exp(t1[11:12, 1]))
t2 <- solEM$tab[, 1:2]
t2[12, 2] <- t2[11, 2]
t3 <- solEMIRLS$tab[, 1:2]
t3[12, 2] <- t3[11, 2]
tab <- round(cbind(t1, t2, t3), 5)
colnames(tab) <- c("parameters", "sd", "parameters", "sd", "parameters", "sd")

kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
  kable_styling() %>%
  add_header_above(c("SolL" = 2, "SolEM" = 2, "SolEMIRLS" = 2)) %>%
  pack_rows("Beta 1", 1, 3) %>%
  pack_rows("Beta 2", 4, 6) %>%
  pack_rows("Nu 1", 7, 8) %>%
  pack_rows("Nu 2", 9, 10) %>%
  pack_rows("Pi", 11, 12)

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

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

t1 <- solLRisk$tab[, 1:2]
t1 <- round(t1, 5)
t2 <- solEMRisk$tab[, 1:2]
t2[, 1] <- round(t2[, 1], 5)
t3 <- round(solEMIRLSRisk$tab[, 1:2], 5)
tab <- cbind(t1, t2, t3)
tab[11:14, 1] <- tab[11:12, 1] - tab[11:14, 1]
tab[11:14, 5] <- tab[11:12, 5] - tab[11:14, 5]
#tab=tab[-11,]
colnames(tab) <- c("parameters", "sd", "parameters", "sd", "parameters", "sd")

kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
  kable_styling() %>%
  add_header_above(c("SolLRisk" = 2, "SolEMRisk" = 2, "SolEMIRLSRisk" = 2)) %>%
  pack_rows("Beta 1", 1, 3) %>%
  pack_rows("Beta 2", 4, 6) %>%
  pack_rows("Nu 1", 7, 8) %>%
  pack_rows("Nu 2", 9, 10) %>%
  pack_rows("Theta", 11, 14)

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

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

t1 <- solLTCOV$tab[, 1]
t1[13:14] <- exp(t1[13:14]) / sum(exp(t1[13:14]))
t2 <- solEMTCOV$tab[, 1]
t3 <- solEMIRLSTCOV$tab[, 1]
tab <- round(cbind(t1, t2, t3), 5)
colnames(tab) <- c("SolLTCOV", "SolEMTCOV", "SoEMIRLSTCOV")

kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
  kable_styling() %>%
  pack_rows("Beta 1", 1, 3) %>%
  pack_rows("Beta 2", 4, 6) %>%
  pack_rows("Nu 1", 7, 8) %>%
  pack_rows("Nu 2", 9, 10) %>%
  pack_rows("Pi", 13, 14) %>%
  pack_rows("TCOV", 11, 12)

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
)