--- title: "Zero Inflated Poisson model" author: "Cédric NOEL - Jang SHILTZ" email: "cedric.noel@univ-lorraine.fr" date: "`r Sys.Date()`" output: rmarkdown::html_vignette bibliography: - biblio.bib vignette: > %\VignetteIndexEntry{Zero Inflated Poisson model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- 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. ```{r} 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"̀ . ```{r, message = FALSE} # 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. ```{r, message = FALSE} # 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 ) ``` ```{r, echo = FALSE} 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. ```{r} 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. ```{r} # 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 =`. ```{r, message = FALSE} 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. ```{r, message = FALSE} 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 ```{r, echo = FALSE} 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. ```{r, message = FALSE} 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. ```{r, message = FALSE} 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 ```{r, echo = FALSE} 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. ```{r, message = FALSE} 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 ) ```