This function fit a Zero Inflated Bivariate Poisson model under Holgate (1964) and Geoffroy et. al (2021).
ZIBP_HOL(
l1.fo,
l2.fo,
psi.fo,
data,
initial.values = NULL,
optimizer = "nlminb"
)a formula object to explain the parameter \(\lambda_1\), the first response \(X_1\) on the left of an ~ operator, and the terms, separated by operators, on the right.
a formula object to explain the parameter \(\lambda_2\), the second response \(X_2\) on the left of an ~ operator, and the terms, separated by operators, on the right.
a formula object to explain the \(\psi\) parameter. An example could be: ~ age + salary.
a data frame containing the variables occurring in the formula.
a vector with possible initial values for the parameters.
the optimizer to be used: nlminb, optim.
Returns a object with class "ZIBP_HOL".
Kouakou, K. J. G., Hili, O., & Dupuy, J. F. (2021). Estimation in the zero-inflated bivariate Poisson model with an application to health-care utilization data. Afrika Statistika, 16(2), 2767-2788.
# Example 1 ---------------------------------------------------------------
# Estimating parameters as a regression model
l1 <- 5
l2 <- 3
psi <- 0.20
l0 <- 1.20
set.seed(12345)
data1 <- rZIBP_HOL(n=500, l1=l1, l2=l2, l0=l0, psi=psi)
data1 <- as.data.frame(data1)
colnames(data1) <- c("y1", "y2")
mod1 <- ZIBP_HOL(l1.fo=y1~1,
l2.fo=y2~1,
psi.fo=~1,
data=data1)
summary(mod1)
#> ---------------------------------------------------------------
#> Fixed effects for log(l1)
#> ---------------------------------------------------------------
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.552910 0.054075 28.718 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Fixed effects for log(l2)
#> ---------------------------------------------------------------
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.043192 0.086603 12.046 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Fixed effects for logit(psi)
#> ---------------------------------------------------------------
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -1.26625 0.10801 -11.724 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Estimation for l0
#> ---------------------------------------------------------------
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.39196 0.23802 5.8482 4.969e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
# To obtain the estimated parameters
c(exp(mod1$par[1:3]), mod1$par[4])
#> l1_hat l2_hat l0_hat psi_hat
#> 4.7251993 2.8382636 0.2818859 1.3919595
moments_estim_ZIBP_HOL(data1)
#> l1_hat l2_hat l0_hat psi_hat
#> 4.8858732 3.3560129 1.1656262 0.2114351
# Example 2 ---------------------------------------------------------------
# Replicating application from section 5.1 from Geoffroy et. al (2021).
if (FALSE) { # \dontrun{
library(AER)
data("NMES1988")
# To obtain some statistics as in the paper
y1 <- NMES1988$nvisits
y2 <- NMES1988$novisits
cor(y1, y2)
cor.test(y1, y2)
cov(y1, y2)
mean(y1==0 & y2==0)
# Model as in section 5.1 from Geoffroy et. al (2021).
mod2 <- NULL
mod2 <- ZIBP_HOL(l1.fo=nvisits ~health+chronic+age+gender+
married+school+income+medicaid,
l2.fo=novisits~health+chronic+age+gender+
married+school+income+medicaid,
psi.fo=~chronic+gender+school+medicaid,
data=NMES1988)
summary(mod2)
} # }